Numerical calculation model, system and method for hydraulic transition process of pumped storage power station

Abstract

The invention discloses a pumped storage power station hydraulic transition process numerical calculation model, system and method, and belongs to the field of water conservancy and hydropower engineering. The model comprises a runner characteristic curve conversion sub-model, a pipeline system parameter calculation sub-model and a transition process main calculation sub-model, and the runner characteristic curve conversion model adopts double algorithms of a Suter method and a length method to solve the problem of connection of small-opening curves; the pipeline system parameter calculation model outputs basic parameters such as water hammer wave velocity; and the transition process main calculation model integrates a plurality of sub-models to realize all-condition dynamic solution. The system comprises a pre-processing module, a main calculation module and a post-processing module and is used for completing data standardization, calculation and result analysis. The method is executed according to data preparation and pretreatment, runner characteristic curve conversion, pipeline system parameter calculation, transition process main calculation, post-treatment and result analysis. The method adapts to a complex system and multiple working conditions, is high in calculation precision and high in intelligence, and can provide technical support for power station design optimization, equipment type selection and safe operation.

Description

Numerical Calculation Model, System and Method for Hydraulic Transient Process in Pumped Storage Power Stations Technical Field

[0001] This invention relates to the field of water conservancy and hydropower engineering technology, and in particular to a numerical calculation model, system, and method for the hydraulic transient process of pumped storage power stations based on the Suter method. It is applicable to the simulation calculation of multi-condition transient processes in mixed-flow reversible pump-turbine pumped storage power stations and conventional hydropower stations, providing technical support for power station design optimization, equipment selection, and safe operation. Background Technology

[0002] With the adjustment of energy structure and the increasing demand for power system flexibility, pumped storage power stations are rapidly developing towards large capacity, high parameters, and complex systems. The hydraulic transient process, as a dynamic process involving the coupling of multiple systems (water flow, machinery, and electrical systems) during power station operation, directly determines the stability design, equipment selection rationality, and operational economy of the power station based on the accuracy of its calculation results. It is a core technical aspect of power station construction and operation and maintenance.

[0003] Existing transient process calculation programs suffer from the following technical defects, making it difficult to meet the engineering requirements of complex pumped storage power stations: Limited applicability: Most programs are only developed for single-type power stations such as conventional hydropower stations, pumped storage power stations with specific structures, or simple pipeline systems. They cannot adapt to complex scenarios involving multiple units operating in parallel or multiple surge tanks, resulting in poor versatility. Insufficient calculation accuracy: The transformation of the runner's full characteristic curve is the core foundation of transient process calculation. Existing programs mostly use a single algorithm to process the characteristic curve. Under small opening conditions, especially guide vane opening ≤0.1 or 0, problems such as curve connection breaks and excessive interpolation errors easily occur, leading to calculation results that are inconsistent with actual operating conditions. Significant deviations; low data processing efficiency: input parameters involve multiple dimensions such as impeller characteristics, pipeline system, and operating condition settings, resulting in a large amount of data with complex formats. Existing programs lack effective automated error detection mechanisms, relying on manual checks for missing data and logical contradictions, which is inefficient and prone to omissions; insufficient intelligence in result analysis: calculation results require professionals to extract key information through complex data processing, lacking automatic stability assessment, extreme value parameter summarization, and visualization functions, making engineering applications difficult; poor ease of operation: it requires users to have professional knowledge of hydraulic machinery, numerical simulation, etc., and the process is cumbersome, making it difficult to meet the needs of efficient iteration in engineering design.

[0004] Compared to conventional hydropower stations, pumped storage power stations must simultaneously consider combined operating conditions such as turbine operation (e.g., power generation), pump operation (e.g., pumping, start-up, shutdown, and load shedding). The transformation of the turbine's full characteristic curves is more complex, and the variety of surge tank types (e.g., overflow, impedance) leads to more prominent coupling characteristics in the hydraulic transient process and greater computational difficulty. Therefore, there is an urgent need to develop a numerical calculation model and system that is widely applicable, highly accurate, intelligent, and easy to operate, to address the pain points of existing technologies. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the prior art and provide a numerical calculation model, system and calculation method for the hydraulic transient process of pumped storage power stations based on the Suter method. Through multi-module collaborative design and dual-algorithm fusion architecture, it can achieve accurate calculation of the transient process under multiple working conditions and complex systems, reduce the operation threshold and improve the practicality of engineering.

[0006] To achieve the above objectives, the technical solution adopted by this invention is as follows: I. Numerical Calculation Model This model is based on the Suter method, integrating the advantages of the length method and the characteristic line method to construct a multi-module collaborative calculation system, specifically including three core sub-models. Each sub-model is functionally independent but logically closely related: 1. Runner Characteristic Curve Conversion Model The model adopts a dual-algorithm architecture combining the Suter method and the length method to achieve accurate conversion and connection of the runner's full characteristic curves: The Suter method is the TRN program module: it converts the runner experimental data, including guide vane opening, unit speed N, etc. 11 Unit flow rate Q 11 Unit torque M 11 The data is converted to standardized data in WH(x,y) and WB(x,y) forms. Equal spacing is achieved through parabolic and linear interpolation, converting the experimental data into standardized data in WH(x,y) and WB(x,y) forms, generating an MPSI1.D data file for the main computational system, adapting to the computational needs of guide vane large opening conditions such as opening > 0.1.

[0007] The length method, or NQM program module, uses the length Lq of the flow characteristic curve as a reference to transform the impeller characteristic curve into Lq-N. 11 Lq-Q 11 Lq-M 11 The associated curve is used to generate the MPSI5.D data file, which is specifically designed to solve the problem of unsmooth connection of characteristic curves under small opening conditions such as guide vane opening ≤0.1 and 0 opening conditions, and ensures that the interpolation error in the full opening range is ≤0.5%.

[0008] The typical characteristic curve is α-N 11 α-Q 11 α-M 11That is, the guide vane opening α is used as the reference, but when the opening is small, α ≤ 0.1, α changes sensitively and the curve is easily broken; while the flow characteristic curve length Lq and Q 11 Directly related, Q 11 The larger the value, the longer Lq becomes. Within the small opening range, Lq changes more smoothly. Therefore, Lq is used instead of α as a unified reference, transforming the three independent curves into a curve with Lq as the abscissa and N as the ordinate. 11 / Q 11 / M 11 The MPSI5.D data file uses the length of the flow characteristic curve Lq as the uniform horizontal axis, storing Lq-N as the correlation curve on the vertical axis to avoid breaks in the connection at small openings. 11 Lq-Q 11 Lq-M 11 Standardized data for three sets of correlation curves.

[0009] 2. The pipeline system parameter calculation model realizes the standardized processing and calculation of pipeline system parameters through the READY program module: Input parameters: equivalent pipe diameter, pipe length, pipe wall thickness, friction loss coefficient, local loss coefficient, elastic modulus of pipe material, elastic modulus of surrounding rock, water density and other basic physical parameters.

[0010] Core calculations: The water hammer wave velocity is adjusted by considering the coupling effect of the elastic modulus of the pipeline material, surrounding rock, and water body. An adaptive segmentation algorithm is used to divide the pipe segments and the calculation nodes. Based on the Courant condition, Δt≤Δx / a, where Δx is the length of the pipe segment and a is the water hammer wave velocity, the calculation time step Δt is determined to avoid numerical calculation oscillations.

[0011] Output data: Generates data files MPSI2.D and MPSI4.D. Data file MPSI2.D includes water hammer wave velocity of pipe section and calculation time step, while data file MPSI4.D includes elevation of each calculation node and node number. It also supports automatic generation of pipeline layout diagrams and error detection function for reasonable input parameters.

[0012] 3. The main calculation model for the transient process is based on solving the transient flow control equations of the pipeline using the method of characteristics. It integrates a multi-physics coupled sub-model to realize numerical calculation of the transient process under all operating conditions. Specifically, it includes: Governor model: compatible with PI and PID governors, supports guide vane break line closing / opening law settings, and allows customization of key parameters such as permanent droop rate, transient droop rate, buffer time constant, and acceleration time constant to adapt to the dynamic response characteristics of different types of speed control systems.

[0013] Transient model of surge tank: compatible with eight types of surge tanks, including overflow type, cylindrical type, cylindrical type with upper and lower chambers, impedance type, differential type, double chamber type, impedance type with overflow weir, and cylindrical type with riser pipe. By establishing the control equation for surge tank water level fluctuation, the model accurately calculates the surge tank water level change curve and the dynamic change of bottom pressure during the transition process.

[0014] Unit boundary model: Boundary control equations are established for turbine and pump operating conditions respectively. The dynamic coupling relationship between guide vane opening, unit speed and flow rate is considered. Influencing factors such as unit rotational inertia and electromagnetic torque are introduced to ensure that the boundary conditions are consistent with the actual operating conditions.

[0015] Operating conditions applicable: Fully covers turbine operating conditions, pump operating conditions, and combined operating conditions. Turbine operating conditions include sudden partial load shedding, sudden 100% rated load shedding, hydraulic interference, guide vane failure to move, and emergency pressure regulating valve operation; pump operating conditions include power failure shutdown, hydraulic interference, and guide vane failure to move; combined operating conditions include start-up-load shedding and start-up-power failure.

[0016] II. Numerical Calculation System The numerical calculation system of the present invention is the hardware carrier and software implementation of the above-mentioned numerical calculation model, including: hardware layer: at least one processor; and memory, input device and output device communicatively connected to the at least one processor.

[0017] Software layer: Computer programs stored in memory and executable by processor. When executed by processor, the computer programs implement the computational logic of the above numerical calculation model. The software layer is divided into three major modules according to function, and the logical connection between the modules is as follows: 1. Preprocessing module: including data formatting program (MMM) and intelligent error detection program (TYJ). MMM program is a data format mapping program that converts input parameters from different sources and in different formats, such as discrete tables of experimental data and design values ​​of engineering parameters, into text files with a unified structure, i.e., the data format commonly used in Fortran programs (.D file). The preprocessing module takes as input experimental data of the turbine runner characteristics, pipeline system parameters, unit parameters, and operating condition setting parameters. The MMM program converts the input data into a standardized format and generates intermediate data files such as TRNI1.D, TRNI2.D, NQMI.D, READYI.D, and MPSI3.D. The TYJ program then performs integrity checks on the intermediate data files, such as detecting missing parameters, performing logical checks, such as checking the monotonicity of the surge tank elevation data, and judging the rationality and validity of pipeline parameters, such as determining whether the initial operating condition parameters such as flow rate and guide vane opening are within a reasonable range, and outputs abnormal prompt information.

[0018] 2. Main Calculation Module: Integrating the above numerical calculation model, the MPS main program calls the TRN, NQM, READY, and subroutines such as READ.FOR, INTE.FOR, PUMP.FOR, and SURGE.FOR to read intermediate data files generated by the preprocessing module. It completes the conversion of the runner characteristic curve, the calculation of pipeline system parameters, and the main calculation of the transient process, outputting MPSO1.D (time series data) containing parameters such as pressure, flow rate, rotational speed, and opening degree at each time point; and MPSO3.D (extreme value parameter data) and other result files. The main calculation module also has a built-in abnormal operating condition handling mechanism. When the surge tank experiences roof collapse (i.e., water level exceeding the design upper limit) or bottom leakage (i.e., water level falling below the design lower limit), insufficient or abnormal characteristic curve data, or calculation non-convergence, it automatically outputs warning information and prompts adjustment schemes.

[0019] Among them, the TRN program, as a module for converting the runner characteristic curve model, uses the Suter method to convert runner experimental data such as guide vane opening and unit speed N. 11 Unit flow rate Q 11 Unit torque M 11 Convert the data to standardized forms WH(x,y) and WB(x,y) to adapt to large opening conditions (opening > 0.1), and output the file MPSI1.D.

[0020] NQM program: Employs the length method, using the length Lq of the flow characteristic curve as a reference, to transform the characteristic curve into Lq-N. 11 Lq-Q 11 Lq-M 11 The associated curve solves the curve connection problem for small opening conditions, i.e., opening ≤0.1 and 0 opening conditions. The output file is MPSI5.D.

[0021] The READY program serves as the implementation vehicle for the pipeline system parameter calculation model. It explicitly inputs basic parameters such as pipe segment diameter, length, wall thickness, elastic modulus, and water density. Through wave velocity adjustment, pipe segment segmentation, and Courant condition time step calculation, it outputs pipe segment wave velocity MPSI2.D data files and node elevation MPSI4.D data files, and supports pipeline layout diagram generation and parameter error detection.

[0022] The READ.FOR subroutine reads standardized data files such as MPSI1.D and MPSI2.D and passes the data into the main calculation process. It is a regular data interface module of the numerical calculation program.

[0023] The INTE.FOR subroutine is used to perform numerical integration of the transient flow equation and calculate the dynamic changes in pressure and flow at pipe nodes. It is a standard numerical implementation module of the method of characteristics.

[0024] The PUMP.FOR subroutine is used to establish the unit boundary equations for the pump operating condition, calculate the dynamic response of speed and torque under the pump operating condition, and form a full operating condition coverage with the corresponding module for the turbine operating condition. PUMP.FOR calls the small opening characteristic data of MPSI5.D.

[0025] The SURGE.FOR subroutine is used to solve the control equation for water level fluctuations in surge tanks, calculate the extreme values ​​and dynamic curves of surge waves in different types of surge tanks such as impedance-type and overflow-type. SURGE.FOR calls the node elevation data of MPSI4.D. 3. Post-processing module: including SS program, SPE program (i.e., the program for calculating speed decay characteristics and quantification index), LAT program (i.e., the program for summarizing extreme parameters under large fluctuation conditions), and DRAW program (i.e., the program for generating visualization curves). The post-processing module reads the result file output by the main calculation module. The SS program is used for stability analysis under small fluctuation conditions to determine whether the system meets the small disturbance stability condition. The SPE program is used to calculate quantitative indicators such as speed transition time, speed fluctuation attenuation coefficient, and pressure pulsation frequency. The LAT program is used to summarize the extreme parameters of large fluctuation conditions, including the maximum or minimum pressure and corresponding time of each pipeline node, the highest or minimum speed of the unit and corresponding time, and the highest or minimum surge value and corresponding time of the surge tank. TECPLOT is a professional drawing software widely used in the fields of water conservancy and hydropower and numerical simulation. Its data format is an industry standard structured text format, which contains core information such as coordinates and parameter values. The DRAW program is used to generate drawing data in TECPLOT format and output visual curves such as pressure-time, flow-time, speed-time, and opening-time.

[0026] III. Calculation Steps The calculation steps for the hydraulic transient process of a pumped storage power station using the above numerical calculation model and system are as follows, with each step executed in logical order: 1. Data Preparation and Preprocessing Step S11: Organize the experimental data on the runner characteristics, including the unit rotational speed N under different guide vane openings. 11 Unit flow rate Q 11 Unit torque M 11S11: Enter the data into the system according to the prescribed format to generate TRNI1.D, TRNI2.D, and NQMI.D files; S12: Enter the pipeline system parameters, unit parameters, governor parameters, and operating conditions to generate READYI.D and MPSI3.D files; Pipeline system parameters include pipe section diameter, length, wall thickness, etc.; Unit parameters include single unit capacity, rated speed, moment of inertia, etc.; Governor parameters include droop rate, time constant, etc.; Operating conditions include calculated operating condition type, initial opening / flow rate, etc.; S13: Run the data formatting program, i.e., the MMM program, to format the files generated in steps S11 and S12 into standardized data that the system can recognize; S14: Run the intelligent error detection program, i.e., the TYJ program, to verify the completeness, logic, and validity of the standardized data. If there are problems such as missing data, logical contradictions, or unreasonable parameters, output prompt information and return to S11 or S12 for correction; If the verification passes, proceed to the next step 2.

[0027] 2. Rotor Characteristic Curve Conversion Steps: S21: Run the TRN program, read the TRNI1.D and TRNI2.D files, and use the Suter method to convert the rotor experimental data into WH (x,y) and WB (x,y) standardized data, generating the MPSI1.D file; S22: Run the NQM program, read the NQMI.D file, and use the length method to convert the rotor characteristic curve into Lq-N... 11 Lq-Q 11 Lq-M 11 Associate the curve data and generate the MPSI5.D file; S23: The system automatically determines the guide vane opening range. For large openings (e.g., opening greater than 0.1), the MPSI1.D file data is used first. For small openings (e.g., opening less than or equal to 0.1) and 0 openings, the MPSI5.D file data is used first, ensuring smooth connection of the characteristic curves across the entire opening range.

[0028] 3. Pipeline system parameter calculation steps: S31: Run the READY program and read the basic data such as pipeline system parameters, elastic modulus, and water density from the READYI.D file; S32: Consider the coupling effect of the elastic modulus of pipeline material, surrounding rock, and water body, and calculate the water hammer wave velocity of the pipe section; S33: Use an adaptive segmentation algorithm to segment the pipe section, divide the calculation nodes, and determine the elevation of each node; S34: Calculate a reasonable time step Δt based on the Courant condition, and generate MPSI2.D (wave velocity and time step file) and MPSI4.D (node ​​elevation file).

[0029] 4. Transient Process Main Calculation Steps: S41: Run the MPS main program and read the MPSI1.D, MPSI5.D, MPSI2.D, MPSI4.D, and MPSI3.D files; S42: Call the governor model, surge tank transient model, and unit boundary model, solve the pipeline transient flow equation based on the method of characteristics, and perform numerical calculations of the transient process; S43: Monitor abnormal operating conditions in real time during the calculation process, including surge tank roof fall / bottom leakage, data anomalies, and calculation non-convergence. If an anomaly occurs, output a warning message and interrupt the calculation. After correcting the parameters, re-execute S41; If the calculation is normal, output MPSO1.D (time series data) and MPSO3.D (extreme parameter data) files, and display key parameters such as time, surge tank water depth, and unit speed in real time.

[0030] 5. Post-processing and result analysis: S51: Run the SS program to analyze the stability of small fluctuation conditions based on the MPSO1.D file and output the stability judgment results; S52: Run the SPE program to calculate quantitative indicators such as speed transition time and attenuation coefficient and output an indicator report; S53: Run the LAT program to extract extreme parameters of large fluctuation conditions from the MPSO3.D file and summarize them into an extreme parameter table of large fluctuation conditions. The extreme parameter table includes maximum / minimum pressure, maximum speed, and surge extreme value;

[0031] S54: Run the DRAW program to generate visualized curve data in TECPLOT format, supporting subsequent plotting and analysis; S55: Combine the stability judgment results of step S51, the quantitative indicators of step S52, the extreme value parameter table of step S53, and the visualized curves of step S54 to complete the comprehensive analysis of the transient process, providing a basis for power plant design or operation optimization.

[0032] IV. Preferred Scheme 1. When converting the runner characteristic curve, for the medium opening range, i.e., 0.3≤guide vane opening≤0.8, a combination of linear interpolation and parabolic interpolation is used to improve the conversion accuracy of the medium opening range; for the small opening range, i.e., guide vane opening≤0.1, length method data is preferred to ensure smooth connection of characteristic curves.

[0033] 2. In the calculation of pipeline system parameters, the segment length Δx is adaptively adjusted according to the pipeline diameter. The larger the diameter, the larger the segment length, but the maximum segment length does not exceed 100m to ensure a balance between calculation accuracy and efficiency. The time step Δt is controlled by the Courant condition to satisfy Δt≤Δx / a, where Δx is the segment length and a is the water hammer wave velocity, to avoid numerical oscillation.

[0034] 3. In the main calculation module, for special working conditions such as pump start-up and guide vane failure to move, a variable step size calculation method is adopted (the initial step size is 1 / 2 of the normal step size, and the normal step size is restored after the calculation stabilizes) to optimize the calculation convergence.

[0035] 4. In the post-processing module, users can customize extreme value parameter filtering conditions, such as filtering the pressure extreme value of a certain pipe section or the speed extreme value within a certain time period, which improves the flexibility of result analysis.

[0036] The beneficial effects of this invention are as follows: 1. Wide range of applications: It supports mixed-flow reversible pumped storage power stations and conventional hydropower stations, is compatible with parallel operation of 1-4 units, various pipeline system layouts and eight types of surge tanks, and covers the calculation of the transition process of turbines, pumps and combinations under all working conditions, thus solving the problem of the limited applicable scenarios of existing programs.

[0037] 2. High calculation accuracy: Adopting a dual algorithm architecture of "Suter method + length method", it overcomes the problem of characteristic curve connection under small opening and 0 opening conditions. The interpolation error in the full opening range is ≤0.5%, and the error between the calculation results and the field test data is ≤3%, which meets the high precision requirements of engineering design.

[0038] 3. High level of intelligence: It integrates automatic data error detection, abnormal working condition prompts, and intelligent result analysis functions, eliminating the need for manual data error checking and automatically outputting stability judgments, quantitative indicators, and extreme value parameters, reducing the reliance on the user's professional knowledge.

[0039] 4. Convenient and efficient operation: Standardized input data, automated calculation process, and intuitive output results. Output results include quantifiable indicators, extreme value tables, and visualization curves. The entire process is integrated into modules, eliminating the need for complex manual intervention and significantly improving engineering design efficiency.

[0040] 5. High engineering applicability: It can be directly applied to power plant design optimization, such as surge tank selection, pipeline diameter determination, equipment selection such as governor parameter matching, and operation safety assessment such as safety verification of sudden load shedding conditions, providing reliable data support for engineering decisions and effectively saving the total investment of the power plant. Attached Figure Description

[0041] Figure 1 is a module architecture diagram of the numerical calculation system of the present invention. Detailed Implementation

[0042] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific examples described herein are only some embodiments of this invention, not all embodiments, and are not intended to limit the invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without inventive effort are within the scope of protection of this invention.

[0043] In this invention: in WH(x,y), x is the guide vane opening α; y is the unit rotational speed N. 11WH(x,y) corresponds to a guide vane opening α(x) and a unit rotational speed N. 11 (y) represents the unit energy parameter of the impeller, reflecting the energy conversion efficiency and head characteristics of water flowing through the impeller.

[0044] The independent variables of WB(x,y) and WH(x,y) are the same: x is the guide vane opening α; y is the unit rotational speed N. 11 Ensure that the parameters of both are matched and can be called synchronously; WB(x,y) corresponds to a certain guide vane opening α(x) and unit rotational speed N. 11 (y) The unit torque parameter of the impeller reflects the torque transmission characteristics between the impeller and the water flow, that is, the torque exerted by the water flow on the impeller, or the torque applied by the impeller to the water flow.

[0045] The MPSI1.D data file contains standardized data obtained using the Suter method. Specifically, it achieves equal-interval processing through parabolic interpolation and linear interpolation, including the runner experimental data such as guide vane opening and unit rotational speed N. 11 Unit flow rate Q 11 Unit torque M 11 Convert the data into standardized data in the forms WH(x,y) and WB(x,y).

[0046] The MPSI2.D data file includes the water hammer wave velocity of the pipe section and the calculation time step.

[0047] The MPSI3.D data file is a standardized data file for operating conditions and unit parameters in the numerical calculation system for the hydraulic transient process of a pumped storage power station. It stores key non-runner, non-pipeline parameters required for calculation, such as unit parameters (single unit capacity, rated speed, moment of inertia, electromagnetic torque, etc.); governor parameters (governor type, permanent droop rate, transient droop rate, buffer time constant, acceleration time constant, etc.); and operating condition parameters (calculation operating condition type, initial guide vane opening, initial flow rate, etc.). The MMM program in the preprocessing module formats the raw data, converting it into a structured format recognizable by the system. It uses a structured text format commonly used in Fortran programs, storing parameters in a fixed row and column order, such as parameter name-value-unit or a pure numerical matrix.

[0048] The MPSI4.D data file includes the elevation of each calculation node and the node number.

[0049] The MPSI5.D file is a standardized data file of impeller characteristics generated based on the length-based NQM program module. The MPSI5.D format uses the length Lq of the flow characteristic curve as the unified horizontal axis and stores Lq-N. 11 Lq-Q 11 Lq-M 11 Standardized data for three sets of correlation curves, where N 11For unit speed, Q 11 Unit flow rate, M 11 The torque is expressed in units.

[0050] The TRNI1.D data file stores the raw experimental data of discrete rotor characteristics, providing an input data source for the TRN program module to perform curve transformation using the Suter method. Its stored parameters include four categories: guide vane opening α, a dimensionless quantity typically ranging from 0 to 1.0, corresponding to fully closed to fully open; unit rotational speed N; and other parameters. 11 The unit is rpm, which represents the standardized rotational speed of the impeller and reflects its rotational characteristics; the unit flow rate is Q. 11 Unit m 3 / s, the standardized flow rate of the impeller, reflecting its flow capacity; unit torque M. 11 The unit is kN·m, representing the standardized torque of the runner, reflecting its torque transmission characteristics; users measure the N of the runner at different guide vane openings through physical experiments. 11 Q 11 M 11 Data; Enter data according to the basic format specified by the system, such as row and column correspondence parameters and units, to generate a TRNI1.D file, which is stored together with TRNI2.D and NQMI.D files of the same type in the intermediate data directory of the preprocessing module.

[0051] The TRNI2.D data file is a supplementary intermediate input file for runner characteristic experimental data. Its function is to supplement and store the raw discrete data of runner characteristic experiments. In conjunction with TRNI1.D, it covers a more comprehensive range of guide vane opening or parameter ranges. It is generated by the user after organizing the runner experimental data and entering it into the system. The parameter types stored in TRNI2.D are completely consistent with those in TRNI1.D. Both are discrete raw data measured by runner physical experiments. The only difference between TRNI2.D and TRNI1.D is the data coverage range: TRNI1.D usually stores the main opening range, such as the core data of large opening range of 0.2~1.0; TRNI2.D supplements and covers other key ranges, such as medium opening data of 0.1~0.3, or partial data split when the experimental data volume is large. After the two are combined, the data coverage of the guide vane large opening range α>0.1 is achieved.

[0052] The NQMI.D data file is a dedicated intermediate input file for turbine characteristic test data. It provides raw data for the NQM program module (i.e., the length method) and serves as the data source for converting turbine characteristic curves under small opening and 0 opening conditions. Its function is to store raw turbine test data adapted for small opening conditions, providing input for the NQM program to convert characteristic curves using the length method, thus solving the problem of broken connections in conventional curves under small opening conditions. It stores parameters of the same type as TRNI1.D and TRNI2.D, all being discrete raw data from turbine physical test measurements. Users select discrete data within the small opening range α≤0.1 from the turbine physical test data; data is entered according to the system's basic format, such as parameter row and column correspondence and unified units, ensuring the integrity of key small opening node data; after being formatted by the MMM program in the preprocessing module, the NQMI.D file is generated and stored together with TRNI1.D, TRNI2.D, etc., as intermediate data. NQMI.D is the data source file for MPSI5.D. After reading the small-aperture raw data from NQMI.D, the NQM program converts it into Lq-N using the length method. 11 Lq-Q 11 Lq-M 11 The associated curve data is used to generate the MPSI5.D file for the main calculation module to use.

[0053] TRNI1.D / TRNI2.D are provided for the TRN program (Suter method) to adapt to large opening conditions, and NQMI.D is provided for the NQM program (length method) to adapt to small opening conditions, together covering the data requirements of the full opening range.

[0054] The READYI.D data file is the original input intermediate file for the piping system and related basic parameters. Its function is to store basic physical parameters of the piping system, surrounding rock, and water body. It is generated by the system after the user compiles the engineering design parameters and enters them into the system. After processing by the READY program, it outputs standardized MPSI2.D files including wave velocity and time step, and MPSI4.D files including node elevations, for use by the main calculation module. The user extracts original parameters such as piping system, material properties, and structure dimensions from the power plant engineering design data; enters the data according to the system's specified format, ensuring consistency in parameter name, value, and unit; and generates the READYI.D file after formatting by the MMM program in the preprocessing module.

[0055] The MPSO1.D data file includes time-series data, containing parameters such as pressure, flow rate, rotational speed, and opening degree at various times.

[0056] The MPSO3.D data file is a summary data file of extreme value parameters for the transient process. It extracts extreme values ​​and their corresponding occurrence times from time series data. MPSO3.D is the "extreme value extraction file" of MPSO1.D. MPSO1.D records the dynamic data throughout the process, such as parameter values ​​every 0.08 seconds, while MPSO3.D only extracts the extreme values.

[0057] The core of the transient process calculation in this invention is to output three types of indicators: hydraulic, mechanical, and stability. The purpose is to provide accurate data support for the design optimization, equipment selection, and safe operation of pumped storage power stations, and to solve the stability and safety problems under complex operating conditions.

[0058] Hydraulic system indicators reflect the dynamic characteristics of water flow. These indicators include pressure, surge, and flow rate indicators. Pressure indicators include the maximum / minimum pressure and occurrence time of each pipeline node, and the maximum vacuum degree of the tailrace pipe. Surge indicators include the highest / lowest surge values ​​and occurrence time of the surge tank. Flow rate indicators include the time series changes in the flow rate within the pipeline, such as the flow rate decay process after a sudden load shedding.

[0059] Mechanical system indicators reflect the unit's operating status. These indicators include speed-related, opening-related, and torque-related indicators. Speed-related indicators include the unit's highest / lowest speeds and their occurrence times, as well as speed fluctuation curves. Opening-related indicators include the dynamic changes in guide vane opening, such as the opening timing data during the broken-line closure process. Torque indicators include the unit torque (M) of the runner. 11 The dynamic changes are adapted to the torque transmission characteristics under different working conditions.

[0060] Stability quantification indicators are used to assess the safety of system operation. Stability quantification indicators include attenuation characteristics, fluctuation parameters, and stability criteria. Attenuation characteristics include speed fluctuation attenuation coefficient and transition time; fluctuation parameters include pressure pulsation frequency and speed oscillation frequency; stability criteria include small fluctuation stability judgment results, including stable / unstable / critically stable.

[0061] The specific implementation of this invention is as follows: I. Hardware and software environment Hardware environment: Processor Pentium II 266 or above; memory ≥ 64M; hard disk space 200M-400M.

[0062] Software environment: The development platform is FortranPowerStation 4.0 or above; the plotting software is Tecplot 7.0 (or above).

[0063] II. Parameter Settings for an Example A pumped storage power station with four generating units is used as an example. The specific parameters are as follows: 1. Basic parameters of the power station: It adopts an impedance-type surge tank, with a design head of 90m, a rated speed of 300rpm, a single unit capacity of 300MW, and a unit rotational inertia of 1200t・m. 2 .

[0064] 2. Pipeline system parameters: Water diversion tunnel diameter 8.5m, length 1200m, wall thickness 0.8m, friction loss coefficient 0.012; tailrace tunnel diameter 9.0m, length 1000m, friction loss coefficient 0.013; pipeline material elastic modulus 206GPa, surrounding rock elastic modulus 30GPa, water density 1000kg / m³ 3 This generates the READYI.D file.

[0065] 3. Speed ​​controller parameters: PI type speed controller, permanent droop rate 5%, buffer time constant 8s, guide vane broken line closing law: first stage closing time 8s, second stage closing time 12s, inflection point opening 50%.

[0066] 4. Calculation condition: The turbine suddenly sheds 100% of its rated load.

[0067] III. Specific Calculation Process Figure 1 shows the modular architecture of the numerical calculation system of this invention, illustrating the relationship between the three major modules—preprocessing, main calculation, and post-processing—and their subroutines. The specific calculation steps of this invention are as follows: 1. Data Preparation and Preprocessing: Organize the runner experimental data (N) in the guide vane opening range of 0.2-1.0. 11 Q 11 M 11 Generate TRNI1.D, TRNI2.D, and NQMI.D files; input the basic parameters of the power station, pipeline system parameters, governor parameters, and operating conditions, and generate READYI.D and MPSI3.D files; run the MMM program to format the data, run the TYJ program to verify, the monotonicity of the surge tank elevation data (480m-520m) is normal, and there are no missing data or abnormal prompts, the preprocessing is complete.

[0068] 2. Rectifier Characteristic Curve Conversion: Run the TRN program to generate the MPSI1.D file, which contains standardized data obtained using the Suter method. Run the NQM program to generate the MPSI5.D file, which converts the rectifier characteristic curves to Lq-N using the length method. 11 Lq-Q 11 Lq-M 11 The system automatically adapts to the opening range, and for small openings, i.e. guide vane openings, it calls MPSI5.D data. The characteristic curves are smoothly connected without any breaks.

[0069] The guide vane opening interval determination involves real-time reading of the calculated guide vane opening value, comparison with a threshold of 0.1, and adaptive switching of data files. Specifically, the guide vane opening α is a dynamic variable in the transient process calculation, not a fixed input value. It iteratively updates with the time step, e.g., calculating the current opening every 0.08 seconds. The opening value is calculated by the governor model based on the unit speed deviation, such as the guide vane action command output by the PI / PID control law, and fed back to the main calculation module through the unit boundary model as the input parameter for interval determination. The guide vane opening interval determination is not an independent process but is embedded in each time step of the main transient process calculation, e.g., Δt=0.08s in the example, and is executed synchronously with the calculation of parameters such as pressure and speed. Within each time step, the main calculation module MPS program first calculates the current guide vane opening α, then performs the interval determination, and finally calls the corresponding data file MPSI1.D or MPSI5.D for subsequent calculations based on the determination result. A dimensionless guide vane opening threshold α0=0.1 is set, corresponding to 10% of the guide vane opening in the project, as the basis for interval division. If the currently calculated guide vane opening α > 0.1: it is determined to be in the large opening range, and the MPSI1.D file generated by the TRN program is called; if the currently calculated guide vane opening α ≤ 0.1 (including 0 opening): it is determined to be in the small opening range, and the MPSI5.D file generated by the NQM program is called; when α fluctuates around 0.1, if there is a small jump caused by calculation error, the system automatically adds hysteresis judgment to avoid calculation oscillation caused by frequent switching of data files. For example, when α drops from 0.105 to 0.1, it switches to MPSI5.D, and only switches back to MPSI1.D when α rises back to above 0.105.

[0070] Equal-spacing processing using parabolic and linear interpolation is a preprocessing operation for turbine runner characteristic experimental data. Its purpose is to transform discrete, non-equal-spacing raw experimental data into uniformly spaced, continuous, standardized data in WH(x,y) and WB(x,y) formats. This ensures that subsequent hydraulic transient process calculations can quickly and accurately retrieve runner characteristic parameters such as unit flow rate and unit torque at any opening / speed, avoiding calculation errors caused by uneven data spacing. Turbine runner experimental data is obtained through physical experiments and is typically non-equally spaced; for example, guide vane opening measurements are only taken at discrete values ​​such as 0.2, 0.4, 0.6, and 1.0, and unit speed measurements are only taken at discrete points such as 250, 300, and 350 rpm. Numerical calculations require uniformly spaced data grids, such as opening intervals of 0.01 from 0 to 1.0 and speed intervals of 5 rpm from 200 to 400 rpm, to solve for the continuous dynamic process using a model. Using linear or parabolic interpolation, data is filled in the gaps between the original discrete data to form a continuous data matrix with equal spacing. This ensures smooth connection of characteristic curves across the full opening and full speed range, with an interpolation error ≤0.5%. Specific processing methods for achieving equal spacing through parabolic and linear interpolation include: Data preprocessing: Defining the original data and target format. The original data consists of experimental measurements, representing four discrete sets of core parameters, as shown in the table below:

[0071] The target data is the result of equal spacing: the guide vane opening α is divided into 101 points from 0 to 1.0 at intervals of 0.01, with a unit rotational speed N. 11 Divide the system into 21 points at intervals of 5 rpm from 250 to 350 rpm, and complete each (α, N) point. 11 The corresponding Q 11 and M 11 This forms equally spaced two-dimensional data matrices, namely WH(x,y) and WB(x,y) files.

[0072] Regarding the specific application of the two interpolation methods in different scenarios, interpolation involves knowing two or three adjacent original data points and using mathematical formulas to calculate the value of a certain intermediate point, ensuring that the result is consistent with the trend of the original data and the curve is smooth.

[0073] (1) Linear interpolation: suitable for intervals with close data intervals and gentle trends. Assume that the characteristic parameter between two adjacent original data points is Q. 11 Depending on the independent variable such as α or N 11 It changes in a straight line; the value of the midpoint can be calculated using the equation of the straight line.

[0074] Using the guide vane opening α as the independent variable, complete Q. 11 For example: Given the original point A(α1, Q) 111) and point B(α2, Q) 112 The Q values ​​of the equidistant points α0 (α1 < α0 < α2) need to be completed. 110 : For example, when α = 0.2, Q 11 When Q = 1.25 and α = 0.4 11 =1.56, Q needs to be completed for equidistant points α=0.3. 11 : ,

[0075] Linear interpolation is suitable for large opening ranges α>0.1 or regions where the rotational speed changes gradually. It is simple to calculate, efficient, and can meet the basic smoothness requirements.

[0076] (2) Parabolic interpolation: It is suitable for intervals with large data intervals and non-linear trends. It assumes that the characteristic parameters between three adjacent original data points change in a parabolic quadratic curve. The value of the intermediate point is calculated by fitting a quadratic function. It is more accurate than linear interpolation and can capture the non-linear trend of the data.

[0077] Using the guide vane opening α as the independent variable, complete Q. 11 For example: Given the original point A(α1, Q) 111 Point B (α2, Q) 112 Point C (α3, Q) 113 (α1<α2<α3), Q needs to be completed at the intermediate point α0. 110 First, construct a quadratic function: Q 11=a α 2 +bα+c, substitute the three original points to solve for the coefficients a, b, and c, then substitute α0 to calculate Q. 110 : After solving for a, b, and c, Q 110 = aα0 2 Given α = 0.2 (Q = 1.25), α = 0.4 (Q = 1.56), and α = 0.6 (Q = 1.82), we need to complete the Q for α = 0.5. 11 Substituting these values ​​into the equation, we get a = -0.125, b = 0.875, and c = 0.95. Therefore, Q 110 = -0.125 x 0.5 2 + 0.875 x 0.5 + 0.95 = 1.6975.

[0078] If linear interpolation is used to calculate the intermediate value of α from 0.4 to 0.6, the result is 1.69. Parabolic interpolation is closer to the actual nonlinear trend.

[0079] Parabolic interpolation is suitable for medium opening ranges of 0.3≤α≤0.8 or regions with large speed variations. It can improve interpolation accuracy and avoid curve bends caused by linear interpolation.

[0080] The complete process for equal spacing processing: Data partitioning: The original impeller experimental data are divided according to the guide vane opening α and the unit rotational speed N. 11 Two-dimensional classification determines the equidistant grid that needs to be filled, such as α: 0 to 1.0, step size 0.01; N 11 : 200 to 400 rpm, in 5 rpm increments.

[0081] Interpolation selection: For small opening intervals α≤0.1: the length method will be used separately later, and interpolation here is only a transition; For medium opening intervals 0.3≤α≤0.8: a combination of linear interpolation and parabolic interpolation is used, with parabolic interpolation used for key nodes and linear interpolation used for ordinary nodes; For large opening intervals α>0.1: linear interpolation is used as the main method to ensure efficiency.

[0082] Data validation: After completing all equally spaced points, check the characteristic curve such as Q. 11 -α curve, M 11 -N 11 The curve should be smooth, without breaks or abrupt changes, and the interpolation error should be ≤0.5%.

[0083] Formatted output: Q after equal spacing 11 M 11 The data is organized into WH(x,y) files (unit energy characteristics) and WB(x,y) files (unit torque characteristics) to generate an MPSI1.D file for use by the main calculation module.

[0084] The purpose of interpolation in this invention is to fill in the gaps in the data, rather than to modify the original experimental data. Ultimately, it is necessary to ensure that the interpolation result is consistent with the trend of the original data, and the error is controlled within the engineering allowable range, i.e., less than or equal to 0.5%. Interpolation is the core function of the turbine characteristic curve conversion model, i.e., the TRN program module. Subsequently, combined with the length method, i.e., the NQM program module, the curve connection problem of small opening α≤0.1 and 0 opening conditions can be further solved.

[0085] The MPSI5.D file of this invention is a standardized data file of turbine runner characteristics generated by the NQM program module based on the length method. It is specifically adapted for hydraulic transient process calculations in pumped storage power stations under small-opening conditions (guide vane opening α ≤ 0.1 and 0 opening). It serves as the data carrier connecting the characteristic curves in the small-opening interval within the numerical calculation model. It complements the MPSI1.D file generated by the Suter method for large-opening conditions (α > 0.1), jointly ensuring a smooth connection of the entire turbine runner characteristic curve. The MPSI5.D file is generated by the NQM program module, with input data including turbine runner experimental data, guide vane opening α, and unit rotational speed N. 11 Unit flow rate Q 11 Unit torque M 11The MPSI5.D format stores Lq-N with the length Lq of the flow characteristic curve as the uniform horizontal axis. 11 Lq-Q 11 Lq-M 11 Standardized data for three sets of correlation curves.

[0086] Lq-N 11 Lq-Q 11 Lq-M 11 The generation of the three sets of correlation curves is explained below:

[0087] The typical characteristic curve is α-N 11 α-Q 11 α-M 11 That is, the guide vane opening α is used as the reference, but when the opening α ≤ 0.1, α changes sensitively and the curve is prone to breakage; while Lq and Q 11 Directly related, Q 11 The larger the value, the longer Lq becomes. Within the small opening range, Lq changes more smoothly. Therefore, Lq is used instead of α as a unified reference, transforming the three independent curves into a curve with Lq as the abscissa and N as the ordinate. 11 / Q 11 / M 11 To avoid breaks in the connection between small openings, the specific transformation method for the correlation curve on the ordinate consists of four steps: first, calculate Lq; then, establish the relationship between Lq and N. 11 / Q 11 / M 11 The mapping relationship is implemented entirely through the NQM program module. The steps are as follows: Step 1: Organize the original experimental data as the input basis. Extract Q under different guide vane openings α from the runner experimental data. 11 N 11 M 11 The resulting discrete data sets are as follows, with the following being the key data for the small opening interval:

[0088] Step 2: Calculate the length Lq of the flow characteristic curve. Starting from 0 opening (α=0), calculate Lq, i.e., Q, corresponding to each α. 11 -The cumulative length of the -α curve from α=0 to the current α.

[0089] Since the experimental data are discrete, a piecewise piecewise linear approximation method is used to calculate Lq, as shown in the following formula: For the i-th opening degree Its corresponding cumulative length L qi = Cumulative length L of the previous opening qi-1 + The curve length ΔL from segment i-1 to segment i qi-1→i ; Single segment length ΔL qi-1→i The calculation is as follows: Single segment length ΔL qi-1→i Essentially, it's about the length of the hypotenuse of a right triangle: the difference in the horizontal coordinates is α, and the difference in the vertical coordinates is Q. 11 Change.

[0090] The calculation process based on the data from step 1 is as follows: Sequence 1 (α=0): Lq1=0, starting point, cumulative length is 0; Sequence 2 (α=0.02): ΔLq1→2= ≈0.351m, then Lq2=0+0.351=0.351m; Serial number 3 (α=0.05): ΔLq2→3= ≈ 0.431m, then Lq3 = 0.351 + 0.431 = 0.782m; Serial number 4 (α = 0.08): ΔLq3 → 4 = ≈0.421m, then Lq4 = 0.782 + 0.421 = 1.203m; Serial number 5 (α = 0.10): ΔLq4 → 5 = ≈ 0.361m, then Lq5 = 1.203 + 0.361 = 1.564m; finally, the α-Lq correspondence is obtained and added to the original data:

[0091] Step 3: Establish Lq and N 11 / Q 11 / M 11 The mapping relationship replaces the x-coordinate α of the original data with Lq, forming three sets of discrete mapping relationships: Lq-N 11 Mapping: (Lq1,N) 111 (Lq2,N) 112 ), ..., (Lqᵢ,N 11 ᵢ);Lq-Q 11 Mapping: (Lq1,Q) 111 (Lq2,Q) 112 ), ..., (Lqᵢ,Q 11 ᵢ); Lq-M 11 Mapping: (Lq1,M) 111 (Lq2,M) 112 ), ..., (Lqᵢ,M 11 ᵢ); Step 4: Equal-interval interpolation to complete and generate correlation curves. To meet the continuous calling requirements of numerical calculation, equal-interval interpolation is performed on the three sets of mapping relationships to complete all data points in the Lq interval, and finally a smooth correlation curve is formed: Determine the equal-interval range of Lq: Take the maximum Lq calculated in step 2, such as 1.564m in the example, as the upper limit, and divide it according to a fixed step size, such as 0.01m, to obtain equal-interval Lq values: 0, 0.01, 0.02, ..., 1.56, 1.564m.

[0092] Interpolation parameter completion: The trend is gentle in the small opening range, and linear interpolation is accurate enough. Linear interpolation is used for each equally spaced Lq value to complete the corresponding N. 11 Q 11 M 11 .

[0093] To complete the parameters when Lq=0.5m: find the interval where Lq=0.5m is located: between Lq=0.351m in step 1 (number 2) and Lq=0.782m in step 1 (number 3).

[0094] Interpolation calculation N 11 :N 11 =220+(0.5 - 0.351) / (0.782 - 0.351)×(250 - 220)≈220+0.346×30≈230.4rpm. Similarly, interpolate Q. 11 M 11 This yields a complete set of equally spaced Lq-parameter data.

[0095] Generate correlation curves: Plot the completed data into three continuous curves, Lq-N 11 Curve, Lq-Q 11 Curve, Lq-M 11 The curve is ultimately saved as an MPSI5.D file.

[0096] The length method, i.e., the connection between NQM and Suter methods: When the guide vane opening α > 0.1, the system switches to the WH(x,y) / WB(x,y) data generated by the Suter method, i.e., the MPSI1.D file; when α ≤ 0.1, it switches to the Lq correlation curve data, i.e., the MPSI5.D file, ensuring seamless connection of characteristic curves across the full opening range. Lq-N 11 Curve, Lq-Q 11 Curve, Lq-M 11 The curve transformation process is completed automatically by the NQM program. Users do not need to manually calculate Lq; they only need to input the original experimental data, which greatly reduces the complexity of operation.

[0097] 3. Pipeline system parameter calculation: Run the READY program, considering the coupling effect of the elastic modulus of the pipeline, surrounding rock, and water body, and calculate the water hammer wave velocity as 1200m / s; divide the water diversion tunnel into 12 segments of 100m each and the tailrace tunnel into 10 segments of 100m each according to the adaptive segmentation algorithm; determine the calculation time step Δt=0.08s based on the Courant condition; generate the MPSI2.D file containing the water hammer wave velocity of 1200m / s and the step size of 0.08s, and the MPSI4.D file containing the elevation data of each node.

[0098] Regarding the calculation of water hammer wave velocity: If the pipeline is not constrained by surrounding rock, such as an overhead pipeline, the formula for water hammer wave velocity is: Where: a0 represents the water hammer wave velocity without surrounding rock constraint, in m / s; K represents the bulk elastic modulus of water (Pa), which is approximately 2.1 × 10⁻⁶ for water at room temperature. 9 Pa is a fixed constant; ρ represents the density of water, which is 1000 kg / m³ at room temperature. 3 , a fixed constant; D represents the equivalent pipe diameter, i.e., the inner diameter of the pipe. In this embodiment of the invention, the water diversion tunnel D=8.5m and the tailrace tunnel D=9.0m; E represents the elastic modulus (Pa) of the pipe material. In this embodiment of the invention, the pipe material is steel, E=206GPa=206×10 9 Pa; e represents the pipe wall thickness. In this embodiment of the invention, the water diversion tunnel has an e value of 0.8m.

[0099] When considering the constraints of surrounding rock, the formula for water hammer wave velocity needs to be modified: In actual pumped storage power stations, pipelines are mostly underground tunnels, such as diversion / tailrace tunnels. The constraint effect of surrounding rock on the pipeline needs to be considered. The surrounding rock restricts the radial deformation of the pipeline, increasing the wave velocity. Therefore, a surrounding rock constraint correction factor C needs to be introduced. The final practical engineering formula is: C represents the surrounding rock constraint coefficient, a dimensionless quantity that depends on the elastic modulus Erock of the surrounding rock and the contact conditions between the pipeline and the surrounding rock. The formula is: E rock This represents the elastic modulus of the surrounding rock, expressed in Pa.

[0100] The calculated water hammer wave velocity 'a' will be directly used to: determine the calculation time step Δt, which must satisfy the Courant condition Δt ≤ Δx / a. The Courant condition, also known as the CFL condition (Courant-Friedrichs-Lewy condition), is a core criterion for judging the stability of solving the transient flow equation in a pipe during numerical calculations. Its function is to avoid numerical oscillations and ensure that the calculation results converge and closely match the actual physical process. The Courant condition requires that, in numerical calculations, the distance the water hammer wave travels within one time step Δt cannot exceed the segment length Δx of the pipe section. Δx is the segment length of the pipe section. In this embodiment of the invention, Δx = 100m, and Δt = 0.08s ≤ 100 / 1200 ≈ 0.083s, which meets the requirements; solve the transient flow equation of the pipe. The transient flow equation of the pipe is a core parameter of the method of characteristics and affects the accuracy of pressure fluctuation calculations; generate the MPSI2.D data file for the main calculation module to call.

[0101] The following explains the MPSI4.D node elevation data file generated by this invention: MPSI4.D is a pipeline system node elevation data file generated by the READY program module. It automatically outputs standardized node elevation data through the process of inputting basic parameters → pipe segmentation and node division → elevation calculation and verification. The MPSI4.D file contains structured data of calculation node numbers and corresponding elevation values, covering the segmentation nodes of all pipe segments in the water intake system and tailrace system, as well as key nodes such as surge tanks and unit interfaces. It is used to determine the gravity head of each node in subsequent transition process calculations, and the gravity head affects pressure calculation.

[0102] The complete steps for generating MPSI4.D are as follows: Step 1: Prepare input data. This is the input basis for the READY program. The basic parameters of the pipeline system need to be organized, entered into the specified format, and a `READYI.D file is generated. This is the intermediate file output by the preprocessing module. Input parameters include: Pipeline segment geometric parameters: the starting elevation and ending elevation of the water diversion tunnel and tailrace tunnel, and the length of the pipe segment. For example, if the length of the water diversion tunnel is 1200m, the starting elevation is assumed to be 520m and the ending elevation is assumed to be 480m; if the length of the tailrace tunnel is 1000m, the starting elevation is assumed to be 480m and the ending elevation is assumed to be 450m; Pipeline segment layout parameters: pipeline slope or horizontal / vertical layout, elbow / branch position; Key structure parameters: bottom elevation and top elevation of the surge tank, unit spiral casing inlet elevation, and tailrace pipe outlet elevation; Calculation control parameters: The maximum segment length of the pipe is limited to ≤100m in this embodiment of the invention, which is determined by the preferred scheme.

[0103] Step 2: Run the READY program. After starting the pipe segmentation and node division READY program and reading the READYI.D file, process it according to the following logic: Pipe segmentation: Adaptive segmentation algorithm is used to divide the pipe into segments according to the pipe length and the maximum segmentation limit, ensuring that the length of each segment is ≤100m. For example, if the length of the water diversion tunnel is 1200m, it is divided into 12 segments at 100m / segment; if the length of the tailrace tunnel is 1000m, it is divided into 10 segments at 100m / segment.

[0104] Calculation node generation: Calculation nodes are set at the beginning and end of each pipe section. At the same time, key locations such as the bottom of the surge tank and the unit interface are set as special nodes and uniformly numbered from upstream to downstream according to the water flow direction, such as nodes 1~12 for the water diversion tunnel, node 13 for the surge tank, node 14 for the unit, and nodes 15~24 for the tailrace tunnel.

[0105] Step 3: Calculate the elevation of each node. The READY program calculates the elevation of each node based on the starting-end elevation and segment length of the pipe section, using linear interpolation or the actual slope of the pipe, as follows: If the pipe has the most common linear slope in the project, the node elevation calculation formula is as follows: , where: H i H represents the elevation (m) of the i-th node.start H represents the starting elevation of the pipe section (m). end Indicates the elevation (m) of the end point of the pipe section; L total L represents the total length of the pipe section (m). cumulative This represents the cumulative length (m) from the starting point of the pipe segment to the i-th node.

[0106] For example, in this embodiment of the invention, the water diversion tunnel starts at an elevation of 520m, ends at an elevation of 480m, and has a total length of 1200m, with each segment being 100m long; the starting elevation of node 1 is 520m; the elevation of node 2 is 520 - (520-480) / 1200 × 100 = 516.67m; the end point of node 12 of the water diversion tunnel, the bottom elevation of the surge tank is 480m; the elevation of node 13 of the surge tank is 480m, level with the bottom.

[0107] Special node handling: The top elevation of the surge tank and the unit interface elevation are directly taken from the design values ​​in the input parameters and are not involved in the interpolation calculation.

[0108] Step 4: Data Verification and Formatted Output Elevation Verification: The READY program automatically verifies the monotonicity of the elevation data, such as the gradual decrease in the elevation of the water diversion tunnel nodes from upstream to downstream without reverse fluctuations. If abrupt elevation changes occur, such as a jump in node elevation, the output parameters will show an abnormality and the program needs to return to correct the starting / ending elevation of the pipe section. Formatted File Generation: After the verification passes, the program generates a structured text file MPSI4.D in a fixed format of node number plus elevation value.

[0109] When performing transient process calculations, the MPS main program reads the node elevation data from MPSI4.D and, in conjunction with parameters such as water hammer velocity and flow rate, solves for the pressure at each node using the method of characteristics. Pressure = head + gravity head corresponding to the elevation.

[0110] 4. Transient process main calculation: Run the MPS main program, call each subroutine to read the above data files, and perform transient process calculation. The calculation time is 45 seconds. There are no abnormal operating condition prompts during the calculation. Output MPSO1.D time series data and MPSO3.D extreme value parameter data files.

[0111] The core computational logic of the MPS main program is based on solving the transient flow equations of the pipeline using the method of characteristics. It couples the governor, surge tank, and unit boundary multiphysics model. Through a process of data reading, model initialization, step-by-step iterative calculation, and result output, it ultimately obtains the time series data MPSO1.D and the extreme parameter data MPSO3.D. The computation time of 45 seconds is derived from iterative accumulation based on the complexity of the operating conditions and the time step. The following is a detailed description of the computational method and result generation process: Preparation before computation: Data reading and initialization. After the MPS main program starts, it first completes data loading and parameter initialization, laying the foundation for the computation.

[0112] Read all input data files: Runner characteristic data: MPSI1.D uses the Suter method and large opening; MPSI5.D uses the length method and small opening, including WH(x,y), WB(x,y), and Lq correlation curve data; Piping system data: MPSI2.D includes water hammer wave velocity of 1200m / s and time step Δt=0.08s; MPSI4.D includes calculation node elevation and number; Operating conditions and unit parameters: MPSI3.D includes PI type governor parameters, sudden 100% load shedding condition, and unit rotational inertia of 1200t·m. 2 wait.

[0113] Automatic adaptation calculation includes: operating condition identification: confirming the current turbine operating condition - sudden shedding of 100% rated load, and calling the corresponding turbine boundary equation; opening range judgment: if the initial guide vane opening is >0.1, MPSI1.D data is called first; if the guide vane is closed to ≤0.1 during the calculation, it automatically switches to MPSI5.D data; surge tank model adaptation: loading the water level fluctuation control equation of the impedance-type surge tank.

[0114] Initial steady-state parameters: Hydraulic system: steady-state pressure at each node of the pipeline, corresponding to the pressure at a design head of 90m; steady-state flow rate, corresponding to the rated load; steady-state water level of the surge tank, set to the range of 480m-520m; Mechanical system: steady-state speed of the unit is 300rpm; steady-state opening of the guide vanes is 50% of the inflection point opening; Control parameters: initial output of the governor, control signal adapted to the steady-state opening.

[0115] The core calculation method for the transient process employs the method of characteristics combined with multi-model coupling iteration. It calculates the transient process in step-by-step iterations with a time step of Δt = 0.08 s. Each step couples the dynamic changes in pipeline flow, unit speed, and governor operation. The specific steps are as follows: Solving the transient flow in the pipeline using the method of characteristics. The transient flow in the pipeline, including pressure and flow rate changes, is the core of the transient process. The method of characteristics is used to transform the partial differential equations into ordinary differential equations for easier numerical solution. Equation discretization: The water diversion tunnel is divided into 12 segments, each 100 m long, and the tailrace tunnel is divided into 10 segments, each 100 m long. The transient flow control equations for the water diversion tunnel and tailrace tunnel, including the continuity equation and momentum equation, are discretized into two characteristic line equations for forward propagation and reverse propagation using the method of characteristics. Where H is the node head, Q is the flow rate, a is the water hammer wave velocity (1200 m / s in this embodiment), g is the gravitational acceleration, A is the cross-sectional area of ​​the pipe, f is the friction loss coefficient (0.012 for the water diversion tunnel), and D is the pipe diameter.

[0116] Node parameter solution: For each upstream and downstream node of a pipe segment, at a time step Δt=0.08s, substitute the pressure and flow data from the previous moment to solve for the node pressure H and flow rate Q at the current moment. For example: In the first time step t=0.08s: based on the steady-state pressure 90m and flow rate, calculate the initial pressure fluctuation and flow rate decay in the pipeline after a sudden load shear; in each subsequent step: using the calculation results from the previous step as input, iteratively update the pressure and flow rate of all nodes to ensure the continuity of the transient process of water flow.

[0117] Unit dynamics calculation describes the dynamic change of rotational speed: combining pipeline flow changes, it solves for the dynamic response of the unit's rotational speed. The core of the calculation is the torque balance equation, which is the balance relationship between the unit's moment of inertia J, electromagnetic torque Me, and turbine torque Mt, reflecting the rate of change of rotational speed dn / dt. Turbine torque Mt: WB(x,y) or Lq-M from runner characteristic data MPSI1.D / MPSI5.D 11 In the middle, based on the current guide vane opening and unit rotational speed N 11 The query results are then standardized to the actual torque; Electromagnetic torque Me: After a sudden 100% load shedding, Me quickly drops to 0 and the load is completely cut off, and subsequently fluctuates slightly with changes in rotational speed.

[0118] Iterative calculation of rotational speed: Within each time step Δt, the change in rotational speed Δn is calculated using the above equation, resulting in the current rotational speed n = rotational speed at the previous time step + Δn. For example: 8 seconds after load shedding, i.e., the 100th time step, 8 / 0.08 = 100, the rotational speed rises to a maximum of 375 rpm, 1.25 times the rated speed, and then gradually decreases.

[0119] The governor control logic calculation describes the guide vane action response. The governor automatically adjusts the guide vane opening based on the unit speed deviation to achieve closed-loop speed control. The core is the PI regulation law: PI regulation equation: The change in guide vane opening Δα is jointly determined by the proportional element (P) and the integral element (I). Where Δα is the guide vane opening deviation, which is the current opening minus the steady-state opening; Kp is the proportional gain; Ki is the integral gain, which are derived from the permanent droop rate of 5% and the buffer time constant of 8s; and n0 is the rated speed of 300rpm.

[0120] Guide vane closing pattern execution: The guide vane closes according to a preset pattern, that is, within the first segment of 8 seconds, it closes from the initial opening to 50%, with 50% being the inflection point opening of the guide vane. Within the second segment of 12 seconds, it closes from 50% to the minimum opening. At each time step, the current opening α is calculated as the previous opening plus Δα to ensure a smooth closing process.

[0121] Transient calculation of surge tanks describes water level fluctuations. Based on the control equations of impedance-type surge tanks, it calculates the water level change at each step and reflects the system's pressure buffering effect.

[0122] The transient calculation equation for well control, also known as the control equation for water level fluctuations in surge tanks, states that the rate of change of water level h in the surge tank is related to the flow difference between the inflow and outflow from the surge tank. , where A t For the cross-sectional area of ​​the surge tank, Q in The flow rate Q from the water diversion tunnel into the surge tank out This refers to the flow rate from the surge tank to the generating unit.

[0123] The relationship between the bottom pressure and the water level is: P = gh + P0, where ρ is the water density, g is the acceleration due to gravity, h is the water level, and P0 is the atmospheric pressure.

[0124] Water level iteration: At each time step, Q is calculated based on the change in pipeline flow rate. in and Q out The difference is used to obtain the water level change Δh, and the current water level h is updated. For example: after load shedding, the pipeline flow rate decreases sharply, the surge tank water level rises, and the maximum surge wave is 3.2m; then the water level falls back, and the minimum surge wave is -2.8m.

[0125] The relative vacuum degree in this invention The calculation method is as follows: In the formula: P0 is atmospheric pressure, which is often taken as the height of a water column corresponding to 10m water column or standard atmospheric pressure in engineering; ρ is water density; g is gravitational acceleration; H 节点压力 The piezometric head of the tailrace pipe node is obtained by the method of characteristics, and z is the elevation of the node, which comes from the MPSI4.D file.

[0126] After calculation, the vacuum data is synchronously stored in the MPSO1.D time series file and correlated with parameters such as pressure and flow rate.

[0127] Real-time monitoring of abnormal operating conditions includes: synchronously monitoring three types of abnormalities during the calculation process to ensure calculation stability: surge tank roof fall / bottom leakage: monitoring whether the water level exceeds the design upper limit such as 520m or falls below the lower limit such as 480m; data validity: checking whether the impeller characteristic data is sufficient and whether the pipeline parameters are reasonable and without abnormalities; calculation convergence: performing convergence verification on the pressure, speed and water level calculation results of each step to avoid numerical oscillation.

[0128] The calculation time of 45s in this invention is the cumulative result of time step multiplied by the number of iteration steps, not a preset value: Time step Δt = 0.08s, determined by the Courant condition Δt ≤ Δx / a: Δx = 100m, a = 1200m / s, Δt ≤ 0.083s, taking 0.08s ensures stability; Number of iteration steps: 45s ÷ 0.08s / step ≈ 562.5 steps, actually taking 563 steps; Termination condition: When the unit speed fluctuation drops to within 5% of the rated speed, approximately 315rpm~285rpm, and the pressure and water level tend to stabilize, the calculation automatically terminates, and this invention satisfies this condition in 45s.

[0129] Result File Generation: MPSO1.D and MPSO3.D data files are generated. The specific contents of MPSO1.D and MPSO3.D data files are as follows: MPSO1.D is time series data: It records the entire process dynamically, storing the core parameters at each time step to form a continuous dynamic data sequence, facilitating subsequent analysis. Time column: 0s, 0.08s, 0.16s, ..., 45s, a total of 563 time points; Parameter column includes: Hydraulic parameters: pressure at each pipeline node, such as the maximum pressure of the water intake pipeline (117m), flow rate; surge tank level, such as the maximum surge (3.2m), minimum (-2.8m); tailrace vacuum degree, such as -4.2m; Mechanical parameters: unit speed, such as 375rpm at 8s, guide vane opening, such as 50% at 8s, minimum opening at 20s; MPSO1.D format: structured text format, each line corresponds to one time point, and each column corresponds to one parameter, which can be directly read by SS, SPE, and DRAW programs.

[0130] MPSO3.D contains extreme value parameter data: key extreme values ​​are extracted. From the time series of MPSO1.D, the maximum / minimum values ​​and corresponding occurrence times of each parameter are filtered and summarized into an extreme value table for engineering use: extreme values ​​of rotational speed: maximum 375 rpm, 8s; minimum stable rotational speed, such as 290 rpm, 45s; extreme values ​​of pressure: maximum 117m in the water intake pipe, 10s; minimum pressure, such as 85m, 30s; maximum vacuum in the tailrace pipe, -4.2m, 12s; extreme values ​​of surge: maximum 3.2m in the surge tank, 9s; minimum -2.8m, 15s; MPSO3.D format: tabular text containing three columns: parameter name, extreme value, and occurrence time, which can be directly used for engineering evaluation, such as determining whether the rotational speed or pressure exceeds the limit.

[0131] In this invention, multiple models are coupled: the pipeline transient flow characteristic line method, the unit dynamic torque balance, the governor control PI regulation, and the surge tank transient water level equation are simultaneously iterated. Each step is interconnected; for example, flow rate changes affect rotational speed, rotational speed changes trigger guide vane action, and guide vane action changes flow rate. This invention adopts step-by-step iteration, with 0.08s as the minimum time unit, to progressively advance the calculation, ensuring the continuity and accuracy of the dynamic process. The results of this invention are output in layers: MPSO1.D records the entire dynamic process, and MPSO3.D extracts the core extreme values, which not only meets the needs of detailed analysis but also facilitates rapid engineering application. The entire calculation process does not require manual intervention; the MPS main program automatically calls subroutines and iteratively solves the problem, outputting the results after 45s. The error between the results and the field test data is ≤3%, ensuring the reliability of the results.

[0132] 5. Post-processing and result analysis: Running the SS program, the system's stability under small fluctuations was determined to be good; running the SPE program, the calculated speed transition time was 45s, and the speed fluctuation attenuation coefficient was 0.05; running the LAT program, the extreme parameters were summarized: the unit's highest speed was 375rpm, which was 1.25 times the rated speed, and it occurred 8s after load shedding; the maximum surge in the surge tank was 3.2m, and the minimum surge was -2.8m; the maximum pressure in the water intake pipeline was 117m, which was 1.3 times the design head, and the maximum vacuum in the tailrace pipe was -4.2m; running the DRAW program, the speed-time, pressure-time, and surge tank water level-time visualization curves were generated.

[0133] The table below shows the correspondence between the functions and calculation methods of each program in the post-processing module:

[0134] This invention utilizes the SS program to determine the good stability of a system under small disturbances. The SS program, or Small Disturbance Stability Analysis Program, is a specialized tool that uses linearization analysis combined with eigenvalue analysis to determine whether a system can recover stable operation under small disturbances. Its judgment method and standards follow the stability theory of hydraulic-mechanical coupled systems and are closely integrated with actual engineering needs, as detailed below: SS Program Judgment Method: Small Disturbance Linearization and Eigenvalue Analysis. The essence of small disturbance stability is whether the system, after being subjected to small disturbances such as small guide vane oscillations, small load fluctuations, or small flow pulsations at its rated operating conditions or steady-state operating point, can automatically attenuate the disturbance and return to steady state. The SS program's judgment method is based on linearization modeling and eigenvalue solving. The steps are as follows: Step 1: Determine the steady-state operating point as the calculation benchmark. The SS program first extracts the steady-state operating parameters of the system from the MPSO1.D time series data output by the main calculation module, i.e., the steady state before the occurrence of small disturbances, as the benchmark point for linearization analysis. The parameters include: Hydraulic system: steady-state pressure and steady-state flow rate of each pipeline node; steady-state water level of the surge tank; Mechanical system: steady-state speed of the unit, such as the rated speed of 300 rpm; steady-state opening of the guide vanes; Control parameters: steady-state output of the governor, such as the steady-state gain of the PI / PID controller.

[0135] Step 2: Establish small-disturbance linearized equations. The nonlinear coupled equations of the system—the transient flow equation, the unit dynamics equation, and the governor control equation—are expanded using Taylor at the steady-state operating point. Higher-order terms are ignored, transforming them into linearized state-space equations: ,in The variables are the rate of change of state variables, such as the rate of change of rotational speed, the rate of change of pressure, and the rate of change of opening; A represents the system matrix, the dimension of which depends on the number of state variables and reflects the coupling relationship of each parameter; X is the vector of state variables, including: rotational speed deviation Δn, guide vane opening deviation Δα, pipeline pressure deviation ΔP, surge tank water level deviation Δh, etc.

[0136] The elements of system matrix A are coefficients that reflect the coupling relationship between various state variables. Each element corresponds to the correlation coefficient between the rate of change of one state variable and the deviation of another state variable. It is the derivative result of the nonlinear equation of the system after linearization at the steady-state operating point, and embodies the dynamic coupling logic of multiple systems including water flow, mechanics, and control.

[0137] The number of columns / rows in matrix A is determined by the state variables: Dimension of matrix A = Number of state variables × Number of state variables (taking four variables as an example): X1: Unit speed deviation Δn, actual speed - rated speed, unit: rpm; X2: Guide vane opening deviation Δα, actual opening - steady-state opening, dimensionless; X3: Pipeline pressure deviation ΔP, actual pressure - steady-state pressure, unit: mH2O; X4: Surge well water level deviation Δh, actual water level - steady-state water level, unit: m.

[0138] At this point, matrix A is a 4×4 square matrix, in the following form: , element a ij The j-th element in the i-th row represents the j-th state variable X. j For every unit change, the rate of change of the i-th state variable is... The change in value, expressed in units of variables per second: Based on the multi-model coupling relationships in the paper, including unit dynamics, governor control, pipeline transient flow, and surge tank transient flow, the elements of the 4×4 matrix A have the following meanings: the first row corresponds to the rate of change of rotational speed. = n, a 11 The coefficient of influence of rotational speed deviation on its own rate of change. D is the mechanical damping coefficient of the unit, and J is the moment of inertia. The larger the speed deviation, the faster the speed decays, and it is usually a negative number.

[0139] a 12 : The coefficient of influence of opening deviation on the rate of change of rotational speed , J is the partial derivative of the torque with respect to the opening degree, and J is the moment of inertia. As the opening degree of the guide vane increases, the turbine torque increases, and the speed increases. It is usually a positive number.

[0140] a 13 The coefficient of influence of pressure deviation on the rate of change of rotational speed. This stems from the coupling between the turbine torque and the pipeline pressure. The partial derivative of torque with respect to pressure means that as the pressure in the pipeline increases, the torque of the turbine increases, and thus the rotational speed increases; it is usually a positive number.

[0141] a 14 The coefficient of influence of water level deviation on the rate of change of rotational speed originates from the indirect effect of the surge tank water level on the pipeline pressure. A rise in water level leads to a rise in pressure, which in turn increases the torque. The formula is: , This is the partial derivative of the torque with respect to the water level, which is usually a positive number. The water level deviation indirectly affects the rotational speed through the pressure. When the water level rises, the pressure rises, which in turn increases the torque and thus the rotational speed.

[0142] The second row corresponds to the rate of change in opening degree. = a 21 The coefficient of influence of rotational speed deviation on the rate of change of opening is given by the formula: a 21 = -Kp, where Kp is the proportional gain. Physically, it means that when the speed increases, the governor commands the guide vanes to close, resulting in a negative rate of change in the opening, which is usually a negative number.

[0143] a 22 The influence coefficient of the opening deviation on its own rate of change is given by the formula a. 22 = - , Let a be the time constant of the speed controller. 22 The physical meaning is that the larger the opening deviation, the slower the speed regulator adjusts the speed due to the damping effect, which is usually a negative number.

[0144] a 23 The influence coefficient of pressure deviation on the opening change rate, derived from pressure feedback control, is given by the formula: a 23 = -KpP, where KpP is the pressure proportional gain. Its physical meaning is that if the speed controller introduces pressure compensation, the speed controller will close the guide vanes when the pressure increases. It is usually a negative number, and it is 0 when there is no pressure feedback.

[0145] a 24 The influence coefficient of water level deviation on the rate of change of opening originates from the indirect effect of water level on pressure. A rise in water level leads to a rise in pressure, which in turn causes the governor to adjust. The formula is: , This is the partial derivative of pressure with respect to water level, which is usually negative. As water level rises, pressure rises, which closes the guide vane. It is 0 when there is no pressure feedback.

[0146] Third row: Pressure change rate = a 31 The coefficient of influence of rotational speed deviation on the rate of pressure change originates from the effect of rotational speed on flow rate. Increased rotational speed leads to increased flow rate and consequently, a change in pipeline pressure. It is typically a positive number, and the formula is: , Let be the partial derivative of pressure with respect to flow rate. Let be the partial derivative of flow rate with respect to rotational speed.

[0147] a 32 The coefficient of influence of opening deviation on the pressure change rate stems from the direct impact of opening degree on flow rate. Increased opening degree leads to increased flow rate and thus increased pressure. The formula is: , The partial derivative symbol is used to represent the rate of change of one variable when the influence of other variables is ignored; P is the actual pressure at the pipe node; Q is the flow rate in the pipe. This refers to the guide vane opening.

[0148] a 33 The influence coefficient of pressure deviation on its own rate of change originates from the hydraulic damping characteristics of the pipeline. The larger the pressure deviation, the faster the pressure decreases due to hydraulic damping. It is usually a negative number, and the formula is: Where: g is the acceleration due to gravity, A is the cross-sectional area of ​​the pipe, L is the pipe length, and f is the loss coefficient.

[0149] a 34The influence coefficient of water level deviation on the rate of pressure change stems from the direct relationship between the water level in the surge tank and the pipeline pressure. As the water level rises, the pressure at the bottom of the well rises; it is usually a positive number, and the formula is: a 34 = ,in This is the density of water.

[0150] The fourth row corresponds to the rate of change of water level. = a 41 The coefficient of influence of rotational speed deviation on the rate of change of water level originates from the effect of rotational speed on flow rate. Increased rotational speed leads to increased flow rate and consequently, a decrease in the water level of the surge tank. It is usually a negative number, and the formula is: , where: A t This refers to the cross-sectional area of ​​the pressure regulating well.

[0151] a 42 The influence coefficient of the opening deviation on the rate of change of water level stems from the direct impact of the opening on the flow rate. An increase in opening leads to an increase in flow rate and a decrease in water level. It is usually a negative number. The formula is: , where: A t This refers to the cross-sectional area of ​​the pressure regulating well. The sign of the partial derivative; Q is the flow rate in the pipe; This refers to the guide vane opening.

[0152] a 43 The influence coefficient of pressure deviation on the rate of change of water level stems from the coupling of pressure and flow rate. Increased pressure leads to decreased flow rate and thus a rise in water level. It is usually a positive number, and the formula is: Where: P is the actual pressure of the pipeline node.

[0153] a 44 The influence coefficient of water level deviation on its own rate of change originates from the damping characteristics of the surge tank. The larger the water level deviation, the slower the water exchange. It is usually a negative number, and the formula is: Where: a is the water hammer wave velocity, L t This is the equivalent length of a pressure regulating well.

[0154] The number of elements in matrix A increases or decreases with the state variables: If more variables are considered, such as multi-pipeline node pressure and multi-unit coupling, the matrix dimension will increase, but the physical meaning of the elements remains the same, which are all coupling coefficients between variables. In engineering applications, no manual calculation is required: The element values ​​are calculated by the SS program by reading the steady-state parameters and sub-model equations in the MPSO1.D file, calculating the partial derivatives and filling matrix A.

[0155] Step 3: Solve for eigenvalues ​​and determine intrinsic stability. Solve for the eigenvalues ​​λ of the system matrix A, in the form λ = σ + jω, where σ is the real part, ω is the imaginary part, and j is the imaginary unit. Determine the intrinsic stability of the system by the sign of the real part of the eigenvalues.

[0156] Solve for the eigenvalues ​​λ of the system matrix A of the linearized equation. The eigenvalues ​​satisfy the equation |A - λI| = 0, where I is the identity matrix, and λI is formed by replacing each diagonal element 1 of the identity matrix I with λ, creating a diagonal matrix. For example, if it is a 4x4 diagonal matrix: The system matrix A is solved by finding the eigenvalues ​​λ. The sign of the real part of the eigenvalues ​​is the basis for judging the stability of small fluctuations. If the real part of all eigenvalues ​​is less than 0, the system is stable. The sign of the real part of the eigenvalues ​​directly determines the stability of the system. The SS program solves all eigenvalues ​​using numerical algorithms such as QR decomposition and then judges the stability based on the distribution of eigenvalues. QR decomposition is a common algorithm in the field of numerical linear algebra and is common knowledge for numerical calculation personnel in water conservancy and hydropower engineering. Technical personnel in this field can implement QR decomposition by calling common numerical calculation libraries in Fortran such as the LAPACK library, which will not be elaborated here.

[0157] Criterion: If the real part of all eigenvalues ​​σ < 0: After a small disturbance, the parameter deviations such as speed and pressure will return to steady state in the form of damped oscillation or monotonically decaying, and the system is judged to be stable; If there is a real part of eigenvalues ​​σ > 0: The deviation will continue to amplify, and the system is judged to be unstable; If the real part of some eigenvalues ​​σ = 0: The deviation oscillates with equal amplitude, and the system is judged to be critically stable, but is considered unstable in engineering.

[0158] Conclusion of the embodiments of the present invention: The eigenvalues ​​are all negative real parts, such as λ1=-0.08+j1.2 and λ2=-0.12+j0.9, and are essentially stable.

[0159] Step 4: To better reflect actual operation, the final conclusion can be revised based on engineering constraints. Governor constraints: The guide vane opening must be within the range of 0~1.0 fully closed to fully open, with no risk of overtravel; Pressure constraints: The pipeline pressure must not exceed 1.3 times the design pressure, and there must be no vacuum damage in the tailrace pipe; Speed ​​constraints: The speed fluctuation range must not exceed ±5% of the rated speed.

[0160] The specific calculation of the attenuation coefficient in this invention is explained as follows: The attenuation coefficient is a quantitative indicator reflecting the rate of oscillation attenuation of parameters after a small disturbance in the system. The larger the value, the faster the disturbance attenuates, and the better the system stability. For example, the speed fluctuation attenuation coefficient of 0.05 mentioned in the embodiment indicates that the amplitude attenuates by 5% for each oscillation. The attenuation coefficient is calculated by the speed attenuation characteristic and quantitative indicator calculation program of the SPE program in the post-processing module. The attenuation coefficient is usually calculated in the following ways: Based on eigenvalue calculation: If the eigenvalue of the linearized system is λ=σ+jω (σ is the real part and ω is the imaginary part), the attenuation coefficient can be expressed as |σ|. The larger the absolute value of σ, the faster the attenuation; Based on time series data fitting: Through the time series curves of parameters such as speed and pressure in the MPSO1.D file, the attenuation law of oscillation amplitude is extracted and calculated according to the exponential attenuation formula. An Let A0 be the amplitude of the nth oscillation, A0 be the initial amplitude, and n be the number of oscillations. The attenuation coefficient is fitted and solved.

[0161] IV. Result Verification: The calculation results were compared with the field test data. The errors of key indicators such as unit speed increase, maximum pipeline pressure, and surge in the surge tank were ≤3%, which verified the accuracy of the calculation model of the present invention.

[0162] The results of the transient process calculation in this invention can be directly applied to four major engineering scenarios: design optimization, equipment selection, operation control, and safety assessment of pumped storage power stations, providing accurate data support for decision-making throughout the entire life cycle of pumped storage power stations.

[0163] During the power plant design phase, structural and system parameters are optimized, specifically including: surge tank design and selection: Based on the calculated maximum / minimum surge values ​​of the surge tank, the surge tank height is determined to avoid top or bottom leakage, and the cross-sectional area and type (e.g., impedance type, overflow type) are considered. For example, if calculations show a maximum surge of 3.2m under a certain operating condition, sufficient safety margin must be reserved for the top elevation of the surge tank. Combining the calculation results of multiple surge tank combinations, the spacing of the surge tank group is optimized to mitigate the pressure superposition risk caused by the coupling effect of the group of surge tanks.

[0164] Pipeline system parameter optimization: Adjust the pipe diameter, wall thickness and material according to the maximum / minimum pressure of each node. If the maximum pressure of the water intake pipe exceeds the design value, the pipe diameter can be increased or a material with a higher elastic modulus can be selected to reduce the impact of water hammer pressure.

[0165] Based on the segmented calculation results, the segment length of the water diversion / tailrace tunnel is optimized to a maximum of 100m, balancing calculation accuracy and engineering cost.

[0166] Unit layout and operating condition adaptation: Based on the calculation results of multiple units, such as 1-4 units operating in parallel, optimize the unit spacing and water diversion branch pipe structure to reduce the adverse effects of hydraulic interference on the transition process. Verify the feasibility of combined operating conditions such as start-up-sudden load shedding, adjust the unit start-up sequence and load switching timing to avoid superimposed risks.

[0167] During the equipment selection phase, key equipment parameters are matched, specifically including: governor selection and parameter calibration: based on indicators such as speed fluctuation attenuation coefficient and transition time, a PI or PID governor is selected, and core parameters such as permanent droop rate and buffer time constant are determined. For example, if the speed attenuation is too slow, the proportional gain Kp can be increased to optimize the guide vane break-line closing law. Combined with calculation results of guide vane failure and emergency pressure regulating valve operation conditions, the dynamic response capability of the governor is verified to ensure the effectiveness of regulation under extreme conditions.

[0168] Runner selection: Utilizing the calculation accuracy of the "runner characteristic curve conversion model" with an interpolation error ≤0.5%, compare the unit energy WH(x,y) data and unit torque WB(x,y) data of different runners to select the runner type suitable for all operating conditions of the turbine / pump / combined operating conditions. For curve connection requirements under small opening conditions, prioritize runners calculated using the appropriate length method to avoid sudden torque changes and flow fluctuations during operation.

[0169] Other auxiliary equipment selection: Based on the maximum vacuum degree of the tailrace pipe, select pipe valves and seals with cavitation resistance performance; based on the water hammer wave velocity calculation results, select pressure gauges and sensors that meet the pressure tolerance requirements to ensure accurate monitoring data.

[0170] During the operation and control phase, dynamic operation plans are formulated. For example, when calculations indicate that the surge tank has collapsed or leaked, or when calculations fail to converge, parameters are adjusted according to system prompts, such as increasing pipeline wave velocity or optimizing pipe segmentation, before restarting operation. Operating parameters are dynamically optimized by combining real-time monitoring data and calculation results to iteratively adjust the governor's broken-line closing time and inflection point opening, and to optimize speed decay characteristics, such as controlling the transition time within the engineering allowable range. For special operating conditions such as pump startup and guide vane malfunction, a variable-step-size calculation-optimized operating procedure is adopted to improve operational stability.

[0171] During the safety assessment and operation and maintenance phase, risks are anticipated in advance. This includes: pre-construction risk simulations, simulating dangerous conditions such as sudden 100% load shedding and power outages; predicting whether the unit's maximum speed and pipeline pressure exceed limits; and developing targeted prevention and control measures such as setting speed protection thresholds and pressure relief valves. As a basis for operation and maintenance, indicators such as pressure pulsation frequency and speed oscillation frequency are used to assess equipment fatigue wear risk. If the pressure pulsation frequency of a certain pipe section is close to the resonant frequency, priority should be given to maintenance and reinforcement. The error between long-term operational calculations and measured data is ≤3%, and model parameters are calibrated to improve the accuracy of subsequent operation and maintenance predictions. As a standard for project acceptance, the calculated extreme parameters, such as maximum / minimum pressure and surge extreme values, are used as indicators for project acceptance, ensuring that all power plant performance standards are met before formal operation.

[0172] The above embodiments are merely illustrative examples for clear explanation and are not intended to limit the implementation. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A numerical calculation model for the hydraulic transient process of a pumped storage power station, characterized in that, include: The system comprises a runner characteristic curve conversion model, a pipeline system parameter calculation model, and a transient process main calculation model. The runner characteristic curve conversion model includes a TRN program module and an NQM program module, used to achieve accurate conversion of the runner's full characteristic curve and smooth transition between the guide vane's full opening range. The TRN program module uses the Suter method to convert runner experimental data into standardized data and generate an MPSI1.D data file adapted to the guide vane's large opening. The NQM program module uses the length method to convert runner experimental data into standardized data and generate an MPSI5.D data file adapted to the guide vane's small and zero opening. The pipeline system parameter calculation model... The READY program module provides standardized basic parameters and boundary conditions for pipeline transient flow calculations. The main transient process calculation model is used to solve the dynamic parameters of the transient process under all operating conditions. The impeller characteristic curve conversion model and the pipeline system parameter calculation model are executed in parallel, respectively completing the standardized conversion of impeller characteristic data and the calculation of basic pipeline system parameters, providing input data for the main transient process calculation model. The main transient process calculation model adaptively switches between MPSI1.D or MPSI5.D data files according to the real-time guide vane opening, solves the transient flow and unit dynamic equations, and completes the dynamic parameter solution of the transient process under all operating conditions.

2. The numerical calculation model for the hydraulic transient process of a pumped storage power station according to claim 1, characterized in that: The TRN program module uses the Suter method to process the rotor experimental data, including the guide vane opening. Unit rotational speed N 11 Unit flow rate Q 11 Unit torque M 11 Through equal-interval processing using parabolic and linear interpolation, the data are converted into standardized data in the forms WH(x,y) and WB(x,y), where x is the guide vane opening α and y is the unit rotational speed N. 11 WH(x,y) corresponds to the guide vane opening α(x) and unit rotational speed N. 11 (y) represents the unit energy parameter of the runner; WB(x,y) corresponds to a guide vane opening α(x) and unit rotational speed N. 11 (y) refers to the unit torque parameter of the runner; the MPSI1.D data file includes standardized data in the form of WH(x,y) and WB(x,y), and the MPSI1.D data file is used for calculations under the condition that the guide vane opening is greater than 0.

1.

3. The numerical calculation model for the hydraulic transient process of a pumped storage power station according to claim 1, characterized in that: The NQM program module uses the length method, taking the flow characteristic curve length Lq as a unified benchmark, to transform the impeller characteristic curve into Lq-N. 11 Lq-Q 11 Lq-M 11 The correlation curve form is used to generate Lq-N. 11 Lq-Q 11 Lq-M 11 The MPSI5.D data file contains the correlation curve data. The MPSI5.D data file is adapted for calculations under conditions where the guide vane opening is less than or equal to 0.1 and zero opening; where N... 11 For unit speed, Q 11 Unit flow rate, M 11 The torque is expressed in units.

4. The numerical calculation model for the hydraulic transient process of a pumped storage power station according to claim 1, characterized in that: The pipeline system parameter calculation model takes as input the basic physical parameters of the pipeline and environment, uses an adaptive segmentation algorithm to divide the pipeline into segments and calculation nodes, and determines the calculation time step Δt based on the Courant condition including Δt≤Δx / a, where Δx is the segment length and a is the water hammer velocity. The output is an MPSI2.D data file including the water hammer velocity of the pipeline segment and the calculation time step, and an MPSI4.D data file including the elevation and node number of each calculation node.

5. The numerical calculation model for the hydraulic transient process of a pumped storage power station according to claim 1, characterized in that: The main computational model for the transient process includes a governor model, a surge tank transient model, and a unit boundary model. The governor model regulates the guide vane opening. The surge tank transient model simulates the water level fluctuations and dynamic changes in bottom pressure during the transient process. The unit boundary model establishes the unit operating boundary conditions for the turbine, pumps, and combined operating conditions, coupling the dynamic interaction between the water flow and mechanical systems. When the operating conditions change, and the unit speed or pipeline flow deviates from the steady state, the unit boundary model receives the dynamic changes in the operating conditions. Based on the speed deviation output by the unit boundary model, the governor model initiates the regulation logic to adjust the guide vane opening and change the flow rate. After the flow rate changes, the surge tank transient model calculates the water level surge and bottom pressure, and feeds the pressure change back to the unit boundary model. The unit boundary model receives new guide vane opening and pressure data, calculates speed and torque parameters, and transmits them again to the governor model and the surge tank transient model, iterating until the system stabilizes.

6. A numerical calculation system for the hydraulic transient process of a pumped storage power station, characterized in that: The numerical calculation system for the hydraulic transient process of the pumped storage power station integrates the numerical calculation model for the hydraulic transient process of the pumped storage power station as described in any one of claims 1-5.

7. The numerical calculation system for the hydraulic transient process of a pumped storage power station according to claim 6, characterized in that, The system includes: a preprocessing module, comprising a data formatting output program MMM and a surge tank data error detection program TYJ, used to adjust the format of input data, verify its integrity, and detect logical errors; a main calculation module, integrating a numerical calculation model of the hydraulic transient process of a pumped storage power station, which calls various subroutines through the MPS main program to complete the transient process calculation and output time series data and extreme value parameters; and a post-processing module, comprising SS, SPE, LAT, and DRAW programs, for performing small fluctuation stability analysis, speed decay characteristic evaluation, large fluctuation result summarization, and plotting data generation.

8. The numerical calculation system for the hydraulic transient process of a pumped storage power station according to claim 7, characterized in that: The MMM program includes data format mapping, which converts input parameters from different sources and in different formats, such as discrete tables of experimental data and design values ​​of engineering parameters, into text files with a unified structure, i.e., the data format commonly used in Fortran programs. The input to the preprocessing module is experimental data of turbine characteristics, pipeline system parameters, unit parameters, and operating condition setting parameters. The MMM program converts the input data into a standardized format and generates an intermediate data file. Then, the TYJ program performs integrity verification, logical verification, and judgment on the rationality and validity of pipeline parameters on the intermediate data file, and outputs abnormal prompt information.

9. The numerical calculation system for the hydraulic transient process of a pumped storage power station according to claim 7, characterized in that: The SS program is used for stability analysis under small fluctuation conditions to determine whether the system meets the small disturbance stability condition. The SPE program is used to calculate the speed transition time, speed fluctuation attenuation coefficient, and pressure pulsation frequency. The LAT program is used to summarize the extreme parameters of large fluctuation conditions, including the maximum or minimum pressure and corresponding time of each pipeline node, the highest or minimum speed of the unit and corresponding time, and the highest or minimum surge value and corresponding time of the surge tank. The DRAW program is used to generate plotting data in TECPLOT format and output visual curves of pressure-time, flow-time, speed-time, and opening-time.

10. A numerical calculation method for the hydraulic transient process of a pumped storage power station, characterized in that: Using any of the numerical calculation models described in claims 1-5 or any of the numerical calculation systems described in claims 6-9, Includes the following steps: S1: Data preparation and preprocessing: Organize and input the test data of the runner characteristics, pipeline system parameters, unit parameters and operating condition parameters. After formatting and error checking, generate a standardized data file; S2: Runner characteristic curve conversion: Run the TRN program module and the NQM program module to generate MPSI1.D data file and MPSI5.D data file respectively. Call the MPSI1.D data file or MPSI5.D data file according to the guide vane opening range; S3: Pipeline system parameter calculation. Run the READY program module to calculate the water hammer wave velocity and time step of the pipe section, divide the pipe section and calculation nodes, and generate MPSI2.D and MPSI4.D data files; S4: Transient process main calculation. Run the MPS main program, call each subroutine and the data files generated in steps S2 and S3, complete the transient process calculation, and output time series data and extreme value parameter data; S5: Post-processing and result analysis. Use the SS program, SPE program, LAT program and DRAW program to complete stability analysis, quantitative index calculation, extreme value parameter summary and visualization curve generation.

Images