A method for predicting flow field parameter distribution of rotating detonation engine based on space propulsion algorithm

Abstract

The present invention discloses a method for predicting the distribution of flow field parameters of a rotating detonation engine based on a space propulsion algorithm, comprising the steps of: estimating the flow parameters before the detonation wave; determining the flow parameters after the detonation wave according to the CJ detonation model of a single specific heat ratio; assigning the flow parameters after the detonation wave and the far downstream flow parameters to the initial value line; determining the flow parameters in the shock wave coordinate system by the space propulsion algorithm; determining the starting position and the ending position of the injection of the combustible mixture; iteratively converging the flow field parameters of the rotating detonation engine, and obtaining the flow parameters in the laboratory coordinate system by coordinate transformation. The present invention couples the space propulsion algorithm and the CJ detonation wave model, obtains the flow parameters in the laboratory coordinate system by coordinate transformation, can quickly obtain the flow field structure and outlet flow parameter distribution of the rotating detonation engine, overcomes the shortcomings of the existing algorithm that it takes a long time to obtain the flow field of the rotating detonation engine, and reduces the design cycle of the rotating detonation engine configuration.

Description

Technical Field

[0001] The invention belongs to the field of numerical calculation of flow field of rotating detonation engine, and in particular relates to a method for predicting flow field parameter distribution of rotating detonation engine based on space propulsion algorithm. Background Art

[0002] As a new type of propulsion device, the rotating detonation engine Figure 2 As shown in the figure, the common structure is a concentric annular column with a combustible mixture injection slot at the bottom of the annular column. There are single or multiple detonation waves in the annular cavity that are inclined toward the injection plane and propagate circumferentially. The gas immediately behind the detonation wave releases heat rapidly, and the temperature and pressure are significantly increased, and then a series of expansion waves accelerate the expansion of the airflow. When the pressure of the combustion products close to the injection plane is lower than the total injection pressure of the combustible mixture, the combustible mixture can be injected into the annular cavity. The detonation wave and the combustible mixture / combustion product discontinuity surface form a triangular area filled with combustible mixture to ensure the sustainability of detonation combustion. In terms of the flow field structure, the detonation wave, shear layer and induced shock wave intersect at the upper vertex of the triangular area. Most of the airflow is discharged along the axial direction of the rotating detonation engine outlet, but the flow parameters still vary significantly along the circumferential direction.

[0003] Compared with gas turbine engines and ramjets, rotating detonation engines have less entropy increase and higher thermal efficiency during operation, and have received extensive attention from various research institutions. In order to obtain complete information about the flow field of a rotating detonation engine, high-precision numerical simulation technology is essential. During the detonation combustion process, the moving shock wave is strongly coupled with the chemical exothermic reaction, resulting in the traditional time-advanced numerical simulation grid scale and discrete time step being too small. In order to obtain the aerodynamic performance of a rotating detonation engine, it often takes a lot of time and computing resources. Summary of the invention

[0004] Technical purpose: In order to solve the above technical problems, the present invention provides a method for predicting the flow field parameter distribution of a rotating detonation engine based on a space propulsion algorithm. By coupling the space propulsion algorithm and the CJ detonation wave model of a single γ, and utilizing the flow characteristics of the streamline grid corresponding to the space propulsion algorithm, the flow parameters after the detonation wave and the flow parameters far downstream of the detonation wave in the shock wave coordinate system are determined, and the flow parameters in the laboratory coordinate system are obtained through coordinate transformation, which effectively reduces the time consumption of solving the flow field of the rotating detonation engine and realizes the rapid prediction of the flow parameter distribution of the rotating detonation engine.

[0005] Technical solution: To achieve the above purpose, the technical solution adopted by the present invention is:

[0006] A method for predicting flow field parameter distribution of a rotating detonation engine based on a space propulsion algorithm, characterized in that it comprises the following steps:

[0007] (1) Estimate the average velocity at the injection surface inlet in advance and determine the detonation wave front flow parameters based on the compressible flow relationship;

[0008] (2) Substituting the flow parameters before the detonation wave into the CJ detonation model with a single specific heat ratio, the flow parameters after the detonation wave are determined;

[0009] (3) Assign the flow parameters after the detonation wave and the flow parameters far downstream to the initial value line; in the first iteration, the flow parameters far downstream of the detonation wave are obtained by expanding the pressure after the detonation wave to the pressure before the detonation wave, and the flow parameters far downstream of the detonation wave in subsequent iterations are based on the previous calculation results;

[0010] (4) Substituting the initial value line into the space propulsion algorithm, the flow field structure and flow parameter distribution of the rotating detonation engine in the shock wave coordinate system are determined;

[0011] (5) Compare the pressure value of the solution point on the injection surface with the total injection pressure. When the pressure value of the solution point is less than the total injection pressure, the injection flow field and the flow field downstream of the detonation wave are coupled and solved. When the pressure value of the solution point is greater than the total injection pressure, only the flow field downstream of the detonation wave is solved.

[0012] (6) Determine the area average of the flow parameters in front of the detonation wave at the far downstream position of the detonation wave, return to step (2) to step (5) to iterate until the average airflow angle far downstream is 0, exit the loop, and execute step (7);

[0013] (7) According to the inclination angle of the detonation wave, the rotation angle of the shock wave coordinate system is determined, and the flow parameter distribution of the rotating detonation engine in the laboratory coordinate system is determined according to the velocity triangle.

[0014] Furthermore, the step (1) comprises the following steps:

[0015] (11) According to the flight conditions of the rotating detonation engine, determine the total pressure Pt and total temperature Tt of the rotating detonation inlet combustible mixture, and estimate the average speed of the rotating detonation engine inlet

[0016] (12) According to the compressible flow relationship, the average static temperature of the detonation wave front is determined Average static pressure And the average density

[0017] Furthermore, the step (2) comprises the following steps:

[0018] (21) Construct a CJ knock model with a single specific heat ratio. The model equation is as follows:

[0019]

[0020]

[0021] Among them, M CJ represents the detonation wave Mach number, H represents the dimensionless heat release, Q represents the chemical reaction heat, R g represents the gas constant, γ is the specific heat ratio of the combustible mixture and the combustion products;

[0022] (22) According to the law of conservation of mass, momentum and energy before and after the detonation wave, the flow parameters after the detonation wave are determined, including the static pressure after the detonation wave. Post-blast density and post-blast temperature

[0023] Furthermore, the step (4) comprises the following steps:

[0024] (41) The intermediate variable χ is determined by the inviscid flux E in the mainstream direction, and then the state quantity Q on the initial value line is obtained. The expression of the intermediate variable χ is as follows:

[0025]

[0026] (42) According to the initial value of the linear velocity distribution, determine the coordinates of the downstream solution point Create a streamline mesh, represents the x-coordinate of the jth downstream solution point in the x-direction, Represents the x-coordinate of the jth initial value line point in the x-direction;

[0027] (43) According to the intrinsically non-oscillatory interpolation method and the limit function, we obtain the half-step along the x direction. and The numerical solution of the step length is 1 / 2, and then the state quantity is obtained. and

[0028] (44) According to the shock wave polar curve theory, the pressure and airflow angle of the solution point at the position where the step length is 1 / 2 of the set value corresponding to the spatial advancement algorithm are determined through two adjacent points in the spatial x direction, and then the Riemann self-similar solution is obtained. The expression of the state curve of the two adjacent points is as follows:

[0029]

[0030]

[0031] Among them, Φ B The state curve representing the points below the initial value line, θ B represents the airflow angle at the lower point of the initial value line, f(M B ,αB ) represents the Rankine-Hsuehgon relationship for the points below the initial value line, α B It represents the ratio of the pressure at any point in the downstream space to the pressure at the point below the initial value line, ν(M B ) represents the Prandtl-Meyer function of the lower point of the initial value line, Φ T Represents the state curve of the upper point of the initial value line, θ T represents the airflow angle of the upper point on the initial value line, f(M T ,α T ) represents the Rankine-Hsuehgon relationship of the side points on the initial value line, α T It represents the ratio of the pressure at any point in the downstream space to the pressure at the upper point on the initial value line, ν(M T ) represents the Prandtl-Meyer function of the upper point of the initial value line;

[0032] (45) By differentially discretizing the compressible Euler equation, the flow parameters under a single spatial progressive length are obtained. The differential equation along the x direction of the flow downstream is as follows:

[0033]

[0034] in, Indicates the y-coordinate difference between two adjacent points on the initial value line. represents the y-coordinate difference between two adjacent solution points downstream, represents the inviscid flux along the y direction of the upper grid point at the position of 1 / 2 length of the spatial step, Indicates the slope of the upper grid edge at the 1 / 2 length position of the space step, represents the inviscid flux along the x direction at the upper grid point at the position of 1 / 2 length of the spatial step, represents the inviscid flux along the y direction at the lower grid point at the position of 1 / 2 length of the spatial step, Indicates the slope of the lower grid edge at the 1 / 2 length position of the space step, Represents the inviscid flux along the x-direction at the lower grid point at the position of 1 / 2 length of the spatial step.

[0035] Furthermore, the step (6) comprises the following steps:

[0036] (61) According to the flow parameter distribution along the y direction far downstream, the interval where the combustible mixture before the detonation wave is located is identified;

[0037] (62) Determine the area average of the flow parameters in the interval where the combustible mixture before the detonation wave is located. The area average formula is as follows:

[0038]

[0039]

[0040]

[0041]

[0042] in, represents the x-component of the average velocity of the detonation wave front, represents the y component of the average velocity of the detonation wave front, n x represents the x component of the detonation wave surface normal vector, n y represents the y component of the detonation wave surface normal vector, Represents the internal energy per unit mass. According to the Newton-Raphson multivariate iteration method, the four equations of the area average formula are solved simultaneously to obtain the detonation wave front flow parameters as well as

[0043] Furthermore, the step (7) comprises the following steps:

[0044] (71) According to the airflow angle in front of the detonation wave, the inclination angle of the detonation wave relative to the injection plane is corrected. The correction formula is as follows:

[0045]

[0046] in, represents the inclination angle of the modified detonation wave relative to the injection plane, represents the inclination angle of the uncorrected detonation wave relative to the injection plane;

[0047] (72) According to the inclination angle of the detonation wave relative to the injection plane after multiple iterations, the rotation angle of the shock wave coordinate system is determined, and the outlet direction is turned to the horizontal through the two-dimensional rotation matrix. The rotation angle and the two-dimensional rotation matrix are in the following form:

[0048]

[0049]

[0050]

[0051] in, represents the rotation angle of the shock wave coordinate system, x new represents the x-coordinate and y-coordinate of any point in the flow field after the shock wave coordinate system is rotated. new represents the y coordinate of any point in the flow field after the shock wave coordinate system is rotated, and x old represents the x-coordinate and y-coordinate of any point in the flow field without rotation of the shock wave coordinate system. old represents the y coordinate of any point in the flow field without rotation of the shock wave coordinate system, u newrepresents the velocity x component of any point in the flow field after the shock wave coordinate system is rotated, v new represents the velocity y component of any point in the flow field after the shock wave coordinate system is rotated, u old represents the x-component of the velocity at any point in the flow field without rotation of the shock wave coordinate system, v old represents the y component of the velocity at any point in the flow field when the shock wave coordinate system is not rotated;

[0052] (73) According to the principle of vector addition and subtraction, the velocity in the shock wave coordinate system is converted into the velocity in the laboratory coordinate system by the velocity triangle method. The conversion relationship is as follows:

[0053]

[0054] v lab =v new

[0055] Among them, u lab represents the x component of the velocity at any point in the laboratory coordinate system, v lab It represents the y component of the velocity at any point in the laboratory coordinate system.

[0056] Furthermore, the airflow angle of the upper section of the section where the combustible mixture before the detonation wave is located in step (6) is corrected to 0 in the shock wave coordinate system of the new round of iteration.

[0057] Beneficial effects: Due to the adoption of the above technical solution, the present invention has the following technical effects:

[0058] The present invention provides a method for predicting the flow field parameter distribution of a rotating detonation engine based on a space propulsion algorithm. By coupling the space propulsion algorithm and a single γ CJ detonation wave model, the flow parameters after the detonation wave and the flow parameters far downstream of the detonation wave in the shock wave coordinate system are determined using the flow characteristics of the streamline grid corresponding to the space propulsion algorithm, and the flow parameters in the laboratory coordinate system are obtained through coordinate transformation. The flow field of the rotating detonation engine can be quickly solved, the aerodynamic performance of the rotating detonation engine can be quickly evaluated, and the flow structures such as the detonation wave, shear layer and induced shock wave in the flow field can be accurately captured, thereby greatly shortening the configuration design cycle of the rotating detonation engine. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 The calculation domain in the shock wave coordinate system determined in steps (1) to (6) of the method for predicting the flow field parameter distribution of a rotating detonation engine based on a space propulsion algorithm of the present invention;

[0060] Figure 2 It is a schematic diagram of the physical structure of a rotating detonation engine;

[0061] Figure 3It is the intersection process of the shock polar curves of two adjacent points on the initial value line;

[0062] Figure 4 It is the stepping process of the spatial advancement algorithm from the initial value line to the solution point, EP represents the grid edge point, CP represents the grid center point, Const X represents the initial value line and the straight line where the solution point is located, PM represents the Prandtl-Meyer expansion wave, Sh represents the shock wave, and Sl represents the inviscid slip line;

[0063] Figure 5 It is the convergence process of detonation wave parameters in the flow field of rotating detonation engine;

[0064] Figure 6 is the streamline grid of the rotating detonation engine in the shock wave coordinate system;

[0065] Figure 7 The temperature field of the rotating detonation engine obtained for the space propulsion algorithm;

[0066] Figure 8 The velocity field of the rotating detonation engine in the laboratory coordinate system obtained for the space propulsion algorithm;

[0067] In the figure, 11-the initial value line corresponding to the upstream detonation wave, 12-the initial value line corresponding to the far downstream flow field of the detonation wave, 13-the initial value line corresponding to the redundant section of the far downstream flow field of the detonation wave, 14-the imaginary non-stick wall surface with an inclination angle consistent with the airflow turning angle of the redundant section, 15-the injection plane blocking section with an inclination angle equal to the inclination angle of the detonation wave, 16-the simulated intake section of the combustible mixture, 17-the straight line corresponding to the solution point of the far downstream flow field of the detonation wave, 18-the location of the solution point corresponding to the downstream detonation wave, 19-the imaginary non-stick wall surface with an inclination angle consistent with the airflow turning angle of the simulated intake section of the combustible mixture. DETAILED DESCRIPTION

[0068] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0069] The present invention discloses a method for predicting the flow field parameter distribution of a rotating detonation engine based on a space propulsion algorithm, comprising the following steps: pre-estimating the flow parameters before the detonation wave; determining the flow parameters after the detonation wave according to a CJ detonation model of a single specific heat ratio; assigning the flow parameters after the detonation wave and the far downstream flow parameters to an initial value line; determining the flow parameters in a shock wave coordinate system by means of a space propulsion algorithm; determining the starting position and the ending position of the combustible mixture injection; iteratively converging the flow field parameters of the rotating detonation engine, and obtaining the flow parameters in a laboratory coordinate system by means of coordinate transformation.

[0070] Combination Figure 1As shown, 11 is the initial value line corresponding to the upstream detonation wave, 12 is the initial value line corresponding to the far downstream flow field of the detonation wave, and 13 is the initial value line corresponding to the redundant section of the far downstream flow field of the detonation wave. 14 corresponds to the imaginary non-stick wall surface with an inclination angle consistent with the airflow turning angle of the redundant section, 15 corresponds to the injection plane blocking section with an inclination angle equal to the inclination angle of the detonation wave, 16 corresponds to the simulated intake section of the combustible mixture, 17 corresponds to the straight line where the solution point of the far downstream flow field of the detonation wave is located, 18 corresponds to the location of the solution point of the downstream detonation wave, and 19 corresponds to the imaginary non-stick wall surface with an inclination angle consistent with the airflow turning angle of the simulated intake section of the combustible mixture.

[0071] The shock wave coordinate system is a coordinate system formed by assuming that the moving detonation wave line segment AB remains stationary, with the direction of line segment AB as the y coordinate direction and the direction perpendicular to line segment AB as the x coordinate direction. The laboratory coordinate system is a coordinate system formed by taking the direction group of line segment AM of the injection plane as the x coordinate direction and the direction perpendicular to line segment AM as the y coordinate direction.

[0072] θ ave Angle is the average airflow angle, located slightly upstream of line segment IJ, and angle θ is the complementary angle between the injection plane (the inclined line segment AJ in the shock wave coordinate system) and the detonation wave AB.

[0073] The prediction method specifically includes the following steps:

[0074] Step (1), pre-estimating the average velocity at the injection surface inlet, and determining the detonation wave front flow parameters according to the compressible flow relationship;

[0075] Step (2), substituting the flow parameters before the detonation wave into the CJ detonation model with a single specific heat ratio γ to determine the flow parameters after the detonation wave;

[0076] Step (3), the flow parameters far downstream of the detonation wave in the first iteration are obtained by expanding the pressure after the detonation wave to the pressure before the detonation wave. The flow parameters far downstream of the detonation wave in subsequent iterations are assigned to the initial value line based on the previous calculation results. Figure 1 As shown, the initial value line is composed of line segments AB and BD, and the flow parameters after the detonation wave AB are consistent with the flow parameters after the detonation wave of line segment IJ calculated last time. At the same time, the flow parameter distribution of the line segment BD above the detonation wave is consistent with the flow parameter distribution of the line segment IE calculated last time, and the specific values ​​are obtained by linear interpolation.

[0077] Step (4), substituting the initial value line into the space propulsion algorithm to determine the flow field structure and flow parameter distribution of the rotating detonation engine in the shock wave coordinate system;

[0078] Step (5), compare the pressure value of the solution point of the injection surface with the total injection pressure. When the pressure value of the solution point is less than the total injection pressure, the injection flow field and the flow field downstream of the detonation wave are coupled and solved; if the pressure value of the solution point is greater than or equal to the total pressure, the injection flow field does not exist, and at this time only the flow field downstream of the detonation wave is solved.

[0079] Step (6), determining the area average of the flow parameters in front of the detonation wave at the far downstream position of the detonation wave, returning to step (2) to iterate until the average airflow angle at the far downstream position is 0, and then exiting the loop;

[0080] Step (7), determining the rotation angle of the shock wave coordinate system according to the inclination angle of the detonation wave, and determining the flow parameter distribution of the rotating detonation engine in the laboratory coordinate system according to the velocity triangle;

[0081] Step (8), outputting the flow parameter distribution of the rotating detonation engine in the laboratory coordinate system into a data file.

[0082] Combination Figure 3 In the present invention, the intersection parameters of the state curves of two adjacent points are obtained by the Newton-Raphson iterative method. Figure 4 In the present invention, in the process of advancing from the initial value line to the downstream solution point, the flow parameters are obtained by directly calculating the interaction between the shock wave and the expansion wave, without the need for an additional unit point construction process, and a streamline grid is formed by adaptively changing the slope of the grid edge.

[0083] Step (1) of the present invention specifically comprises the following steps:

[0084] (11) According to the flight conditions of the rotating detonation engine, determine the total pressure Pt and total temperature Tt of the rotating detonation inlet combustible mixture, and estimate the average speed of the rotating detonation engine inlet

[0085] (12) According to the compressible flow relationship, the detonation wave front flow parameters are determined, including the average static temperature Average static pressure And the average density

[0086] The following relationship exists:

[0087]

[0088]

[0089]

[0090] Among them, the C p and γ are the constant pressure specific heat and specific heat ratio of the combustible mixture and the combustion products, respectively. g represents the gas constant.

[0091] Step (2) specifically includes the following steps:

[0092] (21) According to the CJ knock model of double γ, the CJ knock model of single specific heat ratio γ can be simplified:

[0093]

[0094]

[0095] (22) According to the law of conservation of mass, momentum and energy before and after the detonation wave, the expression of the flow parameters after the detonation wave is determined as follows:

[0096]

[0097]

[0098]

[0099] Among them, M CJ represents the detonation wave Mach number, H represents the dimensionless heat release, Q represents the chemical reaction heat, R g is the gas constant, represents the average static temperature before the detonation wave, represents the average density of the detonation wave front, represents the average static pressure before the detonation wave, represents the static pressure after the detonation wave, represents the density after the detonation wave, Represents the temperature after the detonation wave.

[0100] Step (4) specifically includes the following steps:

[0101] (41) The intermediate variable χ is determined by the inviscid flux E in the mainstream direction, and then the state quantity Q on the initial value line is obtained. The expression of the intermediate variable χ is as follows:

[0102]

[0103] E1, E2, E3, and E4 refer to the inviscid flux in the mainstream direction;

[0104] (42) According to the initial value of the linear velocity distribution, the coordinates of the downstream solution point are determined and the streamline grid is established. The expression of the downstream solution point is as follows:

[0105]

[0106]

[0107] represents the x-coordinate of the jth downstream solution point in the x-direction, represents the x-coordinate of the jth initial value line point in the x-direction, Δx represents the step length of the grid along the x-direction, represents the y-coordinate of the jth downstream solution point in the y direction, represents the x-coordinate of the jth initial value line point in the y direction, represents the y-direction velocity component of the j-th initial value line point in the y direction, It represents the x-direction velocity component of the j-th initial value line point in the y direction;

[0108] (43) According to the intrinsically non-oscillatory interpolation method and the limit function, we obtain the half-step along the x direction. and The numerical solution of the step length is 1 / 2, and then the state quantity is obtained. and

[0109] (44) According to the shock wave extreme curve theory, the pressure and airflow angle of the solution point at the position where the step length is 1 / 2 of the space advancement algorithm corresponding to the step length is 1 / 2 are determined through two adjacent points in the space x direction with a step length of 1 / 2 of the set value, and then the Riemann self-similar solution is obtained. The state curves of the two adjacent points used include the following shock wave extreme curve equation:

[0110]

[0111]

[0112] Among them, Φ B The state curve representing the points below the initial value line, θ B represents the airflow angle at the lower point of the initial value line, f(M B ,α B ) represents the Rankine-Hsuehgon relationship for the points below the initial value line, α B It represents the ratio of the pressure at any point in the downstream space to the pressure at the point below the initial value line, ν(M B ) represents the Prandtl-Meyer function of the lower point of the initial value line, Φ T Represents the state curve of the upper point of the initial value line, θ T represents the airflow angle of the upper point on the initial value line, f(M T ,α T ) represents the Rankine-Hsuehgon relationship of the side points on the initial value line, α T It represents the ratio of the pressure at any point in the downstream space to the pressure at the upper point on the initial value line, ν(M T ) represents the Prandtl-Meyer function of the upper point of the initial value line.

[0113] (45) By differentially discretizing the compressible Euler equation, the flow parameters under a single spatial progressive length are obtained. The differential equation along the x direction of the flow downstream is as follows:

[0114]

[0115] in, Indicates the y-coordinate difference between two adjacent points on the initial value line. represents the y-coordinate difference between two adjacent solution points downstream, represents the inviscid flux along the y direction of the upper grid point at the position of 1 / 2 length of the spatial step, Indicates the slope of the upper grid edge at the 1 / 2 length position of the space step, represents the inviscid flux along the x direction at the upper grid point at the position of 1 / 2 length of the spatial step, represents the inviscid flux along the y direction at the lower grid point at the position of 1 / 2 length of the spatial step, Indicates the slope of the lower grid edge at the 1 / 2 length position of the space step, Represents the inviscid flux along the x-direction at the lower grid point at the position of 1 / 2 length of the spatial step.

[0116] In step (5), the pressure value of the solution point of the injection surface is compared with the total injection pressure. The solution point is the downstream point of the n+1th ConstX line obtained by steps (44) and (45) from the upstream point on the nth ConstX line in the flow field region. Where ConstX represents Figure 1 A straight line parallel to the blast wave AB.

[0117] The step (6) comprises the following steps:

[0118] Step (61), identifying the interval where the combustible mixture before the detonation wave is located according to the flow parameter distribution along the y direction far downstream;

[0119] Step (62), determining the area average value of the flow parameters of the interval where the combustible mixture before the detonation wave is located, the area average formula is as follows:

[0120]

[0121]

[0122]

[0123]

[0124] in, represents the x-component of the average velocity of the detonation wave front, represents the y component of the average velocity of the detonation wave front, n xrepresents the x component of the detonation wave surface normal vector, n y represents the y component of the detonation wave surface normal vector, represents the internal energy per unit mass, and A represents the area occupied by the detonation wave. According to the Newton-Raphson multivariate iteration method, the four equations of the area average formula are solved simultaneously to obtain the detonation wave front flow parameters. as well as

[0125] During iteration, the previous calculation result used in assigning values ​​in step (3) refers to the post-detonation wave flow parameters obtained by further calculating the pre-detonation wave flow parameters obtained in step (62) through step (2).

[0126] The step (7) comprises the following steps:

[0127] (71) According to the airflow angle in front of the detonation wave, the inclination angle of the detonation wave relative to the injection plane is corrected. The correction formula is as follows:

[0128]

[0129] in, represents the inclination angle of the modified detonation wave relative to the injection plane, represents the inclination angle of the uncorrected detonation wave relative to the injection plane;

[0130] (72) According to the inclination angle of the detonation wave relative to the injection plane after multiple iterations, the rotation angle of the shock wave coordinate system is determined, and the outlet direction is turned to the horizontal through the two-dimensional rotation matrix. The rotation angle and the two-dimensional rotation matrix are in the following form:

[0131]

[0132]

[0133]

[0134] in, represents the rotation angle of the shock wave coordinate system, x new represents the x-coordinate and y-coordinate of any point in the flow field after the shock wave coordinate system is rotated. new represents the y coordinate of any point in the flow field after the shock wave coordinate system is rotated, and x old represents the x-coordinate and y-coordinate of any point in the flow field without rotation of the shock wave coordinate system. old represents the y coordinate of any point in the flow field without rotation of the shock wave coordinate system, u new represents the velocity x component of any point in the flow field after the shock wave coordinate system is rotated, v new represents the velocity y component of any point in the flow field after the shock wave coordinate system is rotated, u oldrepresents the x-component of the velocity at any point in the flow field without rotation of the shock wave coordinate system, v old represents the y component of the velocity at any point in the flow field when the shock wave coordinate system is not rotated;

[0135] (73) According to the principle of vector addition and subtraction, the velocity in the shock wave coordinate system is converted into the velocity in the laboratory coordinate system by the velocity triangle method. The conversion relationship is as follows:

[0136]

[0137] v lab =v new

[0138] Among them, u lab represents the x component of the velocity at any point in the laboratory coordinate system, v lab represents the y component of the velocity at any point in the laboratory coordinate system. The velocity triangle is as follows Figure 1 As shown, It represents the flow field velocity in the AMLKED region in the shock wave coordinate system, and the conversion between velocities is obtained by vector addition and subtraction. Since the spatial position of the shock wave is fixed in the shock wave coordinate system, and the flow structure such as the shock wave moves at high speed along the tangent direction of the injection plane in the actual flow (laboratory coordinate system), it is necessary to obtain the flow field of the actual rotating detonation engine through coordinate transformation according to the transformation rules of relative motion.

[0139] Furthermore, the airflow angle of the upper section of the section where the combustible mixture before the detonation wave is located in step (6) is corrected to 0 in the shock wave coordinate system of the new round of iteration.

[0140] The method of the present invention can be used to quickly evaluate the aerodynamic performance of a rotating detonation engine, that is, by establishing a mathematical model of the detonation process and its downstream flow field structure, the time-consuming and computationally resource-intensive characteristics of a traditional time-advancing algorithm can be avoided.

[0141] In order to better illustrate the present invention and facilitate understanding of the technical solution of the present invention, typical but non-limiting embodiments of the present invention are as follows:

[0142] The total injection pressure of the combustible mixture at the inlet of the rotating detonation engine is 5×10 5 Pa, the total injection temperature is 300K, the initial value of the average inlet velocity is 450m / s, the specific heat ratio of the combustible premixed gas is 1.3961, and the gas constant of the combustible premixed gas is 395.75J·kg -1 ·K -1 The specific heat ratio of the combustion products is 1.1653, and the gas constant of the combustion products is 346.20 J·kg -1 ·K -1 , the heat release per unit mass is 5.4704×10 6J·kg -1 The combustible mixture activation temperature is 15100K, and the single-step chemical reaction pre-exponential factor is 1.0×10 9 s -1 The rotating detonation engine center ring has an unfolded length of 314mm and an axial length of 75mm.

[0143] Figure 5 is the convergence process of the detonation wave parameters during the fast solution iteration of the present invention, θ inc represents the detonation wave inclination, H D represents the length of the detonation wave, P1 / P2 represents the static pressure ratio before and after the detonation wave, U lab Represents the detonation wave velocity. Figure 6 It is the streamline grid in the shock wave coordinate system during the rapid solution of the present invention. Figure 7 It is the temperature field of the rotating detonation engine obtained by the present invention, which accurately captures the shock wave supercharging process and the airflow expansion acceleration process after the detonation wave. Figure 8 It is the velocity field of the rotating detonation engine obtained by the present invention in the laboratory coordinate system, which accurately captures the airflow direction change process of the rotating detonation engine. The laboratory coordinate system is defined as a coordinate system formed by taking the direction group of the line segment AM of the injection plane as the x-coordinate direction and the direction perpendicular to the line segment AM as the y-coordinate direction.

[0144] Table 1 is a comparison of the detonation wave parameters obtained by the rapid solution of the present invention and the detonation wave parameters obtained by the traditional time-marching algorithm.

[0145] Table 1

[0146]

[0147] Compared with the detonation wave parameters of the rotating detonation engine obtained by the traditional time-marching algorithm, the relative deviation of the detonation wave height obtained by the present invention is 1.38%, the relative deviation of the detonation wave inclination is 1.35%, the relative deviation of the static pressure before and after the detonation wave is 1.05%, and the relative deviation of the detonation wave velocity is 0.97%. It can be seen that the space-marching algorithm has the same accuracy as the time-marching algorithm in solving the rotating detonation flow field.

[0148] Table 2 is the comparison data of the time consumption of solving the flow field of the rotating detonation engine by the present invention and the time consumption of solving by the traditional time-marching algorithm under the environment of i7 CPU with 4 cores and 8 threads.

[0149] Table 2

[0150] algorithm CPU number and model time consuming Space Advancement Algorithm 4 cores and 8 threads, i7 5h Traditional time marching algorithm 4 cores and 8 threads, i7 48h

[0151] The present invention is a method for predicting the distribution of flow field parameters of a rotating detonation engine based on a space-advancing algorithm. The method uses a 4-core 8-thread i7 CPU to solve the calculation, and the total time consumption is 5 hours. The traditional time-advancing algorithm uses 4 cores and 8 threads to solve the calculation, and the total time consumption is 48 hours. The present invention is 8.6 times faster than the traditional time-advancing algorithm, and can effectively reduce the design and selection cycle of the rotating detonation engine configuration.

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

Claims

1. A method for predicting flow field parameter distribution of a rotating detonation engine based on a space propulsion algorithm, characterized in that: The following steps are involved: (1) Pre-estimating the average velocity at the injection surface inlet and determining the detonation wave front flow parameters according to the compressible flow relationship, including the following steps: (11) According to the flight conditions of the rotating detonation engine, the total pressure P of the combustible mixture at the inlet of the rotating detonation engine is determined. t and total temperature T t , and estimate the average speed at the inlet of the rotating detonation engine (12) According to the compressible flow relationship, the average static temperature of the detonation wave front is determined Average static pressure P1 and average density C p and γ are the constant pressure specific heat and specific heat ratio of the combustible mixture and the combustion products, respectively. g represents the gas constant; (2) Substituting the flow parameters before the detonation wave into the CJ detonation model of a single specific heat ratio to determine the flow parameters after the detonation wave, including the steps of: (21) Construct a CJ knock model with a single specific heat ratio. The model equation is as follows: M CJ represents the detonation wave Mach number, H represents the dimensionless heat release, and Q represents the chemical reaction heat; (22) According to the law of conservation of mass, momentum and energy before and after the detonation wave, the flow parameters after the detonation wave are determined, including the static pressure after the detonation wave. Post-blast density and post-blast temperature (3) Assign the flow parameters after the detonation wave and the flow parameters far downstream to the initial value line; in the first iteration, the flow parameters far downstream of the detonation wave are obtained by expanding the pressure after the detonation wave to the pressure before the detonation wave, and the flow parameters far downstream of the detonation wave in subsequent iterations are based on the previous calculation results; (4) Substituting the initial value line into the space propulsion algorithm to determine the flow field structure and flow parameter distribution of the rotating detonation engine in the shock wave coordinate system, including the steps of: (41) The intermediate variable χ is determined by the inviscid flux E in the mainstream direction, and then the state quantity Q on the initial value line is obtained. The expression of the intermediate variable χ is as follows: E1, E2, E3, and E4 refer to the inviscid flux in the mainstream direction; (42) According to the initial value of the linear velocity distribution, determine the coordinates of the downstream solution point Create a streamline mesh: represents the x-coordinate of the jth downstream solution point in the x-direction, represents the x-coordinate of the jth initial value line point in the x-direction, Δx represents the step length of the grid along the x-direction, represents the y coordinate of the jth downstream solution point in the y direction, represents the x-coordinate of the jth initial value line point in the y direction, represents the y-direction velocity component of the j-th initial value line point in the y direction, It represents the x-direction velocity component of the j-th initial value line point in the y direction; (43) According to the intrinsically non-oscillatory interpolation method and the limit function, we obtain the half-step along the x direction. and The numerical solution of the step length is 1 / 2, and then the state quantity is obtained. and (44) According to the shock wave polar curve theory, the pressure and airflow angle of the solution point at the position where the step length is 1 / 2 of the set value corresponding to the spatial advancement algorithm are determined through two adjacent points in the spatial x direction, and then the Riemann self-similar solution is obtained. The expression of the state curve of the two adjacent points is as follows: Among them, Φ B The state curve representing the points below the initial value line, θ B represents the airflow angle at the lower point of the initial value line, f(M B ,α B ) represents the Rankine-Hsuehgon relationship for the points below the initial value line, α B It represents the ratio of the pressure at any point in the downstream space to the pressure at the point below the initial value line, ν(M B ) represents the Prandtl-Meyer function of the lower point of the initial value line, Φ T Represents the state curve of the upper point of the initial value line, θ T represents the airflow angle of the upper point on the initial value line, f(M T ,α T ) represents the Rankine-Hsuehgon relationship of the side points on the initial value line, α T It represents the ratio of the pressure at any point in the downstream space to the pressure at the upper point on the initial value line, ν(M T ) represents the Prandtl-Meyer function of the upper point of the initial value line; (45) By differentially discretizing the compressible Euler equation, the flow parameters under a single spatial progressive length are obtained. The differential equation along the x direction of the flow downstream is as follows: in, Indicates the y-coordinate difference between two adjacent points on the initial value line. represents the y-coordinate difference between two adjacent solution points downstream, represents the inviscid flux along the y direction of the upper grid point at the position of 1 / 2 length of the spatial step, Indicates the slope of the upper grid edge at the 1 / 2 length position of the space step, represents the inviscid flux along the x direction at the upper grid point at the position of 1 / 2 length of the spatial step, represents the inviscid flux along the y direction at the lower grid point at the position of 1 / 2 length of the spatial step, Indicates the slope of the lower grid edge at the 1 / 2 length position of the space step, represents the inviscid flux along the x direction at the lower grid point at the position of 1 / 2 length of the spatial step; (5) Compare the pressure value of the solution point on the injection surface with the total injection pressure. When the pressure value of the solution point is less than the total injection pressure, the injection flow field and the flow field downstream of the detonation wave are coupled and solved. When the pressure value of the solution point is greater than or equal to the total injection pressure, only the flow field downstream of the detonation wave is solved. (6) Determine the area average of the flow parameters in front of the detonation wave at the far downstream position of the detonation wave, return to step (2) to step (5) to iterate until the average airflow angle far downstream is 0, exit the loop, and execute step (7); The step of determining the area average of the flow parameters in front of the detonation wave at a position far downstream of the detonation wave comprises the following steps: (61) According to the flow parameter distribution along the y direction far downstream, the interval where the combustible mixture before the detonation wave is located is identified; (62) Determine the area average of the flow parameters in the interval where the combustible mixture before the detonation wave is located. The area average formula is as follows: in, represents the x-component of the average velocity of the detonation wave front, represents the y component of the average velocity of the detonation wave front, n x represents the x component of the detonation wave surface normal vector, n y represents the y component of the detonation wave surface normal vector, Represents the internal energy per unit mass. According to the Newton-Raphson multivariate iteration method, the four equations of the area average formula are solved simultaneously to obtain the detonation wave front flow parameters as well as A represents the area occupied by the blast wave; (7) According to the inclination angle of the detonation wave, the rotation angle of the shock wave coordinate system is determined, and according to the velocity triangle, the flow parameter distribution of the rotating detonation engine in the laboratory coordinate system is determined, including the steps of: (71) According to the airflow angle in front of the detonation wave, the inclination angle of the detonation wave relative to the injection plane is corrected. The correction formula is as follows: in, represents the inclination angle of the modified detonation wave relative to the injection plane, represents the inclination angle of the uncorrected detonation wave relative to the injection plane; (72) According to the inclination angle of the detonation wave relative to the injection plane after multiple iterations, the rotation angle of the shock wave coordinate system is determined, and the outlet direction is turned to the horizontal through the two-dimensional rotation matrix. The rotation angle and the two-dimensional rotation matrix are in the following form: in, represents the rotation angle of the shock wave coordinate system, xnew represents the x coordinate of any point in the flow field after the shock wave coordinate system is rotated, and ynew represents the y coordinate of any point in the flow field after the shock wave coordinate system is rotated. old represents the x-coordinate and y-coordinate of any point in the flow field without rotation of the shock wave coordinate system. old represents the y coordinate of any point in the flow field without the rotation of the shock wave coordinate system, unew represents the x component of the velocity of any point in the flow field after the rotation of the shock wave coordinate system, and v new represents the velocity y component of any point in the flow field after the shock wave coordinate system is rotated, u old represents the x-component of the velocity at any point in the flow field without rotation of the shock wave coordinate system, v old represents the y component of the velocity at any point in the flow field when the shock wave coordinate system is not rotated; (73) According to the principle of vector addition and subtraction, the velocity in the shock wave coordinate system is converted into the velocity in the laboratory coordinate system by the velocity triangle method. The conversion relationship is as follows: v lab =v new Among them, u lab represents the x component of the velocity at any point in the laboratory coordinate system, v lab It represents the y component of the velocity at any point in the laboratory coordinate system.

2. The method for predicting flow field parameter distribution of a rotating detonation engine based on a space propulsion algorithm according to claim 1, characterized in that: In the step (6), the airflow angle of the upper section of the section where the combustible mixture before the detonation wave is located is corrected to 0 in the shock wave coordinate system of the new round of iteration.

Images