Simulation and calculation method of aircraft engine surge dynamic process

Abstract

The present invention discloses a simulation calculation method for the dynamic process of surge of an aero-engine as a whole. The functions of each component of the engine are simplified and modeled, and the complex system of the engine is decomposed into a system composed of several simple functional units. By combining the units and adjusting the parameters of each unit, aero-engine systems of different configurations can be obtained. The control equation of each unit is a two-dimensional Euler equation with a source term. The parameters added to the control equation are determined according to the function of each unit. The four-part Runge-Kutta time marching is performed by determining the solution conditions of the inlet and outlet boundary parameters to obtain the characteristics of the engine under various working states. This method can take into account the retention of the main characteristics of the flow field in the dynamic process of the engine and the limitation of computing resources. It can perform quantitative analysis on typical working conditions in a short time, and has strong adaptability for studying different types of aero-engines. It has universal applicability for studying the dynamic process of engine surge.

Description

Technical Field

[0001] The present invention relates to the field of aero-engines, and in particular to a method for simulating and calculating the dynamic process of surge of an aero-engine as a whole. Background Art

[0002] As the heart of an aircraft, ensuring stable operation of the aircraft under varying flight conditions is essential for achieving high safety and reliability. Since the advent of gas turbine aircraft engines, research on their aerodynamic stability has continued to evolve alongside improvements in their performance. In modern times, the pursuit of high thrust-to-weight ratios and low fuel consumption rates has led to demands for high pressure ratios, high speeds, and high combustion chamber outlet temperatures, further increasing the importance of aerodynamic stability research. Modern aircraft engines, especially military aircraft engines, often operate under multiple extreme operating conditions within their flight envelope. These conditions are accompanied by various intake conditions caused by aircraft maneuvers and missile launches, which collectively affect the engine's stable operating range. This not only conflicts with the high thrust and efficiency performance requirements of aircraft engines, but in severe cases can even cause irreversible damage to the aircraft engine itself. Due to the operating characteristics of a compressor, it exerts a large axial pressure rise on the airflow within a short axial distance, which in turn generates a large axial load on the blades. Furthermore, since the compressor operates at high speed, the blades also experience large centrifugal loads. Once instability occurs, engine performance may degrade or even stall, or even worse, blades may break and fly out, causing casing containment failure and seriously impacting the safety of the aircraft and personnel. In the modern aeroengine design cycle, the study of the aerodynamic stability of the entire engine is playing an increasingly important role. Conducting aerodynamic stability studies during the design phase can quantitatively analyze the impact of various components on the stability of the entire engine, guide the optimization of each component, reduce the cost of subsequent improvements, and mitigate the risk of engine instability. Most current research on aeroengine surge and rotating stall is based on a simple one-dimensional model of the compression system. However, this approach cannot account for the nonuniform flow pattern along the circumferential dimension of the engine, and its mathematical model often fails to accurately describe the effect of aeroengine geometry on its post-stall characteristics. Therefore, it is necessary to develop a dynamic theoretical model that can describe the influence of intake distortion and engine geometry on the post-stall characteristics of the entire engine. Summary of the Invention

[0003] The goal to be achieved by the present invention is to perform simulation calculations on the surge dynamic process of an aircraft engine as a whole, which requires formulating a set of calculation processes and providing key algorithms in the calculation processes.

[0004] The present invention adopts the following technical solutions to realize the simulation calculation of the dynamic process of the whole aircraft engine surge:

[0005] The simulation calculation method of the dynamic process of the whole aircraft engine surge includes the following steps:

[0006] Step 1) Based on the engine configuration and the functions of each component, the engine is simplified into two-dimensional grid control volume units arranged along the axial and circumferential directions. The airflow parameters of each two-dimensional grid control volume unit are assumed to be uniformly distributed along the radial direction and only the changes in the airflow parameters along the axial and circumferential directions are considered. The airflow parameters stored in the two-dimensional grid control volume unit are determined, and the two-dimensional Euler equation with source term is used to describe the airflow of the two-dimensional grid control volume unit;

[0007] Step 2) Initialize the model flow field, calculate the parameters of each control unit through the given working point parameters and boundary conditions, and assign values ​​to them to establish the initial calculation domain field of the steady-state working point;

[0008] In step 3), based on the initial field established in step 2, the engine model is subjected to unsteady time-marching calculations. By changing the boundary parameters or adjusting the engine control laws, the engine performance under the target operating conditions is studied, especially the simulation calculation of the whole machine surge dynamic process is performed.

[0009] Preferably, in step 1), the two-dimensional Euler equations with source terms are used to describe the airflow of the two-dimensional grid control volume unit, which includes the continuity equation, the axial momentum equation, the circumferential momentum equation, and the energy equation. These four equations are used to describe the two-dimensional inviscid flow of the airflow. At the same time, by adding volume force to the axial momentum equation and the circumferential momentum equation, the force of the two-dimensional grid control volume unit on the airflow is simulated, reflecting the parameter changes of the airflow when the compressor works on the airflow and the airflow works on the turbine; by adding heat to the energy equation, the heating effect of the combustion chamber on the airflow is reflected, wherein:

[0010] Continuity equation:

[0011]

[0012] Axial momentum equation:

[0013]

[0014] Energy equation:

[0015]

[0016] Circumferential momentum equation:

[0017]

[0018] In the above control equations: ρ is the internal airflow density of the two-dimensional grid control unit; V is the volume of the two-dimensional grid control unit; is the surface airflow velocity of the two-dimensional grid control unit, with the outflow of the two-dimensional grid control unit as positive; S is the surface area of ​​a single two-dimensional grid control unit; C a and C u are the axial velocity and the circumferential velocity respectively; P is the surface pressure of the two-dimensional grid control volume unit; and are the axial and circumferential components of the surface pressure respectively; F a and F u are the axial and circumferential body forces respectively; E is the total energy of the gas per unit mass within the two-dimensional grid control unit; H is the total enthalpy per unit mass at the current total temperature of the airflow flowing into the boundary of the two-dimensional grid control unit; Q is the heat provided by the combustion chamber per unit time; W is the power provided or consumed by the compressor and turbine.

[0019] Preferably,

[0020] make Among them, ad1-ad4 represent the rate of change of airflow mass, axial momentum, energy and circumferential momentum in the two-dimensional grid control volume unit, and also reflect the flux of the above parameters at the boundary of the two-dimensional grid control volume; take static pressure P, static temperature T, mass flow rate G, circumferential velocity C u As the control parameter in the two-dimensional mesh control volume element, the time derivative of the grid center parameter of each two-dimensional mesh control volume element is obtained:

[0021]

[0022] Where R is the gas constant; C f is the mixing coefficient of the two-dimensional grid control volume unit, reflecting the ratio of the length of the leafless area in the grid to the total length of the grid; L is the grid length; α is the average airflow angle in the two-dimensional grid;

[0023] For each two-dimensional grid control body unit, the circumferential sides are connected by side boundaries, and the axial sides are connected by inlet and outlet boundaries; a y array is defined to store the grid center parameters of each two-dimensional grid control body unit, and the y array stores the static pressure, static temperature, mass flow rate and circumferential velocity of each two-dimensional grid control body unit in the order of axial direction first and circumferential direction later; if there are n two-dimensional grid control body units in the axial direction and mf two-dimensional grid control body units in the circumferential direction, then each axial sequence of the y array will store ns=4*n data, and the entire two-dimensional computational domain grid will store a total of nr=mf*ns data; in addition, four additional storage spaces are reserved for storing the relative speeds of the low, medium and high pressure rotors and the relative fuel quantities, so the final y array stores a total of nrh=nr+4 data.

[0024] Preferably, the implementation process of step 2) is:

[0025] According to the engine inlet condition, the total inlet pressure P of the two-dimensional computational domain grid * , total temperature T * , dimensionless density flow q(λ) and the geometric interface area S of each two-dimensional grid control unit, calculate the incoming flow mass flow rate G; through the law of conservation of mass, the mass flow rate of the outlet interface of each two-dimensional grid control unit is obtained, and the total pressure P of the outlet interface is obtained by * and total temperature T * The dimensionless dense flow q(λ) and velocity coefficient λ of the interface are obtained, and then the static pressure P, static temperature T, and density ρ parameters are obtained;

[0026] For the outlet interface total pressure P of the two-dimensional grid control volume element * and total temperature T * The calculation needs to be selected according to the type of two-dimensional mesh control volume element:

[0027] (1) The type of the two-dimensional grid control volume unit is a general unit: the general unit is the casing pipe unit, and the total temperature and total pressure of the inlet and outlet interfaces of the casing pipe unit are the same;

[0028] (2) The type of the two-dimensional grid control unit is a diversion unit: the diversion unit will distribute the airflow to the inner duct and the outer duct of the turbofan engine. The total temperature and total pressure of its outlet interface are the same as the inlet interface of the first unit of the inner and outer ducts. The distribution of its mass flow is achieved by artificially setting the bypass ratio of the current working point.

[0029] (3) The type of the two-dimensional grid control volume unit is a hybrid unit: the mass flow rate at the outlet of the hybrid unit is the sum of the mass flow rates at the interface of the inner and outer culverts: G3 = G1 + G2; the specific heat ratio at the interface of the hybrid unit outlet is: The gas constant at the outlet of the mixing unit is: The total temperature of the interface at the outlet of the mixing unit is: The total pressure at the mixing unit outlet interface is:

[0030] Where subscript 1 represents the interface parameters of the inner channel outlet, 2 represents the interface parameters of the outer channel outlet, and 3 represents the interface parameters of the mixing unit outlet; G is the mass flow rate; k is the specific heat ratio; R is the gas constant; T * is the total pressure; P * is the total pressure; σ is the total pressure recovery coefficient;

[0031] (4) The type of the two-dimensional grid control unit is the combustion unit: during the flow field initialization process, the fuel quantity is calculated by the formula given by the artificial combustion unit outlet temperature, and the outlet interface gas mass flow rate is obtained by adding the combustion unit inlet mass flow rate;

[0032] (5) The type of the two-dimensional grid control volume unit is the nozzle unit: if the tail nozzle operates in a supercritical state, the dimensionless density flow of its throat interface is artificially set to 1, and then the area and geometric parameters of the interface are recalculated;

[0033] (6) The type of the two-dimensional grid control unit is a compression unit: the compression unit interpolates the steady-state total pressure ratio π in its pressure ratio characteristic diagram through the dimensionless density flow q(λ) at its unit inlet and the current inlet interface temperature and rotor speed. * and total temperature ratio θ * , and then calculate the outlet interface total pressure P according to the inlet interface parameters * and total temperature T * Due to the complexity of the real engine system, it is necessary to give the steady-state operating point parameters of the compressor under the real engine working conditions and modify the mathematical model:

[0034] First, the inlet area of ​​the compression unit needs to be corrected:

[0035]

[0036] Where S0 is the original inlet area of ​​the compression unit, q(λ0) is the dimensionless dense flow result of the mathematical model parameters recursively transferred to the interface, and q(λ given ) is the artificially given compressor operating point inlet flow rate. Under the condition of keeping the mass flow rate unchanged, the current compression unit inlet area S' is corrected to meet the given value of the inlet dimensionless density flow;

[0037] According to the compression unit inlet temperature T1 and the current physical speed n, as well as the artificially given maximum reduced speed U 100 , and the current relative reduced speed is:

[0038]

[0039] However, the relative reduced speed of the compressor under actual engine operating conditions still needs to be given Therefore, the maximum reduced speed needs to be corrected:

[0040] After the correction is completed, even if the current relative reduced speed is equal to the given value under the actual engine working condition:

[0041]

[0042] (7) The type of the two-dimensional grid control volume unit is a turbine unit: the inlet area and reduced speed of the turbine unit are corrected in the same way as the compressor unit. In addition, the turbine characteristics must be additionally corrected to meet the balance of compressor and turbine power:

[0043] Total temperature ratio correction factor:

[0044]

[0045] Total pressure ratio correction factor:

[0046]

[0047]

[0048] Among them, W k and W t are the compressor work and turbine work before correction, θ * and π * are the turbine total temperature ratio and total pressure ratio before correction, and is the correction coefficient calculated in the last recursive process, and its initial value is 1; according to the above formula, the total temperature ratio correction coefficient C ft and total pressure ratio correction factor C fp Then the turbine characteristics are updated in the next round of parameter iteration:

[0049]

[0050] Through the above unit equations and the inlet and outlet boundary conditions of the entire engine system, all the airflow parameters of the control center and boundary of each two-dimensional grid unit are calculated, thereby completing the initialization of the entire engine two-dimensional grid calculation domain.

[0051] According to the initial field established in step 2), the engine model is calculated in an unsteady time-advancing manner; step 1 gives the calculation formula of the center parameters of each control body and its derivative, and stores the parameters of all units through the y array, so at the nth time layer, there will be a set of definite y n An array that records the parameters of each unit in the time layer; define a g array whose data format corresponds one-to-one with the y array to record the derivative of each parameter of the y array at the current time step;

[0052] The time marching term of the governing equation adopts the four-step Runge-Kutta time marching method, which is formulated as follows:

[0053]

[0054] Where h is the time step of single-step calculation, and f is the calculation function for the derivative of each parameter of the y array.

[0055] The calculation is performed in the same order as the y array, starting from the first axial 2D mesh control volume unit of each circumferential 2D mesh control volume unit to the last axial 2D mesh control volume unit. The inlet and outlet axial boundaries and left and right circumferential boundary control volume parameters of each 2D mesh control volume unit are found. The boundary flux is calculated using these parameters, and then the partial derivatives of the centroid parameters of each 2D mesh control volume unit are calculated and stored in the g array. Since the time marching term of the control equation is calculated using the four-step Runge-Kutta marching method, it is used to calculate the y array of the next time layer.

[0056] In the unsteady calculation, since the steady-state characteristics are used when interpolating the compressor characteristics, a good calculation effect can only be obtained at the steady-state operating point, which cannot meet the requirements of the post-stall characteristic research. Therefore, the first-order lag model of the compressor characteristics is adopted, and its formula is as follows:

[0057]

[0058] Among them, C is the actual pressure ratio of the compressor, C ss The interpolated pressure ratio of the steady-state characteristic is the larger the difference between the two, the greater the rate of change of the pressure ratio; the hysteresis factor τ is defined as the ratio of the axial length L of the compressor unit to the axial velocity C of the airflow. a The ratio of pressure ratio to pressure ratio represents the time required for the disturbance to be transmitted from the compressor inlet to the outlet. Its effect is that when the compressor flow rate is large, the response lag of the pressure ratio is small, and when the compressor flow rate is small, the response lag of the pressure ratio is also large:

[0059]

[0060] For the four-step Runge-Kutta time marching process, each step in the inner iteration process will re-interpolate the parameters of the previous time step to obtain a new steady-state characteristic pressure ratio. The actual pressure ratio should also be corrected accordingly, and the final pressure ratio is obtained after a complete time step. The formula of the iterative process imitates the iteration of the y array, and the formula is as follows:

[0061]

[0062] Compared with the prior art, the present invention adopts the above technical solution and has the following technical effects:

[0063] 1. Compared with the one-dimensional engine model and the three-dimensional unsteady calculation model, this two-dimensional model calculation method can preserve the main characteristics of the flow field during the engine dynamic process and limit the computing resources when studying the dynamic process of the entire aircraft engine surge. It can perform quantitative analysis of typical operating conditions in a short time, which is beneficial to the overall design of the engine.

[0064] 2. This two-dimensional model can simplify the mathematical models of various engine components. At the same time, by superimposing different engine control laws and changing the size and arrangement order of each component, this method has strong adaptability for studying different types of aircraft engines and is universal for studying the dynamic process of engine surge. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 It is a schematic diagram of the unit grid airflow parameters;

[0066] Figure 2 It is the data structure of the engine model;

[0067] Figure 3 It is the flow chart of initialization of flow field at steady-state working point;

[0068] Figure 4 This is the flow chart for calculating the parameters of the unsteady control body. Specific implementation methods

[0069] The technical solution of the present invention is further described in detail below with reference to the accompanying drawings:

[0070] In step 1), the airflow parameters of each unit are assumed to be uniformly distributed in the radial direction, and only the changes of the airflow parameters in the axial and circumferential directions are considered. According to the engine configuration and the functions of each component, the engine is simplified into two-dimensional grid control body units arranged in the axial and circumferential directions. Figure 1 The parameter relationship between the center and boundary of each control volume of the two-dimensional grid is shown, mainly including static pressure P, static temperature T, mass flow G, tangential velocity C u As the basic calculation variables, subscript 1 represents the axial inlet boundary, subscript 2 represents the axial outlet boundary, subscript L represents the circumferential left boundary, and subscript R represents the circumferential right boundary. Figure 1 The relationship between the grid center parameters and the boundary parameters is given: the grid center static pressure is equal to the axial inlet boundary parameter; the grid center static temperature, mass flow rate, and circumferential velocity are equal to the axial outlet boundary parameter; at the same time Figure 1 The calculation formulas for its left and right boundary parameters are also given.

[0071] As the basis for the simulation of the dynamic process of aircraft engine surge, the two-dimensional Euler equation with source term is used to describe the airflow within the computational domain grid. The governing equations used are:

[0072] Continuity equation:

[0073]

[0074] Axial momentum equation:

[0075]

[0076] Energy equation:

[0077]

[0078] Circumferential momentum equation:

[0079]

[0080] In the above control equation: ρ is the airflow density inside the control volume; V is the unit volume; is the airflow velocity on the unit surface, with the velocity outflowing from the control body as positive; C a and C u are the axial velocity and the circumferential velocity respectively; P is the surface pressure of the control body; F a and F u are the axial and circumferential body forces, respectively; E is the total energy per unit mass of the gas in the control volume; H is the total enthalpy per unit mass of the airflow flowing into the control volume boundary at the current total temperature; Q is the heat provided by the combustion chamber per unit time; W is the power provided or consumed by the compressor and turbine;

[0081] make The time derivative of each control body center parameter can be obtained:

[0082]

[0083] For each basic calculation unit, the units are connected by side boundaries in the circumferential direction and by inlet and outlet boundaries in the axial direction. Figure 2 The y array is defined to store the grid center parameters of each cell. It stores the static pressure, static temperature, mass flow rate, and circumferential velocity of each cell, first in the axial direction and then in the circumferential direction. If there are n cells in the axial direction and mf cells in the circumferential direction, the y array will store ns = 4 * n data points for each axial sequence, and nr = mf * ns data points for the entire two-dimensional computational domain grid. In addition, four additional storage spaces are reserved for storing the relative speeds of the low, medium, and high pressure rotors and the relative fuel quantities, so the y array will ultimately store nrh = nr + 4 data points.

[0084] Step 2) Initialize the model flow field, calculate the parameters of each control unit through the given working point parameters and boundary conditions, and assign values ​​to them to establish the initial calculation domain field of the steady-state working point. Figure 3 This paper demonstrates the initialization process for the computational domain flow field at a steady-state operating point. Its purpose is to establish the initial computational domain field at a steady-state operating point. Using the given operating point parameters and boundary conditions, the parameters of each control unit are calculated and assigned. By recursively deducing the mesh inlet boundary parameters of the entire machine model, the mesh flow field parameters for the entire computational domain can be obtained.

[0085] According to the total inlet pressure P * , total temperature T* , dimensionless density flow q(λ) and the geometric interface area S of each unit, the incoming mass flow rate G can be obtained. Through the law of conservation of mass, the mass flow rate of each unit outlet interface is obtained, and the total pressure P of the interface is obtained. * and total temperature T * The dimensionless density flow q(λ) and velocity coefficient λ of the interface can be obtained, and then the static pressure P, static temperature T, density ρ and other parameters can be obtained for the total pressure P of the unit outlet interface * and total temperature T * The calculation needs to be selected according to the unit type:

[0086] (1) General unit: This is usually a casing pipe unit. It is generally assumed that the flow inside the casing is isentropic, and the total temperature and total pressure at the inlet and outlet interfaces are the same. If there is a total pressure loss due to friction or blockage, the outlet total pressure can be calculated based on the total pressure recovery coefficient and the inlet total pressure.

[0087] (2) Diversion unit: The total temperature and total pressure at the outlet of the diversion unit, i.e., the interface between the inner and outer ducts, are equal to the diversion unit inlet parameters. The distribution of its mass flow is generally achieved by artificially setting the bypass ratio of the current working point.

[0088] (3) Mixing unit: The mass flow rate at the outlet of the mixing unit is the sum of the mass flow rates at the outlet interface of the inner and outer culverts. The remaining parameters of the outlet interface have been given in Section 2.1.5.

[0089] (4) Combustion unit: Section 2.1.6 gives the formula for calculating the fuel quantity by artificially setting the outlet temperature of the combustion unit during the flow field initialization process. The outlet interface gas mass flow rate can be obtained by adding the inlet mass flow rate of the combustion unit.

[0090] (5) Nozzle unit: Different from other units, the nozzle unit does not obtain the dimensionless dense flow by calculating the mass flow rate and the geometric interface area, but artificially makes the dimensionless dense flow of the interface equal to 1 and recalculates the interface area and geometric parameters.

[0091] (6) Compression unit: The compression unit needs to interpolate the steady-state total pressure ratio π from its pressure ratio characteristic diagram through the unit inlet dimensionless density flow q(λ) and the current inlet interface temperature and rotor speed. * and total temperature ratio θ * , and then calculate the outlet interface total pressure P according to the inlet interface parameters * and total temperature T *However, it's important to note that the mathematical model developed for numerical calculations cannot fully capture the actual engine operating point. Due to the complexity of the actual engine system, the compressor inlet airflow parameters, under identical inlet conditions and internal flow path geometry, often deviate from the numerically derived results. Therefore, it's necessary to determine the steady-state compressor operating point parameters under actual engine operating conditions so that the mathematical model can be modified to meet the mission requirements.

[0092] First, the compression unit inlet area needs to be corrected:

[0093]

[0094] Where S0 is the original inlet area of ​​the compression unit, q(λ0) is the dimensionless dense flow result of the mathematical model parameters recursively transferred to the interface, and q(λ given ) is the artificially given compressor operating point inlet flow rate. Under the condition of ensuring that the mass flow rate remains unchanged, the current compression unit inlet area S' is corrected to meet the given value of the inlet dimensionless density flow.

[0095] According to the compression unit inlet temperature T1 and the current physical speed n, as well as the artificially given maximum reduced speed U 100 , the current relative reduced speed can be obtained:

[0096]

[0097] However, the relative reduced speed of the compressor under actual engine operating conditions still needs to be given Therefore, the maximum reduced speed needs to be corrected:

[0098] After the correction is completed, the current relative reduced speed can be equal to the given value under the actual working conditions of the engine:

[0099]

[0100] (7) Turbine unit: The inlet area and reduced speed of the turbine unit are corrected in the same way as the compressor unit. In addition, the turbine characteristics must be additionally corrected to achieve the balance between compressor and turbine power:

[0101] Total temperature ratio correction factor:

[0102]

[0103] Total pressure ratio correction factor:

[0104]

[0105] Among them, W k and W t are the compressor work and turbine work before correction, θ* and π * are the turbine total temperature ratio and total pressure ratio before correction, and is the correction coefficient calculated in the last recursive process, and its initial value is 1. According to the above formula, the total temperature ratio correction coefficient C ft and total pressure ratio correction factor C fp Then the turbine characteristics can be updated in the next round of parameter iteration:

[0106]

[0107] Taking the dual-rotor mixed-displacement turbofan example, three flow field initialization recursions are required. The first is to calculate the flow field based on the given turbine characteristics, and the second and third are to perform correction recursions based on the number of turbine rotors. In addition, an additional final recursion is required to ensure that the pressure of the inner and outer duct outlet units is the same.

[0108] Step 3) Based on the initial field established in step 2, the engine model is subjected to unsteady time marching calculation. Step 1 gives the calculation formula of each control body center parameter and its derivative, and stores the parameters of all units through the y array. Therefore, at the nth time layer, there will be a set of definite y n An array that records the parameters of each unit in the time layer. Define a g array whose data format corresponds one-to-one with the y array to record the derivative of each parameter in the y array at the current time step.

[0109] The time marching term of the governing equation adopts the four-step Runge-Kutta time marching method, which is formulated as follows:

[0110]

[0111] Where h is the time step of single-step calculation, and f is the calculation function for the derivative of each parameter of the y array.

[0112] Figure 4 The process of differentiating each parameter in the y array is shown. The calculation is performed in the same control volume arrangement order as the y array, and the parameters of the inlet and outlet axial boundaries and the left and right circumferential boundaries of each control volume are found. The boundary flux is calculated using the above parameters, and then the partial derivatives of the center parameters of each control volume are calculated and stored in the g array. Since the time advancement term of the control equation is calculated using the four-step Runge-Kutta advancement method, Figure 4 The calculation process shown is used four times to calculate the y array of the next time layer.

[0113] In the unsteady calculation, since the steady-state characteristics are used when interpolating the compressor characteristics, a good calculation effect can only be obtained at the steady-state operating point, which cannot meet the requirements of the post-stall characteristic research. Therefore, the first-order lag model of the compressor characteristics is adopted, and its formula is as follows:

[0114]

[0115] Among them, C is the actual pressure ratio of the compressor, C ss The pressure ratio is interpolated for the steady-state characteristic. The greater the difference between the two, the greater the rate of change of the pressure ratio. The hysteresis factor τ is defined as the ratio of the axial length of the compressor unit to the axial velocity of the airflow. It represents the time required for the disturbance to propagate from the compressor inlet to the outlet. Its effect is that when the compressor flow rate is large, the response lag of the pressure ratio is small, and when the compressor flow rate is small, the response lag of the pressure ratio is also large:

[0116]

[0117] For the four-step Runge-Kutta time marching process, each step in the inner iteration process will re-interpolate the parameters of the previous time step to obtain a new steady-state characteristic pressure ratio. The actual pressure ratio should also be corrected accordingly, and the final pressure ratio is obtained after a complete time step. The formula of the iterative process imitates the iteration of the y array, and the formula is as follows:

[0118]

[0119] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. A method for simulating and calculating the dynamic process of surge of an aircraft engine, characterized in that: The following steps are involved: Step 1) Based on the engine configuration and the functions of each component, the engine is simplified into two-dimensional grid control volume units arranged along the axial and circumferential directions. The airflow parameters of each two-dimensional grid control volume unit are assumed to be uniformly distributed along the radial direction and only the changes in the airflow parameters along the axial and circumferential directions are considered. The airflow parameters stored in the two-dimensional grid control volume unit are determined, and the two-dimensional Euler equation with source term is used to describe the airflow of the two-dimensional grid control volume unit; Step 2) Initialize the model flow field, calculate the parameters of each control unit through the given working point parameters and boundary conditions, and assign values ​​to them to establish the initial calculation domain field of the steady-state working point; Step 3) Based on the initial field established in step 2, an unsteady time marching calculation is performed on the engine model. By changing boundary parameters or adjusting the engine control law, the engine performance under the target operating condition is studied, and the surge dynamic process of the whole engine is simulated. In step 1), the two-dimensional Euler equations with source terms are used to describe the airflow of the two-dimensional grid control volume unit, including the continuity equation, axial momentum equation, circumferential momentum equation, and energy equation. These four equations are used to describe the two-dimensional inviscid flow of the airflow. At the same time, by adding volume force to the axial momentum equation and the circumferential momentum equation to simulate the force of the two-dimensional grid control volume unit on the airflow, the parameter changes of the airflow when the compressor works on the airflow and the airflow works on the turbine are reflected; by adding heat to the energy equation, the heating effect of the combustion chamber on the airflow is reflected, where: Continuity equation: Axial momentum equation: Energy equation: Circumferential momentum equation: In the above control equations: ρ is the airflow density inside the two-dimensional grid control volume unit; V is the volume of the two-dimensional grid control volume unit; is the airflow velocity on the surface of the two-dimensional grid control unit, with the velocity outflowing from the two-dimensional grid control unit as positive; S is the surface area of ​​a single two-dimensional grid control unit; C a and C u are the axial velocity and the circumferential velocity respectively; P is the surface pressure of the two-dimensional grid control volume unit; and are the axial and circumferential components of the surface pressure respectively; F a and F u are the axial and circumferential body forces respectively; E is the total energy of the gas per unit mass within the two-dimensional grid control unit; H is the total enthalpy per unit mass at the current total temperature of the airflow flowing into the boundary of the two-dimensional grid control unit; Q is the heat provided by the combustion chamber per unit time; W is the power provided or consumed by the compressor and turbine.

2. The method for simulating and calculating the dynamic process of surge of an aircraft engine according to claim 1, characterized in that: make Among them, ad1-ad4 represent the rate of change of airflow mass, axial momentum, energy and circumferential momentum in the two-dimensional grid control volume unit, and also reflect the flux of the above parameters at the boundary of the two-dimensional grid control volume; take static pressure P, static temperature T, mass flow rate G, circumferential velocity C u As the control parameter in the two-dimensional mesh control volume element, the time derivative of the grid center parameter of each two-dimensional mesh control volume element is obtained: Where R is the gas constant; C f is the mixing coefficient of the two-dimensional grid control volume unit, reflecting the ratio of the length of the leafless area in the grid to the total length of the grid; L is the grid length; α is the average airflow angle in the two-dimensional grid; For each two-dimensional grid control body unit, the circumferential sides are connected by side boundaries, and the axial sides are connected by inlet and outlet boundaries; a y array is defined to store the grid center parameters of each two-dimensional grid control body unit, and the y array stores the static pressure, static temperature, mass flow rate and circumferential velocity of each two-dimensional grid control body unit in the order of axial direction first and circumferential direction later; if there are n two-dimensional grid control body units in the axial direction and mf two-dimensional grid control body units in the circumferential direction, then each axial sequence of the y array will store ns=4*n data, and the entire two-dimensional computational domain grid will store a total of nr=mf*ns data; in addition, four additional storage spaces are reserved for storing the relative speeds of the low, medium and high pressure rotors and the relative fuel quantities, so the final y array stores a total of nrh=nr+4 data.

3. The method for simulating and calculating the dynamic process of surge of an aircraft engine according to claim 1, characterized in that: The implementation process of step 2) is: According to the engine inlet condition, the total inlet pressure P of the two-dimensional computational domain grid * , total temperature T * , dimensionless density flow q(λ) and the geometric interface area S of each two-dimensional grid control unit, calculate the incoming flow mass flow rate G; through the law of conservation of mass, the mass flow rate of the outlet interface of each two-dimensional grid control unit is obtained, and the total pressure P of the outlet interface is obtained by * and total temperature T * The dimensionless dense flow q(λ) and velocity coefficient λ of the interface are obtained, and then the static pressure P, static temperature T, and density ρ parameters are obtained; For the outlet interface total pressure P of the two-dimensional grid control volume element * and total temperature T * The calculation needs to be selected according to the type of two-dimensional mesh control volume element: (1) The type of the two-dimensional grid control volume unit is a general unit: the general unit is the casing pipe unit, and the total temperature and total pressure of the inlet and outlet interfaces of the casing pipe unit are the same; (2) The type of the two-dimensional grid control unit is a diversion unit: the diversion unit will distribute the airflow to the inner duct and the outer duct of the turbofan engine. The total temperature and total pressure of its outlet interface are the same as the inlet interface of the first unit of the inner and outer ducts. The distribution of its mass flow is achieved by artificially setting the bypass ratio of the current working point. (3) The type of the two-dimensional grid control volume unit is a hybrid unit: the mass flow rate at the outlet of the hybrid unit is the sum of the mass flow rates at the interface of the inner and outer culverts: G3 = G1 + G2; the specific heat ratio at the interface of the hybrid unit outlet is: The gas constant at the outlet of the mixing unit is: The total temperature of the interface at the outlet of the mixing unit is: The total pressure at the mixing unit outlet interface is: Where subscript 1 represents the interface parameters of the inner channel outlet, 2 represents the interface parameters of the outer channel outlet, and 3 represents the interface parameters of the mixing unit outlet; G is the mass flow rate; k is the specific heat ratio; R is the gas constant; T * is the total pressure; P * is the total pressure; σ is the total pressure recovery coefficient; (4) The type of the two-dimensional grid control unit is the combustion unit: during the flow field initialization process, the fuel quantity is calculated by the formula given by the artificial combustion unit outlet temperature, and the outlet interface gas mass flow rate is obtained by adding the combustion unit inlet mass flow rate; (5) The type of the two-dimensional grid control volume unit is the nozzle unit: if the tail nozzle operates in a supercritical state, the dimensionless density flow of its throat interface is artificially set to 1, and then the area and geometric parameters of the interface are recalculated; (6) The type of the two-dimensional grid control unit is a compression unit: the compression unit interpolates the steady-state total pressure ratio π in its pressure ratio characteristic diagram through the dimensionless density flow q(λ) at its unit inlet and the current inlet interface temperature and rotor speed. * and total temperature ratio θ * , and then calculate the outlet interface total pressure P according to the inlet interface parameters * and total temperature T * Due to the complexity of the real engine system, it is necessary to give the steady-state operating point parameters of the compressor under the real engine working conditions and modify the mathematical model: First, the inlet area of ​​the compression unit needs to be corrected: Where S0 is the original inlet area of ​​the compression unit, q(λ0) is the dimensionless dense flow result of the mathematical model parameters recursively transferred to the interface, and q(λ given ) is the artificially given compressor operating point inlet flow rate. Under the condition of keeping the mass flow rate unchanged, the current compression unit inlet area S' is corrected to meet the given value of the inlet dimensionless density flow; According to the compression unit inlet temperature T1 and the current physical speed n, as well as the artificially given maximum reduced speed U 100 , and the current relative reduced speed is: However, the relative reduced speed of the compressor under actual engine operating conditions still needs to be given Therefore, the maximum reduced speed needs to be corrected: After the correction is completed, even if the current relative reduced speed is equal to the given value under the actual engine working condition: (7) The type of the two-dimensional grid control volume unit is a turbine unit: the inlet area and reduced speed of the turbine unit are corrected in the same way as the compressor unit. In addition, the turbine characteristics must be additionally corrected to meet the balance of compressor and turbine power: Total temperature ratio correction factor: Total pressure ratio correction factor: Among them, W k and W t are the compressor work and turbine work before correction, θ * and π * are the turbine total temperature ratio and total pressure ratio before correction, and is the correction coefficient calculated in the last recursive process, and its initial value is 1; according to the above formula, the total temperature ratio correction coefficient C ft and total pressure ratio correction factor C fp Then the turbine characteristics are updated in the next round of parameter iteration: Through the above unit equations and the inlet and outlet boundary conditions of the entire engine system, all the airflow parameters of the control center and boundary of each two-dimensional grid unit are calculated, thereby completing the initialization of the entire engine two-dimensional grid calculation domain.

4. The method for simulating and calculating the dynamic process of surge of an aircraft engine as claimed in claim 2, wherein: According to the initial field established in step 2), the engine model is calculated in an unsteady time-advancing manner; step 1 gives the calculation formula of the center parameters of each control body and its derivative, and stores the parameters of all units through the y array, so at the nth time layer, there will be a set of definite y n An array that records the parameters of each unit in the time layer; define a g array whose data format corresponds one-to-one with the y array to record the derivative of each parameter of the y array at the current time step; The time marching term of the governing equation adopts the four-step Runge-Kutta time marching method, which is formulated as follows: Where h is the single-step calculation time step, and f is the calculation function for the derivative of each parameter of the y array; The calculation is performed in the same order as the y array, starting from the first axial 2D mesh control volume unit of each circumferential 2D mesh control volume unit to the last axial 2D mesh control volume unit. The inlet and outlet axial boundaries and left and right circumferential boundary control volume parameters of each 2D mesh control volume unit are found. The boundary flux is calculated using these parameters, and then the partial derivatives of the centroid parameters of each 2D mesh control volume unit are calculated and stored in the g array. Since the time marching term of the control equation is calculated using the four-step Runge-Kutta marching method, it is used to calculate the y array of the next time layer. In the unsteady calculation, since the steady-state characteristics are used when interpolating the compressor characteristics, a good calculation effect can only be obtained at the steady-state operating point, which cannot meet the requirements of the post-stall characteristic research. Therefore, the first-order lag model of the compressor characteristics is adopted, and its formula is as follows: Among them, C is the actual pressure ratio of the compressor, C ss The interpolated pressure ratio of the steady-state characteristic is the larger the difference between the two, the greater the rate of change of the pressure ratio; the hysteresis factor τ is defined as the ratio of the axial length L of the compressor unit to the axial velocity C of the airflow. a The ratio of , which represents the time required for the disturbance to be transmitted from the compressor inlet to the outlet, has the effect that when the compressor flow rate is large, the response lag of the pressure ratio is small, and when the compressor flow rate is small, the response lag of the pressure ratio is also large: For the four-step Runge-Kutta time marching process, each step in the inner iteration process will re-interpolate the parameters of the previous time step to obtain a new steady-state characteristic pressure ratio. The actual pressure ratio should also be corrected accordingly, and the final pressure ratio is obtained after a complete time step. The formula of the iterative process imitates the iteration of the y array, and the formula is as follows:

Images