Overlapping grid shock wave assembly disturbance domain propulsion method for aerodynamic simulation of blunt-nosed aircraft

Abstract

The present invention relates to an overlapping grid shock wave assembly disturbance domain propulsion method for blunt-nosed aircraft aerodynamic simulation, belonging to the field of computational fluid dynamics. The method addresses the adaptability issues of existing shock wave assembly methods for large deformation problems, as well as the difficulties encountered in disturbance domain propulsion methods in achieving slow / difficult convergence of strong detached shock waves in high-speed blunt-nosed aircraft aerodynamic simulation, and the easy divergence of simulated strong detached shock waves. The overlapping grid shock wave assembly disturbance domain propulsion method of the present invention effectively addresses the poor applicability of shock wave assembly methods by efficiently coupling three methods: disturbance domain propulsion, overlapping grids, and shock wave assembly, and significantly improves computational efficiency.

Description

Technical Field

[0001] The present invention relates to the technical field of computational fluid dynamics, and in particular to an overlapping grid shock wave assembly disturbance domain propulsion method for aerodynamic simulation of a blunt-nosed aircraft. Background Art

[0002] High-speed blunt-nosed aircraft hold broad prospects for military and civilian applications, and have become a key area of ​​competition among major aerospace powers. To ensure the flight performance and safety of high-speed aircraft, efficient and accurate prediction of their aerodynamic and thermal loads is essential during the design phase. Computational fluid dynamics (CFD) numerical simulation is one of the primary predictive tools.

[0003] However, there are strong detached shock waves in the flow field of high-speed blunt-nosed aircraft, and the flow parameter gradients near them are extremely large, which poses a huge challenge to the numerical accuracy and stability of CFD methods. For shock wave simulation, the existing methods mainly use the shock wave capture method. This method has strong universality and can simulate any unknown situation. However, for strong detached shock waves, the shock wave capture method is very prone to numerical oscillation and produces non-physical solutions. In contrast, the shock wave assembly method is more suitable for the simulation of strong detached shock waves. This method treats the shock wave as a mathematical discontinuity and no longer calculates the flow parameter gradients at the shock wave. It has the advantages of high accuracy in the shock wave area, good numerical stability, and low computational complexity. However, the shock wave assembly method must know the approximate shape of the shock wave in advance, which reduces its universality and limits its development. In the past decade, in order to improve the universality of the shock wave assembly method, domestic and foreign scholars have developed a new shock wave assembly method based on unstructured dynamic grids.

[0004] Significantly improving computational efficiency while maintaining solution accuracy is a key development requirement for CFD methods in aircraft design. Invention patent application number CN201810250654.8 discloses a novel acceleration technique called "Disturbance Region Update Method for Steady Compressible Flow." This method establishes a dynamic computational domain that expands as disturbances propagate and shrinks as the solution converges. By avoiding the inefficient computations associated with traditional methods, it significantly improves computational efficiency while maintaining solution accuracy. Summary of the Invention

[0005] In view of the above problems, the present invention provides an overlapping grid shock wave assembly disturbance domain propulsion method for aerodynamic simulation of blunt-nosed aircraft, which solves the adaptability problem of the shock wave assembly method in the prior art in large deformation problems and the two difficulties of the disturbance domain propulsion method in aerodynamic simulation of high-speed blunt-nosed aircraft, namely, slow / difficult convergence of strong detached shock waves and easy divergence of simulated strong detached shock waves.

[0006] The present invention provides an overlapping grid shock wave assembly disturbance domain propulsion method for aerodynamic simulation of a blunt-nosed aircraft, which specifically comprises the following steps:

[0007] Step S1: Obtain a fixed grid of the blunt-nosed vehicle flow field;

[0008] Step S2: using the perturbation domain advancing method to determine the inviscid initial flow field of the fixed grid;

[0009] Step S3: obtaining the initial position of the shock wave based on the density of each grid cell in the inviscid initial flow field of step S2;

[0010] Step S4: generating a moving grid based on the current shock wave boundary position;

[0011] Step S5: Using the fixed grid as the background grid and the moving grid as the subgrid to obtain an overlapping grid;

[0012] Step S6: obtaining the inherited flow field and the inherited viscous dynamic domain;

[0013] Step S7: Overlapping grid boundary condition processing;

[0014] Step S8: performing residual estimation of the Navier-Stokes equation in the inherited viscous dynamic domain to obtain the residual of the current conserved quantity of the overlapped grid;

[0015] Step S9: In the inherited viscous dynamic domain, time-integrating the time derivative of the residual of the current conserved quantity of the overlapping grid obtained in step S8 to obtain the updated amount of the current conserved quantity of the overlapping grid, and determining the current conserved quantity of the overlapping grid based on the updated amount of the current conserved quantity of the overlapping grid;

[0016] Step S10: Determine whether the current conservation value update value of the overlapping grid has converged. If not, proceed to step S11; if converged, proceed to step S13;

[0017] Step S11: Determine whether the shock wave boundary position on the shock wave curved surface has moved. If not, proceed to step S12; if moved, return to step S4;

[0018] Step S12: reducing the inherited sticky dynamic domain;

[0019] Step S13: Obtain the final conservation quantity of the flow field.

[0020] Optionally, step S2 specifically includes the following steps:

[0021] Step S201: obtaining an initialized fixed grid flow field;

[0022] Step S202: Under the initialization of the fixed grid flow field, the multi-layer grid cells adjacent to the wall boundary in the fixed grid are selected as the initial grid cells of the convective dynamic domain;

[0023] Step S203: Processing fixed grid boundary conditions, estimating the residual of the Euler equation in the convection calculation domain, performing time integration on the time derivative, and determining the current conservation quantity update amount;

[0024] Step S204: Determine to increase the convective dynamic domain based on the current conservation quantity update amount;

[0025] Step S205: If the convective dynamic domain is increasing, return to step S203 and perform the operation of the next grid unit; if the convective dynamic domain is not increasing, obtain the inviscid initial flow field of the fixed grid and enter step S3.

[0026] Optionally, the initial position of the shock wave is determined by the density gradient in the flow velocity direction of each grid cell in the inviscid initial flow field of the fixed grid.

[0027] Optionally, the moving grid is an interpolation unit; all interpolation units are traversed, and each interpolation unit and the grid unit with the closest grid center in the overlapping fixed grid form an overlapping grid.

[0028] Optionally, the specific steps of step S6 are:

[0029] If step S6 is entered for the first time: Inheriting the flow field means interpolating the flow field of the fixed grid at the same position for the moving grid area on the overlapping grid; Inheriting the viscous dynamic domain means copying the range covered by the convective dynamic domain on the fixed grid in step S2 to the viscous dynamic domain. The viscous dynamic domain is the area where the flow field is actually calculated when solving the Navier-Stokes equations in subsequent steps;

[0030] If it is not the first entry: inherit the inherited flow field and inherited viscous dynamic domain range of the overlapping grid in the previous step.

[0031] Optionally, the specific steps of step S12 are: traversing the boundary cells of the inherited sticky dynamic domain, and for any boundary cell, if it meets the judgment condition, removing the boundary cell from the inherited sticky dynamic domain.

[0032] Optionally, the judgment conditions are: the boundary unit solution has converged, the boundary unit is located at the upstream, the boundary unit no longer affects the solution of other units in the inherited viscous dynamic domain, and the boundary unit is no longer affected by other units in the inherited viscous dynamic domain.

[0033] Optionally, the final conserved quantities of the flow field include flow field velocity, density, temperature and pressure.

[0034] Compared with the prior art, the present invention has at least the following beneficial effects:

[0035] (1) The overlapping grid shock wave assembly disturbance domain propulsion method of the present invention effectively solves the problem of poor applicability of the shock wave assembly method by efficiently coupling the three methods of disturbance domain propulsion, overlapping grid and shock wave assembly, and can significantly improve the computational efficiency.

[0036] (2) Steps S2 and S3 of the overlapping grid shock wave assembly disturbance domain propulsion method of the present invention are responsible for predicting the initial position of the shock wave using the shock wave capture method, establishing a convective dynamic domain, and solving the Euler equation only in the convective dynamic domain. Since the convective dynamic domain only contains the wave-post region in the fixed grid, the present invention can predict the initial position of the shock wave faster than the traditional method. Steps S4 to S13 are responsible for accurately solving the flow field using the shock wave assembly method. Step S4 can simply and efficiently generate a high-quality grid in the shock wave region without the need for special encryption, thereby not only reducing the amount of calculation, but also avoiding the slow convergence caused by the over-density of the shock wave region grid in the traditional method. Steps S8 and S9 solve the Navier-Stokes equations only in the inherited viscous dynamic domain. Step S12 will reduce the inherited viscous dynamic domain as the calculation converges, thereby effectively reducing the amount of calculation for solving the flow control equations, and can significantly improve the solution accuracy and calculation efficiency of the detached strong shock wave flow field simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] The drawings are only for purposes of illustrating particular embodiments and are not to be considered limiting of the invention.

[0038] Figure 1 A flow chart of the overlapping grid shock wave assembly disturbance domain propulsion method for aerodynamic simulation of a blunt-nosed aircraft according to the present invention;

[0039] Figure 2 This is a schematic diagram of the grid and flow field for solving the supersonic blunt-nosed body flow problem in the present invention. DETAILED DESCRIPTION

[0040] In order to more clearly understand the above-mentioned objects, features and advantages of the present invention, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments can be combined with each other. In addition, the present invention can also be implemented in other ways different from those described herein. Therefore, the scope of protection of the present invention is not limited by the specific embodiments disclosed below.

[0041] A specific embodiment of the present invention, as Figure 1-Figure 2 The present invention provides an overlapping grid shock wave assembly disturbance domain propulsion method for aerodynamic simulation of a blunt-nosed aircraft, which specifically includes the following steps:

[0042] Step S1: Obtain the fixed grid and calculation conditions of the blunt-nosed vehicle flow field.

[0043] Specifically, the fixed grid includes the grid coordinates and boundary conditions of the grid cells. The computational conditions include the spatial discretization format, the time marching format, and the Courant-Friedrichs-Lewy number (CFL number).

[0044] It can be understood that the blunt-nose aircraft is a blunt-nose hypersonic aircraft.

[0045] Step S2: using the perturbation domain advancing method to determine the inviscid initial flow field of the fixed grid;

[0046] Step S201: obtaining an initialized fixed grid flow field;

[0047] Specifically, the conservation quantities of all grid cells of the fixed grid in step S1 are assigned as the incoming flow values.

[0048] Step S202: Under the initialization of the fixed grid flow field, the multi-layer grid units adjacent to the wall boundary in the fixed grid are taken as the initial grid units of the convective dynamic domain.

[0049] Preferably, the number of layers is 10 layers.

[0050] Step S203: Process the fixed grid boundary conditions, estimate the residual of the Euler equation in the convection calculation domain, perform time integration on the time derivative, and determine the current conservation quantity update amount. The specific steps are:

[0051] Step S203 - 1 : Determine fixed grid boundary conditions using a virtual grid method.

[0052] Boundary conditions include the far field / inflow boundary, supersonic outflow boundary, internal boundary, symmetric boundary, and wall boundary (including inviscid and viscous wall boundaries). The imaginary grids are labeled "-1" and "0," and the real grids are labeled "1" and "2." "-1" to "2" are four adjacent cells. The boundary conditions are assigned as follows:

[0053] 1) Far field / inflow boundary; the virtual grid of the far field / inflow boundary satisfies the incoming flow state, that is, W -1 = W 0= W ∞， in, W -1 represents the conserved quantity on the second layer of imaginary grid, W 0 represents the conserved quantity on the first layer of imaginary grid, W ∞ represents the conservation quantity corresponding to the far field / inflow boundary, and the incoming flow value of the dimensionless turbulence variable is: ∞ =0.1, k ∞ =10 -5 、ω ∞ =10k∞ , where ν ∞ represents the turbulent fluctuation field variable corresponding to the far field / inflow boundary, k ∞ represents the turbulent kinetic energy corresponding to the far field / inflow boundary, ω ∞ represents the specific dissipation rate of the far field / inflow boundary.

[0054] 2) Supersonic outflow boundary: Any point in the supersonic flow is only affected by the flow in its upstream dependent domain. The virtual grid of the supersonic outflow boundary takes the value of the real grid adjacent to the boundary, that is, W -1 = W 0= W 1, among which, W 1 indicates the conserved quantity on the first grid layer.

[0055] 3) Inner boundary: There is an inner boundary for multiple grid blocks. The virtual grid of the inner boundary takes the value of the real grid of the corresponding grid block.

[0056] 4) Symmetrical boundary and inviscid wall boundary: The symmetric boundary and inviscid wall boundary have the same boundary conditions and meet the normal non-penetration condition, that is, V=0, then the virtual grid assignment method is ρ0=ρ1, u 0= u 1-2V n w , E0=E1, the corresponding relationship between the virtual grid "-1" and the real grid "2" is the same, where V is the inversion speed, ρ0 and ρ1 represent the density values ​​corresponding to the first layer of virtual grid and the first layer of grid respectively, u 0 and u 1 represents the velocity corresponding to the first layer of virtual grid and the first layer of grid, n w represents the unit normal vector of the boundary, E0 and E1 represent the specific total energy corresponding to the first layer of imaginary grid and the first layer of grid, respectively.

[0057] 5) Viscous wall boundary: The viscous wall boundary condition satisfies the no-slip condition. Among them, the adiabatic wall boundary also satisfies the zero temperature gradient at the wall, and the isothermal wall boundary also satisfies the given temperature value T at the wall. w Therefore, the assignment method of the adiabatic wall is: ρ0=ρ1, u 0=- u 1. E0=E1, the assignment method of isothermal wall is P0=P1, u 0=- u 1. T0=2T w -T1, where P is the pressure. The correspondence between the imaginary grid "-1" and the real grid "2" is the same. P0 and P1 represent the pressure corresponding to the first layer of imaginary grid and the first layer of grid respectively. T0, T wand T1 represent the temperatures corresponding to the first virtual grid, the wall, and the first grid, respectively. In addition, for turbulent flow, since the eddy viscosity at the wall is 0, ν0=-ν1, k0=-k1; ω0=2ω w -ω1, where ν0 and ν1 represent the turbulent fluctuation field variables corresponding to the first imaginary grid and the first grid, respectively; k0 and k1 represent the turbulent kinetic energy corresponding to the first imaginary grid and the first grid, respectively; ω0, ω w and ω1 represent the specific dissipation rates corresponding to the first layer of imaginary grid, wall and first layer of grid, respectively.

[0058] Furthermore, the expression of the dimensionless specific dissipation rate ω at the wall is:

[0059] (1)

[0060] in, represents the dynamic viscosity coefficient; represents the k-ω turbulence model correction coefficient; It represents the distance between the cell center of the adjacent wall unit and the wall; represents the Reynolds number; represents the density at the center of the grid cell.

[0061] Step S203-2: On the current grid cell, use the discrete Euler equation to estimate the residual, i.e. the right side of equation (2). The expression of the discrete Euler equation is:

[0062] (2)

[0063] in, represents the current conserved quantity of the fixed grid; t Indicates time; Indicates the volume of the current grid cell; Represents the Euler equation convection flux on the current grid cell; 、 Respectively represent the number and area of ​​the cell face of the current grid cell; represents the source term of the turbulence model equation; Indicates the first unit faces.

[0064] Furthermore, the convective flux on the current grid cell The expression is:

[0065] (3)

[0066] in, p 、 H denote the pressure and total enthalpy at the center of the grid cell, respectively; , , represents the velocity component, , , Represents the out-of-plane normal component of the mesh element, where the subscripts 1, 2, and 3 correspond to the x, y, and z axes of the Cartesian coordinates, respectively; represents the inverse velocity in the normal direction outside the grid cell, ;

[0067] Step S203-3: performing time integration on the time derivative term of the conservation quantity of the discrete Euler equation to determine the current conservation quantity update amount of the fixed grid;

[0068] In the convective dynamic domain, the time derivative term in the discrete Euler equation in step S203-2 is discretized using the time marching format , determines the current conservation update amount of the fixed grid .

[0069] Step S204: Determine to increase the convective dynamic domain based on the current conservation quantity update amount;

[0070] Specifically, the current conserved quantity is updated by Determine whether all boundary grid cells on the boundary of the convective dynamic domain are disturbed; and add the possibly disturbed adjacent cells of the disturbed boundary grid cells to the convective dynamic domain.

[0071] Furthermore, the current conservation quantity update amount is used to determine whether the boundary grid unit is disturbed according to the disturbance condition;

[0072] Specifically, the perturbation condition is: if the current conservation update amount is greater than the convection threshold, the boundary grid cell is perturbed, and the expression is:

[0073] (4)

[0074] in, Indicates the newly added convection threshold.

[0075] Furthermore, the perturbation is performed along the unit vector of the boundary grid cell center pointing to the boundary grid cell point Directional propagation satisfies the unit vector When the velocity component in the direction is positive, the adjacent cells containing the grid point are added to the convective dynamic domain. The expression is:

[0076] (5)

[0077] in, represents the flow velocity vector, is the speed of sound.

[0078] S205: If the convective dynamic domain is increasing, return to step S203 and perform the operation of the next grid unit; if the convective dynamic domain is not increasing, obtain the inviscid initial flow field of the fixed grid and enter step S3.

[0079] Step S3: Based on the density of each grid unit in the inviscid initial flow field of step S2, the initial position of the shock wave is obtained.

[0080] Furthermore, the initial position of the shock wave is determined along the density gradient of each grid cell in the flow velocity direction in the inviscid initial flow field of the aircraft. The specific steps are as follows:

[0081] Calculate the density gradient of a grid cell The dot product of the flow velocity is:

[0082] (6)

[0083] in, is the 2-norm of the flow velocity vector; Represents the dot product of density gradient and flow velocity.

[0084] Furthermore, according to the density gradient of the grid cells in the inviscid initial flow field The density gradient along the flow direction is obtained by multiplying the density gradient along the flow direction by the flow velocity. , the expression is:

[0085] (7)

[0086] Furthermore, the grid cells whose density gradients meet the shock wave screening conditions are screened out to obtain the shock wave point cloud, and the shock wave point cloud is used to generate a smooth surface, which is used as the shock wave surface. The initial position of the shock wave is obtained based on the shock wave surface.

[0087] Furthermore, the shock wave screening conditions are:

[0088] (8)

[0089] in, represents a constant value greater than zero, preferably, The value of is 0.1.

[0090] Step S4: Generate a moving grid based on the current shock wave boundary position.

[0091] Specifically, the moving grid is a single-block structure grid that moves as the method steps progress. It contains one shock wave surface normal direction and two shock wave surface tangential directions. Each direction has two end faces, for a total of six end faces. In the shock wave surface normal direction, end face 1 is the shock wave surface, and end face 2 is at a distance of d The isosurface ofd Take 1 / 2 of the minimum distance between the shock wave and the wall. Tangentially to the shock wave, the two end faces are perpendicular planes at the shock wave endpoints. The number of mesh elements in each direction is calculated based on the fixed mesh element size; to ensure convergence, the number of elements in the same direction must not be less than 10.

[0092] Based on the six normal and tangential end faces of the shock wave at the shock wave boundary, an algebraic method is used to generate the internal computational domain of the moving mesh. The specific implementation can be achieved by using the transfinite interpolation method. The one-way interpolation mapping equation between every two opposite faces in the physical area of ​​the moving mesh is:

[0093] (9)

[0094] in, express The physical coordinates obtained by interpolation between the two opposite end faces of 0 and 1; express The physical coordinates obtained by interpolation between the two opposite end faces of 0 and 1; express The physical coordinates obtained by interpolation between the two opposite end faces of 0 and 1; Represent the three coordinates of the moving grid calculation space; Indicates the physical coordinates on the end face.

[0095] It can be understood that the coordinates in the computational space are equivalent to the index labels of each grid, and the corresponding physical coordinates are the real coordinates in the real flow field. There is a mapping relationship between the computational coordinates and the physical coordinates; the starting point and end point in each coordinate direction each correspond to an end face, and there are a total of 6 end faces.

[0096] The tensor product of three bilinear functions is expressed as:

[0097] (10)

[0098] in, Represents interpolation and The intersection of Represents interpolation and The intersection of Represents interpolation and The intersection of .

[0099] The cubic linear transformation is:

[0100] (11)

[0101] in, Represents interpolation , and The intersection of .

[0102] The final interpolation formula of three-dimensional transfinite interpolation is:

[0103] (12).

[0104] Substituting equations (9), (10) and (11) into equation (12), the physical coordinates of the grid points in the calculation domain inside the moving grid are obtained by calculating the coordinates based on the shock wave surface normal and tangential six end faces at the shock wave position.

[0105] Step S5: Using the fixed grid as the background grid and the moving grid as the subgrid to obtain an overlapping grid.

[0106] Specifically, all moving grids are interpolation units; all interpolation units are traversed, the grid unit number of each interpolation unit that is closest to the grid center in the fixed grid that overlaps with it is determined and recorded, and the overlapping grid is obtained according to the closest grid unit number.

[0107] When searching for the corresponding interpolation relationship of the first cell in the moving mesh, since it is difficult to estimate its corresponding cell, the cell in the fixed mesh with the smallest wall distance to the first cell in the moving mesh is used as the starting cell. The template walk method is used to search for the corresponding cell of the first cell in the moving mesh. For the remaining cells in the moving mesh, the corresponding cell of the previous cell is used as the starting cell, and the template walk method is used to search for their corresponding cells.

[0108] Step S6: obtaining the inherited flow field and the inherited viscous dynamic domain;

[0109] If this step is the first time, inheriting the flow field interpolates the flow field from the fixed grid at the same location on the moving grid area on the overlapping grid. Inheriting the viscous dynamic domain copies the range covered by the convective dynamic domain on the fixed grid in step 2 to the viscous dynamic domain. The viscous dynamic domain is the area where the flow field is actually calculated in subsequent steps when solving the Navier-Stokes equations.

[0110] If it is not the first entry: inherit the inherited flow field and inherited viscous dynamic domain range of the overlapping grid in the previous step.

[0111] Step S7: Overlapping grid boundary condition processing;

[0112] The virtual grid method is used to handle the physical boundary conditions such as the far field, inflow and outflow of the overlapping grids, as well as the artificial boundary conditions such as symmetry, internal boundary and periodic boundary. The shock wave boundary of the shock wave surface is treated with the shock wave assembly method, using the variable values ​​such as the local conservation quantity of the fixed grid as the upstream variable value, that is, the shock wave front variable value, and the shock wave RH relationship is used to calculate the variable value of the flow field downstream of the shock wave, that is, the shock wave post variable value. The overlapping boundaries are treated with conservative interpolation.

[0113] It can be understood that the variable values ​​such as the local conservation variables of the fixed grid are the local flow field variable values ​​of the fixed grid, and each fixed grid point stores the local flow field variable value, such as the conservation variable or velocity density.

[0114] For example, the variable values ​​of the flow field downstream of the shock wave are the Mach number and speed of the local flow field of the aircraft, and the density and pressure of the environment in which the aircraft is located.

[0115] Step S8: Perform Navier-Stokes equation residual estimation in the inherited viscous dynamic domain to obtain the current conservation quantity of the overlapped grid W v The residual R ( W v );

[0116] Specifically, in the inherited viscous dynamic domain obtained in step S6, the current conservation quantity of the overlapped grid is estimated by the Navier-Stokes equation (13): W v The residual R ( W v ).

[0117] Among them, the current conservation quantity of the overlapping grid is W v The residual R ( W v ) is:

[0118] (13)

[0119] in, represents the viscous flux; represents the convective flux on the grid for the Navier-Stokes equations.

[0120] Step S9: In the inherited viscous dynamic domain, the current conservation quantity of the overlapping grid obtained in step 8 is W v The time derivative of the residual is integrated to obtain the current conservation update Δ of the overlapping grid Wv , based on the current conservation quantity update Δ of the overlapping grid W v Determine the current conservation amount for overlapping meshes W v .

[0121] Specifically, in the inherited viscous dynamic domain, the residual obtained in step S8 is R ( W v ), substitute the time derivative term and the residual in Equation (14) to obtain the time derivative term, and further use the time marching format to discretize the time derivative term of Equation (14) to obtain the current conservation quantity update Δ of the overlapping grid W v , and finally the current conservation quantity update amount Δ W v Add the conservation amount of the overlapping grid in the previous iteration step to obtain the conservation amount of the current overlapping grid W v .

[0122] Furthermore, the relationship between the time derivative and the residual is expressed as:

[0123] (14)

[0124] Step S10: Determine the current conservation update value Δ of the overlapping grid W v Whether it has converged, if not, go to step S11; if it has converged, go to step S13.

[0125] Furthermore, the amount Δ is updated by the current conservation amount of the overlapping grid W v Whether the calculation is converged is determined by whether it approaches 0 or whether the residual value of the Navier-Stokes equation is small enough. If not, the process goes to step S11; if converged, the process goes to step S13.

[0126] Step S11: Determine whether the shock wave boundary position on the shock wave curved surface has moved. If not, proceed to step S12; if moved, return to step S4;

[0127] Calculate the movement speed and direction of the shock wave node located on the shock wave surface, and combine it with the time step to obtain the position of the shock wave node in the next time step. Use the displacement modulus of all shock wave surface nodes to determine whether the shock wave boundary position has moved. If not, go to step S12; if so, return to step S4.

[0128] Step S12: reducing the inherited sticky dynamic domain;

[0129] Traverse the boundary cells of the inherited sticky dynamic domain. For any boundary cell, if it meets the following four judgment conditions, the boundary cell will be removed from the inherited sticky dynamic domain. The specific judgment conditions are as follows:

[0130] (1) The boundary element solution has converged;

[0131] make Indicates a given deletion threshold, which is 10 -7 , then the boundary element solution has converged and can be described as .

[0132] (2) The boundary unit is located at the upstreammost position;

[0133] If the boundary cell is located at the most upstream, all its adjacent cells in the convective dynamic domain should satisfy:

[0134] (15)

[0135] Where, Indicates the upstream element tolerance angle of the boundary element, which is 10° for supersonic flow and 45° for subsonic flow.

[0136] (3) Boundary elements no longer affect the solution of other elements in the inherited viscous dynamic domain;

[0137] Specifically, within incompressible flows with Mach numbers less than 0.3, there is no order of convergence between upstream and downstream flows. Therefore, compressible flows with Mach numbers greater than 0.3 placed at the far downstream end of the inherited viscous dynamic domain no longer affect the solution of other elements within the inherited viscous dynamic domain. However, incompressible flows with Mach numbers less than 0.3 cannot be removed.

[0138] (4) The boundary element is no longer affected by other elements in the inherited viscous dynamic domain;

[0139] If the convergence condition is still satisfied after considering the influence of the adjacent cells on the updated conservation quantity of the boundary cell, then it can be considered that the boundary cell is no longer affected by other cells, that is, it satisfies:

[0140] (16)

[0141] (17)

[0142] Where, represents the iteration step size; The CFL number representing the time-marching format, I , J , K Respectively represent the grid directions of the boundary cells; The residual of the adjacent cells in the inherited viscous dynamic domain to the boundary cells of the inherited viscous dynamic domain is expressed as i The influence of direction, represents the change in convective flux, that is, the difference between the current step and the previous step; Indicates that the boundary cells of the inherited viscous dynamic domain are along the positive and negative Direction of the adjacent unit; subscript Indicates that the boundary cells of the inherited viscous dynamic domain are along the positive and negative The unit surface of The Jacobian matrix represents the convective flux along spectral radius of the direction; Indicates that the coordinate number is between With number The area of ​​the unit faces between them.

[0143] Return to step S7.

[0144] Step S13: Obtain the final conservation quantity of the flow field.

[0145] It can be understood that the final conserved quantity of the flow field can be converted into flow field physical quantities, such as flow field velocity, density, temperature, pressure, etc.

[0146] See also Figure 2 First, based on the fixed grid, the perturbation domain propulsion method is used to calculate the low-precision supersonic blunt-head body flow field of the shock wave capture method, and determine the initial position of the shock wave. Secondly, based on the initial position of the shock wave, a moving grid is generated and the overlapping grid assembly with the background grid is completed. The above-mentioned fixed grid is used to complete the interpolation of the flow field data at the same position of the overlapping grid, and the range covered by the convective dynamic domain in the fixed grid is copied to the inherited viscous dynamic domain of the overlapping grid. For the moving grid, its shock wave boundary is processed using the RH shock wave relationship, and the overlapping boundary is interpolated through the background grid. Subsequently, the shock wave assembly method is used to iteratively solve on the overlapping grid. As shown Figure 2 As shown, the present invention uses the shock wave assembly method to process the shock wave into a strong discontinuity. Compared with the calculation results of the traditional shock wave capture method where the shock wave thickness is much larger than the actual thickness, the final shock wave flow field obtained by the present invention is more accurate.

[0147] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by any technician familiar with this technical field within the technical scope disclosed by the present invention should be covered by the scope of protection of the present invention.

Claims

1. An overlapping grid shock wave assembly disturbance domain propulsion method for aerodynamic simulation of a blunt-nosed aircraft, characterized in that: The specific steps include: Step S1: Obtain a fixed grid of the blunt-nosed vehicle flow field; Step S2: Use the perturbation domain advancing method to determine the inviscid initial flow field of the fixed grid. The specific steps are as follows: Step S201: obtaining an initialized fixed grid flow field; Step S202: Under the initialization of the fixed grid flow field, the multi-layer grid cells adjacent to the wall boundary in the fixed grid are selected as the initial grid cells of the convective dynamic domain; Step S203: Processing fixed grid boundary conditions, estimating the residual of the Euler equation in the convection calculation domain, performing time integration on the time derivative, and determining the current conservation quantity update amount; Step S204: determining whether to increase the convective dynamic domain based on the current conservation quantity update amount; Step S205: If the convective dynamic domain is increasing, return to step S203 and perform the operation on the next grid unit; if the convective dynamic domain is not increasing, obtain the inviscid initial flow field of the fixed grid and proceed to step S3; Step S3: Based on the density of each grid cell in the inviscid initial flow field of step S2, the initial position of the shock wave is obtained. The specific steps are as follows: The density gradient of each grid cell in the flow velocity direction of the inviscid initial flow field of the aircraft is used to determine the grid cells that meet the shock wave screening conditions, and a shock wave point cloud is obtained. The shock wave point cloud is then used to generate a smooth surface, which serves as the shock wave surface. The initial shock wave position is then obtained based on the shock wave surface. Step S4: generating a moving grid based on the current shock wave boundary position; Step S5: Using the fixed grid as the background grid and the moving grid as the subgrid to obtain an overlapping grid; Step S6: obtaining the inherited flow field and the inherited viscous dynamic domain; Step S7: Overlapping grid boundary condition processing; The virtual grid method is used to handle the physical boundary conditions of the far field, inflow and outflow of the overlapping grid, as well as the artificial boundary conditions of the symmetric, internal and periodic boundaries. The shock wave boundary of the shock wave surface is handled by the shock wave assembly method, using the local conservation quantity of the fixed grid as the upstream variable value, and the shock wave RH relationship to calculate the flow field variable value downstream of the shock wave. Conservative interpolation is used to handle the overlapping boundaries. Step S8: performing residual estimation of the Navier-Stokes equation in the inherited viscous dynamic domain to obtain the residual of the current conserved quantity of the overlapped grid; Step S9: In the inherited viscous dynamic domain, time-integrating the time derivative of the residual of the current conserved quantity of the overlapping grid obtained in step S8 to obtain the updated amount of the current conserved quantity of the overlapping grid, and determining the current conserved quantity of the overlapping grid based on the updated amount of the current conserved quantity of the overlapping grid; Step S10: Determine whether the current conservation value update amount of the overlapping grid has converged. If not, proceed to step S11; if converged, proceed to step S13; Step S11: Determine whether the shock wave boundary position on the shock wave curved surface has moved. If not, proceed to step S12; if moved, return to step S4; Step S12: reducing the inherited sticky dynamic domain; Step S13: Obtain the final conservation quantity of the flow field.

2. The overlapping grid shock wave assembly disturbance domain propulsion method according to claim 1, characterized in that: The initial position of the shock wave is determined by the density gradient in the flow velocity direction of each grid cell in the inviscid initial flow field along a fixed grid.

3. The overlapping grid shock wave assembly disturbance domain propulsion method according to claim 1, characterized in that: The moving grid is the interpolation unit; all interpolation units are traversed, and each interpolation unit and the grid unit with the closest grid center in the overlapping fixed grid form an overlapping grid.

4. The overlapping grid shock wave assembly disturbance domain propulsion method according to claim 1, characterized in that: The specific steps of step S6 are: If step S6 is entered for the first time: Inheriting the flow field means interpolating the flow field of the fixed grid at the same position for the moving grid area on the overlapping grid; Inheriting the viscous dynamic domain means copying the range covered by the convective dynamic domain on the fixed grid in step S2 to the viscous dynamic domain. The viscous dynamic domain is the area where the flow field is actually calculated when solving the Navier-Stokes equations in subsequent steps; If it is not the first entry: inherit the inherited flow field and inherited viscous dynamic domain range of the overlapping grid in the previous step.

5. The overlapping grid shock wave assembly disturbance domain propulsion method according to claim 1, characterized in that: The specific steps of step S12 are: traversing the boundary cells of the inherited sticky dynamic domain, and for any boundary cell, if it meets the judgment condition, removing the boundary cell from the inherited sticky dynamic domain.

6. The overlapping grid shock wave assembly disturbance domain propulsion method according to claim 5, characterized in that: The judgment conditions are: the boundary element solution has converged, the boundary element is located at the upstream, the boundary element no longer affects the solution of other elements in the inherited viscous dynamic domain, and the boundary element is no longer affected by other elements in the inherited viscous dynamic domain.

7. The overlapping grid shock wave assembly disturbance domain propulsion method according to any one of claims 1 to 6, characterized in that: The final conserved quantities of the flow field include flow velocity, density, temperature and pressure.

Images