Engine variable-dimension numerical simulation method and device

Abstract

The invention belongs to the technical field of engine simulation, and relates to an engine variable-dimension numerical simulation method and device, and the method comprises the steps: calculating the parameters of a contra-rotating propeller fan based on a contra-rotating propeller fan three-dimensional lift line model; determining the axial average speed of air flow at the inlet of the air inlet; determining section parameters at an inlet of the air inlet channel; calculating section parameters at an outlet of the air inlet; the air flow at the outlet of the air inlet channel is divided into two parts, one part of the air flow is distorted, the other part of the air flow is distortionless, and section parameters at the outlet of the low-pressure compressor are determined through mixing at the outlet; calculating inlet and outlet section parameters of each component according to the zero-dimensional model of each component of the engine; a balance equation is constructed according to circulation of airflow in all parts of the engine; constructing a balance equation according to the balance relation between the shaft power and the torque output by the free turbine and the power and the torque calculated by the contra-rotating propeller fan three-dimensional lift line model; the balance equation is iteratively solved to determine an initial guess value, and when iteration converges, engine performance is calculated. The method is high in simulation precision and high in speed.

Description

Technical Field

[0001] The present application belongs to the technical field of engine simulation, and in particular relates to a method and device for numerical simulation of engine variable dimensions. Background Art

[0002] Aircraft engine whole machine simulation technology is an important means in the process of aircraft engine conceptual design, research and development, manufacturing, testing and use. With the help of more accurate aircraft engine whole machine simulation tools, the design and test cycle can be greatly reduced and the cost can be reduced. However, the traditional zero-dimensional whole machine model of aircraft engines based on component characteristics is highly dependent on component characteristics, and as the engine configuration becomes more complex, the zero-dimensional simulation model has convergence problems in the application process. Quasi-three-dimensional or three-dimensional aircraft engine whole machine simulation requires a lot of computing resources and has slow calculation speed. At the same time, the large-scale complex regional flow field simulation with multi-field coupling makes it difficult to guarantee its convergence and stability. Therefore, it is necessary to develop a high-precision whole machine simulation method that can achieve fast calculation and good model stability and convergence.

[0003] As a new generation of efficient and energy-saving green power, the open rotor engine has been a hot topic in recent years. At present, the simulation method of the open rotor engine is mainly based on the component-level zero-dimensional model, and the simulation calculation accuracy is limited.

[0004] In the prior art, in terms of variable dimension simulation of engines, the Chinese invention patent with publication number CN114781153A simulated the fan components in three dimensions and proposed a whole machine variable dimension simulation performance simulation process control method. The Chinese invention patent with publication number CN118114595A simulated the compressor and turbine in three dimensions and provided a coupling calculation method for variable dimension simulation of combustion engines. However, there is no research on variable dimension simulation of open rotor engines.

[0005] In addition, for a gear-driven open rotor engine with a traction structure, due to the interference of the counter-rotating propeller fan wake vortex, a certain degree of intake distortion will occur at the engine inlet, affecting the stability of the compressor. However, in the prior art, there is no research on the intake distortion of an open rotor engine. Summary of the invention

[0006] In order to solve the above problems, the present application provides a method and device for numerical simulation of engine variable dimensions.

[0007] The first aspect of the present application provides an engine variable dimension numerical simulation method, which mainly includes: Step S1, given an initial guess value of the front propeller speed, the front propeller pitch angle, the rear propeller speed, and the rear propeller pitch angle of the counter-rotating propeller fan, the total thrust, the front propeller power, the rear propeller power, the front propeller torque, and the rear propeller torque of the counter-rotating propeller fan are calculated based on a three-dimensional lift line model of the counter-rotating propeller fan; Step S2, determining the average axial velocity of the airflow at the inlet of the air inlet behind the counter-rotating propeller fan according to the total thrust of the counter-rotating propeller fan; Step S3, determining the cross-sectional parameters at the inlet of the air intake duct according to the average axial velocity of the airflow at the inlet of the air intake duct; Step S4, calculating the cross-sectional parameters at the inlet duct outlet according to the zero-dimensional model of the inlet duct; Step S5, dividing the airflow at the outlet of the inlet duct into two streams, one stream flows to the distorted low-pressure compressor, and the other stream flows to the undistorted low-pressure compressor, determining the cross-sectional parameters at the inlet of the distorted low-pressure compressor according to the given pressure distortion index, and then determining the cross-sectional parameters at the outlet of the distorted low-pressure compressor according to the given initial guess value of the pressure ratio of the distorted low-pressure compressor and the initial guess value of the low-pressure compressor speed, and at the same time, determining the cross-sectional parameters at the outlet of the undistorted low-pressure compressor according to the given initial guess value of the pressure ratio of the undistorted low-pressure compressor and the initial guess value of the low-pressure compressor speed, mixing the airflow at the outlet of the distorted low-pressure compressor with the airflow at the outlet of the undistorted low-pressure compressor, calculating the cross-sectional parameters at the outlet of the low-pressure compressor through a mixing model, and constructing a first equilibrium equation according to the static pressure equilibrium conditions of the outlet of the distorted low-pressure compressor and the outlet of the undistorted low-pressure compressor; Step S6, given the initial guess value of the high-pressure compressor pressure ratio, the initial guess value of the high-pressure compressor speed, the initial guess value of the combustion chamber fuel supply, the initial guess value of the high-pressure turbine pressure ratio, the initial guess value of the low-pressure turbine pressure ratio, the initial guess value of the free turbine pressure ratio and the initial guess value of the free turbine speed, based on the zero-dimensional model of the high-pressure compressor, the combustion chamber, the high-pressure turbine, the low-pressure turbine, the free turbine, and the tail nozzle components connected in sequence behind the low-pressure compressor, calculate the inlet and outlet cross-sectional parameters of each component; Step S7, constructing the second to eighth balance equations according to the seven balance conditions of flow continuity between the high-pressure compressor and the low-pressure compressor, flow continuity between the combustion chamber and the high-pressure turbine, power balance on the high-pressure shaft, flow continuity between the high-pressure turbine and the low-pressure turbine, power balance on the low-pressure shaft, flow continuity between the low-pressure turbine and the free turbine, and flow continuity between the free turbine and the tail nozzle; Step S8, constructing a ninth balance equation according to the power balance relationship that needs to be satisfied between the shaft power output by the free turbine and the sum of the front propeller power and the rear propeller power calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; Step S9, calculating the sun shaft torque of the planetary gear reducer connected to the free turbine according to the shaft power and speed outputted by the free turbine, distributing the sun shaft torque to the front propeller fan and the rear propeller fan according to the connection relationship between the planetary gear reducer and the front propeller fan and the rear propeller fan, and constructing the tenth balance equation and the eleventh balance equation according to the relationship between the distributed torque and the front propeller torque and the rear propeller torque calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; Step S10, selecting three as control parameters from the fourteen independent variables that need to be given initial guess values, solving the common working equation group consisting of the eleven equilibrium equations constructed above by iterating the initial guess values ​​of the remaining eleven independent variables, and calculating the engine performance when the iteration converges; The cross-sectional parameters include one or more of flow rate, area, total temperature, total pressure and total enthalpy.

[0008] Preferably, step S2 further comprises: Step S21, determining the average induced axial velocity at the counter-rotating propeller according to the total thrust of the counter-rotating propeller; Step S22, calculating the average axial velocity of the airflow at the inlet of the air inlet according to the distance between the counter-rotating propeller fan and the inlet of the air inlet and the average axial induced velocity.

[0009] Preferably, in step S5, the outlet static pressures of the distorted low-pressure compressor and the non-distorted low-pressure compressor are calculated according to their respective inlet cross-sectional parameters, based on the assumption that the inlet and outlet areas of the low-pressure compressors remain unchanged, and by the inlet areas in their respective inlet cross-sectional parameters.

[0010] Preferably, in step S10, the three selected control parameters include the front propeller pitch angle, the rear propeller pitch angle and the fuel supply to the combustion chamber.

[0011] Preferably, in step S10, the common working equation group consisting of the eleven equilibrium equations constructed above is solved by Newton's iteration method, and the Jacobian matrix constructed by the Newton's iteration method is calculated in parallel by multi-threading.

[0012] The second aspect of the present application provides an engine variable dimension numerical simulation device, which mainly includes: A contra-rotating propeller fan parameter calculation module is used to calculate the total thrust, front propeller power, rear propeller power, front propeller torque and rear propeller torque of the contra-rotating propeller fan based on a three-dimensional lift line model of the contra-rotating propeller fan given an initial guess of the front propeller speed, the front propeller pitch angle, the rear propeller speed and the rear propeller pitch angle; The contra-rotating propeller fan and the air inlet aerodynamic coupling module is used to determine the axial average velocity of the airflow at the inlet inlet behind the contra-rotating propeller fan according to the total thrust of the contra-rotating propeller fan, and then determine the cross-sectional parameters at the inlet inlet according to the axial average velocity of the airflow at the inlet inlet; An intake duct parameter calculation module, used to calculate the cross-sectional parameters at the intake duct outlet according to the intake duct zero-dimensional model; A low-pressure compressor parameter calculation module is used to divide the airflow at the inlet outlet into two streams, one stream flows to the distorted low-pressure compressor, and the other stream flows to the undistorted low-pressure compressor, determine the cross-sectional parameters at the inlet of the distorted low-pressure compressor according to a given pressure distortion index, and then determine the cross-sectional parameters at the outlet of the distorted low-pressure compressor according to a given initial guess value of the pressure ratio of the distorted low-pressure compressor and an initial guess value of the low-pressure compressor speed, and at the same time, determine the cross-sectional parameters at the outlet of the undistorted low-pressure compressor according to a given initial guess value of the pressure ratio of the undistorted low-pressure compressor and an initial guess value of the low-pressure compressor speed, mix the airflow at the outlet of the distorted low-pressure compressor with the airflow at the outlet of the undistorted low-pressure compressor, calculate the cross-sectional parameters at the outlet of the low-pressure compressor through a mixing model, and construct a first balance equation according to the static pressure balance conditions of the outlet of the distorted low-pressure compressor and the outlet of the undistorted low-pressure compressor; The module for calculating parameters of other components is used to calculate the inlet and outlet cross-sectional parameters of each component based on the zero-dimensional model of the high-pressure compressor, combustion chamber, high-pressure turbine, low-pressure turbine, free turbine, and tail nozzle components connected in sequence behind the low-pressure compressor, given the initial guess values ​​of the high-pressure compressor pressure ratio, the initial guess value of the high-pressure compressor speed, the initial guess value of the combustion chamber fuel supply, the initial guess value of the high-pressure turbine pressure ratio, the initial guess value of the low-pressure turbine pressure ratio, and the initial guess value of the free turbine speed; The engine working together balance equation construction module is used to construct the second to eighth balance equations according to the seven balance conditions of flow continuity between the high-pressure compressor and the low-pressure compressor, flow continuity between the combustion chamber and the high-pressure turbine, power balance on the high-pressure shaft, flow continuity between the high-pressure turbine and the low-pressure turbine, power balance on the low-pressure shaft, flow continuity between the low-pressure turbine and the free turbine, and flow continuity between the free turbine and the tail nozzle; The gear reducer power transfer calculation module is used to construct the ninth balance equation according to the power balance relationship that needs to be satisfied between the shaft power output by the free turbine and the sum of the front propeller power and the rear propeller power calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; a gear reducer torque transmission calculation module, for calculating the sun shaft torque of the planetary gear reducer connected thereto according to the shaft power and speed output by the free turbine, distributing the sun shaft torque to the front propellers and the rear propellers according to the connection relationship between the planetary gear reducer and the front propellers and the rear propellers, and constructing the tenth balance equation and the eleventh balance equation according to the relationship between the distributed torque and the front propeller torque and the rear propeller torque calculated based on the three-dimensional lift line model of the counter-rotating propellers; A common working equation group solving module is used to select three as control parameters from the above fourteen independent variables that need to be given initial guess values, solve the common working equation group consisting of the above eleven equilibrium equations by iterating the initial guess values ​​of the remaining eleven independent variables, and calculate the engine performance when the iteration converges; The cross-sectional parameters include one or more of flow rate, area, total temperature, total pressure and total enthalpy.

[0013] Preferably, the counter-rotating propeller fan and air inlet aerodynamic coupling module comprises: A counter-rotating propeller axial average induced speed calculation unit, used to determine the counter-rotating propeller axial average induced speed at the counter-rotating propeller according to the total thrust of the counter-rotating propeller; The airflow axial average velocity calculation unit is used to calculate the axial average velocity of the airflow at the inlet of the air inlet according to the distance between the counter-rotating propeller fan and the inlet of the air inlet and the axial average induced velocity.

[0014] Preferably, in the low-pressure compressor parameter calculation module, the outlet static pressure of the distorted low-pressure compressor and the non-distorted low-pressure compressor is calculated according to their respective inlet cross-sectional parameters, based on the assumption that the inlet and outlet areas of the low-pressure compressor are constant, and the respective outlet static pressures are calculated from the inlet areas in the respective inlet cross-sectional parameters.

[0015] Preferably, in the module for solving the common working equation group, the three selected control parameters include the front propeller pitch angle, the rear propeller pitch angle and the fuel supply to the combustion chamber.

[0016] Preferably, in the common working equation group solving module, the common working equation group consisting of the eleven equilibrium equations constructed above is solved by Newton iteration method, and the Jacobian matrix constructed by Newton iteration method is calculated in parallel by multi-threading.

[0017] This application achieves high-fidelity calculation of the overall performance of an open rotor engine with high calculation accuracy and high speed. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 This is a flow chart of a preferred embodiment of the engine variable dimension numerical simulation method of the present application.

[0019] Figure 2 This is a schematic diagram of the position relationship and cross-sectional numbering of the components of a traction-type open rotor engine. DETAILED DESCRIPTION

[0020] In order to make the purpose, technical scheme and advantages of the implementation of this application clearer, the technical scheme in the implementation of this application will be described in more detail in combination with the drawings in the implementation of this application. In the drawings, the same or similar reference numerals throughout represent the same or similar elements or elements with the same or similar functions. The described implementation is a part of the implementation of this application, not all of the implementations. The implementation described below with reference to the drawings is exemplary and is intended to be used to explain this application, and cannot be understood as a limitation on this application. Based on the implementation in this application, all other implementations obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application. The implementation of this application is described in detail below in combination with the drawings.

[0021] The present application provides a method and device for numerical simulation of engine variable dimensions. Taking into account the shortcomings of the prior art in that the three-dimensional model of the whole machine has a large amount of calculation and the zero-dimensional model of the whole machine has an inaccurate calculation, the present application combines the three-dimensional model with the zero-dimensional model. In the following step S1, the counter-rotating propellers in the engine structure use a three-dimensional or quasi-three-dimensional lift line model, which can directly obtain the true characteristics of the components, avoid the dependence on the characteristic diagram, thereby avoiding the error caused by characteristic scaling and improving the simulation accuracy of the component. In addition to the counter-rotating propellers, other engine components, including the zero-dimensional model of the air inlet introduced in step S4, the zero-dimensional model of the distorted low-pressure compressor and the zero-dimensional model of the undistorted low-pressure compressor introduced in step S5, and the zero-dimensional models of the high-pressure compressor, combustion chamber, high-pressure turbine, low-pressure turbine, free turbine, and tail nozzle components introduced in step S6, still use the zero-dimensional or quasi-one-dimensional models of the prior art, while in the following steps S2-S3, aerodynamic coupling is performed between the three-dimensional lift line model of the counter-rotating propeller and the zero-dimensional model of the engine air inlet, and in the following steps S8-S9, structural coupling is performed between the zero-dimensional model of the engine free turbine and the three-dimensional lift line model of the counter-rotating propeller. These two coupling methods realize the variable-dimensional numerical simulation of the entire open rotor engine. The following is a detailed description.

[0022] The first aspect of the present application provides a variable dimension numerical simulation method for an engine, such as Figure 1 As shown, it mainly includes: Step S1, given an initial guess value of the front propeller speed, the front propeller pitch angle, the rear propeller speed, and the rear propeller pitch angle of the counter-rotating propeller fan, the total thrust, the front propeller power, the rear propeller power, the front propeller torque, and the rear propeller torque of the counter-rotating propeller fan are calculated based on a three-dimensional lift line model of the counter-rotating propeller fan; Step S2, determining the average axial velocity of the airflow at the inlet of the air inlet behind the counter-rotating propeller fan according to the total thrust of the counter-rotating propeller fan; Step S3, determining the cross-sectional parameters at the inlet of the air intake duct according to the average axial velocity of the airflow at the inlet of the air intake duct; Step S4, calculating the cross-sectional parameters at the inlet duct outlet according to the zero-dimensional model of the inlet duct; Step S5, dividing the airflow at the outlet of the inlet duct into two streams, one stream flows to the distorted low-pressure compressor, and the other stream flows to the undistorted low-pressure compressor, determining the cross-sectional parameters at the inlet of the distorted low-pressure compressor according to the given pressure distortion index, and then determining the cross-sectional parameters at the outlet of the distorted low-pressure compressor according to the given initial guess value of the pressure ratio of the distorted low-pressure compressor and the initial guess value of the low-pressure compressor speed, and at the same time, determining the cross-sectional parameters at the outlet of the undistorted low-pressure compressor according to the given initial guess value of the pressure ratio of the undistorted low-pressure compressor and the initial guess value of the low-pressure compressor speed, mixing the airflow at the outlet of the distorted low-pressure compressor with the airflow at the outlet of the undistorted low-pressure compressor, calculating the cross-sectional parameters at the outlet of the low-pressure compressor through a mixing model, and constructing a first equilibrium equation according to the static pressure equilibrium conditions of the outlet of the distorted low-pressure compressor and the outlet of the undistorted low-pressure compressor; Step S6, given the initial guess value of the high-pressure compressor pressure ratio, the initial guess value of the high-pressure compressor speed, the initial guess value of the combustion chamber fuel supply, the initial guess value of the high-pressure turbine pressure ratio, the initial guess value of the low-pressure turbine pressure ratio, the initial guess value of the free turbine pressure ratio and the initial guess value of the free turbine speed, based on the zero-dimensional model of the high-pressure compressor, the combustion chamber, the high-pressure turbine, the low-pressure turbine, the free turbine, and the tail nozzle components connected in sequence behind the low-pressure compressor, calculate the inlet and outlet cross-sectional parameters of each component; Step S7, constructing the second to eighth balance equations according to the seven balance conditions of flow continuity between the high-pressure compressor and the low-pressure compressor, flow continuity between the combustion chamber and the high-pressure turbine, power balance on the high-pressure shaft, flow continuity between the high-pressure turbine and the low-pressure turbine, power balance on the low-pressure shaft, flow continuity between the low-pressure turbine and the free turbine, and flow continuity between the free turbine and the tail nozzle; Step S8, constructing a ninth balance equation according to the power balance relationship that needs to be satisfied between the shaft power output by the free turbine and the sum of the front propeller power and the rear propeller power calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; Step S9, calculating the sun shaft torque of the planetary gear reducer connected to the free turbine according to the shaft power and speed outputted by the free turbine, distributing the sun shaft torque to the front propeller fan and the rear propeller fan according to the connection relationship between the planetary gear reducer and the front propeller fan and the rear propeller fan, and constructing the tenth balance equation and the eleventh balance equation according to the relationship between the distributed torque and the front propeller torque and the rear propeller torque calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; Step S10, selecting three as control parameters from the fourteen independent variables that need to be given initial guess values, solving the common working equation group consisting of the eleven equilibrium equations constructed above by iterating the initial guess values ​​of the remaining eleven independent variables, and calculating the engine performance when the iteration converges; The cross-sectional parameters include one or more of flow rate, area, total temperature, total pressure and total enthalpy.

[0023] First, before step S1, the airflow velocity at the inlet of the counter-rotating propeller fan is calculated according to the given atmospheric conditions (altitude, flight Mach number, temperature, etc.). ,temperature ,pressure Then, in step S1, given the initial guess values ​​of the front propeller speed n1, the front propeller pitch angle β1, the rear propeller speed n2 and the rear propeller pitch angle β2, the total thrust F is calculated through the established three-dimensional lift line model of the counter-rotating propeller fan. CRP , total power P CRP , efficiency η CRP , front propeller power P1, rear propeller power P2, front propeller torque Q1 and rear propeller torque Q2.

[0024] Then, in step S2, the average axial velocity of the airflow at the inlet of the air duct is calculated.

[0025] In some optional implementations, step S2 further includes: Step S21, determining the average induced axial velocity at the counter-rotating propeller according to the total thrust of the counter-rotating propeller; Step S22, calculating the average axial velocity of the airflow at the inlet of the air inlet according to the distance between the counter-rotating propeller fan and the inlet of the air inlet and the average axial induced velocity.

[0026] In step S21, the axial average induced velocity of the counter-rotating propeller position is calculated according to the momentum theorem: : .

[0027] Where ρ is the atmospheric density, A disk is the fan area of ​​the counter-rotating propeller fan.

[0028] Then, in step S22, the average axial velocity of the airflow at the inlet position of the engine inlet under the interference of the counter-rotating propeller is calculated according to the flow channel contraction principle. : ; Wherein, z refers to the distance between the downstream of the counter-rotating propeller fan and the inlet of the air inlet, and R is the fan radius of the counter-rotating propeller fan.

[0029] At the same time, the three-dimensional distribution of the airflow velocity at the inlet section of the air inlet can also be obtained based on the three-dimensional lift line model of the counter-rotating propeller fan.

[0030] Then in step S3, the average velocity at the inlet of the air duct is , the inlet cross-sectional parameters after considering the interference of the counter-rotating propellers can be obtained. The cross-sectional parameters include one or more of flow rate, area, total temperature, total pressure and total enthalpy. Here, the total temperature T of the inlet can be mainly included. t1 , total pressure p t1 , total enthalpy h t1 . refer to Figure 2 , the front section of the inlet is numbered 1, so the inlet section parameters are usually represented by the subscript t1. For example, the low-pressure compressor outlet section is numbered 25, and the total enthalpy of the low-pressure compressor outlet is h t25 , the total outlet temperature is T t25 .

[0031] With the above inlet duct inlet cross-sectional parameters, refer to Figure 2 , the existing zero-dimensional model of each engine component can be used for simulation calculation later, wherein, in step S4, the average total temperature of the inlet duct outlet is calculated by the zero-dimensional model of the inlet duct. , average total pressure , average total enthalpy .

[0032] Afterwards, in step S5, due to the influence of the contra-rotating propellers' wake, the actual inlet airflow is uneven, and there is an intake distortion problem. Therefore, the present application considers the intake distortion problem by establishing a parallel compressor model for the low-pressure compressor.

[0033] In some optional implementations, in step S5, the outlet static pressures of the distorted low-pressure compressor and the non-distorted low-pressure compressor are calculated based on their respective inlet cross-sectional parameters and the assumption that the inlet and outlet areas of the low-pressure compressors are constant, and the inlet areas in the inlet cross-sectional parameters are used to calculate their respective outlet static pressures. The parallel compressor model is described in detail below.

[0034] (1) The airflow from the inlet outlet is divided into two streams, one stream flows to the distorted low-pressure compressor a, and the other stream flows to the undistorted low-pressure compressor b. The distorted low-pressure compressor a and the undistorted low-pressure compressor b only have different inlet conditions, but the speed is the same. The two compressors work in parallel and then mix at the outlet. The pressure distortion index of a given distortion sector is and temperature distortion index , so the total inlet temperature of the distorted low-pressure compressor a can be calculated and total pressure .

[0035] (2) For the low-pressure compressor a, based on the given distorted low-pressure pressure ratio Z LCa The initial guess value (the "pressure ratio" here is a variable z, and its expression is z=(π-π min ) / (π max -πmin ), where π is the pressure ratio, π max is the maximum value of the pressure ratio on the equal speed characteristic line of the characteristic diagram, π min is the minimum value of the pressure ratio on the characteristic line of the characteristic diagram) and the low-pressure compressor speed n LC The initial guess value of the low-pressure compressor flow W is obtained by interpolating the characteristic diagram 2a,cor and efficiency η LCa , using the zero-dimensional model of the distorted low-pressure compressor, the outlet flow rate W of the distorted low-pressure compressor a can be further calculated 2a , total temperature , total pressure , total enthalpy and power P LCa By W 2a and The inlet and outlet cross-sectional area A of the distorted low-pressure compressor a can be calculated 2a Assuming that the inlet and outlet areas of the low-pressure compressor remain unchanged, according to the area A 2a And the total parameters of the outlet of the distorted low-pressure compressor a, the outlet static pressure of the distorted low-pressure compressor a can be iteratively calculated .

[0036] (3) For the distortion-free low-pressure compressor b, based on the given distortion-free low-pressure pressure ratio Z LCb The initial guess value and the low-pressure compressor speed n given in the previous step LC Initial guess, through the characteristic diagram interpolation to obtain the converted flow rate W of the low pressure compressor 2b,cor and efficiency η LCb Using the zero-dimensional model of the distortion-free low-pressure compressor, the outlet flow rate W of the distortion-free low-pressure compressor b can be further calculated. 2b , total temperature , total pressure , total enthalpy and power P LCb By W 2b and The inlet and outlet cross-sectional area A of the distortion-free low-pressure compressor b can be calculated 2b Assuming that the inlet and outlet areas of the distortion-free low-pressure compressor b remain unchanged, according to the area A 2b The total outlet parameters of the distortion-free low-pressure compressor b can be used to iteratively calculate the outlet static pressure of the distortion-free low-pressure compressor b. .

[0037] (4) Assume that the two air flows are mixed after coming out, and the outlet flow rate of the low-pressure compressor after mixing is W 25 =W 2a +W 2b , the total outlet enthalpy is h t25 =(W 2a h t25,a+W 2b h t25,b ) / W 25 , by h t25 Iterative calculation of the total temperature T at the outlet of the low-pressure compressor t25 By P LC =P LCa +P LCb Calculate the total power P of the low pressure compressor LC .

[0038] The mixing process of the two airflows needs to meet the static pressure balance condition at the low-pressure compressor outlet section, namely: , so we get the equilibrium equation that needs to be satisfied: y1=(p s25,a -p s25,b ) / p s25,b .

[0039] Then, in step S6, refer to Figure 1 and Figure 2 The mixed airflow enters the high-pressure compressor, combustion chamber, high-pressure turbine, low-pressure turbine, free turbine, and tail nozzle in turn. These components all use the traditional zero-dimensional model. By transferring the inlet and outlet parameters of each component between different components, a quasi-one-dimensional simulation model of the engine can be obtained.

[0040] In the quasi-one-dimensional engine model from the inlet to the tail nozzle, in addition to the three initial guesses required for the parallel low-pressure compressor model in step S5 (the pressure ratio of the undistorted low pressure, the ratio Z LCb , the pressure ratio of the distorted low pressure is Z LCa and low pressure compressor speed n LC ), the remaining components also need to give the following 7 initial guess values: High pressure ratio ratio Z HC , high pressure compressor speed n HC , fuel supply to combustion chamber W fb 、High-pressure turbine pressure drop ratio π HT 、Low-pressure turbine pressure drop ratio π LT 、Free turbine pressure drop ratio π FT , free turbine speed n FT .

[0041] In step S7, according to the balance conditions of the engine (continuous flow, equal power, equal speed, and balanced pressure), the following balance equations need to be satisfied in the quasi-one-dimensional engine model from the air inlet to the tail nozzle: .

[0042] In the formula, In the parallel compressor model of the low-pressure compressor, the outlet static pressure balance parameter of the two parallel low-pressure compressors a and b is 0, indicating balance. It is the flow continuity parameter between the high-pressure compressor and the low-pressure compressor. When it is 0, it means the flow is continuous. The following is equivalent and will not be repeated. is the low pressure compressor outlet flow rate, is the high pressure compressor inlet flow rate, is the flow continuity parameter between the combustion chamber and the high-pressure turbine, is the combustion chamber outlet flow rate, is the high pressure turbine inlet flow rate, is the power balance parameter on the high-voltage shaft, is the high pressure turbine output power, is the high pressure compressor power, High voltage power transfer efficiency, is the flow continuity parameter between the high-pressure turbine and the low-pressure turbine, is the high pressure turbine outlet flow rate, is the low-pressure turbine inlet flow rate, is the power balance parameter on the low pressure axis, is the output power of the low-pressure turbine, is the low voltage power transfer efficiency, is the total power of the low-pressure compressor, is the flow continuity parameter between the low-pressure turbine and the free turbine, is the low-pressure turbine outlet flow rate, is the free turbine inlet flow, is the flow continuity parameter between the free turbine and the tail nozzle, is the free turbine outlet flow, Tail nozzle inlet flow rate.

[0043] Then in step S9, refer to Figure 2 , the shaft power P output from the free turbine FT and speed n FT , transmitted to the sun shaft of the planetary gear reducer through the mechanical shaft. At this time, the torque of the sun shaft is Q sun . Shaft power P FT The power is distributed by a planetary gear reducer, the front propellers are driven by the planetary carrier of the reducer, and the rear propellers are driven by the outer ring gear of the reducer.

[0044] The radius ratio R of the planetary gear and sun gear of the planetary gear reducer is given by p / R s The torque Q of the planet carrier can be calculated according to the following formula star The torque of the outer ring gear is Q ring : .

[0045] in, For transmission efficiency, the power transmitted to the front propeller fan through the planetary gear reducer is P1, and the power transmitted to the rear propeller fan is P2. The shaft power output from the free turbine is P FT , the total power P calculated in step S1 is required CRP Equal, the planetary carrier torque Q of the planetary gear reducer star It needs to be equal to the front propeller torque Q1 of the front propeller fan, and the outer ring gear torque Q ring It needs to be equal to the rear propeller torque Q2 of the rear propeller fan. Therefore, the following three balance conditions need to be met between the free turbine, planetary gear reducer and counter-rotating propeller fan: .

[0046] Where ΔP CRP is the balance parameter of the output power of the counter-rotating propeller fan and the free turbine, ΔQ1 is the torque balance parameter of the front propeller and the planetary gear carrier; ΔQ2 is the torque balance parameter of the rear propeller and the annular gear ring.

[0047] Finally, in step S10, according to the previous steps, the open rotor engine has 14 independent variables that need to be given initial guess values, namely: β1, n1, β2, n2, Z LCa , Z LCb , n LC , Z HC , n HC , W fb , π HT , π LT , π FT , n FT However, the engine has 11 equilibrium equations. In order to close the equations, it is necessary to select 3 variables from the above independent variables as control parameters through the control strategy, and the remaining 11 independent variables are used as free variables x required by the engine's common working equation group y.

[0048] In some optional embodiments, in step S10, the three control parameters selected include the front propeller pitch angle β1, the rear propeller pitch angle β2 and the combustion chamber fuel supply W fb .

[0049] Based on this, the common working equations of the open rotor engine considering intake distortion can be expressed as: .

[0050] In some optional implementations, in step S10, the common working equation group consisting of the eleven equilibrium equations constructed above is solved by Newton's iteration method, and the Jacobian matrix constructed by the Newton's iteration method is calculated in parallel by multi-threading.

[0051] The Newton iteration method is used to solve the 11-dimensional engine common working nonlinear equations of the open rotor engine variable dimension simulation process in step 10, so that y=y(x)=0. Take any initial value x0 of x and construct the following iterative format for iterative calculation: ; Where k is the number of iterations, J -1 is the inverse matrix of the Jacobian matrix J, and its specific expression is as follows: .

[0052] When the residual vector of the common working equations ||y||<10 -4 , the calculation is considered to converge, and the working state of the engine can be obtained, otherwise the value of the independent variable x is updated through the Jacobian matrix J. When the program converges, the performance of the engine is calculated, including thrust, fuel consumption, power, propulsion efficiency and other performance parameters.

[0053] In the variable-dimensional numerical simulation program for open rotor engines, since the 11-dimensional Jacobian matrix needs to be calculated, the quasi-3D lift line model of the counter-rotating propeller fan needs to be called 12 times for each iteration, 11 of which are for calculating the 11-dimensional Jacobian matrix and 1 for judging the convergence state. Since the counter-rotating propeller fan module takes about 20 seconds to calculate a working point, if a single thread is used for calculation, each iterative calculation takes 4 minutes. If the iteration converges after 5 times, it will take at least 20 minutes to obtain the performance of a non-design point, which cannot meet the purpose of rapid design and calculation.

[0054] To this end, this application adopts a multi-threaded parallel computing method when calculating the Jacobian matrix. Because the partial derivatives of the Jacobian matrix are independent of each other, by using 11 threads to simultaneously carry out parallel calculations of small perturbations of 11 independent variables, one iteration only takes 40 seconds, and 5 convergences only take 3.3 minutes to obtain a non-design point performance, which greatly saves calculation time.

[0055] This application can realize the rapid calculation of engine performance. The calculation time of the design point is about 20s. For the calculation process of the non-design point, the convergence steps of Newton iteration are 4-6 steps and the calculation time is about 2.5min-4min. The program has good convergence and stability.

[0056] The second aspect of the present application provides an engine variable dimension numerical simulation device corresponding to the above method, mainly comprising: A contra-rotating propeller fan parameter calculation module is used to calculate the total thrust, front propeller power, rear propeller power, front propeller torque and rear propeller torque of the contra-rotating propeller fan based on a three-dimensional lift line model of the contra-rotating propeller fan given an initial guess of the front propeller speed, the front propeller pitch angle, the rear propeller speed and the rear propeller pitch angle; The contra-rotating propeller fan and the air inlet aerodynamic coupling module is used to determine the axial average velocity of the airflow at the inlet inlet behind the contra-rotating propeller fan according to the total thrust of the contra-rotating propeller fan, and then determine the cross-sectional parameters at the inlet inlet according to the axial average velocity of the airflow at the inlet inlet; An intake duct parameter calculation module, used to calculate the cross-sectional parameters at the intake duct outlet according to the intake duct zero-dimensional model; A low-pressure compressor parameter calculation module is used to divide the airflow at the inlet outlet into two streams, one stream flows to the distorted low-pressure compressor, and the other stream flows to the undistorted low-pressure compressor, determine the cross-sectional parameters at the inlet of the distorted low-pressure compressor according to a given pressure distortion index, and then determine the cross-sectional parameters at the outlet of the distorted low-pressure compressor according to a given initial guess value of the pressure ratio of the distorted low-pressure compressor and an initial guess value of the low-pressure compressor speed, and at the same time, determine the cross-sectional parameters at the outlet of the undistorted low-pressure compressor according to a given initial guess value of the pressure ratio of the undistorted low-pressure compressor and an initial guess value of the low-pressure compressor speed, mix the airflow at the outlet of the distorted low-pressure compressor with the airflow at the outlet of the undistorted low-pressure compressor, calculate the cross-sectional parameters at the outlet of the low-pressure compressor through a mixing model, and construct a first balance equation according to the static pressure balance conditions of the outlet of the distorted low-pressure compressor and the outlet of the undistorted low-pressure compressor; The module for calculating parameters of other components is used to calculate the inlet and outlet cross-sectional parameters of each component based on the zero-dimensional model of the high-pressure compressor, combustion chamber, high-pressure turbine, low-pressure turbine, free turbine, and tail nozzle components connected in sequence behind the low-pressure compressor, given the initial guess values ​​of the high-pressure compressor pressure ratio, the initial guess value of the high-pressure compressor speed, the initial guess value of the combustion chamber fuel supply, the initial guess value of the high-pressure turbine pressure ratio, the initial guess value of the low-pressure turbine pressure ratio, and the initial guess value of the free turbine speed; The engine working together balance equation construction module is used to construct the second to eighth balance equations according to the seven balance conditions of flow continuity between the high-pressure compressor and the low-pressure compressor, flow continuity between the combustion chamber and the high-pressure turbine, power balance on the high-pressure shaft, flow continuity between the high-pressure turbine and the low-pressure turbine, power balance on the low-pressure shaft, flow continuity between the low-pressure turbine and the free turbine, and flow continuity between the free turbine and the tail nozzle; The gear reducer power transfer calculation module is used to construct the ninth balance equation according to the power balance relationship that needs to be satisfied between the shaft power output by the free turbine and the sum of the front propeller power and the rear propeller power calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; a gear reducer torque transmission calculation module, for calculating the sun shaft torque of the planetary gear reducer connected thereto according to the shaft power and speed output by the free turbine, distributing the sun shaft torque to the front propellers and the rear propellers according to the connection relationship between the planetary gear reducer and the front propellers and the rear propellers, and constructing the tenth balance equation and the eleventh balance equation according to the relationship between the distributed torque and the front propeller torque and the rear propeller torque calculated based on the three-dimensional lift line model of the counter-rotating propellers; A common working equation group solving module is used to select three as control parameters from the above fourteen independent variables that need to be given initial guess values, solve the common working equation group consisting of the above eleven equilibrium equations by iterating the initial guess values ​​of the remaining eleven independent variables, and calculate the engine performance when the iteration converges; The cross-sectional parameters include one or more of flow rate, area, total temperature, total pressure and total enthalpy.

[0057] In some optional embodiments, the counter-rotating propeller fan and the air inlet aerodynamic coupling module includes: A counter-rotating propeller axial average induced speed calculation unit, used to determine the counter-rotating propeller axial average induced speed at the counter-rotating propeller according to the total thrust of the counter-rotating propeller; The airflow axial average velocity calculation unit is used to calculate the axial average velocity of the airflow at the inlet of the air inlet according to the distance between the counter-rotating propeller fan and the inlet of the air inlet and the axial average induced velocity.

[0058] In some optional embodiments, in the low-pressure compressor parameter calculation module, the outlet static pressure of the distorted low-pressure compressor and the non-distorted low-pressure compressor is calculated based on their respective inlet cross-sectional parameters and based on the assumption that the inlet and outlet areas of the low-pressure compressor are constant, and the respective outlet static pressures are calculated from the inlet areas in the respective inlet cross-sectional parameters.

[0059] In some optional embodiments, in the module for solving the common working equation group, the three control parameters selected include the front propeller pitch angle, the rear propeller pitch angle and the fuel supply to the combustion chamber.

[0060] In some optional implementations, in the common working equation group solving module, the common working equation group consisting of the eleven equilibrium equations constructed above is solved by Newton's iteration method, and the Jacobian matrix constructed by the Newton's iteration method is calculated in parallel by multi-threading.

[0061] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application shall be based on the protection scope of the claims.

Claims

1. A variable dimension numerical simulation method for an engine, characterized in that: include: Step S1, given an initial guess value of the front propeller speed, the front propeller pitch angle, the rear propeller speed, and the rear propeller pitch angle of the counter-rotating propeller fan, the total thrust, the front propeller power, the rear propeller power, the front propeller torque, and the rear propeller torque of the counter-rotating propeller fan are calculated based on a three-dimensional lift line model of the counter-rotating propeller fan; Step S2, determining the average axial velocity of the airflow at the inlet of the air inlet behind the counter-rotating propeller fan according to the total thrust of the counter-rotating propeller fan; Step S3, determining the cross-sectional parameters at the inlet of the air intake duct according to the average axial velocity of the airflow at the inlet of the air intake duct; Step S4, calculating the cross-sectional parameters at the inlet duct outlet according to the zero-dimensional model of the inlet duct; Step S5, dividing the airflow at the outlet of the inlet duct into two streams, one stream flows to the distorted low-pressure compressor, and the other stream flows to the undistorted low-pressure compressor, determining the cross-sectional parameters at the inlet of the distorted low-pressure compressor according to the given pressure distortion index, and then determining the cross-sectional parameters at the outlet of the distorted low-pressure compressor according to the given initial guess value of the pressure ratio of the distorted low-pressure compressor and the initial guess value of the low-pressure compressor speed, and at the same time, determining the cross-sectional parameters at the outlet of the undistorted low-pressure compressor according to the given initial guess value of the pressure ratio of the undistorted low-pressure compressor and the initial guess value of the low-pressure compressor speed, mixing the airflow at the outlet of the distorted low-pressure compressor with the airflow at the outlet of the undistorted low-pressure compressor, calculating the cross-sectional parameters at the outlet of the low-pressure compressor through a mixing model, and constructing a first equilibrium equation according to the static pressure equilibrium conditions of the outlet of the distorted low-pressure compressor and the outlet of the undistorted low-pressure compressor; Step S6, given the initial guess value of the high-pressure compressor pressure ratio, the initial guess value of the high-pressure compressor speed, the initial guess value of the combustion chamber fuel supply, the initial guess value of the high-pressure turbine pressure ratio, the initial guess value of the low-pressure turbine pressure ratio, the initial guess value of the free turbine pressure ratio and the initial guess value of the free turbine speed, based on the zero-dimensional model of the high-pressure compressor, the combustion chamber, the high-pressure turbine, the low-pressure turbine, the free turbine, and the tail nozzle components connected in sequence behind the low-pressure compressor, calculate the inlet and outlet cross-sectional parameters of each component; Step S7, constructing the second to eighth balance equations according to the seven balance conditions of flow continuity between the high-pressure compressor and the low-pressure compressor, flow continuity between the combustion chamber and the high-pressure turbine, power balance on the high-pressure shaft, flow continuity between the high-pressure turbine and the low-pressure turbine, power balance on the low-pressure shaft, flow continuity between the low-pressure turbine and the free turbine, and flow continuity between the free turbine and the tail nozzle; Step S8, constructing a ninth balance equation according to the power balance relationship that needs to be satisfied between the shaft power output by the free turbine and the sum of the front propeller power and the rear propeller power calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; Step S9, calculating the sun shaft torque of the planetary gear reducer connected to the free turbine according to the shaft power and speed outputted by the free turbine, distributing the sun shaft torque to the front propeller fan and the rear propeller fan according to the connection relationship between the planetary gear reducer and the front propeller fan and the rear propeller fan, and constructing the tenth balance equation and the eleventh balance equation according to the relationship between the distributed torque and the front propeller torque and the rear propeller torque calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; Step S10, selecting three as control parameters from the fourteen independent variables that need to be given initial guess values, solving the common working equation group consisting of the eleven equilibrium equations constructed above by iterating the initial guess values ​​of the remaining eleven independent variables, and calculating the engine performance when the iteration converges; The cross-sectional parameters include one or more of flow rate, area, total temperature, total pressure and total enthalpy.

2. The engine variable dimension numerical simulation method according to claim 1, characterized in that: Step S2 further comprises: Step S21, determining the average induced axial velocity at the counter-rotating propeller according to the total thrust of the counter-rotating propeller; Step S22, calculating the average axial velocity of the airflow at the inlet of the air inlet according to the distance between the counter-rotating propeller fan and the inlet of the air inlet and the average axial induced velocity.

3. The engine variable dimension numerical simulation method according to claim 1, characterized in that: In step S5, the outlet static pressures of the distorted low-pressure compressor and the non-distorted low-pressure compressor are calculated according to their respective inlet cross-sectional parameters and based on the assumption that the inlet and outlet areas of the low-pressure compressors remain unchanged, and according to the inlet areas in their respective inlet cross-sectional parameters.

4. The engine variable dimension numerical simulation method according to claim 1, characterized in that: In step S10, the three control parameters selected include the front propeller pitch angle, the rear propeller pitch angle and the fuel supply to the combustion chamber.

5. The engine variable dimension numerical simulation method according to claim 1, characterized in that: In step S10, the common working equation group consisting of the eleven equilibrium equations constructed above is solved by Newton's iteration method, and the Jacobian matrix constructed by the Newton's iteration method is calculated in parallel by multi-threading.

6. An engine variable dimension numerical simulation device, characterized in that: include: A contra-rotating propeller fan parameter calculation module is used to calculate the total thrust, front propeller power, rear propeller power, front propeller torque and rear propeller torque of the contra-rotating propeller fan based on a three-dimensional lift line model of the contra-rotating propeller fan given an initial guess of the front propeller speed, the front propeller pitch angle, the rear propeller speed and the rear propeller pitch angle; The contra-rotating propeller fan and the air inlet aerodynamic coupling module is used to determine the axial average velocity of the airflow at the inlet inlet behind the contra-rotating propeller fan according to the total thrust of the contra-rotating propeller fan, and then determine the cross-sectional parameters at the inlet inlet according to the axial average velocity of the airflow at the inlet inlet; An intake duct parameter calculation module, used to calculate the cross-sectional parameters at the intake duct outlet according to the intake duct zero-dimensional model; A low-pressure compressor parameter calculation module is used to divide the airflow at the inlet outlet into two streams, one stream flows to the distorted low-pressure compressor, and the other stream flows to the undistorted low-pressure compressor, determine the cross-sectional parameters at the inlet of the distorted low-pressure compressor according to a given pressure distortion index, and then determine the cross-sectional parameters at the outlet of the distorted low-pressure compressor according to a given initial guess value of the pressure ratio of the distorted low-pressure compressor and an initial guess value of the low-pressure compressor speed, and at the same time, determine the cross-sectional parameters at the outlet of the undistorted low-pressure compressor according to a given initial guess value of the pressure ratio of the undistorted low-pressure compressor and an initial guess value of the low-pressure compressor speed, mix the airflow at the outlet of the distorted low-pressure compressor with the airflow at the outlet of the undistorted low-pressure compressor, calculate the cross-sectional parameters at the outlet of the low-pressure compressor through a mixing model, and construct a first balance equation according to the static pressure balance conditions of the outlet of the distorted low-pressure compressor and the outlet of the undistorted low-pressure compressor; The module for calculating parameters of other components is used to calculate the inlet and outlet cross-sectional parameters of each component based on the zero-dimensional model of the high-pressure compressor, combustion chamber, high-pressure turbine, low-pressure turbine, free turbine, and tail nozzle components connected in sequence behind the low-pressure compressor, given the initial guess values ​​of the high-pressure compressor pressure ratio, the initial guess value of the high-pressure compressor speed, the initial guess value of the combustion chamber fuel supply, the initial guess value of the high-pressure turbine pressure ratio, the initial guess value of the low-pressure turbine pressure ratio, and the initial guess value of the free turbine speed; The engine working together balance equation construction module is used to construct the second to eighth balance equations according to the seven balance conditions of flow continuity between the high-pressure compressor and the low-pressure compressor, flow continuity between the combustion chamber and the high-pressure turbine, power balance on the high-pressure shaft, flow continuity between the high-pressure turbine and the low-pressure turbine, power balance on the low-pressure shaft, flow continuity between the low-pressure turbine and the free turbine, and flow continuity between the free turbine and the tail nozzle; The gear reducer power transfer calculation module is used to construct the ninth balance equation according to the power balance relationship that needs to be satisfied between the shaft power output by the free turbine and the sum of the front propeller power and the rear propeller power calculated based on the three-dimensional lift line model of the counter-rotating propeller fan; a gear reducer torque transmission calculation module, for calculating the sun shaft torque of the planetary gear reducer connected thereto according to the shaft power and speed output by the free turbine, distributing the sun shaft torque to the front propellers and the rear propellers according to the connection relationship between the planetary gear reducer and the front propellers and the rear propellers, and constructing the tenth balance equation and the eleventh balance equation according to the relationship between the distributed torque and the front propeller torque and the rear propeller torque calculated based on the three-dimensional lift line model of the counter-rotating propellers; A common working equation group solving module is used to select three as control parameters from the above fourteen independent variables that need to be given initial guess values, solve the common working equation group consisting of the above eleven equilibrium equations by iterating the initial guess values ​​of the remaining eleven independent variables, and calculate the engine performance when the iteration converges; The cross-sectional parameters include one or more of flow rate, area, total temperature, total pressure and total enthalpy.

7. The engine variable dimension numerical simulation device according to claim 6, characterized in that: The counter-rotating propeller fan and the air inlet aerodynamic coupling module comprises: A counter-rotating propeller axial average induced speed calculation unit, used to determine the counter-rotating propeller axial average induced speed at the counter-rotating propeller according to the total thrust of the counter-rotating propeller; The airflow axial average velocity calculation unit is used to calculate the axial average velocity of the airflow at the inlet of the air inlet according to the distance between the counter-rotating propeller fan and the inlet of the air inlet and the axial average induced velocity.

8. The engine variable dimension numerical simulation device according to claim 6, characterized in that: In the low-pressure compressor parameter calculation module, the outlet static pressure of the distorted low-pressure compressor and the non-distorted low-pressure compressor is calculated based on their respective inlet cross-sectional parameters and the assumption that the inlet and outlet areas of the low-pressure compressor remain unchanged, and the respective outlet static pressures are calculated based on the inlet areas in their respective inlet cross-sectional parameters.

9. The engine variable dimension numerical simulation device according to claim 6, characterized in that: In the joint working equation group solving module, the three selected control parameters include the front propeller pitch angle, the rear propeller pitch angle and the combustion chamber fuel supply.

10. The engine variable dimension numerical simulation device according to claim 6, characterized in that: In the common working equation group solving module, the common working equation group consisting of the eleven equilibrium equations constructed above is solved by Newton iteration method, and the Jacobian matrix constructed by Newton iteration method is calculated in parallel by multi-threading.

Images