A method for dynamic modeling of aeroengine based on pseudo-Jacobian matrix

Abstract

The present invention discloses a method for dynamic modeling of an aeroengine based on a pseudo-Jacobian matrix. Aiming at the problem that the component-level model of the Jacobian matrix calculated by bidirectional difference has poor real-time performance, the present invention linearizes the common equation of the turboshaft engine based on the nonlinear model dynamic linearization strategy, constructs a pseudo-Jacobian matrix evaluation function, calculates the Jacobian matrix by optimal estimation approximation, obtains a recursive calculation method of the pseudo-Jacobian matrix, avoids repeated calls to the component model, and greatly improves the real-time performance of the dynamic model.

Description

Technical Field

[0001] The invention belongs to the field of system control and simulation in aerospace propulsion theory and engineering, and specifically relates to a method for modeling a dynamic real-time model of an aero-engine. Background Art

[0002] Numerical simulation technology of aircraft engines plays an extremely important role in the development of aircraft engines and the design and verification of control systems. Using mathematical models instead of real engines to carry out numerical simulation or semi-physical simulation can conduct preliminary verification of control algorithms, fault diagnosis technologies, health management technologies, etc., saving experimental costs and reducing experimental risks. Since the late 1980s, Western countries have successively carried out research on aircraft engine simulation technology and developed a variety of numerical simulation systems for aircraft engines. After years of research, development and application, the relevant theoretical technology has matured, and the simulation confidence and accuracy have reached a very high level. With the development of model-based control technology, the requirements for model accuracy and real-time performance have been further improved, and a lot of research work has also been carried out.

[0003] The mathematical model of an aero-engine introduces the basic parameters of each component of a real engine, establishes aerodynamic thermodynamic models of each component along the engine gas path and according to the characteristics of the components, and obtains the parameters of each component that matches the current working state by solving the common working equations between the components. The convergence speed and accuracy of the common working equations are directly related to the accuracy and real-time performance of the mathematical model of the engine.

[0004] The common methods for solving the common working equations of aircraft engines are the Newton-Raphson (NR) method and the Broyden quasi-Newton method. The NR method needs to repeatedly call the engine component model to calculate the Jacobian matrix, which seriously affects the real-time performance of the model; while the quasi-Newton method approximates the Jacobian matrix and obtains the approximate Jacobian matrix through recursion, which reduces the computational complexity. However, the algorithm iterative convergence requires a smaller deviation of the initial solution and is affected by the calculation step size, so the convergence ability is weaker than the NR method.

[0005] Based on the above two methods, Chen Yuchun et al. from Northwestern Polytechnical University [1] The convergence problem of the solution is improved by changing the calculation step size of the NR method and limiting the initial value, but the Jacobian matrix still needs to be calculated by difference. Liao Guanghuang et al. from Nanjing University of Aeronautics and Astronautics [2]A neural network-based method for solving the engine joint working equation was proposed. Offline training data calculated based on the NR method was collected, and the residual was used as input. The correction value of the joint working equation guess was used as output. The neural network was used for training, which avoided the Jacobian matrix calculation and the iterative solution process of the equation, and effectively improved the real-time performance of the model. However, the neural network adopted an offline training method, and the model accuracy was affected by the generalization ability of the network, and the adaptability to uncertainty was weak. Wang Yuan, Wu Qian, et al. from Nanjing University of Aeronautics and Astronautics [3-4] Combining the square convergence of the NR method with the superlinear convergence characteristics of the quasi-Newton method, a self-correcting Broyden quasi-Newton method is proposed. By judging the convergence trend of the model and adjusting the calculation step size, the convergence accuracy and speed are improved. However, at working points with poor convergence trends, the correction matrix still needs to be initialized by difference, resulting in low real-time performance at individual working points. Pang Shuwei, Nanjing University of Aeronautics and Astronautics [5] A method for establishing engine component-level models based on precise partial derivatives is proposed. The Jacobian matrix is ​​constructed by replacing differential calculations with precise partial derivatives, which improves the real-time performance of the model. However, the algorithm calculation process is complex and difficult to implement. Summary of the invention

[0006] The technical problem to be solved by the present invention is to overcome the deficiencies of the prior art and provide an aero-engine dynamic modeling method based on a pseudo-Jacobian matrix, which can avoid repeated calls to the model by a differential algorithm and effectively improve the real-time performance of the model.

[0007] The present invention specifically adopts the following technical solutions to solve the above technical problems:

[0008] A method for dynamic modeling of an aeroengine based on a pseudo-Jacobian matrix is ​​proposed. The method constructs and solves the dynamic co-working balance equation of the aeroengine during the working process of the aeroengine to obtain the parameters of each component of the aeroengine that matches the current working state. When solving the dynamic co-working balance equation of the aeroengine, the guessed value of the dynamic co-working balance equation is iteratively corrected by the following method so that the current dynamic co-working balance equation meets the convergence condition:

[0009] Step 1, linearize the dynamic joint working equilibrium equation at time k;

[0010] Step 2: Based on the following optimization objective function, calculate the optimal estimate of the Jacobian matrix J(k) of the linearized dynamic joint working equilibrium equation: (pseudo-Jacobian matrix):

[0011]

[0012] Among them, ε(k) and ε(k-1) represent the residuals of the dynamic joint working balance equation at time k and time k-1 respectively, Δx(k)=x(k)-x(k-1) represents the difference between the guessed values ​​x(k) and x(k-1) of the dynamic joint working balance equation at time k and time k-1, μ is the penalty factor, The optimal estimate of the Jacobian matrix representing the linearized dynamic co-working equilibrium equation at time k-1;

[0013] Step 3: Update the dynamic joint working balance equation guess value x(k+1) at time k+1 through iterative calculation:

[0014]

[0015] Among them, λ∈(0,1] is the iteration step size.

[0016] Preferably, the optimal estimate Calculated by the following formula:

[0017]

[0018] Among them, η∈(0,1] is the step size factor, Δε(k)=ε(k)-ε(k-1) represents the difference between the residuals ε(k) and ε(k-1) of the dynamic joint working equilibrium equation at time k and time k-1, and the superscript “T” represents the matrix transpose, “|| || 2 ” represents the square of the vector 2 norm.

[0019] Compared with the prior art, the technical solution of the present invention has the following beneficial effects:

[0020] (1) High real-time performance: The present invention replaces the Jacobian matrix calculation in differential form with the recursive pseudo-Jacobian matrix calculation, which greatly reduces the number of component-level model calls and greatly improves the real-time performance of dynamic model calculation.

[0021] (2) Strong versatility: The calculation of the pseudo-Jacobian matrix is ​​based on the dynamic linearization method, which is itself a data-driven method with strong adaptability to changes in engine characteristics. It can be used for the calculation of dynamic real-time models of various types of engines. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 It is a structural schematic diagram of a twin-rotor turboshaft engine;

[0023] Figure 2 This is a flow chart of dynamic modeling of turboshaft engine based on pseudo-Jacobian matrix;

[0024] Figure 3 Comparison of key engine output parameters between the bidirectional difference algorithm and the pseudo-Jacobian matrix algorithm of the present invention (1km altitude);

[0025] Figure 4 Comparison of key engine output parameters between the bidirectional difference algorithm and the pseudo-Jacobian matrix algorithm of the present invention (3km altitude). DETAILED DESCRIPTION

[0026] In view of the problem that the existing Jacobian matrix calculation based on the difference method has poor real-time performance, the present invention combines the nonlinear model dynamic linearization method and uses optimal estimation to calculate the pseudo-Jacobian matrix, thereby improving the real-time performance of the dynamic model.

[0027] The present invention specifically adopts the following technical solutions to solve the above technical problems:

[0028] A method for dynamic modeling of an aeroengine based on a pseudo-Jacobian matrix is ​​proposed. The method constructs and solves the dynamic co-working balance equation of the aeroengine during the working process of the aeroengine to obtain the parameters of each component of the aeroengine that matches the current working state. When solving the dynamic co-working balance equation of the aeroengine, the guessed value of the dynamic co-working balance equation is iteratively corrected by the following method so that the current dynamic co-working balance equation meets the convergence condition:

[0029] Step 1, linearize the dynamic joint working equilibrium equation at time k;

[0030] Step 2: Based on the following optimization objective function, calculate the optimal estimate of the Jacobian matrix J(k) of the linearized dynamic joint working equilibrium equation: (pseudo-Jacobian matrix):

[0031]

[0032] Among them, ε(k) and ε(k-1) represent the residuals of the dynamic joint working balance equation at time k and time k-1 respectively, Δx(k)=x(k)-x(k-1) represents the difference between the guessed values ​​x(k) and x(k-1) of the dynamic joint working balance equation at time k and time k-1, μ is the penalty factor, The optimal estimate of the Jacobian matrix representing the linearized dynamic co-working equilibrium equation at time k-1;

[0033] Step 3: Update the dynamic joint working balance equation guess value x(k+1) at time k+1 through iterative calculation:

[0034]

[0035] Among them, λ∈(0,1] is the iteration step size.

[0036] The applicable objects of the dynamic modeling method of the present invention include but are not limited to turboshaft engines, turboprop engines, turbofan engines, variable cycle engines, combined engines, etc.

[0037] In order to facilitate public understanding, the technical solution of the present invention is described in detail below through a specific embodiment and in conjunction with the accompanying drawings:

[0038] This embodiment uses Figure 1 The twin-rotor turboshaft engine shown in the figure is the research object. The gas turbine and compressor of the twin-rotor turboshaft engine are coaxial, and the power turbine and the main rotor are coaxial. During the working process, air enters from the air inlet, flows to the combustion chamber through the compressor, mixes with the fuel and burns to produce high-temperature and high-pressure gas. When the gas flows through the gas turbine, a part of the energy is extracted to drive the compressor. The gas leaving the gas turbine flows through the power turbine, and the power turbine extracts the remaining energy to drive the main rotor to meet the working needs of the helicopter.

[0039] for Figure 1 The twin-rotor turboshaft engine shown in the figure, the dynamic modeling method proposed by the present invention is as follows: Figure 2 As shown, t is the dynamic model simulation time, k is the number of iterations for solving the equilibrium equation, and T is the simulation step. The method specifically includes the following steps:

[0040] Step A: Obtain the flight conditions at time t (helicopter flight altitude H, forward flight speed V x ) and fuel flow W f ;

[0041] Step B, calculating the thermal parameters of each component and the rotor load power requirement in sequence from the air inlet to the tail nozzle along the airflow process;

[0042] Step C, constructing a dynamic common working balance equation of the turboshaft engine;

[0043] Under the conditions of known flight conditions and fuel flow, the analysis is carried out according to the flow continuity and pressure balance conditions that must be met during the dynamic operation of the twin-rotor turboshaft engine. Usually three common equations are selected, denoted by ε i , i=1,2,3 are the residuals of the dynamic joint working equation, including

[0044] (1) Gas turbine inlet flow continuity equation φ1(x):

[0045] φ1(x)=(W 41xs -Q 41xs ) / Q 41xs =ε1 (1)

[0046] Where x is the guess value in the dynamic joint working equation, W 41xs is the similar flow rate of the gas turbine inlet under the current working conditions calculated from the gas turbine component characteristic curve, Q 41xs is the similar flow rate from the gas turbine guide vane into the gas turbine obtained by calculating the gas path process, ε1 is the residual of equation φ1(x);

[0047] (2) Continuity equation of power turbine inlet flow φ2(x):

[0048] φ2(x)=(W 44xs -Q 44xs ) / Q 44xs =ε2 (2)

[0049] Among them, W 44xs is the similar flow rate of the power turbine inlet under the current working conditions calculated from the power turbine component characteristic curve, Q 44xs is the similar flow rate from the power turbine guide vane into the power turbine obtained by calculating the gas path process, ε2 is the residual of equation φ2(x);

[0050] (3) Tail nozzle outlet pressure balance equation φ3(x):

[0051] φ3(x)=(p c7 -p7) / p7=ε3 (3)

[0052] Among them, p c7 is the total pressure of the airflow entering the tail nozzle outlet, p7 is the nozzle back pressure, ε3 is the residual of equation φ3(x);

[0053] In the component-level model working balance equation, the guess value x=[Z C Z G Z P ] T ; Among them, Z C is the compressor pressure ratio coefficient, Z G is the gas turbine pressure ratio coefficient, Z P is the power turbine pressure ratio coefficient; Z c For example, the pressure ratio coefficient is defined as:

[0054]

[0055] Among them, π c is the compressor pressure ratio at the current operating point. The subscripts min and max represent the minimum and maximum pressure ratio on the current speed line, respectively.

[0056] Step D: determine whether the dynamic joint working balance equation meets the convergence condition. If flag=1, jump to step F; otherwise, flag=0 and execute step E;

[0057]

[0058] Among them, k max is the maximum number of iterations allowed.

[0059] The convergence conditions in this embodiment are as follows:

[0060] ε i <10 -5 ,i=1,2,3,or k>k max .

[0061] Step E, calculate the pseudo-Jacobian matrix, correct the guessed value, and return to step B; specifically includes the following sub-steps:

[0062] Step E1, linearize the dynamic joint working equilibrium equation at time k:

[0063] ε(k)=φ(x(k)) (6)

[0064] Assume that the partial derivative of φ(x(k)) with respect to x is continuous and satisfies ||ε(k)-ε(k-1)||≤b||x(k)-x(k-1)||. That is, if the guess value changes finitely, the change of the residual of the equation is also finite. This condition is satisfied when the model converges.

[0065] Linearize equation (6) to obtain the following linearized equation:

[0066]

[0067] Among them, Δε(k)=ε(k)-ε(k-1), Δx(k)=x(k)-x(k-1), is the pseudo-Jacobian matrix; Equation (7) is the linearized equilibrium equation;

[0068] Step E2: Calculate the pseudo-Jacobian matrix using the optimal estimation method

[0069] The specific evaluation indicators are as follows:

[0070]

[0071] Among them, μ is the penalty factor, J(k) is the Jacobian matrix;

[0072] Formula (8) shows that when L(J(k)) = 0, ε(k)-ε(k-1) = J(k)Δx(k), That is, the linearized equation (6) holds, and the change of the Jacobian matrix is ​​constrained;

[0073] The Jacobian matrix is ​​shown in formula (9):

[0074]

[0075] Since the working process of an aircraft engine is a highly nonlinear and complex aerodynamic and thermal process, its mathematical modeling relies on the interpolation and iterative calculation of component characteristic diagrams and aerodynamic and thermal parameters.i (x), i = 1, 2, 3 is an implicit equation that contains the calculation process of each component of the engine. It is impossible to directly calculate the partial derivatives, and it can only be obtained through the method of numerical difference. In order to improve the accuracy of the model, the intermediate difference method is usually used to calculate the partial derivatives in the Jacobian matrix, that is:

[0076]

[0077] In the process of calculating the partial derivative of formula (10), it is necessary to guess the value x i Small disturbances upward and downward are performed respectively, and the component-level mathematical model is called to perform differential calculations. There are three guessed values, so a total of 6 calculations are required, which seriously affects the real-time performance of the dynamic model. For this reason, the present invention adopts a recursive method to calculate the pseudo-Jacobian matrix.

[0078] Find the minimum value of J(k) for equation (8) to obtain the optimal estimate of J(k) have:

[0079]

[0080] Move item sorting

[0081]

[0082] Subtract both sides of the equation

[0083]

[0084] So there is

[0085]

[0086] Among them, η∈(0,1] is the step size factor;

[0087] When t = 0, Initialize by difference, otherwise, Here, the subscript t represents the simulation time of the dynamic model, that is, the initial value of the pseudo-Jacobian matrix at the current simulation time is the pseudo-Jacobian matrix value at the end of the iteration of the previous simulation time;

[0088] Formula (14) calculates the pseudo-Jacobian matrix in a recursive manner. During the calculation process, only the increment of the residual of the current equation and the increment of the guessed value are required. There is no need to call the component model for differential calculation, which greatly simplifies the calculation process.

[0089] Step E3: Update the guess value through iterative calculation:

[0090]

[0091] Among them, λ∈(0,1] is the iteration step size.

[0092] Step F, calculating the rotor acceleration and updating the rotation speed;

[0093]

[0094] Among them, n G 、n P is the gas turbine and power turbine speed, J G , J P is the moment of inertia of the gas turbine rotor and the power turbine rotor, η G , η P is the mechanical efficiency of the gas turbine rotor and the power turbine rotor, W GT , W PT is the gas turbine and power turbine power, W C , W R is the compressor and helicopter rotor power.

[0095] In order to verify the advantages of the pseudo-Jacobian matrix algorithm used in the present invention in solving the dynamic joint working equations of the turboshaft engine, a simulation is performed on the component-level model of the T700 turboshaft engine, and the initial flight altitude of the turboshaft engine is set to H = 1 km and the initial forward flight speed V x =20m / s, the dynamic working process of the engine is simulated during the simulation, and the models of the two algorithms are called 10,000 times respectively. The model is developed based on VC++6.0, the computer system is Windows 10 Home Edition, CPU i5-8250U1.60GHz, 8GB RAM. The average calculation time of the pseudo-Jacobian matrix algorithm model for 10,000 calls is 3.19s, while the average calculation time of the bidirectional difference calculation Jacobian matrix model is 6.97s. It can be seen that the time used by the pseudo-Jacobian matrix algorithm is about 3 / 7 of the bidirectional difference algorithm, which verifies the superiority of the algorithm in real-time performance.

[0096] The accuracy of the model is further verified within the envelope range. Two sets of simulation results are given here. Figure 3 , Figure 4 As shown in the figure, the simulation results in other flight areas and engine working conditions are similar. In the figure, PJ (Pseudo-Jacobian) represents the component model simulation effect using the pseudo-Jacobian matrix algorithm; BDD (Bi-directional differential) represents the model simulation effect using the bidirectional differential Jacobian matrix algorithm. The figure shows the change curve of the key parameters of the turboshaft engine, where sfc is the fuel consumption rate, T 41 is the total temperature at the gas turbine inlet, e represents the output deviation between the PPD algorithm and the differential algorithm model, and the deviation calculation formula is:

[0097]

[0098] (1) Simulation test at a flight altitude of 1 km

[0099] Set the helicopter flight altitude H = 1km, forward flight speed V x From 20m / s to 50m / s, the output parameter change curve is as follows Figure 3 shown.

[0100] like Figure 3 It can be seen that compared with the bidirectional difference algorithm, the speed deviation of the component-level model output based on the pseudo-Jacobian matrix is ​​less than 0.15%, the temperature deviation is less than 0.25%, and the maximum deviation of the fuel flow in the closed-loop control is less than 3% due to its own small order of magnitude. The calculation deviation of the fuel consumption rate caused by the fuel flow deviation is less than 1.2%. The detailed error information is shown in Table 1. From Table 1 and Figure 3 It can be seen that the mathematical model of the turboshaft engine based on the pseudo-Jacobian matrix has achieved high accuracy, with an average error of less than 0.14%, which verifies the effectiveness of the algorithm in this paper in terms of dynamic model accuracy.

[0101] Table 1 Errors of the pseudo-Jacobian matrix model compared to the bidirectional difference model

[0102]

[0103] (2) Simulation test at a flight altitude of 3 km

[0104] Set the helicopter flight altitude H = 3km, forward flight speed V x From 20m / s to 50m / s, the output parameter change curve is as follows Figure 4 shown.

[0105] Depend on Figure 4 It can be seen that compared with the bidirectional difference algorithm, the speed deviation output by the component-level model based on the pseudo-Jacobian matrix at an altitude of 3km is less than 0.3%, and the temperature deviation is less than 0.4%. Since the fuel flow rate itself is of a small order of magnitude, the maximum deviation in the closed-loop control is less than 1.1%. The calculation deviation of the fuel consumption rate caused by the fuel flow deviation is also less than 1.1%. The maximum modeling error is better than that at an altitude of 1km, which verifies its applicability within the envelope.

[0106] References:

[0107] [1] Chen Yuchun, Xu Siyuan, Yang Yunkai, et al. Methods to improve the convergence of aeroengine performance calculations [J]. Journal of Aerospace Power, 2008, 23(12): 2242-2248.

[0108] [2] Liao Guanghuang, Jiao Yang, Li Qiuhong, et al. Research on high-precision real-time component-level model of turboshaft engine[J]. Propulsion Technology, 2016, 37(01): 25-33.

[0109] [3] Wang Yuan, Li Qiuhong, Huang Xianghua. Numerical calculation of aeroengine model based on self-correcting Broyden quasi-Newton method[J]. Journal of Aerospace Power, 2016, 31(1): 249-256.

[0110] [4] Wu Qian, Li Qiuhong, Shan Ruibin. Research on the solution method of aerodynamic and thermal system model of aircraft engines[J]. Computer Simulation, 2019, 36(01):76-81.

[0111] [5]Pang, S.; Li, Q.; Ni, B. Improved nonlinear MPC for aircraft gas turbine engine based on semi-alternative optimization strategy [J]. Aerosp.Sci.Technol.2021,118,106983.

Claims

1. A method for dynamic modeling of an aircraft engine based on a pseudo-Jacobian matrix, which obtains parameters of each component of the aircraft engine matching the current working state by constructing and solving the dynamic common working balance equation of the aircraft engine during the working process of the aircraft engine; characterized in that: When solving the dynamic co-working balance equation of an aircraft engine, the guessed value of the dynamic co-working balance equation is iteratively corrected by the following method so that the current dynamic co-working balance equation meets the convergence condition: Step 1, linearize the dynamic joint working equilibrium equation at time k; Step 2: Based on the following optimization objective function, calculate the optimal estimate of the Jacobian matrix J(k) of the linearized dynamic joint working equilibrium equation: Among them, ε(k) and ε(k-1) represent the residuals of the dynamic joint working balance equation at time k and time k-1 respectively, Δx(k)=x(k)-x(k-1) represents the difference between the guessed values ​​x(k) and x(k-1) of the dynamic joint working balance equation at time k and time k-1, μ is the penalty factor, The optimal estimate of the Jacobian matrix representing the linearized dynamic co-working equilibrium equation at time k-1; Step 3: Update the dynamic joint working balance equation guess value x(k+1) at time k+1 through iterative calculation: Among them, λ∈(0,1] is the iteration step size.

2. The method for dynamic modeling of an aircraft engine based on a pseudo-Jacobian matrix as claimed in claim 1, characterized in that: Best estimate Calculated by the following formula: Among them, η∈(0,1] is the step size factor, Δε(k)=ε(k)-ε(k-1) represents the difference between the residuals ε(k) and ε(k-1) of the dynamic joint working equilibrium equation at time k and time k-1, and the superscript "T" represents the matrix transpose, "|||| 2 ” represents the square of the vector 2 norm.

Images