A fast initialization method for hypersonic boundary layer flow

Abstract

The present invention discloses a method for rapid initialization of hypersonic boundary layer flow, which belongs to the field of computational fluid dynamics and includes the following steps: firstly, based on the hypersonic inviscid flow hypothesis and the thin shock layer theory, the wall pressure distribution is estimated by the local tilt angle method, secondly, the boundary layer thickness is estimated based on the boundary layer similarity theory, and then the local shock angle is estimated based on the thin shock layer and the small tilt angle hypothesis, and further based on the adiabatic isentropic flow hypothesis, the flow parameters of the outer edge of the boundary layer are estimated by the positive / oblique shock wave relationship, and finally the flow parameter distribution in the normal direction of the boundary layer is calculated by linear interpolation. The present invention greatly improves the convergence speed of the numerical iteration process by quickly establishing a reasonable initial value of the hypersonic boundary layer flow field, and has a significant effect on improving the performance of the hypersonic CFD software in terms of computational efficiency.

Description

Technical Field

[0001] The present invention relates to the field of computational fluid dynamics, and more specifically, to a method for rapid initialization of hypersonic boundary layer flow. Background Art

[0002] In computational fluid dynamics (CFD), flow simulation is based on various grid technologies (structured grid, unstructured grid, adaptive Cartesian grid, etc.), and various numerical methods (finite volume method, finite difference method, etc.) are used to solve the Navier-Stokes (NS) equations. Regardless of the type of numerical method, the common feature is that the nonlinear NS equations are transformed into compatible linear iterative equations through a quasi-linearization process, and then various relaxation iterative algorithms are used to perform numerical calculations on the quasi-linearized equations. From a mathematical point of view, the convergence of the relaxation iterative algorithm is highly dependent on the initial value of the iteration. The appropriate initial value can reduce the number of iterations and speed up the convergence speed; conversely, the inappropriate initial value will affect the convergence speed, and even fail to obtain a converged result. Therefore, good flow initialization conditions (i.e., the initial value of the flow field) are crucial to ensure the convergence and computational efficiency of the flow field CFD calculation.

[0003] In most CFD software, the initial value of the flow field is usually set uniformly by the incoming flow conditions, that is, in the entire calculation area, the initial values ​​of all grid cells are uniformly assigned by the incoming flow conditions. Due to the complexity of hypersonic flow, its numerical simulation process usually adopts such initialization settings. A large number of numerical studies have found that with the advancement of relaxation iteration time, the shock layer in the hypersonic external flow field structure is usually gradually pushed outward from the boundary layer, that is, the shock wave is generally generated at the wall and gradually pushes outward. This is because the airflow is affected by viscosity, and the kinetic energy in the boundary layer is converted into internal energy, so the temperature and pressure at the wall surface first increase. At the same time, the disturbance caused by the normal velocity gradient causes the flow parameters such as pressure and temperature to gradually develop from the wall to the outside, so the flow field characteristics are manifested as shock waves gradually generated from the wall and pushed outward. As the thickness of the boundary layer develops and evolves, there are differences in the velocity gradient (or disturbance velocity) at different wall positions. The closer to the stagnation point area, the thicker the boundary layer, the smaller the normal disturbance velocity, and the slower the development of the shock wave. Conversely, the closer to the downstream flow area, the thinner the boundary layer, the larger the normal disturbance velocity, and the faster the development of the shock wave. Therefore, the shock wave first appears in the downstream area, and then gradually expands to the stagnation point area with time.

[0004] From the boundary layer similarity theory analysis, it can be known that the thickness of the boundary layer is proportional to the square value of the incoming flow Mach number and inversely proportional to the root mean square of the local Reynolds number. Therefore, the increase of flight speed (Mach number), atmospheric altitude and other parameters will lead to an increase in the thickness of the boundary layer, the normal velocity gradient and disturbance velocity near the wall are small, the shock wave advances slowly in the boundary layer, and the numerical iteration of the flow field converges in a long time. At the same time, if the grid in the boundary layer is dense, the disturbance velocity between the grid levels is further weakened by the reduction of the time advancement step, then the development of the shock layer will be suppressed in the boundary layer for a long time, and the advancement is extremely slow, which will eventually further extend the time required for the calculation to converge. Therefore, under the conditions of high altitude and high Mach number flight and dense grids near the wall, the convergence speed of hypersonic flow simulation is extremely slow, and it is very important to adopt necessary acceleration technology.

[0005] Based on the above analysis, the reasonable setting of the initial value of the boundary layer flow, especially the setting of the velocity gradient (or disturbance velocity), is the "hotspot" to accelerate the convergence speed of hypersonic flow field calculation and improve the efficiency of numerical simulation. It is necessary to develop corresponding flow initialization methods to accelerate the convergence speed of hypersonic flow simulation and improve the computational efficiency of hypersonic CFD software. Summary of the invention

[0006] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a method for rapid initialization of hypersonic boundary layer flow, which can accelerate the convergence speed of the numerical iteration process and improve the performance of hypersonic CFD software in terms of computational efficiency by quickly establishing reasonable initial values ​​of the hypersonic boundary layer flow.

[0007] The object of the present invention is achieved through the following solutions:

[0008] A method for rapid initialization of hypersonic boundary layer flow is provided. For a CFD calculation process of a hypersonic flow, a boundary layer flow field initialization process is added to the flow field initialization module thereof. The boundary layer flow field initialization process specifically includes the following steps:

[0009] S1: Based on the hypersonic inviscid flow hypothesis and thin shock layer theory, the wall pressure distribution is estimated and calculated using the local tilt angle method;

[0010] S2: The boundary layer thickness is estimated and calculated using the boundary layer similarity theory;

[0011] S3: Based on the assumption of thin shock layer and small inclination angle, the local shock angle corresponding to the streamline at the outer edge of the boundary layer is estimated and calculated from the local collision angle;

[0012] S4: Based on the adiabatic isentropic flow assumption, the flow parameters at the outer edge of the boundary layer are estimated and calculated using the normal / oblique shock wave relationship;

[0013] S5: Linear interpolation is used to calculate the flow parameter distribution in the normal direction within the boundary layer.

[0014] Further, in step S1, the local tilt angle method includes Newton iteration method and correction method thereof, wedge cutting / cone cutting method, Daheim-Back method and correction method thereof, shock wave-expansion wave method and Prandtl-Meyer method.

[0015] Furthermore, in step S1, the wall pressure distribution is estimated and calculated using the local tilt angle method, which specifically includes the sub-steps of: for the windward area surface element, predicting the pressure value using the Daheim-Back correction method; for the leeward area surface element, calculating the pressure value using the Prandtl-Meyer method.

[0016] Further, in step S2, the boundary layer thickness is estimated and calculated by using the boundary layer similarity theory, which specifically includes the sub-steps of: the boundary layer thickness adopts a unified value and is estimated and calculated by using the following formula:

[0017] ;

[0018] in, δ is the local boundary layer thickness, is the incoming flow Mach number, is the incoming flow Reynolds number.

[0019] Further, in step S3, the local shock wave corresponding to the outer edge streamline of the boundary layer is estimated and calculated by the local collision angle, which specifically includes the sub-steps: the local shock wave angle is determined by the local collision angle, and the estimation calculation is performed using the following formula:

[0020] ;

[0021] in, β is the shock wave angle corresponding to the streamline at the outer edge of the surface element boundary layer, is the local collision angle of the surface element, is the set minimum shock wave angle, γ is the specific heat ratio of gases.

[0022] Further, in step S4, the method of using the normal / oblique shock wave relationship to estimate and calculate the flow parameters at the outer edge of the boundary layer specifically includes the following sub-steps:

[0023] Step ①, based on the zero pressure gradient condition in the boundary layer, the wall pressure Determine the pressure value at the outer edge of the boundary layer ,Right now ;

[0024] Step ②, based on the one-dimensional adiabatic isentropic flow condition, the total pressure after the shock wave and the pressure value at the outer edge of the boundary layer Calculate the Mach number at the outer edge of the boundary layer ;

[0025] Step ③, based on the one-dimensional adiabatic isentropic flow condition, the total temperature after the shock wave is and the Mach number at the outer edge of the boundary layer Calculate the temperature value at the outer edge of the boundary layer ;

[0026] Step ④, from the Mach number at the outer edge of the boundary layer The temperature at the outer edge of the boundary layer Calculate the velocity value at the outer edge of the boundary layer ;

[0027] Step ⑤: Calculate the density value of the outer edge of the boundary layer by the state equation .

[0028] Further, in step ②, based on the one-dimensional adiabatic isentropic flow condition, the total pressure after the shock wave is and the pressure value at the outer edge of the boundary layer Calculate the Mach number at the outer edge of the boundary layer , specifically including the sub-steps of: using the oblique shock wave relationship, the total pressure after the shock wave The calculation method is:

[0029] ;

[0030] ;

[0031] In the formula, β is the shock wave angle corresponding to the streamline at the outer edge of the surface element boundary layer, is the incoming flow Mach number, γ is the gas specific heat ratio, , are the incoming flow temperature and pressure, is the total pressure before the shock wave, is the total pressure after the shock wave, is the pressure at the outer edge of the boundary layer;

[0032] In step ③, based on the one-dimensional adiabatic isentropic flow condition, the total temperature after the shock wave is and the Mach number at the outer edge of the boundary layer Calculate the temperature value at the outer edge of the boundary layer , specifically including the sub-steps:

[0033] The total temperature after the shock wave and the Mach number at the outer edge of the boundary layer Calculate the temperature value at the outer edge of the boundary layer , using the oblique shock wave relationship, the total temperature after the shock wave The calculation method is:

[0034] ;

[0035] In the formula, is the total temperature before the shock wave;

[0036] Temperature at the outer edge of the boundary layer It is expressed as:

[0037] ;

[0038] In step ④, the Mach number at the outer edge of the boundary layer The temperature at the outer edge of the boundary layer Calculate the velocity value at the outer edge of the boundary layer The specific sub-steps include:

[0039] Mach number at the outer edge of the boundary layer The temperature at the outer edge of the boundary layer Calculate the velocity value at the outer edge of the boundary layer , the specific calculation method is:

[0040] ;

[0041] in, is the air constant;

[0042] The velocity at the outer edge of the boundary layer Calculate the value of each velocity component. Assume that after the incoming flow collides with the object surface, the normal momentum is lost, the airflow moves along the tangential direction of the object surface, and the velocity vector at the outer edge of the boundary layer is calculated by the geometric relationship between the surface element normal vector and the incoming flow direction. The specific calculation method of each velocity component value is:

[0043] ;

[0044] in, , , are the component values ​​of the velocity at the outer edge of the boundary layer in the directions of each coordinate axis, is the normal vector of the surface element, The velocity vector of the incoming flow, is the tangential velocity vector at the outer edge of the boundary layer, , , are the tangential velocity vectors in x , y , z Component values ​​of the coordinate axes;

[0045] In step ⑤, the density value of the outer edge of the boundary layer is calculated by the state equation , specifically including the sub-steps:

[0046] ;

[0047] in, is the air constant, is the pressure at the outer edge of the boundary layer.

[0048] Furthermore, in step S5, the linear interpolation calculation is used to obtain the flow parameter distribution in the normal direction in the boundary layer, which specifically includes the following sub-steps:

[0049] The following formula is used for estimation calculation:

[0050] ;

[0051] in, j represents the grid cell number along the normal direction of the face element, σ is the interpolation coefficient, δ is the local boundary layer thickness, is the normal distance from the center of the grid cell to the face element, Represents a component i The quality score, ns is the total number of components in the mixed gas, Indicates the components in the incoming flow i The quality score, i =1,2,…, ns ; , , is the velocity component value at the outer edge of the boundary layer, , , are the temperature, density and pressure at the outer edge of the boundary layer, are the velocity components at the wall, are the temperature, density and pressure at the wall, are the velocity components of the grid cells, are the temperature, density and pressure of the grid cell, respectively.

[0052] Furthermore, the pressure at the wall The temperature at the wall is calculated by the local tilt angle method in S1. Calculated from the wall temperature condition, the velocity component value at the wall is It is calculated by the flow state conditions, and the density value at the wall is calculated by the state equation.

[0053] Furthermore, the pressure at the wall The temperature at the wall is calculated by the local tilt angle method in S1. Calculated from the wall temperature condition, the velocity component value at the wall is It is calculated by the flow state conditions, and the density value at the wall is calculated by the state equation, which specifically includes the following sub-steps:

[0054] For isothermal wall conditions, the wall temperature ,in is the constant parameter given by the calculation condition;

[0055] For adiabatic wall conditions, the wall temperature ,in is the temperature value of the first layer of grid cells;

[0056] For viscous flow, the wall velocity is zero, that is, ;

[0057] For inviscid flow, ,in , , are the velocity component values ​​of the first layer of grid cells respectively;

[0058] Wall density value From the equation of state Calculated, where is the air constant.

[0059] The beneficial effects of the present invention include:

[0060] The present invention greatly improves the convergence speed of the numerical iteration process by quickly establishing a reasonable initial value of the hypersonic boundary layer flow field, and has a significant effect on improving the performance of hypersonic CFD software in terms of calculation efficiency.

[0061] The method of the present invention is implemented based on a strict theoretical basis, can quickly establish reasonable and good initial values ​​of hypersonic boundary layer flow, can greatly accelerate the convergence speed of the flow field numerical iteration process and improve the comprehensive calculation efficiency of flow simulation.

[0062] The method of the present invention is applicable to general CFD calculation processes of hypersonic flows, and the algorithm is easy to transplant and simple to use. Secondly, the method of the present invention is only used in the flow field initialization process, which will not increase excessive calculation consumption and will not affect the accuracy of the final result. It is significantly helpful to improve the performance of hypersonic CFD software in terms of calculation efficiency.

[0063] The method of the present invention has wide applicability and can be used in CFD solvers such as the complete gas model, the thermodynamic single-temperature model, the thermodynamic two-temperature model, etc., and has a significant accelerated convergence effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0065] Figure 1 CFD calculation process for hypersonic flow;

[0066] Figure 2 Schematic diagram of hypersonic boundary layer flow;

[0067] Figure 3 is the wall pressure distribution curve under the convergence condition of the embodiment of the present invention;

[0068] Figure 4 is the wall heat flux distribution curve under the convergence condition of the embodiment of the present invention;

[0069] Figure 5 It is a comparison of the contour line cloud map of the flow field after 100 iterations in the embodiment of the present invention;

[0070] Figure 6 It is a comparison of the contour line cloud map of the flow field iterated 2500 steps in the embodiment of the present invention;

[0071] Figure 7 is a curve showing the variation of the residual with the number of iteration steps in an embodiment of the present invention;

[0072] Figure 8 This is a curve showing the variation of the residual with the CPU calculation time in an embodiment of the present invention. DETAILED DESCRIPTION

[0073] All features disclosed in all embodiments in this specification, or steps in all methods or processes implicitly disclosed, except for mutually exclusive features and / or steps, can be combined and / or expanded or replaced in any manner.

[0074] Explanation of terms

[0075] In the wedge / cone cutting method, the Daheim-Back method and its modified method, the shock-expansion wave method and the Prandtl-Meyer method, "-" and " / " are commonly used expression symbols in the field and have clear and distinct meanings.

[0076] The specific implementation process of the embodiment of the present invention is as follows:

[0077] like Figure 1As shown in the figure, the CFD calculation process of hypersonic flow generally includes the main steps of reading the grid, setting calculation and control conditions, flow field initialization, relaxation iteration, flow field update, post-processing output, etc. Based on this calculation process, only the flow field initialization module is modified, and the rest of the calculation process remains unchanged. The final required steady-state flow parameter distribution is obtained through time advancement and relaxation iteration. The flow field initialization process first uses unified parameters to set the initial value of the global calculation domain flow field. On this basis, the boundary layer flow field initialization process is added. Hypersonic boundary layer flow diagram Figure 2 As shown, the following steps are used to modify and establish the initial flow value in the boundary layer, specifically including:

[0078] Step S1: Based on the hypersonic inviscid flow assumption and the thin shock layer theory, the wall pressure distribution is estimated using the local tilt angle method.

[0079] The local tilt angle method includes Newton's iteration method and its correction method, wedge cutting / cone cutting method, Daheim-Bark method and its correction method, shock wave-expansion wave method, and other engineering prediction methods. In a preferred embodiment, a better combination method mode is selected here, and the Daheim-Bark correction method is used to predict the pressure value for the windward area surface element; and the Prandtl-Meyer method is used to calculate the pressure value for the leeward area surface element.

[0080] Specifically, the Daheim-Back correction method is described as:

[0081] (1)

[0082] in, is the local collision angle (the angle between the incoming flow velocity and the surface element normal vector ), Cp is the surface element pressure coefficient, is the incoming flow Mach number, γ is the gas specific heat ratio, K is the stagnation pressure coefficient after the normal shock wave, ≠20) / ( =20) is the correction factor, is the pressure coefficient after the normal shock wave (at the stagnation point), is the pressure coefficient calculated using the tangential cone method.

[0083] The correction factor is calculated as:

[0084] (2)

[0085] in is the local collision angle described by an angle value, a , n All are correction factors.

[0086] The Prandtl-Meyer method is described as:

[0087] (3)

[0088] in is the local collision angle, is the surface element pressure coefficient, is the incoming flow Mach number, γ is the specific heat ratio of gases.

[0089] Furthermore, the wall pressure is calculated from the pressure coefficient, i.e.

[0090] (4)

[0091] in, is the surface element pressure coefficient, , are the incoming flow density and velocity, , are the inflow pressure and the wall pressure respectively.

[0092] Step S2: using boundary layer similarity theory to estimate the boundary layer thickness.

[0093] The boundary layer thickness adopts a unified value, and its estimation method is expressed as:

[0094] (5)

[0095] in, δ is the local boundary layer thickness, is the incoming flow Mach number, is the incoming flow Reynolds number.

[0096] Step S3: Based on the assumption of a thin shock layer and a small inclination angle, the local shock angle corresponding to the outer edge streamline of the boundary layer is estimated from the local collision angle.

[0097] The local shock wave angle is determined by the local collision angle using the shock wave relationship, and its estimation method is expressed as:

[0098] (5)

[0099] in, β is the shock wave angle corresponding to the streamline at the outer edge of the surface element boundary layer, is the local collision angle, is the set minimum shock wave angle, γ is the specific heat ratio of gases.

[0100] Step S4: Based on the adiabatic isentropic flow assumption, the flow parameters at the outer edge of the boundary layer are estimated using the normal / oblique shock wave relationship. The method for calculating the flow parameters at the outer edge of the boundary layer specifically includes:

[0101] ① Based on the zero pressure gradient condition in the boundary layer, the wall pressure pw Determine the pressure value at the outer edge of the boundary layer ,Right now .

[0102] ② Based on the one-dimensional adiabatic isentropic flow condition, the total pressure after the shock wave and the pressure value at the outer edge of the boundary layer Calculate the Mach number at the outer edge of the boundary layer Using the oblique shock wave relationship, the total pressure after the shock wave The calculation method is:

[0103] (6)

[0104] In the formula, β is the shock wave angle corresponding to the streamline at the outer edge of the surface element boundary layer, is the incoming flow Mach number, γ is the gas specific heat ratio, , are the incoming flow temperature and pressure, , are the total pressure before and after the shock wave, respectively.

[0105] Furthermore, based on the one-dimensional adiabatic isentropic flow condition, the Mach number at the outer edge of the boundary layer is It is expressed as:

[0106] (7)

[0107] in, γ is the gas specific heat ratio, is the total pressure after the shock wave, is the pressure at the outer edge of the boundary layer.

[0108] ③ Based on the one-dimensional adiabatic isentropic flow condition, the total temperature after the shock wave and the Mach number at the outer edge of the boundary layer Calculate the temperature value at the outer edge of the boundary layer Using the oblique shock wave relationship, the total temperature after the shock wave The calculation method is:

[0109] (8)

[0110] In the formula, is the incoming flow Mach number, γ is the gas specific heat ratio, is the incoming flow temperature, , are the total temperatures before and after the shock wave, respectively.

[0111] Furthermore, based on the one-dimensional adiabatic isentropic flow condition, the temperature at the outer edge of the boundary layer is It is expressed as:

[0112] (9)

[0113] In the formula, Mae is the Mach number at the outer edge of the boundary layer, is the total temperature after the shock wave, γ is the specific heat ratio of gases.

[0114] ④ Mach number at the outer edge of the boundary layer The temperature at the outer edge of the boundary layer Calculate the velocity value at the outer edge of the boundary layer , the specific calculation method is:

[0115] (10)

[0116] in, γ is the gas specific heat ratio, is the air constant.

[0117] Furthermore, it is assumed that after the incoming flow collides with the object surface, the normal momentum is lost and the airflow moves along the tangential direction of the object surface. Therefore, the velocity vector at the outer edge of the boundary layer is estimated by the geometric relationship between the surface element normal vector and the incoming flow direction. The calculation method of each velocity component value is:

[0118] (11)

[0119] in, is the velocity at the outer edge of the boundary layer, , , are the component values ​​of the velocity at the outer edge of the boundary layer in the directions of each coordinate axis, is the normal vector of the surface element, The velocity vector of the incoming flow, is the tangential velocity vector at the outer edge of the boundary layer, , , are the tangential velocity vectors in x , y , z The component values ​​of the axis.

[0120] ⑤Finally, the density value of the outer edge of the boundary layer is calculated by the state equation ,Right now ,in is the air constant, , are the pressure and temperature at the outer edge of the boundary layer, respectively.

[0121] Step S5: linear interpolation is used to calculate the flow parameter distribution in the normal direction of the boundary layer.

[0122] like Figure 2 As shown in the figure, the flow parameters in the boundary layer are obtained by linear interpolation calculation along the normal direction of the wall, and the estimation method is:

[0123] (12)

[0124] in, j represents the grid cell number along the normal direction of the face element, σ is the interpolation coefficient, δ is the local boundary layer thickness, is the normal distance from the center of the grid cell to the face element, Represents a component i The quality score, ns is the total number of components in the mixed gas, Indicates the components in the incoming flow i The quality score, i =1,2,…, ns ; , , is the velocity component value at the outer edge of the boundary layer, , , are the temperature, density and pressure at the outer edge of the boundary layer, are the velocity components at the wall, are the temperature, density and pressure at the wall, are the velocity components of the grid cells, are the temperature, density and pressure of the grid cell, respectively.

[0125] Furthermore, the pressure at the wall The temperature at the wall is calculated by the local tilt angle method in S1. Calculated from the wall temperature condition, the velocity at the wall is The density at the wall is calculated from the flow state conditions, and the density at the wall is calculated from the state equation. Specifically, for isothermal wall conditions, the wall temperature ,in, is a constant parameter given by the calculation conditions; for adiabatic wall conditions, the wall temperature ,in is the temperature value of the first layer of grid cells. For viscous flow, the wall velocity is zero, that is, ; For inviscid flow, ,in , , are the velocity components of the first layer of grid cells respectively. Finally, the wall density is calculated by the state equation ,Right now ,in is the air constant.

[0126] In other embodiments, based on the hyper CFD software, a head spherical model of the ELECTRE blunt cone standard model shape is selected, and the method of the present invention is applied to carry out high-altitude high-speed flow simulation. The computational model and experimental conditions refer to "Muylaert J, Walpot L, Hauser J. Standard Model Testing in the European High Enthalpy Facility F4 and Extrapolation to Flight [J]. AIAA 92-3905, 1992.". The incoming flow velocity is 4230 m / s, and the flight altitude is 71 kilometers. The numerical simulation adopts the thermodynamic single-temperature 7-component Park chemical model, the convective flux selects the mixed format of Steger and AUSMpw+, the CFL number is 100, the wall condition is the isothermal complete catalytic wall condition (wall temperature 343K), and the 4-core parallel iteration step number is 20,000 steps.

[0127] Figure 3 and Figure 4 The wall pressure distribution and wall heat flux distribution obtained by the method of the present invention and the conventional method under the calculation convergence condition are respectively given. It can be seen that the pressure and heat flux curves calculated by the method of the present invention are basically consistent with the calculation results of the conventional method, indicating that the new method will not affect the final result of the flow simulation and can ensure the calculation accuracy.

[0128] Figure 5 and Figure 6 The contour line distribution of flow parameters such as flow field pressure and temperature at 100 and 2500 iteration steps are given respectively. It can be seen that the conventional calculation method needs about 2500 iterations to clearly show the flow shock layer structure, while the method of the present invention can produce an obvious flow shock layer structure within 100 steps. The method of the present invention can achieve the rapid advancement and outward expansion of the flow shock layer structure by quickly establishing the initial value of the flow, which shows that the method has a significant effect of accelerating the convergence of calculations.

[0129] Figure 7 and Figure 8The comparison of the average residual with the number of iterations and the CPU time is given respectively. It can be seen that when the same convergence standard is reached (the average residual decreases by 5 orders of magnitude and the residual value is the same), the method of the present invention is at least 6000 steps ahead, and the calculation time is reduced by about 6500 seconds. Compared with the conventional method, the number of iterations and the calculation time are reduced by about 30%, which further proves that the method of the present invention has a significant effect of accelerating the calculation convergence, and has a positive effect on improving the performance of the super CFD software in terms of calculation efficiency.

[0130] The units involved in the embodiments of the present invention may be implemented by software or hardware, and the units described may also be arranged in a processor. The names of these units do not, in some cases, limit the units themselves.

[0131] According to one aspect of an embodiment of the present invention, a computer program product or a computer program is provided, the computer program product or the computer program includes a computer instruction, and the computer instruction is stored in a computer-readable storage medium. A processor of a computer device reads the computer instruction from the computer-readable storage medium, and the processor executes the computer instruction, so that the computer device executes the method provided in the above various optional implementations.

[0132] As another aspect, an embodiment of the present invention further provides a computer-readable medium, which may be included in the electronic device described in the above embodiment; or may exist independently without being assembled into the electronic device. The above computer-readable medium carries one or more programs, and when the above one or more programs are executed by an electronic device, the electronic device implements the method described in the above embodiment.

Claims

1. A method for rapid initialization of hypersonic boundary layer flow, characterized in that: In the CFD calculation process of hypersonic flow, a boundary layer flow field initialization process is added to the flow field initialization module; the boundary layer flow field initialization process specifically includes the following steps: S1: Based on the hypersonic inviscid flow hypothesis and thin shock layer theory, the wall pressure distribution is estimated and calculated using the local tilt angle method; S2: The boundary layer thickness is estimated and calculated using the boundary layer similarity theory; S3: Based on the assumption of thin shock layer and small inclination angle, the local shock angle corresponding to the streamline at the outer edge of the boundary layer is estimated and calculated from the local collision angle; S4: Based on the adiabatic isentropic flow assumption, the flow parameters at the outer edge of the boundary layer are estimated and calculated using the normal / oblique shock wave relationship. In step S4, the estimated and calculated flow parameters at the outer edge of the boundary layer using the normal / oblique shock wave relationship specifically include the following sub-steps: Step ①, based on the zero pressure gradient condition in the boundary layer, the wall pressure Determine the pressure value at the outer edge of the boundary layer ,Right now ; Step ②, based on the one-dimensional adiabatic isentropic flow condition, the total pressure after the shock wave and the pressure value at the outer edge of the boundary layer Calculate the Mach number at the outer edge of the boundary layer ; Step ③, based on the one-dimensional adiabatic isentropic flow condition, the total temperature after the shock wave is and the Mach number at the outer edge of the boundary layer Calculate the temperature value at the outer edge of the boundary layer ; Step ④, from the Mach number at the outer edge of the boundary layer The temperature at the outer edge of the boundary layer Calculate the velocity value at the outer edge of the boundary layer ; Step ⑤: Calculate the density value of the outer edge of the boundary layer by the state equation ; In step ②, based on the one-dimensional adiabatic isentropic flow condition, the total pressure after the shock wave is and the pressure value at the outer edge of the boundary layer Calculate the Mach number at the outer edge of the boundary layer , specifically including the sub-steps of: using the oblique shock wave relationship, the total pressure after the shock wave The calculation method is: ； ； In the formula, β is the shock wave angle corresponding to the streamline at the outer edge of the surface element boundary layer, is the incoming flow Mach number, γ is the gas specific heat ratio, , are the incoming flow temperature and pressure, is the total pressure before the shock wave, is the total pressure after the shock wave, is the pressure at the outer edge of the boundary layer; In step ③, based on the one-dimensional adiabatic isentropic flow condition, the total temperature after the shock wave is and the Mach number at the outer edge of the boundary layer Calculate the temperature value at the outer edge of the boundary layer , specifically including the sub-steps: The total temperature after the shock wave and the Mach number at the outer edge of the boundary layer Calculate the temperature value at the outer edge of the boundary layer , using the oblique shock wave relationship, the total temperature after the shock wave The calculation method is: ； In the formula, is the total temperature before the shock wave; Temperature at the outer edge of the boundary layer It is expressed as: ； In step ④, the Mach number at the outer edge of the boundary layer The temperature at the outer edge of the boundary layer Calculate the velocity value at the outer edge of the boundary layer The specific sub-steps include: Mach number at the outer edge of the boundary layer The temperature at the outer edge of the boundary layer Calculate the velocity value at the outer edge of the boundary layer , the specific calculation method is: ； in, is the air constant; The velocity at the outer edge of the boundary layer Calculate the value of each velocity component. Assume that after the incoming flow collides with the object surface, the normal momentum is lost, the airflow moves along the tangential direction of the object surface, and the velocity vector at the outer edge of the boundary layer is calculated by the geometric relationship between the surface element normal vector and the incoming flow direction. The specific calculation method of each velocity component value is: ； in, , , are the component values ​​of the velocity at the outer edge of the boundary layer in the directions of each coordinate axis, is the normal vector of the surface element, The velocity vector of the incoming flow, is the tangential velocity vector at the outer edge of the boundary layer, , , are the tangential velocity vectors in x , y , z Component values ​​of the coordinate axes; In step ⑤, the density value of the outer edge of the boundary layer is calculated by the state equation , specifically including the sub-steps: ； in, is the air constant, is the pressure at the outer edge of the boundary layer; S5: Linear interpolation is used to calculate the flow parameter distribution in the normal direction within the boundary layer.

2. The method for rapid initialization of hypersonic boundary layer flow according to claim 1, characterized in that: In step S1, the local tilt angle method includes Newton's iteration method and its correction method, wedge cutting / cone cutting method, Daheim-Back method and its correction method, shock wave-expansion wave method and Prandtl-Meyer method.

3. The method for rapid initialization of hypersonic boundary layer flow according to claim 2, characterized in that: In step S1, the wall pressure distribution is estimated and calculated using the local tilt angle method, which specifically includes the following sub-steps: for the windward area element, the pressure value is predicted using the Daheim-Back correction method; for the leeward area element, the pressure value is calculated using the Prandtl-Meyer method.

4. The method for rapid initialization of hypersonic boundary layer flow according to claim 1, characterized in that: In step S2, the boundary layer thickness is estimated and calculated by using the boundary layer similarity theory, which specifically includes the sub-steps of: the boundary layer thickness adopts a unified value and is estimated and calculated by using the following formula: ； in, δ is the local boundary layer thickness, is the incoming flow Mach number, is the incoming flow Reynolds number.

5. The method for rapid initialization of hypersonic boundary layer flow according to claim 1, characterized in that: In step S3, the local shock wave corresponding to the outer edge streamline of the boundary layer is estimated and calculated by the local collision angle, which specifically includes the sub-steps: the local shock wave angle is determined by the local collision angle, and the estimation calculation is performed using the following formula: ； in, β is the shock wave angle corresponding to the streamline at the outer edge of the surface element boundary layer, is the local collision angle of the surface element, is the set minimum shock wave angle, γ is the specific heat ratio of gases.

6. The method for rapid initialization of hypersonic boundary layer flow according to claim 1, characterized in that: In step S5, the linear interpolation calculation is used to obtain the flow parameter distribution in the normal direction of the boundary layer, which specifically includes the following sub-steps: The following formula is used for estimation calculation: ； in, j represents the grid cell number along the normal direction of the face element, σ is the interpolation coefficient, δ is the local boundary layer thickness, is the normal distance from the center of the grid cell to the face element, Represents a component i The quality score, ns is the total number of components in the mixed gas, Indicates the components in the incoming flow i The quality score, i =1,2,…, ns ; , , is the velocity component value at the outer edge of the boundary layer, , , are the temperature, density and pressure at the outer edge of the boundary layer, are the velocity components at the wall, are the temperature, density and pressure at the wall, are the velocity components of the grid cells, are the temperature, density and pressure of the grid cell, respectively.

7. The method for rapid initialization of hypersonic boundary layer flow according to claim 6, characterized in that: Pressure at the wall The temperature at the wall is calculated by the local tilt angle method in S1. Calculated from the wall temperature condition, the velocity component value at the wall is It is calculated by the flow state conditions, and the density value at the wall is calculated by the state equation.

8. The method for rapid initialization of hypersonic boundary layer flow according to claim 7, characterized in that: The pressure at the wall The temperature at the wall is calculated by the local tilt angle method in S1. Calculated from the wall temperature condition, the velocity component value at the wall is It is calculated by the flow state conditions, and the density value at the wall is calculated by the state equation, which specifically includes the following sub-steps: For isothermal wall conditions, the wall temperature ,in is the constant parameter given by the calculation condition; For adiabatic wall conditions, the wall temperature ,in is the temperature value of the first layer of grid cells; For viscous flow, the wall velocity is zero, that is, ; For inviscid flow, ,in , , are the velocity component values ​​of the first layer of grid cells respectively; Wall density value From the equation of state Calculated, where is the air constant.

Images