CN113987694A - Rotary detonation engine flow field parameter distribution prediction method based on space propulsion algorithm - Google Patents
Rotary detonation engine flow field parameter distribution prediction method based on space propulsion algorithm Download PDFInfo
- Publication number
- CN113987694A CN113987694A CN202111092248.1A CN202111092248A CN113987694A CN 113987694 A CN113987694 A CN 113987694A CN 202111092248 A CN202111092248 A CN 202111092248A CN 113987694 A CN113987694 A CN 113987694A
- Authority
- CN
- China
- Prior art keywords
- representing
- detonation
- point
- flow
- detonation wave
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000005474 detonation Methods 0.000 title claims abstract description 215
- 238000000034 method Methods 0.000 title claims abstract description 47
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 43
- 230000035939 shock Effects 0.000 claims abstract description 65
- 238000002347 injection Methods 0.000 claims abstract description 47
- 239000007924 injection Substances 0.000 claims abstract description 47
- 239000000203 mixture Substances 0.000 claims abstract description 19
- 238000006243 chemical reaction Methods 0.000 claims abstract description 14
- 239000000243 solution Substances 0.000 claims description 40
- 230000004907 flux Effects 0.000 claims description 16
- 230000003068 static effect Effects 0.000 claims description 14
- 238000002485 combustion reaction Methods 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 3
- 239000006185 dispersion Substances 0.000 claims description 3
- 238000004134 energy conservation Methods 0.000 claims description 3
- 230000009191 jumping Effects 0.000 claims description 3
- 230000000750 progressive effect Effects 0.000 claims description 3
- 230000008878 coupling Effects 0.000 abstract description 3
- 238000010168 coupling process Methods 0.000 abstract description 3
- 238000005859 coupling reaction Methods 0.000 abstract description 3
- 238000013461 design Methods 0.000 abstract description 3
- 230000007547 defect Effects 0.000 abstract 1
- 230000008569 process Effects 0.000 description 10
- 239000000853 adhesive Substances 0.000 description 7
- 230000001070 adhesive effect Effects 0.000 description 7
- 230000000903 blocking effect Effects 0.000 description 4
- 238000011144 upstream manufacturing Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 2
- 230000001133 acceleration Effects 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000004323 axial length Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000004880 explosion Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000006049 ring expansion reaction Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Computer Hardware Design (AREA)
- Fluid Mechanics (AREA)
- Mathematical Physics (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Computational Mathematics (AREA)
- Combined Controls Of Internal Combustion Engines (AREA)
Abstract
The invention discloses a method for predicting the distribution of flow field parameters of a rotary detonation engine based on a space propulsion algorithm, which comprises the following steps: pre-estimating the flow parameters of the front of the detonation wave; determining a flowing parameter after a detonation wave according to a CJ detonation model with 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 flow parameters under a shock wave coordinate system through a space propulsion algorithm; determining the initial position and the final position of combustible gas mixture injection; and (3) iteratively converging the flow field parameters of the rotary detonation engine, and obtaining the flow parameters in a laboratory coordinate system through coordinate conversion. According to the method, the flow parameters under a laboratory coordinate system are obtained through a coupling space propulsion algorithm and a CJ detonation wave model through coordinate conversion, the flow field structure and the outlet flow parameter distribution of the rotary detonation engine can be rapidly obtained, the defects that the conventional algorithm is long in time consumption for obtaining the flow field of the rotary detonation engine and the like are overcome, and the design cycle of the rotary detonation engine is shortened.
Description
Technical Field
The invention belongs to the field of numerical calculation of a flow field of a rotary detonation engine, and particularly relates to a rotary detonation engine flow field parameter distribution prediction method based on a space propulsion algorithm.
Background
As a novel propelling device, as shown in FIG. 2, a common structure of a rotary detonation engine is a concentric ring column, and a combustible mixture injection notch is formed in the bottom end of the ring column. Within the annular chamber there are single or multiple detonation waves inclined towards the injection plane and propagating circumferentially. The gas immediately after detonation wave is subjected to rapid heat release, the temperature and pressure are obviously improved, and then a series of expansion waves accelerate the expansion of the gas flow. The combustible mixture is injected into the annular chamber when the pressure of the combustion products close to the injection plane is lower than the total injection pressure of the combustible mixture. The detonation wave and the combustible mixture/combustion product discontinuity enclose a triangular area filled with the combustible mixture, so that the sustainability of detonation combustion is ensured. In the flow field structure, the detonation wave, the shear layer and the induced shock wave intersect at the upper vertex of the triangular region. Most of the airflow is discharged axially along the outlet of the rotary detonation engine, but the variation of the flow parameters along the circumferential direction is still significant.
Compared with a gas turbine engine and a ramjet engine, the rotary detonation engine has the advantages of smaller entropy increase and higher thermal efficiency in work, and is widely concerned by various research organizations. In order to obtain complete information of a flow field of the rotary detonation engine, a high-precision numerical simulation technology is indispensable. The kinetic shock wave and the chemical exothermic reaction are strongly coupled in the detonation combustion process, so that the traditional numerical simulation grid scale and discrete time step based on time advance are too small. In order to obtain the aerodynamic performance of a rotary detonation engine, much time and computing resources are consumed.
Disclosure of Invention
The technical purpose is as follows: in order to solve the technical problems, the invention provides a method for predicting the distribution of flow field parameters of a rotary detonation engine based on a space propulsion algorithm.
The technical scheme is as follows: in order to achieve the purpose, the invention adopts the technical scheme that:
a method for predicting the distribution of flow field parameters of a rotary detonation engine based on a space propulsion algorithm is characterized by comprising the following steps of: the method comprises the following steps:
(1) estimating the average speed of an inlet of an injection surface in advance, and determining the flow parameters of the detonation wave front according to a compressible flow relation;
(2) substituting the flow parameters of the detonation wave front into a CJ detonation model with a single specific heat ratio to determine the flow parameters of the detonation wave front;
(3) assigning the flow parameters after the detonation wave and the flow parameters far downstream to an initial value line; 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, and the flow parameters far downstream of the detonation wave in the later iteration are obtained according to the calculation result of the previous iteration;
(4) substituting the initial value line into a space propulsion algorithm, and determining a flow field structure and flow parameter distribution of the rotary detonation engine under a shock wave coordinate system;
(5) comparing 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 smaller than the total injection pressure, performing coupled solution on the injection flow field and the detonation wave downstream flow field, and when the pressure value of the solution point is larger than the total injection pressure, only solving the detonation wave downstream flow field;
(6) determining the area average value of the flow parameters of the detonation wave front at the far downstream position of the detonation wave, returning to the step (2) to the step (5) for iteration, jumping out of a loop until the average airflow angle at the far downstream is 0, and executing the step (7);
(7) and 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 rotary detonation engine under the laboratory coordinate system according to the velocity triangle.
Further, the step (1) comprises the following steps:
(11) determining the total pressure Pt and the total temperature Tt of the combustible mixture at the rotary detonation inlet according to the flight condition of the rotary detonation engine, and estimating the average speed at the inlet of the rotary detonation engine
(12) Determining the average static temperature of the detonation wave front according to the compressible flow relation
Mean static pressure
And average density
Further, the step (2) comprises the following steps:
(21) constructing a CJ knock model with a single specific heat ratio, wherein a model equation is in the form of:
wherein M isCJRepresenting the Mach number of the detonation wave, H representing the dimensionless heat release amount, Q representing the heat of chemical reaction, RgRepresenting a gas constant, gamma is the specific heat ratio of the combustible mixture and the combustion products;
(22) determining post-detonation wave flow parameters including post-detonation wave static pressure parameters according to the mass conservation law, momentum conservation law and energy conservation law before and after the detonation wave
Density after detonation wave
And post detonation wave temperature
Further, the step (4) comprises the following steps:
(41) determining an intermediate variable chi through the non-viscous flux E in the main flow direction, and further obtaining a state quantity Q on an initial value line, wherein the expression of the intermediate variable chi is in the following form:
(42) determining downstream solution point coordinates according to the velocity distribution on the initial value line
A streamline grid is established, and the streamline grid is established,
representing the x coordinate of the jth downstream solution point in the x direction,
representing the x coordinate of the jth initial value line point in the x direction;
(43) obtaining step by step half-steps along the x-direction according to an essentially oscillationless interpolation method and a limiting function
And
to obtain the state quantity of 1/2 position of step length
And
(44) according to the shock wave polar curve theory, the solution point pressure and the airflow angle of the space propulsion algorithm at the position 1/2 in the stepping length are determined through two adjacent points of which the stepping length in the space x direction is 1/2 set values, and then a Riemann self-similarity solution is obtained, wherein the expression of the state curve of the two adjacent points is as follows:
wherein phiBA state curve representing the lower side point of the initial value line, thetaBAir flow angle, f (M), representing the lower side point of the initial value lineB,αB) Rankine allowances relation, alpha, representing the lower side point of the initial value lineBDenotes the ratio of pressure at any point in the downstream space to pressure at the lower side of the initial value line, v (M)B) Prandtl-Meyer function, phi, representing the lower side point of the initial lineTA state curve representing the upper point of the initial value line, thetaTAirflow angle, f (M), representing the upper point of the initial value lineT,αT) Rankine allowances relation, alpha, representing upper side point of initial value lineTDenotes the ratio of pressure at any point in the downstream space to pressure at a point on the initial line, v (M)T) A planter-meier function representing an upper point of the initial value line;
(45) obtaining flow parameters under a single space progressive step length by carrying out differential dispersion on a compressible Euler equation, wherein the downstream differential equation along the flow x direction is in the following form:
wherein,
represents the difference value of y coordinates of two adjacent points of the initial value line,
the y coordinate difference of two adjacent solution points at the downstream is shown,
the inviscid flux in the y-direction of the upper grid point representing the length position of spatial step 1/2,
the slope of the upper grid edge representing the length position of spatial step 1/2,
the unbounded flux in the x-direction of the upper grid point representing the length position of spatial step 1/2,
the inviscid flux in the y-direction of the lower grid point representing the length position of spatial step 1/2,
the slope of the lower lattice edge representing the length position of spatial step 1/2,
the non-viscous flux in the x-direction of the lower grid point representing the length position of spatial step 1/2.
Further, the step (6) comprises the following steps:
(61) identifying the interval of the combustible mixed gas in the detonation wave front according to the flow parameter distribution of the far downstream along the y direction;
(62) determining the area average value of the flow parameters of the interval where the detonation wave front combustible mixed gas is located, wherein the area average formula is as follows:
wherein,
an x component representing the mean velocity of the detonation wave front,
y component representing mean velocity of detonation wave front, nxX component, n, representing the normal vector of the detonation wave frontyRepresents the y-component of the normal vector of the knock wave front,
expressing the internal energy contained in unit mass, simultaneously solving four equations of an area average formula according to a Newton-Laverson multivariate iterative method to obtain a detonation wave front flow parameter
And
further, the step (7) comprises the following steps:
(71) and correcting the inclination angle of the detonation wave relative to the injection plane according to the airflow angle of the detonation wave front, wherein the correction formula is as follows:
wherein,
indicating the inclination of the modified detonation wave relative to the injection plane,
representing the inclination of the uncorrected detonation wave relative to the injection plane;
(72) determining the rotation angle of a shock wave coordinate system according to the inclination angle of the detonation wave corrected by multiple iterations relative to the injection plane, and converting the outlet direction into the horizontal direction through a two-dimensional rotation matrix, wherein the rotation angle and the two-dimensional rotation matrix are in the following forms:
wherein,
representing the rotation angle, x, of the shock coordinate systemnewRepresenting the x coordinate and y coordinate of any point in the flow field after the rotation of the shock wave coordinate systemnewThe y coordinate and x of any point in the flow field after the shock wave coordinate system rotates are expressedoldRepresenting the x coordinate, y, of any point in the non-rotating flow field of the shock coordinate systemoldY-coordinate, u, representing any point in the non-rotating flow field of the shock coordinate systemnewRepresenting the velocity x component, v, of any point in the flow field after the rotation of the shock wave coordinate systemnewRepresenting the velocity y component u of any point in the flow field after the rotation of the shock wave coordinate systemoldIndicating shock wave seatThe velocity x component, v, of any point in the non-rotating flow fieldoldRepresenting the velocity y component of any point in the non-rotating flow field of the shock wave coordinate system;
(73) converting the speed in the shock wave coordinate system into the speed in the laboratory coordinate system by a speed triangle method according to the vector addition and subtraction principle, wherein the conversion relation is as follows:
vlab＝vnew
wherein u islabX component, v, representing the velocity at any point in the laboratory coordinate systemlabThe y component represents the velocity at any point in the laboratory coordinate system.
Further, in the step (6), the airflow angle in the upper section of the section where the knock wavefront combustible mixture is located is corrected to be 0 in the shock coordinate system of the new iteration.
Has the advantages that: due to the adoption of the technical scheme, the invention has the following technical effects:
according to the method for predicting the distribution of the parameters of the flow field of the rotary detonation engine based on the space propulsion algorithm, through coupling the space propulsion algorithm and a single gamma CJ detonation wave model, by utilizing the flow following characteristics of a streamline grid corresponding to the space propulsion algorithm, the flow parameters after the detonation wave and the far and downstream flow parameters of the detonation wave under a shock wave coordinate system are determined, and the flow parameters under a laboratory coordinate system are obtained through coordinate conversion, so that the flow field of the rotary detonation engine can be rapidly solved, the pneumatic performance of the rotary detonation engine can be rapidly evaluated, meanwhile, the flow structures such as the detonation wave, a shear layer and an induced shock wave in the flow field can be accurately captured, and the configuration design period of the rotary detonation engine is greatly shortened.
Drawings
FIG. 1 is a calculation domain under a shock wave coordinate system determined in steps (1) to (6) in a method for predicting the distribution of flow field parameters of a rotary detonation engine based on a space propulsion algorithm;
FIG. 2 is a physical schematic diagram of a rotary detonation engine;
FIG. 3 is the intersection process of the shock polar curves of two adjacent points on the initial line;
FIG. 4 is a step process from an initial line to a solution point of a space propulsion algorithm, wherein E-P represents a grid edge point, C-P represents a grid center point, Const X represents a straight line where the initial line and the solution point are located, P-M represents a Plantt-Meyer expansion wave, Sh represents a shock wave, and Sl represents a non-stick slip line;
FIG. 5 is a detonation wave parameter convergence process in a flow field of a rotary detonation engine;
FIG. 6 is a streamline grid of a rotary detonation engine under a shock wave coordinate system;
FIG. 7 is a rotational knock engine temperature field obtained by a spatial boosting algorithm;
FIG. 8 is a rotating detonation engine speed field in a laboratory coordinate system obtained by a spatial boosting algorithm;
in the figure, 11-an upstream detonation wave corresponds to an initial value line, 12-a far downstream flow field of the detonation wave corresponds to an initial value line, 13-a redundancy section of the far downstream flow field of the detonation wave corresponds to an initial value line, 14-an imaginary non-adhesive wall surface with an inclination angle consistent with the airflow deflection angle of the redundancy section, 15-an injection plane blocking section with an inclination angle equal to the inclination angle of the detonation wave, 16-a combustible mixed gas modeling air inlet section, 17-a far downstream flow field of the detonation wave corresponds to a straight line where a solution point is located, 18-a downstream detonation wave corresponds to a position where the solution point is located, and 19-an imaginary non-adhesive wall surface with an inclination angle consistent with the airflow deflection angle of the combustible mixed gas modeling air inlet section.
Detailed Description
The present invention will be further described with reference to the accompanying drawings and examples.
The invention discloses a method for predicting the distribution of flow field parameters of a rotary detonation engine based on a space propulsion algorithm, which comprises the following steps: estimating a flow parameter of the detonation wave front in advance; determining a flowing parameter after a detonation wave according to a CJ detonation model with 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 flow parameters under a shock wave coordinate system through a space propulsion algorithm; determining the initial position and the final position of combustible gas mixture injection; and (3) iteratively converging the flow field parameters of the rotary detonation engine, and obtaining the flow parameters in a laboratory coordinate system through coordinate conversion.
Referring to fig. 1, 11 is an initial value line corresponding to an upstream detonation wave, 12 is an initial value line corresponding to a far downstream flow field of the detonation wave, and 13 is an initial value line corresponding to a redundancy section of the far downstream flow field of the detonation wave. The method comprises the following steps that 14 a virtual non-adhesive wall surface with an inclination angle consistent with an airflow folding angle of a redundant section is corresponding to the virtual non-adhesive wall surface, 15 an injection plane blocking section with an inclination angle equal to an inclination angle of a detonation wave is corresponding to the injection plane blocking section, 16 a combustible mixture modeling air inlet section is corresponding to the injection plane blocking section, 17 a straight line with a far downstream flow field corresponding to the detonation wave is corresponding to the point, 18 a position with a downstream detonation wave corresponding to the point, and 19 a virtual non-adhesive wall surface with an inclination angle consistent with an airflow folding angle of a combustible mixture modeling air inlet section is corresponding to the virtual non-adhesive wall surface.
The shock wave coordinate system is a coordinate system formed by assuming that a moving detonation wave line segment AB is kept still, taking the direction of the line segment AB as a y coordinate direction and taking the direction vertical to the line segment AB as an x coordinate direction. The laboratory coordinate system is a coordinate system formed by taking a set of line segment AM directions of the injection plane as an x coordinate direction and taking a direction perpendicular to the line segment AM as a y coordinate direction.
θaveThe angle is the average flow angle and is located slightly upstream of the line segment IJ, and the angle θ is the complement of the angle between the injection plane (line segment AJ inclined in the case of a shock coordinate system) and the detonation wave AB.
The prediction method specifically comprises the following steps:
step (1), estimating the average speed of an injection surface inlet in advance, and determining the flow parameters of the detonation wave front according to a compressible flow relation;
substituting the flow parameters of the detonation wave front into a CJ detonation model with a single specific heat ratio gamma to determine the flow parameters of the detonation wave front;
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, and the flow parameters far downstream of the detonation wave in the later iteration are assigned to an initial value line according to the previous calculation result; as shown in fig. 1, the initial value line is composed of line segments AB and BD, the flow parameter after the detonation wave AB is consistent with the flow parameter after the detonation wave of the line segment IJ calculated in the previous 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 in the previous time, and the specific value is obtained by linear interpolation.
Substituting the initial value line into a space propulsion algorithm to determine a flow field structure and flow parameter distribution of the rotary detonation engine under a shock wave coordinate system;
step (5), comparing the pressure value of the solution point of the injection surface with the injection total pressure, and when the pressure value of the solution point is smaller than the injection total pressure, performing coupling solution on an injection flow field and an explosion wave downstream flow field; and if the solution point pressure value is greater than or equal to the total pressure, the injection flow field does not exist, and only the solution of the detonation wave downstream flow field is carried out at the moment.
Step (6), determining the area average value of the flow parameters of the detonation wave front at the far downstream position of the detonation wave, returning to the step (2) for iteration, and jumping out of the cycle until the average airflow angle of the far downstream position is 0;
determining the rotation angle of a shock wave coordinate system according to the size of the detonation wave inclination angle, and determining the flow parameter distribution of the rotary detonation engine in a laboratory coordinate system according to the speed triangle;
and (8) outputting the flow parameter distribution of the rotary detonation engine in the laboratory coordinate system into a data file.
With reference to fig. 3, in the present invention, the intersection point parameters of the state curves of two adjacent points are obtained by a newton-raphson iteration method. With reference to fig. 4, in the present invention, in the process of advancing from the initial 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 of additional construction of unit point process, and the streamline mesh is formed by adaptively changing the slope of the lattice edge.
The method comprises the following steps of (1):
(11) determining the total pressure Pt and the total temperature Tt of the combustible mixture at the rotary detonation inlet according to the flight condition of the rotary detonation engine, and estimating the average speed at the inlet of the rotary detonation engine
(12) Determining a detonation wave front flow parameter, including an average static temperature, from the compressional flow relationship
Mean static pressure
And average density
The following relationships exist:
wherein, the CpAnd gamma is the constant pressure specific heat and specific heat ratio of the combustible mixture and the combustion product, RgRepresenting the gas constant.
The step (2) specifically comprises the following steps:
(21) according to the CJ knock model with double gamma, the CJ knock model with a single specific heat ratio gamma is simplified and obtained:
(22) determining an expression of a flow parameter after the detonation wave as follows according to a mass conservation law, a momentum conservation law and an energy conservation law before and after the detonation wave is crossed:
wherein M isCJRepresenting the Mach number of the detonation wave, H representing the dimensionless heat release amount, Q representing the heat of chemical reaction, RgWhich represents the constant of the gas,
represents the average static temperature of the detonation wave front,
represents the average density of the detonation wave front,
represents the mean static pressure of the detonation wave front,
which represents the static pressure after the detonation wave,
which represents the density after the detonation wave,
indicating the post-detonation temperature.
The step (4) specifically comprises the following steps:
(41) determining an intermediate variable chi through the non-viscous flux E in the main flow direction, and further obtaining a state quantity Q on an initial value line, wherein the expression of the intermediate variable chi is in the following form:
E1、E2、E3、E4the main flow direction has no adhesive flux;
(42) determining downstream solution point coordinates according to the velocity distribution on the initial value line, and establishing a streamline grid, wherein the expression of the downstream solution point is in the following form:
(43) obtaining step by step half-steps along the x-direction according to an essentially oscillationless interpolation method and a limiting function
And
to obtain the state quantity of 1/2 position of step length
And
(44) according to the shock wave polar curve theory, the solution point pressure and the airflow angle of the space propulsion algorithm at the position with the stepping length of 1/2 are determined through two adjacent points with the stepping length of 1/2 in the x direction of the space, so that the Riemann self-similarity solution is obtained, and the state curves of the two adjacent points comprise the following shock wave polar curve equations:
wherein phiBA state curve representing the lower side point of the initial value line, thetaBAir flow angle, f (M), representing the lower side point of the initial value lineB,αB) Rankine allowances relation, alpha, representing the lower side point of the initial value lineBDenotes the ratio of pressure at any point in the downstream space to pressure at the lower side of the initial value line, v (M)B) Prandtl-Meyer function, phi, representing the lower side point of the initial lineTA state curve representing the upper point of the initial value line, thetaTAirflow angle, f (M), representing the upper point of the initial value lineT,αT) Rankine allowances relation, alpha, representing upper side point of initial value lineTTo representThe ratio of the pressure at any point in the downstream space to the pressure at the side point on the initial value line, v (M)T) A plantt-meier function representing the upper point of the initial line.
(45) Obtaining flow parameters under a single space progressive step length by carrying out differential dispersion on a compressible Euler equation, wherein the downstream differential equation along the flow x direction is in the following form:
wherein,
represents the difference value of y coordinates of two adjacent points of the initial value line,
the y coordinate difference of two adjacent solution points at the downstream is shown,
the inviscid flux in the y-direction of the upper grid point representing the length position of spatial step 1/2,
the slope of the upper grid edge representing the length position of spatial step 1/2,
the unbounded flux in the x-direction of the upper grid point representing the length position of spatial step 1/2,
the inviscid flux in the y-direction of the lower grid point representing the length position of spatial step 1/2,
the slope of the lower lattice edge representing the length position of spatial step 1/2,
the non-viscous flux in the x-direction of the lower grid point representing the length position of spatial step 1/2.
In step (5), the pressure value of the injection surface solution point is compared with the injection total pressure, and the solution point is the downstream point of the n +1 st line obtained by the steps (44) and (45) at the upstream point on the nth ConstX line in the flow field area. Wherein ConstX represents a straight line parallel to the detonation wave AB in fig. 1.
The step (6) comprises the following steps:
step (61), identifying an interval where the detonation wave front combustible mixed gas is located according to the flow parameter distribution of far downstream along the y direction;
step (62), determining the area average value of the flow parameters of the interval where the detonation wave front combustible mixed gas is located, wherein the area average formula is as follows:
wherein,
an x component representing the mean velocity of the detonation wave front,
y component representing mean velocity of detonation wave front, nxX component, n, representing the normal vector of the detonation wave frontyRepresents the y-component of the normal vector of the knock wave front,
the internal energy per unit mass is shown, and A is the area occupied by the detonation wave. According to a Newton-Laverson multivariate iterative method, simultaneously solving four equations of an area average formula to obtain a detonation wave front flow parameter
And
during iteration, the previous calculation result used in the assignment of the step (3) refers to the knock wave front flow parameter obtained by further calculating the knock wave front flow parameter obtained in the step (62) through the step (2).
The step (7) comprises the following steps:
(71) and correcting the inclination angle of the detonation wave relative to the injection plane according to the airflow angle of the detonation wave front, wherein the correction formula is as follows:
wherein,
indicating the inclination of the modified detonation wave relative to the injection plane,
representing the inclination of the uncorrected detonation wave relative to the injection plane;
(72) determining the rotation angle of a shock wave coordinate system according to the inclination angle of the detonation wave corrected by multiple iterations relative to the injection plane, and converting the outlet direction into the horizontal direction through a two-dimensional rotation matrix, wherein the rotation angle and the two-dimensional rotation matrix are in the following forms:
wherein,
representing the rotation angle, x, of the shock coordinate systemnewRepresenting the x coordinate and y coordinate of any point in the flow field after the rotation of the shock wave coordinate systemnewThe y coordinate and x of any point in the flow field after the shock wave coordinate system rotates are expressedoldRepresenting the x coordinate, y, of any point in the non-rotating flow field of the shock coordinate systemoldY-coordinate, u, representing any point in the non-rotating flow field of the shock coordinate systemnewRepresenting the velocity x component, v, of any point in the flow field after the rotation of the shock wave coordinate systemnewRepresenting the velocity y component u of any point in the flow field after the rotation of the shock wave coordinate systemoldRepresenting the x component, v, of the velocity at any point in the non-rotating flow field of the shock coordinate systemoldRepresenting the velocity y component of any point in the non-rotating flow field of the shock wave coordinate system;
(73) converting the speed in the shock wave coordinate system into the speed in the laboratory coordinate system by a speed triangle method according to the vector addition and subtraction principle, wherein the conversion relation is as follows:
vlab＝vnew
wherein u islabX component, v, representing the velocity at any point in the laboratory coordinate systemlabThe y component represents the velocity at any point in the laboratory coordinate system. The velocity triangle is shown in figure 1 and,
and the flow field speed in the AMLKED area under the shock wave coordinate system is represented, and the conversion between the speeds is obtained by vector addition and subtraction. Because the space position of the shock wave in the shock wave coordinate system is fixed, and the flow structures such as the shock wave in the actual flow (laboratory coordinate system) move at high speed along the tangential direction of the injection plane, the flow field of the actual rotary detonation engine is obtained through coordinate conversion according to the transformation rule of relative movement.
Further, in the step (6), the airflow angle in the upper section of the section where the knock wavefront combustible mixture is located is corrected to be 0 in the shock coordinate system of the new iteration.
The method can be used for rapidly evaluating the pneumatic performance of the rotary detonation engine, namely, the characteristics that the traditional time propulsion algorithm consumes time and computing resources are avoided by establishing a mathematical model of a detonation process and a downstream flow field structure of the detonation process.
For better illustration of the invention, and to facilitate understanding of the technical solutions thereof, typical but non-limiting examples of the invention are as follows:
the total injection pressure of combustible mixture at the inlet of the rotary detonation engine is 5 multiplied by 105Pa, total injection temperature of 300K, initial value of inlet average speed of 450m/s, specific heat ratio of the combustible premixed gas of 1.3961, and gas constant of the combustible premixed gas of 395.75J-kg-1·K-1The specific heat ratio of the combustion product is 1.1653, and the gas constant of the combustion product is 346.20J-kg-1·K-1The heat release per unit mass is 5.4704X 106J·kg-1Activation temperature of combustible gas mixture 15100K, single-step chemical reaction front factor 1.0X 109s-1. The ring expansion length of the rotary detonation engine is 314mm, and the axial length of the rotary detonation engine is 75 mm.
FIG. 5 is a graph showing the convergence course, θ, of each parameter of the detonation wave in the fast solution iteration of the present inventionincRepresenting the angle of inclination of the detonation wave, HDIndicating the length of the detonation wave, P1/P2Representing the detonation wave front wave-to-wave static pressure ratio, UlabRepresenting the knock wave velocity. FIG. 6 shows the shock coordinate system during fast resolving of the present inventionA grid of streamlines. FIG. 7 shows the temperature field of the rotary detonation engine obtained by the invention, and accurately captures the shock wave supercharging process and the post-detonation gas flow expansion acceleration process. Fig. 8 is a speed field of the rotary detonation engine obtained by the invention in a laboratory coordinate system, which is defined as a coordinate system formed by taking a line segment AM direction group of a jetting plane as an x coordinate direction and taking a direction perpendicular to the line segment AM as a y coordinate direction, and accurately capturing the airflow direction change process of the rotary detonation engine.
Table 1 shows the comparative data between the detonation wave parameters obtained by the rapid calculation and the detonation wave parameters obtained by the conventional time-marching algorithm.
TABLE 1
Compared with the detonation wave parameters of the rotary detonation engine obtained by the traditional time propulsion algorithm, the relative deviation of the height of the detonation wave obtained by the method is 1.38%, the relative deviation of the inclination angle of the detonation wave is 1.35%, the relative deviation of the front static pressure and the rear static pressure of the detonation wave is 1.05%, the relative deviation of the wave speed of the detonation wave is 0.97%, and the visible space propulsion algorithm has the same precision of resolving the rotary detonation flow field by the time propulsion algorithm.
Table 2 is data comparing the time spent in resolving the flow field of the rotary detonation engine and the time spent in resolving the flow field of the traditional time-marching algorithm in the 4-core 8-thread i7CPU environment.
TABLE 2
| Algorithm | Number and model of CPU | Time consuming |
| Space-marching |
4 core 8 threads, i7 | 5h |
| Traditional |
4 core 8 threads, i7 | 48h |
The invention discloses a method for predicting the distribution of flow field parameters of a rotary detonation engine based on a space propulsion algorithm, which uses a CPU (central processing unit) of a 4-core 8-thread i7 to carry out solving operation, and the total time consumption is 5 h. And the traditional time advance algorithm uses 4 cores and 8 threads to carry out solving operation, and the total time consumption is 48 h. Compared with the traditional time propulsion algorithm, the method has the advantage that the solving speed is 8.6 times faster, and the configuration design and the selection period of the rotary detonation engine can be effectively reduced.
The above description is only of the preferred embodiments of the present invention, and it should be noted that: it will be apparent to those skilled in the art that various modifications and adaptations can be made without departing from the principles of the invention and these are intended to be within the scope of the invention.
Claims (7)
1. A method for predicting the distribution of flow field parameters of a rotary detonation engine based on a space propulsion algorithm is characterized by comprising the following steps of: the method comprises the following steps:
(1) estimating the average speed of an inlet of an injection surface in advance, and determining the flow parameters of the detonation wave front according to a compressible flow relation;
(2) substituting the flow parameters of the detonation wave front into a CJ detonation model with a single specific heat ratio to determine the flow parameters of the detonation wave front;
(3) assigning the flow parameters after the detonation wave and the flow parameters far downstream to an initial value line; 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, and the flow parameters far downstream of the detonation wave in the later iteration are obtained according to the calculation result of the previous iteration;
(4) substituting the initial value line into a space propulsion algorithm, and determining a flow field structure and flow parameter distribution of the rotary detonation engine under a shock wave coordinate system;
(5) comparing 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, performing coupled solution on the injection flow field and the detonation wave downstream flow field, and when the pressure value of the solution point is more than or equal to the total injection pressure, only performing solution on the detonation wave downstream flow field;
(6) determining the area average value of the flow parameters of the detonation wave front at the far downstream position of the detonation wave, returning to the step (2) to the step (5) for iteration, jumping out of a loop until the average airflow angle at the far downstream is 0, and executing the step (7);
(7) and 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 rotary detonation engine under the laboratory coordinate system according to the velocity triangle.
2. A method for predicting distribution of flow field parameters of a rotary detonation engine based on a spatial propelling algorithm according to claim 1, wherein the step (1) comprises the following steps:
(11) determining the total pressure Pt and the total temperature Tt of the combustible mixture at the rotary detonation inlet according to the flight condition of the rotary detonation engine, and estimating the average speed at the inlet of the rotary detonation engine
(12) Determining the average static temperature of the detonation wave front according to the compressible flow relation
Mean static pressure
And average density
3. The method for predicting the distribution of the flow field parameters of the rotary detonation engine based on the spatial propelling algorithm is characterized in that the step (2) comprises the following steps:
(21) constructing a CJ knock model with a single specific heat ratio, wherein a model equation is in the form of:
wherein M isCJRepresenting the Mach number of the detonation wave, H representing the dimensionless heat release amount, Q representing the heat of chemical reaction, RgRepresenting a gas constant, gamma is the specific heat ratio of the combustible mixture and the combustion products;
(22) determining post-detonation wave flow parameters including post-detonation wave static pressure parameters according to the mass conservation law, momentum conservation law and energy conservation law before and after the detonation wave
Density after detonation wave
And post detonation wave temperature
4. A method for predicting distribution of flow field parameters of a rotary detonation engine based on a spatial propelling algorithm according to claim 3, wherein the step (4) comprises the following steps:
(41) determining an intermediate variable chi through the non-viscous flux E in the main flow direction, and further obtaining a state quantity Q on an initial value line, wherein the expression of the intermediate variable chi is in the following form:
(42) determining downstream solution point coordinates according to the velocity distribution on the initial value line
A streamline grid is established, and the streamline grid is established,
representing the x coordinate of the jth downstream solution point in the x direction,
representing the x coordinate of the jth initial value line point in the x direction;
(43) obtaining step by step half-steps along the x-direction according to an essentially oscillationless interpolation method and a limiting function
And
to obtain the state quantity of 1/2 position of step length
And
(44) according to the shock wave polar curve theory, the solution point pressure and the airflow angle of the space propulsion algorithm at the position 1/2 in the stepping length are determined through two adjacent points of which the stepping length in the space x direction is 1/2 set values, and then a Riemann self-similarity solution is obtained, wherein the expression of the state curve of the two adjacent points is as follows:
wherein phiBA state curve representing the lower side point of the initial value line, thetaBAir flow angle, f (M), representing the lower side point of the initial value lineB,αB) Rankine allowances relation, alpha, representing the lower side point of the initial value lineBDenotes the ratio of pressure at any point in the downstream space to pressure at the lower side of the initial value line, v (M)B) Prandtl-Meyer function, phi, representing the lower side point of the initial lineTA state curve representing the upper point of the initial value line, thetaTAirflow angle, f (M), representing the upper point of the initial value lineT,αT) Rankine allowances relation, alpha, representing upper side point of initial value lineTDenotes the ratio of pressure at any point in the downstream space to pressure at a point on the initial line, v (M)T) A planter-meier function representing an upper point of the initial value line;
(45) obtaining flow parameters under a single space progressive step length by carrying out differential dispersion on a compressible Euler equation, wherein the downstream differential equation along the flow x direction is in the following form:
wherein,
represents the difference value of y coordinates of two adjacent points of the initial value line,
the y coordinate difference of two adjacent solution points at the downstream is shown,
the inviscid flux in the y-direction of the upper grid point representing the length position of spatial step 1/2,
the slope of the upper grid edge representing the length position of spatial step 1/2,
the unbounded flux in the x-direction of the upper grid point representing the length position of spatial step 1/2,
the inviscid flux in the y-direction of the lower grid point representing the length position of spatial step 1/2,
the slope of the lower lattice edge representing the length position of spatial step 1/2,
the non-viscous flux in the x-direction of the lower grid point representing the length position of spatial step 1/2.
5. A method for predicting distribution of flow field parameters of a rotary detonation engine based on a spatial propelling algorithm according to claim 4, characterized in that the step (6) comprises the following steps:
(61) identifying the interval of the combustible mixed gas in the detonation wave front according to the flow parameter distribution of the far downstream along the y direction;
(62) determining the area average value of the flow parameters of the interval where the detonation wave front combustible mixed gas is located, wherein the area average formula is as follows:
wherein,
an x component representing the mean velocity of the detonation wave front,
y component representing mean velocity of detonation wave front, nxX component, n, representing the normal vector of the detonation wave frontyRepresents the y-component of the normal vector of the knock wave front,
expressing the internal energy contained in unit mass, simultaneously solving four equations of an area average formula according to a Newton-Laverson multivariate iterative method to obtain a detonation wave front flow parameter
And
a represents the area occupied by the detonation wave.
6. A method for predicting distribution of flow field parameters of a rotary detonation engine based on a spatial propelling algorithm according to claim 5, characterized in that the step (7) comprises the following steps:
(71) and correcting the inclination angle of the detonation wave relative to the injection plane according to the airflow angle of the detonation wave front, wherein the correction formula is as follows:
wherein,
indicating the inclination of the modified detonation wave relative to the injection plane,
representing the inclination of the uncorrected detonation wave relative to the injection plane;
(72) determining the rotation angle of a shock wave coordinate system according to the inclination angle of the detonation wave corrected by multiple iterations relative to the injection plane, and converting the outlet direction into the horizontal direction through a two-dimensional rotation matrix, wherein the rotation angle and the two-dimensional rotation matrix are in the following forms:
wherein,
representing the rotation angle, x, of the shock coordinate systemnewRepresenting the x coordinate and y coordinate of any point in the flow field after the rotation of the shock wave coordinate systemnewThe y coordinate and x of any point in the flow field after the shock wave coordinate system rotates are expressedoldRepresenting the x coordinate, y, of any point in the non-rotating flow field of the shock coordinate systemoldY-coordinate, u, representing any point in the non-rotating flow field of the shock coordinate systemnewRepresenting the velocity x component, v, of any point in the flow field after the rotation of the shock wave coordinate systemnewRepresenting the velocity y component u of any point in the flow field after the rotation of the shock wave coordinate systemoldRepresenting the x component, v, of the velocity at any point in the non-rotating flow field of the shock coordinate systemoldRepresenting the velocity y component of any point in the non-rotating flow field of the shock wave coordinate system;
(73) converting the speed in the shock wave coordinate system into the speed in the laboratory coordinate system by a speed triangle method according to the vector addition and subtraction principle, wherein the conversion relation is as follows:
vlab＝vnew
wherein u islabX component, v, representing the velocity at any point in the laboratory coordinate systemlabThe y component represents the velocity at any point in the laboratory coordinate system.
7. The method for predicting the distribution of the flow field parameters of the rotary detonation engine based on the spatial propelling algorithm according to claim 1, wherein the method comprises the following steps: and (4) correcting the airflow angle of the upper side interval of the detonation wave front combustible mixed gas in the step (6) to be 0 in a new iteration shock wave coordinate system.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111092248.1A CN113987694A (en) | 2021-09-17 | 2021-09-17 | Rotary detonation engine flow field parameter distribution prediction method based on space propulsion algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111092248.1A CN113987694A (en) | 2021-09-17 | 2021-09-17 | Rotary detonation engine flow field parameter distribution prediction method based on space propulsion algorithm |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN113987694A true CN113987694A (en) | 2022-01-28 |
Family
ID=79736023
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111092248.1A Pending CN113987694A (en) | 2021-09-17 | 2021-09-17 | Rotary detonation engine flow field parameter distribution prediction method based on space propulsion algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113987694A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117688697A (en) * | 2024-02-02 | 2024-03-12 | 中国人民解放军空军工程大学 | Design method of rotary detonation engine air inlet channel |
| CN118395639A (en) * | 2024-06-20 | 2024-07-26 | 中国人民解放军空军工程大学 | Design method of rotary detonation engine spray pipe |
-
2021
- 2021-09-17 CN CN202111092248.1A patent/CN113987694A/en active Pending
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117688697A (en) * | 2024-02-02 | 2024-03-12 | 中国人民解放军空军工程大学 | Design method of rotary detonation engine air inlet channel |
| CN117688697B (en) * | 2024-02-02 | 2024-04-26 | 中国人民解放军空军工程大学 | Design method of rotary detonation engine air inlet channel |
| CN118395639A (en) * | 2024-06-20 | 2024-07-26 | 中国人民解放军空军工程大学 | Design method of rotary detonation engine spray pipe |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107679319B (en) | Algebraic modeling method for circumferential pulsation stress term in turbine through-flow model | |
| CN112417596B (en) | Parallel grid simulation method for through-flow model of combustion chamber of aero-engine | |
| CN113987694A (en) | Rotary detonation engine flow field parameter distribution prediction method based on space propulsion algorithm | |
| Guo et al. | Flowfield structure characteristics of the hypersonic flow over a cavity: From the continuum to the transition flow regimes | |
| CN105631125A (en) | Aerodynamic-thermal-structural coupling analysis method based on reduced-order model | |
| Xianhong et al. | Aerodynamic design and numerical simulation of over-under turbine-based combined-cycle (TBCC) inlet mode transition | |
| CN111079232A (en) | Calculation method for predicting influence of rotational flow distortion air intake on performance of aircraft engine | |
| Restemeier et al. | Numerical and experimental analysis of the effect of variable blade row spacing in a subsonic axial turbine | |
| Hoopes et al. | The StreamVane method: a new way to generate swirl distortion for jet engine research | |
| CN105069221A (en) | Critical performance calculation method for supersonic speed air inlet passage optimization design | |
| Lengyel et al. | Design of a counter rotating fan-an aircraft engine technology to reduce noise and CO2-emissions | |
| Van Rooij et al. | Improving aerodynamic matching of axial compressor blading using a three-dimensional multistage inverse design method | |
| CN104091003B (en) | Finite element modeling method of large-deformation responses of flexible shell structures during basic movement | |
| Spinner et al. | A blade element theory based actuator disk methodology for modeling of fan engines in RANS simulations | |
| CN112562793B (en) | Two-step reaction model calculation method for fuel detonation combustion | |
| Hoopes | A new method for generating swirl inlet distortion for jet engine research | |
| Cummings et al. | Supersonic, turbulent flow computation and drag optimization for axisymmetric afterbodies | |
| Yuan et al. | A CFD approach to fluid dynamic optimum design of steam turbine stages with stator and rotor blades | |
| CN111523201A (en) | Internal and external flow field coupling iterative calculation method in engine reverse thrust state | |
| Sanders et al. | CFD performance predictions of a serpentine diffuser configuration in an annular cascade facility | |
| Hale et al. | A numerical simulation capability for analysis of aircraft inlet-engine compatibility | |
| Park et al. | Optimization of a 3-Stage Booster: Part 1—The Axisymmetric Multi-Disciplinary Optimization Approach to Compressor Design | |
| Kannan et al. | Coupling of compressible turbomachinery and incompressible combustor flow solvers for aerothermal applications | |
| Hall et al. | Energy efficient engine low pressure subsystem aerodynamic analysis | |
| CN103745030B (en) | A kind of eccentric drum barrel aerodynamics evaluation method with labyrinth gas seals |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |