A method and device for improving numerical shock stability of hypersonic speed

Abstract

The application relates to a kind of computing method and device for improving the numerical shock stability of hypersonic speed, belonging to the technical field of computational fluid dynamics.The method comprises: obtaining the grid information of numerical simulation of hypersonic vehicle, and initializing flow field information; for each calculation grid unit, two calculation localities are constructed along the x direction and the y direction with two adjacent grids; two localities are respectively equipped with stability matrix, and the instability condition of grid in the calculation process is judged by the positive and negative of the eigenvalue thereof; according to the instability condition of the grid unit, a Riemann solver is selected: for the grid unit not prone to instability, a high-resolution HLLC Riemann solver is used for calculation; for the grid unit prone to shock instability, a high-robustness HLL Riemann solver is used for calculation. By using different Riemann solvers for calculation at different positions, the shock instability position in the flow field can be detected, identified and processed, and the stable capture of shock is realized.

Description

Technical Field

[0001] This application relates to the field of computational fluid dynamics, and in particular to a computational method and apparatus for improving the stability of hypersonic numerical shock waves. Background Technology

[0002] Hypersonic vehicles cruise through the atmosphere for extended periods, frequently encountering complex flow field structures such as high-intensity shock waves, shock wave / shock wave interference, and shock wave / boundary layer interference. These complex flow field structures significantly impact the aerodynamic and thermal characteristics of hypersonic vehicles. Therefore, in numerical simulations of hypersonic vehicle flow fields, the quality of shock wave simulation directly affects the accuracy and reliability of aerodynamic and thermal characteristic calculations.

[0003] Currently, shock assembly and shock capture are two commonly used methods for shock simulation. Shock assembly offers advantages such as high computational accuracy and clear physical concepts, but its computational process is complex and computationally intensive, requiring prior knowledge of the approximate location of the shock. In contrast, shock capture, by applying conservation laws within discrete grid cells, can naturally obtain various discontinuous solutions, including the shock. It boasts advantages such as simplicity, clear concepts, ease of use, and automatic shock capture, making it the more widely used numerical shock processing method.

[0004] However, when applying shock wave capture methods to simulate multidimensional strong shock waves, severe shock wave instability problems often occur, such as the well-known carbuncle phenomenon. Figure 1 As shown, the occurrence of shock wave instability leads to turbulence in the flow field structure, resulting in a large deviation from the actual situation, which in turn has a serious adverse effect on the accuracy of aerodynamic and thermal characteristic calculations.

[0005] The occurrence of shock wave instability is closely related to the dissipative properties of the Riemann approximation solver. Generally, dissipative Riemann solvers (such as the HLL and HLLE Riemann solvers) have better shock wave stability and can stably capture shock waves, while low-dissipation Riemann solvers (such as the HLLC and Roe Riemann solvers) are prone to shock wave instability. However, in practical engineering applications, low-dissipation Riemann solvers are usually used to solve the problem in order to better characterize the details of the flow field. Therefore, the occurrence of shock wave instability severely restricts the application of low-dissipation Riemann solvers in solving hypersonic flow fields.

[0006] Hybrid Riemann solvers are a method that combines the different characteristics of dissipative and low-dissipative Riemann solvers for computation. This approach integrates the high robustness of dissipative Riemann solvers with the high resolution of low-dissipative Riemann solvers, enabling accurate and stable flow field calculations. A commonly used approach for hybrid Riemann solvers is domain decomposition, where the computational domain is divided. A robust dissipative Riemann solver is used near the shock wave, while a high-resolution low-dissipative Riemann solver is used in other regions. Existing techniques include matrix stability analysis, which analyzes the eigenvectors of the stability matrix to determine the location of shock wave instability. A robust Riemann solver is then used at the instability location, while a low-dissipative Riemann solver is used elsewhere to achieve stable shock wave capture. However, this method assembles the entire flow field to obtain the stability matrix and solves for its eigenvectors and eigenvalues. When the mesh size is large, the computational load is high, resulting in low efficiency. Summary of the Invention

[0007] Therefore, it is necessary to provide a calculation method and apparatus for improving the stability of hypersonic numerical shock waves in response to the above-mentioned technical problems.

[0008] A calculation method for improving the stability of hypersonic numerical shock waves, the method comprising:

[0009] Obtain the mesh information for the numerical simulation of the hypersonic vehicle, obtain all mesh elements, and initialize the flow field information to obtain the initialized flow field information.

[0010] The flow field information is reconstructed based on the grid information and the initial flow field information to obtain the state variables of the left and right sides of the interface.

[0011] Select a grid cell and construct two computational localities along the x and y directions. Assemble stability matrices for the two computational localities. Based on the stability matrices of the two computational localities, use matrix stability analysis to determine the instability of the grid cell. Select the next grid cell and continue to determine the instability until all grid cells have been traversed.

[0012] When the instability of the grid cells is not easily unstable, the first flux of the interface is calculated by a high-resolution HLLC Riemann solver based on the state variables of the interface.

[0013] When the instability of the grid cells is easy to become unstable, the second flux of the interface is calculated by the highly robust HLL Riemann solver based on the state variables of the interface.

[0014] The initial flow field information is updated based on the first and second fluxes of the interface to obtain the updated flow field information. The next time step is iterated based on the updated flow field information until the preset conditions are met, and the simulation calculation results of the hypersonic numerical shock wave stability of the hypersonic vehicle are output.

[0015] A computational device for improving the stability of hypersonic numerical shock waves, the device comprising:

[0016] The flow field information initialization module is used to obtain the mesh information of the numerical simulation of the hypersonic vehicle, obtain all mesh elements, and initialize the flow field information to obtain the initialized flow field information.

[0017] The network reconstruction module is used to reconstruct the flow field information based on the grid information and the initial flow field information to obtain the state variables of the left and right sides of the interface.

[0018] The grid cell instability judgment module is used to select grid cells and construct two computational localities along the x and y directions. Stability matrices are assembled for the two computational localities. Based on the stability matrices of the two computational localities, matrix stability analysis is used to determine the instability of the grid cell. The next grid cell is selected to continue the instability judgment until all grid cells are traversed.

[0019] The interface flux calculation module is used to calculate the first flux of the interface based on the state variables of the left and right sides of the interface when the instability of the mesh cells is not easy to instability; and to calculate the second flux of the interface based on the state variables of the left and right sides of the interface using a highly robust HLL Riemann solver when the instability of the mesh cells is easy to instability.

[0020] The calculation result output module is used to update the initial flow field information according to the first flux and the second flux of the interface to obtain the updated flow field information, and to perform the next time step iteration according to the updated flow field information until the preset conditions are met, and output the simulation calculation results of the hypersonic numerical shock wave stability of the hypersonic vehicle.

[0021] The aforementioned calculation method and apparatus for improving the stability of hypersonic numerical shock waves include: acquiring grid information from a numerical simulation of a hypersonic vehicle and initializing the flow field information; constructing two computational localities by taking two adjacent grids along the x and y directions for each computational grid cell; assembling stability matrices for each locality and determining the instability of the grid during the calculation process by the sign of its eigenvalues; when the eigenvalues ​​of the local stability matrices in both directions are less than 0, the grid cell is considered to be less prone to instability during the calculation process, and a high-resolution HLLC Riemann solver is used for calculation; when the eigenvalue of the local stability matrix in one direction is greater than 0, the grid cell is considered to be prone to shock wave instability during the calculation process, and a highly robust HLL Riemann solver is used for calculation. By using different Riemann solvers at different locations for targeted calculation, the instability location of the shock wave in the flow field can be detected, identified, and processed, achieving stable shock wave capture. Attached Figure Description

[0022] Figure 1 The present invention compares the shock wave stability and shock wave instability in the prior art, wherein (a) represents the shock wave stability and (b) represents the shock wave instability.

[0023] Figure 2 This is a flowchart illustrating a calculation method for improving the stability of hypersonic numerical shock waves in one embodiment.

[0024] Figure 3 This is a schematic diagram illustrating the lattice-centered finite volume method in another embodiment;

[0025] Figure 4 This is a schematic diagram of computational domain decomposition in another embodiment;

[0026] Figure 5 This is a schematic diagram of the grid cell instability judgment process in another embodiment;

[0027] Figure 6 This is a schematic diagram of matrix stability analysis in another embodiment.

[0028] Figure 7 Another embodiment provides a method for calculating the stability of hypersonic numerical shock waves;

[0029] Figure 8 This is the process of determining the Riemann solver based on local matrix stability analysis in another embodiment;

[0030] Figure 9 The following is a comparison of the normal shock wave results calculated using the HLLC Riemann solver and the hybrid Riemann solver of the present invention in another embodiment, wherein (a) is the calculation result using the HLLC Riemann solver and (b) is the calculation result using the hybrid Riemann solver of the present invention.

[0031] Figure 10 This is a block diagram of a computing device for improving the stability of hypersonic numerical shock waves, as shown in one embodiment. Detailed Implementation

[0032] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0033] The calculation method for improving the stability of hypersonic numerical shock waves provided in this application can be applied to the following application environments. Specifically, the terminal executes a calculation method for improving the stability of hypersonic numerical shock waves. This involves acquiring grid information from a numerical simulation of a hypersonic vehicle, obtaining all grid cells, and initializing the flow field information to obtain initialized flow field information. Based on the grid information and the initialized flow field information, the flow field information is reconstructed to obtain the state variables on the left and right sides of the interface. A grid cell is selected, and two computational localities are constructed along the x and y directions. Stability matrices are assembled for the two computational localities. Based on the stability matrices of the two computational localities, a matrix stability analysis method is used to determine the instability status of the grid cell. The next grid cell is selected to continue the instability assessment. The process continues until all mesh cells are traversed. When a mesh cell is deemed not easily unstable, the first flux of the interface is calculated using a high-resolution HLLC Riemann solver based on the state variables of the interface. When a mesh cell is deemed easily unstable, the second flux of the interface is calculated using a robust HLLC Riemann solver based on the state variables of the interface. The initial flow field information is updated based on the first and second fluxes of the interface to obtain the updated flow field information. The next time step is iterated based on the updated flow field information until a preset condition is met, and the simulation results of the hypersonic numerical shock wave stability of the hypersonic vehicle are output. The terminal can be, but is not limited to, various personal computers, laptops, and tablets.

[0034] In one embodiment, such as Figure 2 As shown, a calculation method for improving the stability of hypersonic numerical shock waves is provided, which includes the following steps:

[0035] Step 200: Obtain the mesh information of the hypersonic vehicle numerical simulation, obtain all mesh elements, and initialize the flow field information to obtain the initialized flow field information.

[0036] Specifically, such as Figure 2 As shown, the solution domain of the hypersonic vehicle flow field is first initialized to obtain the mesh information, and the flow field state variables such as velocity, temperature, pressure, and density are initialized.

[0037] Step 202: Reconstruct the flow field information based on the grid information and the initialized flow field information to obtain the state variables on the left and right sides of the interface.

[0038] Specifically, this invention employs the finite volume method with a lattice center to solve the Euler equations, such as... Figure 3 As shown, it is assumed that the flow field state variables are stored at the center point of the grid. Therefore, in order to calculate the flux at the grid interface, the flow field information needs to be reconstructed.

[0039] The left and right state variables of the interface refer to U L and U R , Subscripts K = L, R.

[0040] Step 204: Select a grid cell and construct two computational localities along the x and y directions. Assemble stability matrices for the two computational localities. Based on the stability matrices of the two computational localities, use matrix stability analysis to determine the instability of the grid cell. Select the next grid cell and continue to determine the instability until all grid cells have been traversed.

[0041] Specifically, in order to improve computational efficiency and locate the instability of the shock wave, the computational domain is decomposed, and the computational locality composed of each element is analyzed.

[0042] like Figure 4 As shown, the grid cell (i,j) and its four surrounding grid cells are divided into two computational localities along the x and y directions. The two computational localities contain grid cells {(i-1,j),(i,j),(i+1,j)} and {(i,j-1),(i,j),(i,j+1)}, respectively. Stability matrices are assembled for each computational locality. The instability of the grid cell (i,j) during computation is determined by the signs of the eigenvalues ​​of the two stability matrices. When the eigenvalues ​​of the stability matrices in both directions are less than or equal to 0, the grid cell is less prone to instability during computation. When one of the eigenvalues ​​of the stability matrix in one direction is greater than 0, the grid cell is more prone to shock wave instability during computation.

[0043] Step 206: When the instability of the grid cell is not easy to be unstable, the first flux of the interface is calculated by using a high-resolution HLLC Riemann solver based on the state variables of the interface.

[0044] Step 208: When the instability of the grid cell is easy to be unstable, the second flux of the interface is calculated by the highly robust HLL Riemann solver based on the state variables of the interface.

[0045] Step 210: Update the initial flow field information according to the first and second fluxes on the interface to obtain the updated flow field information. Iterate the next time step according to the updated flow field information until the preset conditions are met, and output the simulation calculation results of the hypersonic numerical shock wave stability of the hypersonic vehicle.

[0046] Specifically, after obtaining the interface flux, the ordinary differential equations are solved using a time-space format to update the state variables of the flow field and obtain the flow field at the next time step. The system then checks if the calculation is complete; if so, the result is output; otherwise, the loop continues.

[0047] Stable capture of shock waves is achieved by using different Riemann solvers at different locations.

[0048] The aforementioned calculation method for improving the stability of hypersonic numerical shock waves includes: acquiring grid information from the numerical simulation of the hypersonic vehicle and initializing the flow field information. For each computational grid cell, two adjacent grids along the x and y directions are used to construct two computational localities; stability matrices are assembled for each locality, and the instability of the grid during the calculation is determined by the sign of its eigenvalues; when the eigenvalues ​​of the local stability matrices in both directions are less than 0, the grid cell is considered to be less prone to instability during the calculation, and a high-resolution HLLC Riemann solver is used for calculation; when the eigenvalue of the local stability matrix in one direction is greater than 0, the grid cell is considered to be prone to shock wave instability during the calculation, and a highly robust HLL Riemann solver is used for calculation. By using different Riemann solvers at different locations for targeted calculation, the location of shock wave instability in the flow field can be detected, identified, and processed, achieving stable shock wave capture.

[0049] In one embodiment, step 200 includes: acquiring the mesh information of the hypersonic vehicle numerical simulation, obtaining all mesh cells, and initializing the flow field information to obtain the initialized flow field information, which includes velocity, temperature, pressure, and density.

[0050] In one embodiment, step 202 includes: reconstructing the flow field information using a first-order spatial precision format based on the mesh information and the initialized flow field information to obtain the state variables on the left and right sides of the interface; the expressions for the state variables on the left and right sides of the interface are:

[0051]

[0052] in, These are the state variables on the left and right sides of the i+1 / 2 grid interface, respectively, U i U is the center value of the i-th grid cell. i+1 It is the center value of the (i+1)th grid cell.

[0053] In one embodiment, such as Figure 5 As shown, step 204 specifically includes the following steps:

[0054] Step 500: Select grid cell (i,j) and divide grid cell (i,j) and its four surrounding grid cells into two computational localities according to different x and y directions. The grid cells contained in the two computational localities in the x and y directions are: {(i-1,j),(i,j),(i+1,j)} and {(i,j-1),(i,j),(i,j+1)}, respectively, where i and j are integers greater than 0.

[0055] Step 502: Calculate the local assembly stability matrices for the two matrices and calculate the eigenvalues ​​of the local stability matrices in the x and y directions.

[0056] Step 504: When the eigenvalues ​​of the computational local stability matrix in both the x and y directions are less than or equal to 0, the mesh element is less prone to instability during computation; when any eigenvalue of the computational local stability matrix in the x and y directions is greater than 0, the mesh element is prone to shock wave instability during computation. A schematic diagram of matrix stability analysis is shown below. Figure 6 As shown.

[0057] Step 506: Select the next grid cell to continue the instability assessment until all grid cells have been traversed.

[0058] In one embodiment, step 502 includes: setting all grids outside the two computational localities to be virtual grids; when analyzing the evolution of internal perturbations of grid cells, only the evolution of the initial perturbation is considered, then the stability matrices of the two computational localities are:

[0059]

[0060] Among them, S x With S y These are two matrices representing the local stability of computation; U = [ρ ρu ρv ρe] T U represents a conserved quantity. i,j Let ρ, u, v, and e represent the conserved quantities of grid cell (i,j), where ρ, u, v, and e represent density, velocity in the x-direction, velocity in the y-direction, and internal energy per unit area, respectively; the subscripts 1 to 4 of A and F represent the contact surfaces between grid cell (i,j) and the four grid cells on the right, top, left, and bottom, respectively; |Ω i,j | represents the volume of the mesh cell (i,j), which is the area in the two-dimensional case; A represents the area of ​​the contact surface, which is the length in the two-dimensional case; F is the numerical flux of the interface, F=(1-σ i,j )F HLLC +σ i,jF HLL F HLLC For the first flux, F HLL For the second flux, σ i,j For control coefficients,

[0061] Calculate the eigenvalues ​​of the two computationally local stability matrices.

[0062] Specifically, matrix stability analysis is used in two directions to analyze the evolution of perturbations within grid cell (i,j). During the analysis, all grid cells outside the local area of ​​analysis are considered virtual grids, resulting in the following:

[0063]

[0064] In the above formula, U = [ρ ρu ρv ρe] T Let ρ, u, v, and e represent the conserved quantities, respectively, and let U represent the density, velocity in the x-direction, velocity in the y-direction, and internal energy per unit volume. i,j This refers to the conserved quantity of grid cell (i,j). F is the numerical flux at the interface, calculated by the Riemann solver. Subscripts 1 to 4 represent the contact surfaces between grid cell (i,j) and the four grid cells on the right, top, left, and bottom, respectively. Ω i,j Let A represent the volume of grid cell (i,j), which is the area in the two-dimensional case; let A represent the area of ​​the contact surface, which is the length of the contact interface in the two-dimensional case. For the three grid cells in the two localities, we can obtain expressions similar to equation (3). After simplification, we can obtain the evolution relationship of the perturbation in the two localities:

[0065]

[0066] When only considering the evolution of the initial perturbation, the solution to the above equation can be expressed as:

[0067]

[0068] Where t is time, S x With S y These are two local stability matrices, and their expressions are shown in equation (2).

[0069] From the above analysis, it can be seen that the condition for local stability in the x and y directions is the stability matrix S. x With S y The largest eigenvalue is less than or equal to 0, that is:

[0070]

[0071] Where, Re(λ(S) x S represents the stability matrix in the x-direction.x The real part of the eigenvalues, λ(S x S represents the stability matrix in the x-direction. x eigenvalues; Re(λ(S) y S represents the stability matrix in the y-direction. y The real part of the eigenvalues, λ(S y S represents the stability matrix in the y-direction. y eigenvalues.

[0072] The hybrid Riemann solver for calculating the interface flux can be written as:

[0073] F=(1-σ i,j )F HLLC +σ i,j F HLL (7)

[0074] The above equation shows that for mesh elements prone to shock wave instability, the HLL Riemann solver is used to calculate the interface flux; while for stable mesh elements, the HLLC Riemann solver is used to calculate the interface flux.

[0075] In one embodiment, step 206 includes: when the instability of the mesh cells is not easily unstable, calculating the first flux of the interface using a high-resolution HLLC Riemann solver based on the state variables of the interface, wherein the expression for the first flux is:

[0076]

[0077] Among them, F HLLC The first flux; F L F R The flux is calculated using the state variables on both sides of the interface; S L and S R For the two wave velocities on the left and right sides of the interface; S * The wave speed is in the intermediate state. α L With α R As an intermediate variable, α L =ρ L (S L -q L ), α R =ρ R (S R -q R ), ρ L ρ R The densities on the left and right sides of the interface are q, respectively. L q R These are the speeds perpendicular to the left and right sides of the interface, respectively. This is the flux in the intermediate state. and intermediate state U * The state quantities on the left and right sides. Subscript K = L, R p * p represents the pressure in the intermediate state. L p R The pressure is measured on the left and right sides of the interface, respectively.

[0078] Specifically, F HLLC It was calculated using the HLLC Riemann solver. The basic idea of ​​the HLLC Riemann solver is to process the intermediate state U. * Further divided into two states and There is a contact wave between the two states. Because it uses two intermediate states to reconstruct an approximate solution to the Riemann problem, the HLLC Riemann solver has lower dissipation compared to the HLL Riemann solver. It has the characteristics of high accuracy, poor robustness, and can accurately calculate flow field information, but it is prone to shock wave instability problems.

[0079] In one embodiment, step 208 includes: when the instability of the mesh element is easily unstable, calculating the second flux of the interface using a highly robust HLL Riemann solver based on the state variables of the interface on the left and right sides, wherein the expression for the second flux is:

[0080]

[0081]

[0082] Among them, F HLL For the second flux; F L F R The flux is calculated using the state variables on both sides of the interface; S L and S R For the two wave velocities on the left and right sides of the interface, S L and S R The estimation formula for S: L =min(q) L -c L ,q R -c R ), S R =max(q) L +c L ,q R +c R ), where c L,R The speed of sound on the left and right sides of the interface; U L U RHere, q represents the state variables on both sides of the interface; q is the normal velocity of the control volume interface, q = u·n = un x +vn y ; p represents pressure, (n x ,n y ) represents the unit normal vector of the interface; ρ, u, and v represent density, velocity in the x-direction, and velocity in the y-direction, respectively; H is the total enthalpy per unit mass. The preferred specific heat ratio is γ = 1.4.

[0083] Specifically, F HLL The solution was obtained using the HLL Riemann solver, a dissipative Riemann solver characterized by high numerical dissipation, robustness, and stable shock wave capture. The HLL approximate Riemann solution neglects the contact discontinuity, specifically the two wave velocities S to the left and right of the discontinuity. L and S R Only one intermediate state U is assumed between them. * .

[0084] In one embodiment, the method uses a lattice-centered finite volume method to solve the Euler equations, where the state variables of the flow field are stored at the center point of the grid cells.

[0085] It should be understood that, although Figure 2 , Figure 5 The steps in the flowchart are shown sequentially as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified herein, there is no strict order in which these steps are executed, and they can be performed in other orders. Figure 2 , Figure 5 At least some of the steps in the process may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed in turn or alternately with other steps or at least some of the sub-steps or stages of other steps.

[0086] In another embodiment, such as Figure 7 As shown, a calculation method for improving the stability of hypersonic numerical shock waves is provided, including:

[0087] The process includes initializing flow field information, reconstructing flow field information, selecting a Riemann solver based on matrix stability analysis to calculate cross-sectional flux, updating flow field information, determining whether the calculation is complete based on the specific calculation problem, iterating if not, and outputting the calculation results if complete.

[0088] like Figure 8The process of determining the Riemann solver based on local matrix stability analysis is shown below. This includes: selecting grid cells; constructing two computational localities along the x and y directions; assembling the two computational localities to obtain the stability matrix; solving for eigenvalues; determining whether the grid is prone to instability; if the maximum eigenvalues ​​in both directions are less than or equal to 0, the grid is not prone to instability; if there are eigenvalues ​​greater than 0 in both directions, the grid is prone to instability; selecting the next grid cell for stability assessment; traversing all grids; selecting the Riemann solver based on the instability of all grids; using the KLL Riemann solver to calculate the interface flux for easily unstable grid cells; and using the HLCC Riemann solver to calculate the interface flux for less easily unstable grid cells.

[0089] The beneficial effect of this invention is that the hybrid method based on local matrix stability analysis can stably capture shock waves and avoid the occurrence of shock wave instability problems. Figure 9 The results of normal shock wave calculations obtained using the HLLC Riemann solver and the hybrid Riemann solver of this invention are compared. Figure 9 It can be observed that when using the HLLC Riemann solver, the calculation is prone to shock wave instability, while the hybrid Riemann solver proposed in this invention can stably capture the shock wave.

[0090] In one embodiment, such as Figure 10 As shown, a computational device for improving the stability of hypersonic numerical shock waves is provided, comprising: a flow field information initialization module, a network reconstruction module, a grid cell instability judgment module, an interface flux calculation module, and a calculation result output module, wherein:

[0091] The flow field information initialization module is used to obtain the mesh information of the numerical simulation of the hypersonic vehicle, obtain all mesh elements, and initialize the flow field information to obtain the initialized flow field information.

[0092] The network reconstruction module is used to reconstruct the flow field information based on the grid information and the initial flow field information to obtain the state variables on the left and right sides of the interface.

[0093] The grid cell instability judgment module is used to select grid cells and construct two computational localities along the x and y directions. Stability matrices are assembled for the two computational localities. Based on the stability matrices of the two computational localities, matrix stability analysis is used to determine the instability of the grid cell. The next grid cell is selected to continue the instability judgment until all grid cells are traversed.

[0094] The interface flux calculation module is used to calculate the first flux of the interface based on the state variables of the left and right sides of the interface when the instability of the mesh cells is not easy to instability; and to calculate the second flux of the interface based on the state variables of the left and right sides of the interface using a highly robust HLL Riemann solver when the instability of the mesh cells is easy to instability.

[0095] The calculation result output module is used to update the initial flow field information based on the first and second fluxes on the interface, obtain the updated flow field information, and perform the next time step iteration based on the updated flow field information until the preset conditions are met, and output the simulation calculation results of the hypersonic numerical shock wave stability of the hypersonic vehicle.

[0096] In one embodiment, the flow field information initialization module is also used to acquire the grid information of the hypersonic vehicle numerical simulation, obtain all grid cells, and initialize the flow field information to obtain the initialized flow field information, which includes velocity, temperature, pressure, and density.

[0097] In one embodiment, the network reconstruction module is further configured to reconstruct the flow field information using a spatial first-order precision format based on the grid information and the initialized flow field information, so as to obtain the state variables of the left and right sides of the interface; the expressions of the state variables of the left and right sides of the interface are shown in Equation (1).

[0098] In one embodiment, the grid cell instability judgment module is further configured to select grid cell (i,j), divide grid cell (i,j) and its four surrounding grid cells into two computational localities according to different x and y directions, wherein the grid cells contained in the two computational localities in the x and y directions are {(i-1,j),(i,j),(i+1,j)} and {(i,j-1),(i,j),(i,j+1)}, respectively, where i and j are integers greater than 0; assemble stability matrices for the two computational localities, and calculate the eigenvalues ​​of the stability matrices of the computational localities in the x and y directions; when the eigenvalues ​​of the stability matrices of the computational localities in the x and y directions are both less than or equal to 0, the grid cell is not prone to instability during the calculation process; when there are eigenvalues ​​of the stability matrices of the computational localities in the x and y directions that are greater than 0, the grid cell is prone to shock wave instability during the calculation process; select the next grid cell to continue the instability judgment until all grid cells are traversed.

[0099] In one embodiment, the grid cell instability judgment module is further configured to set all grids outside the two computational localities as virtual grids; when analyzing the evolution relationship of the internal disturbance of the grid cell, only the evolution process of the initial disturbance is considered, and the expression of the stability matrix of the two computational localities is shown in Equation (2). Calculate the eigenvalues ​​of the stability matrix of the two computational localities.

[0100] In one embodiment, the interface flux calculation module is further configured to calculate the first flux of the interface based on the state variables of the left and right sides of the interface using a high-resolution HLLC Riemann solver when the instability of the grid cell is not easily unstable. The expression for the first flux is shown in Equation (8).

[0101] In one embodiment, the interface flux calculation module is further configured to calculate the second flux of the interface based on the state variables of the left and right sides of the interface using a highly robust HLL Riemann solver when the instability of the grid cell is easily unstable. The expression for the second flux is shown in equations (9)-(10).

[0102] In one embodiment, the Euler equations are solved using a lattice-centered finite volume method, where the state variables of the flow field are stored at the center points of the grid cells.

[0103] Specific limitations regarding the computational device for improving the stability of hypersonic numerical shock waves can be found in the limitations of the computational methods for improving the stability of hypersonic numerical shock waves described above, and will not be repeated here. Each module in the aforementioned computational device for improving the stability of hypersonic numerical shock waves can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device in hardware form, or stored in the memory of a computer device in software form, so that the processor can call and execute the operations corresponding to each module.

[0104] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0105] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this patent application should be determined by the appended claims.

Claims

1. A calculation method for improving the stability of hypersonic numerical shock waves, characterized in that, The method includes: Obtain the mesh information for the numerical simulation of the hypersonic vehicle, obtain all mesh elements, and initialize the flow field information to obtain the initialized flow field information; The flow field information is reconstructed based on the grid information and the initial flow field information to obtain the state variables of the left and right sides of the interface; Select a grid cell and construct two computational localities along the x and y directions. Assemble stability matrices for the two computational localities. Based on the stability matrices of the two computational localities, use matrix stability analysis to determine the instability of the grid cell. Select the next grid cell and continue to determine the instability until all grid cells have been traversed. When the instability of the grid cells is not easy to be unstable, the first flux of the interface is calculated by a high-resolution HLLC Riemann solver based on the state variables of the interface. When the instability of the grid cell is easy to be unstable, the second flux of the interface is calculated by the highly robust HLL Riemann solver based on the state variables of the interface. The initial flow field information is updated based on the first and second fluxes of the interface to obtain the updated flow field information. The next time step is iterated based on the updated flow field information until the preset conditions are met, and the simulation calculation results of the hypersonic numerical shock wave stability of the hypersonic vehicle are output.

2. The method according to claim 1, characterized in that, Obtain the mesh information for the numerical simulation of the hypersonic vehicle, obtain all mesh elements, and initialize the flow field information to obtain the initialized flow field information, including: Obtain the mesh information of the numerical simulation of the hypersonic vehicle, obtain all mesh cells, and initialize the flow field information to obtain the initialized flow field information, which includes velocity, temperature, pressure, and density.

3. The method according to claim 1, characterized in that, The flow field information is reconstructed based on the mesh information and the initial flow field information to obtain the state variables on the left and right sides of the interface, including: Based on the grid information and the initialized flow field information, the flow field information is reconstructed using a first-order spatial precision format to obtain the state variables on the left and right sides of the interface; the expressions for the state variables on the left and right sides of the interface are: in, U represents the left and right state variables of the interface of the i-th grid cell, respectively. i U is the center value of the i-th grid cell. i+1 It is the center value of the (i+1)th grid cell.

4. The method according to claim 1, characterized in that, Select a grid cell and construct two computational localities along the x and y directions. Assemble stability matrices for the two computational localities and use matrix stability analysis to determine the instability of the grid cell based on the stability matrices of the two computational localities. Select the next mesh cell to continue the instability assessment until all mesh cells have been traversed, including: Select grid cell (i,j), and divide grid cell (i,j) and its four surrounding grid cells into two computational localities according to different x and y directions. The grid cells contained in the two computational localities in the x and y directions are: {(i-1,j),(i,j),(i+1,j)} and {(i,j-1),(i,j),(i,j+1)}, respectively, where i and j are integers greater than 0. For two computational local assembly stability matrices, calculate the eigenvalues ​​of the computational local stability matrices in the x and y directions; When the eigenvalues ​​of the computational local stability matrix in both the x and y directions are less than or equal to 0, the mesh element is not prone to instability during computation; when the eigenvalues ​​of the computational local stability matrix in both the x and y directions are greater than 0, the mesh element is prone to shock wave instability during computation. Select the next grid cell to continue the instability assessment until all grid cells have been traversed.

5. The method according to claim 4, characterized in that, For two computational local assembly stability matrices, calculate the eigenvalues ​​of the computational local stability matrices in the x and y directions, including: Set all grids outside the two computational localities to be virtual grids; When analyzing the evolution of internal perturbations within a grid cell, considering only the evolution of the initial perturbation, the two computational local stability matrices are: Among them, S x With S y These are two matrices representing the local stability of computation; U = [ρ ρu ρv ρe] T U represents a conserved quantity. i,j Ω represents the conserved quantities of grid cell (i,j), where ρ, u, v, and e represent density, velocity in the x-direction, velocity in the y-direction, and internal energy per unit area, respectively; the subscripts 1 to 4 of A and F represent the contact surfaces between grid cell (i,j) and the four grid cells on the right, top, left, and bottom, respectively; i,j Represents the volume or area of ​​the mesh cell (i,j); A is the area of ​​the contact surface or, in two dimensions, the length; F is the numerical flux of the interface, F=(1-σ i,j )F HLLC +σ i,j F HLL F HLLC For the first flux, F HLL For the second flux, σ i,j For control coefficients, Calculate the eigenvalues ​​of the two computationally local stability matrices.

6. The method according to claim 1, characterized in that, When the mesh element instability condition is not easily unstable, the first flux of the interface is calculated using a high-resolution HLLC Riemann solver based on the state variables on both sides of the interface, including: When the mesh element instability condition is not easily unstable, the first flux of the interface is calculated using a high-resolution HLLC Riemann solver based on the state variables on the left and right sides of the interface. The expression for the first flux is: Among them, F HLLC The first flux; F L F R The flux is calculated using the state variables on both sides of the interface; S L and S R For the two wave velocities on the left and right sides of the interface; S * The wave speed is in the intermediate state. This represents the flux in the intermediate state.

7. The method according to claim 1, characterized in that, When the mesh element instability condition is easily unstable, the second flux of the interface is calculated using a robust HLL Riemann solver based on the state variables on both sides of the interface, including: When the mesh element instability condition is easily unstable, the second flux of the interface is calculated using a robust HLL Riemann solver based on the state variables on the left and right sides of the interface. The expression for the second flux is: Among them, F HLL For the second flux, F L F R S is the flux calculated using the state variables on both sides of the interface. L and S R For the two wave velocities on the left and right sides of the interface, U L U R Let q be the state variables on both sides of the interface, q be the normal velocity of the control volume interface, H be the total enthalpy per unit mass, and p be the pressure. x ,n y ) is the unit normal vector of the interface, and ρ, u, and v represent the density, the velocity in the x-direction, and the velocity in the y-direction, respectively.

8. The method according to any one of claims 1-7, characterized in that, The Euler equations are solved using the lattice-centered finite volume method, where the state variables of the flow field are stored at the center points of the grid cells.

9. A computing device for improving the stability of hypersonic numerical shock waves, characterized in that, The device includes: The flow field information initialization module is used to acquire the mesh information of the numerical simulation of the hypersonic vehicle, obtain all mesh elements, and initialize the flow field information to obtain the initialized flow field information. The network reconstruction module is used to reconstruct the flow field information based on the grid information and the initial flow field information to obtain the state variables of the left and right sides of the interface; The grid cell instability judgment module is used to select grid cells and construct two computational localities along the x and y directions. Stability matrices are assembled for the two computational localities. Based on the stability matrices of the two computational localities, matrix stability analysis is used to determine the instability of the grid cell. The next grid cell is selected to continue the instability judgment until all grid cells are traversed. The interface flux calculation module is used to calculate the first flux of the interface based on the state variables of the left and right sides of the interface when the mesh cell is not easily unstable; and to calculate the second flux of the interface based on the state variables of the left and right sides of the interface using a high-robust HLL Riemann solver when the mesh cell is easily unstable. The calculation result output module is used to update the initial flow field information according to the first flux and the second flux of the interface to obtain the updated flow field information, and to perform the next time step iteration according to the updated flow field information until the preset conditions are met, and output the simulation calculation results of the hypersonic numerical shock wave stability of the hypersonic vehicle.

10. The apparatus according to claim 9, characterized in that, The flow field information initialization module is also used to acquire the grid information of the hypersonic vehicle numerical simulation, obtain all grid cells, and initialize the flow field information to obtain the initialized flow field information, which includes velocity, temperature, pressure, and density.

Images