A high-precision shock wave capture method based on finite difference method

Abstract

The present application relates to a high-precision shock wave capture method based on the finite difference method. The method includes: constructing a hyperbolic conservation law equation of the flow field of the aircraft to be simulated, setting a calculation interval; dividing the calculation interval by a semi-uniform grid division method, determining the solution point and the flux point; constructing a nonlinear interpolation format according to the Lagrangian interpolation format, nonlinear weight and linear weight on multiple mutually nested preset sub-templates; obtaining the left upwind value and the right upwind value at the flux point according to the nonlinear interpolation format and the preset sub-template; calculating the numerical flux of the flux point, and using a preset difference format to approximate the derivative of the flux at the solution point; performing numerical simulation by a time-advancing method based on the approximate derivative of the flux, and completing the shock wave capture of the flow field of the aircraft to be simulated. This method can capture shock waves robustly at the boundary of the computational domain, and can also well solve the problem of the numerical solution converging to an erroneous solution or the computational divergence at the interface.

Description

Technical Field

[0001] The present application relates to the technical field of computational fluid dynamics, and in particular to a high-precision shock wave capturing method based on a finite difference method. Background Art

[0002] Computational fluid dynamics (CFD) is the discipline that numerically solves fundamental equations of fluid and gas dynamics, such as the Euler / Navier-Stokes equations, on computers. Its development is closely tied to computational mathematics and computer technology. The rapid development of CFD in recent decades has fundamentally transformed the aerospace design process, effectively reducing the number of wind tunnel experiments required, overcoming the limitations of experimental techniques while also lowering energy consumption and design costs. Among the many numerical methods available, high-order finite difference methods have garnered widespread attention due to their high accuracy, computational efficiency, and ease of implementation.

[0003] However, a major challenge that has long plagued high-order finite difference methods is how to construct a boundary format capable of capturing shock waves. Specifically, this challenge stems from two requirements:

[0004] First, because fluid dynamics problems often produce discontinuous phenomena called "shock waves" during their development (for example, in scenarios like fighter jets flying at supersonic speeds or spacecraft reentering the atmosphere), numerical simulation methods must be able to accurately capture the flow field, especially shock waves generated near the physical boundaries of the computational domain (such as the surface of an aircraft wing). This requires boundary formats with the ability to capture shock waves.

[0005] On the other hand, because the shapes simulated in engineering applications are often complex (such as large passenger aircraft, spacecraft, and automobiles), multi-block grid technology is required. In a multi-block grid, each grid block only needs to accommodate a portion of the entire complex shape, making it relatively easier to generate. However, a subsequent requirement is that the numerical method must ensure that physical information can be naturally transmitted between these spliced grid blocks, especially the critical phenomenon of shock waves. Therefore, the boundary format applied to the grid block boundaries must be able to capture and transmit shock waves.

[0006] Under current technical conditions, the virtual point method is used to treat boundaries. Its purpose is to avoid boundary formats and use the same format throughout the entire computational domain. Although this method is easy to implement, it can lead to computational distortion, loss of conservation properties in the numerical solution, distortion of shock wave morphology, and phase lag when encountering shock waves. Using the virtual point method on multi-block grids can even lead to intermittent stagnation, where the shock wave cannot cross the interface, resulting in simulation failure. Therefore, it is necessary to construct a shock wave capture boundary format that does not use virtual points.

[0007] This invention designs a boundary format construction technology for the high-order precision finite difference method for the first time. The boundary format constructed using this technology can robustly capture shock waves at the boundary of the aircraft flow field and can also realize the natural transmission of shock waves on multiple grids. Summary of the Invention

[0008] Based on this, it is necessary to provide a high-precision shock wave capture method based on the finite difference method to address the above technical problems.

[0009] A high-precision shock wave capture method based on a finite difference method, the method comprising:

[0010] Construct the hyperbolic conservation law equations that control the flow field of the aircraft to be simulated and set the computational domain;

[0011] Based on the computational domain, a semi-uniform grid distribution method is used to divide the computational domain into interlaced solution points and flux points in each coordinate direction. The semi-uniform grid distribution method is as follows: the computational interval is divided into a number of cells using flux points at equal intervals in each coordinate direction, the solution points are arranged at the center of the cells, and then the two solution points closest to the boundary are replaced with non-uniform points.

[0012] Constructing a nonlinear interpolation format according to Lagrange interpolation functions on a plurality of mutually nested preset sub-templates, linear weights, and a nonlinear weight of each preset sub-template;

[0013] According to the variable value at the solution point and a plurality of preset sub-templates, nonlinear interpolation is performed using the nonlinear interpolation format to obtain the left upwind value and the right upwind value at the flux point;

[0014] According to the left upwind value and the right upwind value at the flux point, calculate the numerical flux at the corresponding flux point;

[0015] Using a preset difference format according to the numerical flux, an approximate derivative of the flux in each coordinate direction at the solution point is obtained;

[0016] According to the approximate derivatives of the flux in each coordinate direction at the solution point, the numerical simulation is carried out using the time-marching method to complete the shock wave capture of the aircraft flow field to be simulated.

[0017] The high-precision shock wave capture method based on the finite difference method constructs the hyperbolic conservation law equations for the aircraft flow field to be simulated and sets a calculation interval. A semi-uniform grid is used to determine the solution points and flux points within the calculation interval. A nonlinear interpolation format is constructed based on the Lagrangian interpolation format, nonlinear weights, and linear weights on multiple preset sub-templates. The left and right upwind values at the flux point are obtained based on the nonlinear interpolation format and multiple preset sub-templates. The numerical flux at the corresponding flux point is calculated based on the left and right upwind values at the flux point. Based on the numerical flux, a preset difference format is used to obtain the approximate derivative of the flux at the solution point. Based on the approximate derivative of the flux at the solution point, a time-marching method is used to perform numerical simulation, completing shock wave capture of the aircraft flow field to be simulated. This method can robustly capture shock waves at the boundary of the aircraft flow field calculation domain and effectively address the problem of numerical solutions converging to incorrect solutions or computational divergence at the interface. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 1 is a flow chart of a high-precision shock wave capture method based on a finite difference method in one embodiment;

[0019] Figure 2 The mutually nested preset sub-templates used for the left windward value of the interior point in one embodiment;

[0020] Figure 3 The mutually nested sub-templates used in the boundary nonlinear interpolation format in another embodiment, where (a) is the left boundary flux point x 3 / 2 The sub-template at (j=1), (b) is the left boundary flux point x 5 / 2 The sub-template at (j=2), (c) is the right boundary flux point x N-1 / 2 The sub-template at (j=N-1), (d) is the right boundary flux point x N+1 / 2 The sub-template at (j=N);

[0021] Figure 4 A schematic diagram of calculation results of a strong shock wave impacting a bubble model in another embodiment;

[0022] Figure 5 Schematic diagram of calculation results for a supersonic pipe flow problem in another embodiment, where (a) is a multi-grid result and (b) is a single-grid calculation result;

[0023] Figure 6 Calculation results of a transonic flow problem around a NACA0012 airfoil in another embodiment, where (a) is a multi-grid result and (b) is a single-grid result.

[0024] Figure 7is the surface pressure coefficient C of the NACA0012 airfoil in another embodiment. p Distribution diagram. DETAILED DESCRIPTION

[0025] In order to make the purpose, technical solutions and advantages of this application more clear, the following further describes this application in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not intended to limit this application.

[0026] Parameter description: Assume the numerical solution u j and f′ j is the exact value of u(x j ,t) and f(u(x j ,t)) are numerically approximate, and similar notations are used for the variables at the flux points.

[0027] In one embodiment, Figure 1 As shown, a high-precision shock wave capture method based on a finite difference method is provided, the method comprising the following steps:

[0028] Step 100: Construct the hyperbolic conservation law equations that control the aircraft flow field to be simulated and set the computational domain.

[0029] Specifically, the aircraft flow field to be simulated may be, but is not limited to, the flow field of a fighter jet flying at supersonic speed or the flow field of a spacecraft return capsule re-entering the atmosphere.

[0030] During supersonic flight, shock waves are generated on the wing surface of an aircraft, which requires the boundary format to have the ability to capture shock waves.

[0031] The one-dimensional hyperbolic conservation law equation in the x-axis direction is:

[0032] u t +f(u) x =0 (1)

[0033] where u = u(x, t) and f(u) represent the conserved variable and flux, respectively.

[0034] Step 102: Based on the computational domain, a semi-uniform grid distribution method is used to divide the computational domain into interlaced solution points and flux points in each coordinate direction. The semi-uniform grid distribution method means that the computational interval is divided into a number of cells using flux points at equal intervals in each coordinate direction, and the solution points are arranged at the center of the cells. Then, the two solution points closest to the boundary are replaced with non-uniform points.

[0035] Specifically, taking the format of 5th order accuracy as an example, the grid division method in the x, y, and z directions is the same. Here, the semi-uniform grid distribution method is introduced in the x-axis direction as an example. For the calculation interval [a, b], first, by using the flux point x j+1 / 2 =a+jΔ, j=0,1,2,…,N, which is divided into N units, where Δx=(ba) / N represents the spatial step in the x-axis direction. The 3rd to N-2nd solution points are usually placed in the middle of each corresponding unit, that is, x j =a+(j-1 / 2)Δx, j=3,4,…,N-2; replace the first two solution points and the last two solution points with non-uniform points.

[0036] Step 104: construct a nonlinear interpolation format according to the Lagrangian interpolation functions on the multiple mutually nested preset sub-templates, the linear weights and the nonlinear weight of each preset sub-template.

[0037] Specifically, a hierarchical structure of nested preset sub-templates is used to construct multi-resolution interpolation. A Lagrangian interpolation format is used to calculate the Lagrangian interpolation polynomial for each preset sub-template. This is then weighted summed using linear weights and the nonlinear weights of each preset sub-template to produce a nonlinear interpolation format.

[0038] Step 106: Based on the variable values at the solution point and a plurality of preset sub-templates, a nonlinear interpolation format is used to perform nonlinear interpolation to obtain the left upwind value and the right upwind value at the flux point.

[0039] Step 108: Calculate the numerical flux at the corresponding flux point according to the left upwind value and the right upwind value at the flux point.

[0040] Specifically, the calculation of the numerical flux can adopt any numerical flux function, such as the common Roe flux, HLL flux, or HLLC flux.

[0041] Step 110: Using a preset difference format according to the numerical flux, obtain the approximate derivatives of the flux in each coordinate direction at the solution point.

[0042] Step 112: Based on the approximate derivatives of the flux in each coordinate direction at the solution point, a time-marching method is used to perform numerical simulation to complete the shock wave capture of the aircraft flow field to be simulated.

[0043] Specifically, time advancement may also use any time format, such as TVD-RK3 format, LUSGS format, etc.

[0044] In the high-precision shock wave capture method based on the finite difference method, the hyperbolic conservation law equations of the aircraft flow field to be simulated are constructed and a calculation interval is set. A semi-uniform grid division method is used according to the calculation interval to determine the solution points and flux points within the calculation interval. A nonlinear interpolation format is constructed based on the Lagrangian interpolation format, nonlinear weights, and linear weights on multiple preset sub-templates. Based on the nonlinear interpolation format and multiple preset sub-templates, the left and right upwind values at the flux point are obtained. Based on the left and right upwind values at the flux point, the numerical flux at the corresponding flux point is calculated. Based on the numerical flux, a preset difference format is used to obtain the approximate derivative of the flux at the solution point. Based on the approximate derivative of the flux at the solution point, a time-marching method is used to perform numerical simulation to complete the shock wave capture of the aircraft flow field to be simulated. This method can robustly capture shock waves at the boundary of the computational domain and effectively solve the problem of numerical solutions converging to incorrect solutions or computational divergence at the interface.

[0045] This method provides a general framework for designing boundary schemes with shock-capturing capabilities for high-order finite-difference schemes. Without using virtual points, the constructed scheme can capture shock waves throughout the entire computational domain, including the boundary region, and can also transfer shock waves across the boundaries between adjacent grid blocks.

[0046] In one embodiment, step 102 includes: in the x-axis direction, setting the number of solution points in the x-axis direction to N, the calculation interval to [a, b], and using N+1 flux points x j+1 / 2 =a+jΔx divides the calculation interval in the x-axis direction into N units; where N is an integer greater than 1, a and b are real numbers and a<b, and Δx=(ba) / N is the spatial step in the x-axis direction; for the 2p-1 order format, the solution points of the 2p-2 units close to the left and right boundaries of the calculation interval are non-uniformly distributed, and the solution points of the other N-2p+2 units are arranged at the center of the unit; where p is an integer greater than or equal to 1; in the y and z axis directions, the same grid division method as that in the x-axis direction is adopted to divide the calculation intervals in the y and z axis directions into interlaced solution points and flux points, respectively.

[0047] Specifically, taking the 5th-order precision interpolation format as an example, in the x-axis direction, the number of solution points in the x-axis direction is set to N, and the calculation interval in the x-axis direction is divided into N units according to the number of solution points using flux points at equal intervals; the solution points of the four units near the left and right boundaries of the calculation interval are non-uniformly distributed, and the solution points of the other N-4 units are arranged at the center of the unit; the solution points in the x-axis direction are:

[0048]

[0049] Among them, [a,b] is the calculation interval in the x-axis direction, j is the unit number of the calculation interval in the x-axis direction, and x j is the coordinate of the solution point in the jth unit in the x-axis direction, and Δx=(ba) / N is the spatial step size in the x-axis direction.

[0050] The grid division method in the y- and z-axis directions is the same as that in the x-axis direction. Therefore, the expressions of the solution points in the y- and z-axis directions are the same as those in the x-axis direction.

[0051] In one embodiment, when an interpolation format with 2p-1 order precision is selected, the number of preset sub-templates is p, and the number of solution points included in the m-th preset sub-template is 2m-1, where m=1, 2, ..., p, and p is an integer greater than or equal to 1; step 104 includes the following steps:

[0052] Step 200: Calculate the Lagrangian interpolation function based on each preset sub-template in the x-axis direction to obtain the Lagrangian interpolation format on the corresponding preset sub-template; the Lagrangian interpolation format on the preset sub-template in the x-axis direction is:

[0053] q k =φ k (x j+1 / 2 ),k=1,2,…,p

[0054] Among them, φ k (x) is the Lagrange interpolation function on the kth preset sub-template in the x-axis direction, φ k (x j+1 / 2 ) is the Lagrange interpolation function φ on the kth preset sub-template in the x-axis direction k (x) in x j+1 / 2 The general interpolation formula on the kth preset sub-template in the x-axis direction is:

[0055]

[0056] In the above formula are the Lagrange interpolation coefficients on the 2nd, 3rd, ..., p preset sub-templates in the x-axis direction, They are the variable values at the solution points contained in the 1st, 2nd, 3rd, ..., pth preset sub-templates in the x-axis direction respectively.

[0057] Step 202: The Lagrangian interpolation on the p preset sub-templates is weightedly combined using the linear weight and nonlinear weight of each preset sub-template in the x-axis direction to obtain a nonlinear interpolation format in the x-axis direction. The nonlinear interpolation format in the x-axis direction is:

[0058]

[0059] in, is the interpolation value at the j-th flux point in the x-axis direction, ω k is the nonlinear weight of the kth preset sub-template in the x-axis direction, γ k is the linear weight of the kth preset sub-template in the x-axis direction, k = 1, 2, ..., p. Linear weight γ k The value of is determined based on actual engineering design problems.

[0060] Step 204: In the y and z axis directions, a nonlinear interpolation format in the y and z axis directions is constructed based on the Lagrange interpolation function, linear weights and nonlinear weights of each preset sub-template on the multiple nested preset sub-templates in the y and z axis directions.

[0061] Steps 200 and 202 are an explanation of the method for constructing the nonlinear interpolation format in the x-axis direction. For the y- and z-axis directions, the construction method of the nonlinear interpolation format is the same as that of the x-axis direction. The expression of the nonlinear interpolation format is the same, except that the preset sub-template in the x-axis direction is replaced by the preset sub-template in the y- and z-axis directions.

[0062] Improvements to the nonlinear weighting strategy have been made to suppress oscillations and enhance shock wave capture capabilities. Furthermore, the aforementioned boundary format naturally provides a conservative and easy-to-implement approach to the interface coupling problem between adjacent grid blocks in multi-block grids. This significantly alleviates the smoothness constraints imposed when generating multi-block grids, thereby significantly improving the practicality of high-precision finite-difference methods in complex engineering applications.

[0063] In one embodiment, the following steps are included before step 200: a method for determining the nonlinear weight of the preset sub-template in the x-, y-, and z-axis directions; wherein the method for determining the nonlinear weight of the preset sub-template in the x-axis direction specifically includes:

[0064] Step 300: Obtain smooth indicators on corresponding preset sub-templates according to the Lagrangian interpolation functions of the second to the p-th preset sub-templates in the x-axis direction.

[0065] Step 302: Determine a smooth indicator on the first preset sub-template in the x-axis direction based on the solution point in the first preset sub-template in the x-axis direction, the variable values at two adjacent solution points, and the spatial step length in the x-axis direction; the smooth indicator on the first preset sub-template in the x-axis direction is:

[0066]

[0067] in,

[0068] Step 304: Obtain an unnormalized nonlinear weight value corresponding to the preset sub-template in the x-axis direction according to the smoothness indicator, linear weight, and global smoothness metric factor on the preset sub-template in the x-axis direction.

[0069] Step 306: Normalize the unnormalized nonlinear weight values corresponding to all preset sub-templates to obtain the nonlinear weight corresponding to each preset sub-template in the x-axis direction; the expression of the nonlinear weight in the x-axis direction is:

[0070]

[0071] in, is the unnormalized nonlinear weight value, γ k is a linear weight, preferably, k=1,2,…,p, τ is the global smoothness factor, ε is a constant, preferably, ε=10 -10 .

[0072] Step 308: Replace the preset sub-template in the x-axis direction with the preset sub-templates in the y- and z-axis directions, and use the smooth indicator determination method on the preset sub-template in the x-axis direction to obtain the nonlinear weights of the preset sub-templates in the y- and z-axis directions.

[0073] In one embodiment, in step 300, the smooth indicator β on the kth preset sub-template in the x-axis direction is k for:

[0074]

[0075] Where λ is a positive constant, The solution point x of the jth unit in the x-axis direction of the Lagrange interpolation function on the kth preset sub-template j The nth derivative at x j is the solution point of the jth element in the x-axis direction.

[0076] In one embodiment, the expression of the global smoothness metric factor in the x-axis direction in step 304 is:

[0077]

[0078] Among them, β a is the smooth indicator on the a-th preset sub-template in the x-axis direction, a=1,2,…,p.

[0079] In one embodiment, step 106 includes: taking the flux points with calculation unit number j satisfying p≤j≤N-p+1 as interior points, and the remaining flux points as boundary points; when p≤j≤N-p+1, setting the first preset sub-template as the single-point sub-template S I ={u j}, set the second preset sub-template to the three-point sub-template S II ={u j-1 ,u j ,u j+1}, set the third preset sub-template to five-point sub-template S III ={u j-2 ,u j-1 ,u j ,u j+1 ,u j+2}, and so on, set the p-th preset submodule to 2p-1 point sub-template S p ={u j-p+1 ,u j-p+2 ,…,u j+p-2 ,u j+p-1}; According to multiple preset sub-templates, a nonlinear interpolation format is used to obtain the left upwind value of the interior point; for the boundary point, when j = 0, the leftmost flux point The value is set to the given boundary condition value, or in the case of multiple grids, the value obtained by interpolation of the adjacent area on the left; when j = 1, 2, ..., p-1, the first preset sub-template is set to the single-point sub-template S I ={u j}, the rest of the preset sub-templates S k ,k=2,3,…,p, when x j When there are enough solution points on the left, S k ={u j-k+1 ,u j-k+2 ,…,u j+k-2 ,u j+k-1}, when the preset sub-template uses the leftmost u 1 / 2 When S k ={u 1 / 2 ,u1,…,u 2k-2}, according to the preset sub-template, a nonlinear interpolation format is used to obtain the left upwind value at the j-th flux point; when j = N-p+2, N-p+3, ..., N, the first preset sub-template is set to the single-point sub-template S I ={u j}, the rest of the preset sub-templates S k ,k=2,3,…,p, when x j When there are enough solution points on the right side, S k ={u j-k+1 ,u j-k+2 ,…,uj+k-2 ,u j+k-1}, when the sub-template is preset to use the rightmost u N When S k ={u N-2p+2 ,u N-2p+3 ,…,u N-1 ,u N}, according to the preset sub-template, a nonlinear interpolation format is used to obtain the left upwind value of the j-th flux point; according to the left upwind value of the flux point, the right upwind value of the flux point is obtained by flipping the interpolation template.

[0080] Specifically, the methods for determining the left upwind value and the right upwind value of the flux point in the three directions of the x, y, and z axes are the same, and only the preset sub-templates used are different.

[0081] In one embodiment, the preset differential format in the x-axis direction in step 110 is:

[0082]

[0083] Among them, θ j,k is the differential coefficient; f′ j is the approximate derivative of the flux at the jth solution point in the x-axis direction, Δx is the spatial step size, N is the total number of units in the calculation space, and f k-1 / 2 is the flux at the flux point with subscript k-1 / 2, f j+k-p-1 / 2 is the flux at the flux point with the subscript j+kp-1 / 2, f k+N-2p+3 / 2 is the flux at the flux point with subscript k+N-2p+3 / 2.

[0084] Specifically, the differential formats preset in the y-axis and z-axis directions are the same as the differential format preset in the x-axis direction.

[0085] In a specific embodiment, the interpolation format is 5th order precision, and the number of preset sub-templates is 3.

[0086] In the x-axis direction, the number of solution points in the x-axis direction is set to N. According to the number of solution points, the calculation interval in the x-axis direction is divided into N units using flux points at equal intervals. The solution points of the four units close to the left and right boundaries of the calculation interval are non-uniformly distributed, and the solution points of the other N-4 units are arranged at the center of the unit. The expression of the solution points in the x-axis direction is shown in formula (1).

[0087] The grid division method in the y- and z-axis directions is the same as that in the x-axis direction. Therefore, the expression of the solution point in the y- and z-axis directions is the same as that in the x-axis direction.

[0088] 2) Boundary interpolation format

[0089] Preselected submodules include: single point template S I ={u j}, three-point template S II ={u j-1 ,u j ,u j+1}, Five-point template S III ={u j-2 ,u j-1 ,u j+1 ,u j+2}, for the interior points with subscript 5≤j≤N-4, the solution value u for the flux point j+1 / 2 The interpolation template has Figure 2 Same structure as shown.

[0090] First, calculate the following values on each sub-template:

[0091] q k =φ k (x j+1 / 2 ),k=1,2,3 (9)

[0092] Among them, φ k (x) Figure 2 The Lagrange interpolation polynomial defined on the sub-template shown in k (x j+1 / 2 ) is φ k (x) in x j+1 / 2 The value at the flux point. Specifically, the above interpolation format can be explicitly written in the following general form:

[0093]

[0094]

[0095]

[0096] Then, through the nonlinear weight ω k The value q on each sub-template k By combining the following formulas, we can get the final nonlinear interpolation u j+1 / 2 :

[0097]

[0098] The nonlinear weight ω in the above formula k Calculated by the following formula:

[0099]

[0100]

[0101] Among them, as a preference, the linear weight γ k The value of Here the global smoothness factor τ is defined as:

[0102]

[0103] In the above formula, constant ε = 10 -10 , and β k Represents the corresponding sub-template S k Smooth indicator of upper interpolation. β2 and β3 are calculated as follows:

[0104]

[0105]

[0106] Among them, φ2(x) and φ3(x) are the Lagrange interpolation polynomials used to interpolate q2 and q3 respectively. The derivatives of φ2(x) and φ3(x) can be explicitly expressed as:

[0107]

[0108]

[0109]

[0110]

[0111]

[0112]

[0113] It is also worth noting that in order to obtain better performance in reducing numerical oscillations at discontinuities, the smoothness indicator β k The newly introduced positive parameter λ = 2 in the definition of acts as an amplification factor, thereby enhancing the shock wave capture capability of the format. When encountering a slowly propagating shock wave, if λ is not included, the shock wave will not be able to pass through the interface.

[0114] In particular, due to the lack of template points, β1 needs to be treated separately. To this end, the following formula is used to define it:

[0115]

[0116] in:

[0117]

[0118]

[0119] When the template SIII When it comes to non-uniform solution points, that is, when j = 3, 4 or j = N-3, N-2, the above formula still holds, but the coefficients in formulas (11), (12) and (19)-(24) are different. More precisely, for j = 3 and j = N-3, the sub-template S II The solution points involved are still uniform, so the coefficients and They are still given by formulas (11) and (19)-(20) respectively; for other cases, the corresponding coefficient values in the formulas are given in Tables 1 to 4.

[0120] Table 1 Interpolation coefficients when j=3 and j=N-2 The specific value of

[0121]

[0122]

[0123] Table 2 Interpolation coefficients when j = 3, 4 and j = N-3, N-2 The specific value of

[0124]

[0125] Table 3 Coefficients in the smooth indicator calculation formula when j = 3 and j = N-2 The specific value of

[0126]

[0127] Table 4 Coefficients in the smooth indicator calculation formula when j = 3, 4 and j = N-3, N-2 Value

[0128]

[0129] 3) Boundary interpolation format

[0130] When j = 0, 1, 2 or j = N-1, N, due to insufficient template points at the boundary of the computational domain, Figure 2 The interior point interpolation template shown can no longer be used to calculate u j+1 / 2 Therefore, a special boundary interpolation format must be constructed. For the leftmost flux point u (i.e., when j = 0) 1 / 2 To determine the value, simply set it directly to the given boundary condition value, or take the value obtained by interpolation from the adjacent area on the left in the case of multiple grids.

[0131] For the interpolation format at points j = 1, 2 or j = N-1, N, the special interpolation template used is as follows Figure 3 As shown, where (a) is the left boundary flux point x3 / 2 The sub-template at (j=1), (b) is the left boundary flux point x 5 / 2 The sub-template at (j=2), (c) is the right boundary flux point x N-1 / 2 The sub-template at (j=N-1), (d) is the right boundary flux point x N+1 / 2 Once these templates are determined, the general interpolation formulas (10)-(12) can be used. Therefore, the coefficients in formulas (11)-(12) and (19)-(24) can be determined by using only the Lagrange interpolation polynomials. and and The corresponding coefficient values are presented in the form of tables in Tables 5 to 8. Thus, the format is completed.

[0132] Table 5 Interpolation coefficients when j = 1, 2 and j = N-1, N Value

[0133]

[0134] Table 6 Interpolation coefficients when j = 1, 2 and j = N-1, N Value

[0135]

[0136] Table 7 Coefficients in the smooth indicator calculation formula when j = 1, 2 and j = N-1, N Value

[0137]

[0138] Table 8 Coefficients in the smooth indicator calculation formula when j = 1, 2 and j = N-1, N The specific value of

[0139]

[0140] 4) The pre-set differential format in the x, y, and z directions is used to obtain the approximate derivative of the flux at the j-th solution point, where the pre-set differential format in the x-axis direction is:

[0141]

[0142] Among them, θ j,k is the differential coefficient; f′ j is the approximate derivative of the flux at the jth solution point in the x-axis direction, Δx is the spatial step size, N is the total number of units in the calculation space, and f k-1 / 2 is the flux at the flux point with subscript k-1 / 2, f j+k-7 / 2is the flux at the flux point with the subscript j+k-7 / 2, f N+k-9 / 2 is the flux at the flux point with subscript N+k-9 / 2.

[0143] Difference coefficient θ j,k The coefficient values of the specific fifth-order precision interpolation format are given in Table 9.

[0144] Table 9 Coefficients θ in the fifth-order global conservative difference scheme j,k The specific value of

[0145]

[0146] In one embodiment, the aircraft flow field to be simulated in step 100 includes at least one quadrilateral computational grid.

[0147] It should be understood that although Figure 1 The steps in the flowchart are shown in sequence as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified in this document, there is no strict order restriction for the execution of these steps, and these steps can be executed in other orders. In addition, Figure 1 At least part of the steps may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed 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 part of the sub-steps or stages of other steps.

[0148] In a verification example, an accuracy test, a strong shock wave test, an internal flow problem test, and an external flow problem test were performed. The results of each test are as follows:

[0149] 1) Accuracy test

[0150] First, we verify the high-order accuracy of the format constructed by the method proposed in this work by solving a smooth linear convection problem. Taking the fifth-order format specifically described in the previous example as an example, the test results are shown in Table 10, showing that the format achieves fifth-order accuracy very well.

[0151] Table 10: Fifth-order format accuracy test results.

[0152]

[0153] 2) Strong shock wave test

[0154] Next, the format's ability to capture shock waves in a challenging model of a strong shock wave impacting a bubble was tested, as was its ability to transfer shock waves between different grid blocks on a multi-block grid. To this end, the entire computational domain was divided into four regions of varying mesh fineness, allowing the disturbances generated by the shock wave and bubble impact to interact and reflect, continuously passing through the interfaces between the grid blocks. Figure 4 The flow field image at time t = 0.4 is given. It can be seen that the shock wave and flow field details are captured very accurately. The grid block boundaries represented by the black lines in the figure do not cause any interference with the shock wave propagation. Figure 4 The black lines represent the boundaries between grid blocks.

[0155] 3) Internal flow problem test

[0156] The format's ability to capture shock waves in internal flow problems and transfer shock waves between different grid blocks on a multi-grid grid was tested by solving a duct flow problem with a protrusion. To this end, the computational domain was divided into three interconnected regions. A supersonic flow with a Mach number of 1.7 was applied to the inlet, allowing the shock waves generated within the duct to reflect off the upper and lower walls, forming a steady-state state across the different grid blocks. Figure 5 A schematic diagram of the calculation results of the supersonic pipe flow problem is given, where (a) is the result of multiple grids and (b) is the result of a single grid. Figure 5 The black lines in the middle represent the interfaces between grid blocks. Figure 5 (a) and (b) show the comparison of the flow field pressure cloud maps calculated using multi-block grids and single-block grids respectively. It can be seen that the shock wave structure near the wall boundary is captured clearly and accurately, and the results calculated on the block grid using the boundary format constructed by the present invention are completely consistent with those of the single-block grid, indicating that the boundary format has good shock wave capture and shock wave transmission capabilities, and is suitable for dealing with the multi-block grid interface problem mentioned above.

[0157] 4) Outflow problem test

[0158] The format's shock wave capture capability in external flow problems and its ability to transfer shock waves between different grid blocks on a multi-block grid were tested by solving a transonic flow problem around a NACA0012 airfoil. To this end, the entire computational domain was divided into three interconnected regions. The free flow Mach number at infinity was Ma = 0.8, and the incoming flow angle of attack was a = 1.25, resulting in two shock waves on each of the upper and lower surfaces of the wing. To increase the difficulty of the test, the grid block interface of the multi-block grid was deliberately set at the shock wave, thereby assessing the boundary format's shock wave capture capability. Figure 6 These are the calculation results of the transonic flow problem around a NACA0012 airfoil, where (a) is the result of multiple grid blocks and (b) is the result of a single grid block. The black lines in the figure represent the interfaces between grid blocks. Figure 6 Figures (a) and (b) show the comparison of the Mach number cloud maps of the flow field obtained by using multi-block grids and single-block grids respectively. It can be seen that the shock wave structure on the wing surface and nearby is clearly and accurately captured, and the results calculated on the block grid using the boundary format constructed by the present invention are completely consistent with those of the single-block grid. Figure 7 The airfoil surface pressure coefficient C is given p The distribution further confirms this conclusion, indicating that the boundary format has good shock wave capture and shock wave transmission capabilities, and is suitable for dealing with the multi-block grid interface problem mentioned above.

[0159] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, 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.

[0160] The above-described embodiments merely represent several implementation methods of the present application. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the present invention. It should be noted that a person skilled in the art could make various modifications and improvements without departing from the spirit of the present application, all of which fall within the scope of protection of the present application. Therefore, the scope of protection of the present patent application shall be determined by the appended claims.

Claims

1. A high-precision shock wave capture method based on the finite difference method, characterized in that: The method comprises: Construct the hyperbolic conservation law equations that control the flow field of the aircraft to be simulated and set the computational domain; Based on the computational domain, a semi-uniform grid distribution method is used to divide the computational domain into interlaced solution points and flux points in each coordinate direction. The semi-uniform grid distribution method is as follows: the computational interval is divided into a number of cells using flux points at equal intervals in each coordinate direction, the solution points are arranged at the center of the cells, and then the two solution points closest to the boundary are replaced with non-uniform points. Constructing a nonlinear interpolation format according to Lagrange interpolation functions on a plurality of mutually nested preset sub-templates, linear weights, and a nonlinear weight of each preset sub-template; According to the variable value at the solution point and a plurality of preset sub-templates, nonlinear interpolation is performed using the nonlinear interpolation format to obtain the left upwind value and the right upwind value at the flux point; According to the left upwind value and the right upwind value at the flux point, calculate the numerical flux at the corresponding flux point; Using a preset difference format according to the numerical flux, an approximate derivative of the flux in each coordinate direction at the solution point is obtained; According to the approximate derivatives of the flux in each coordinate direction at the solution point, the numerical simulation is carried out using the time-marching method to complete the shock wave capture of the aircraft flow field to be simulated.

2. The method according to claim 1, characterized in that According to the computational domain, a semi-uniform grid distribution method is used to divide the computational domain into interlaced solution points and flux points in each coordinate direction, including: In the x-axis direction, set the number of solution points in the x-axis direction to N, the calculation interval to [a, b], and use N+1 flux points x j+1 / 2 =a+jΔx divides the calculation interval in the x-axis direction into N units; where N is an integer greater than 1, a and b are real numbers with a<b, and Δx=(ba) / N is the spatial step size in the x-axis direction; For the interpolation format with 2p-1 order accuracy, the solution points of the 2p-2 cells near the left and right boundaries of the calculation interval are non-uniformly distributed, and the solution points of the other N-2p+2 cells are all arranged at the center of the cell; where p is the number of preset sub-templates and p is an integer greater than or equal to 1; In the y- and z-axis directions, the same grid division method as that in the x-axis direction is used to divide the calculation intervals in the y- and z-axis directions into interlaced solution points and flux points.

3. The method according to claim 1, characterized in that When the interpolation format of 2p-1 order accuracy is selected, the number of preset sub-templates is p, and the number of solution points contained in the m-th preset sub-template is 2m-1, where m = 1, 2, ..., p, and p is an integer greater than or equal to 1; According to the Lagrange interpolation function on multiple nested preset sub-templates, the linear weight and the nonlinear weight of each preset sub-template, a nonlinear interpolation format is constructed, including: The Lagrangian interpolation function is calculated according to each preset sub-template in the x-axis direction to obtain the Lagrangian interpolation format on the corresponding preset sub-template; the Lagrangian interpolation format on the preset sub-template in the x-axis direction is: q k =φ k (x j+1 / 2 ),k=1,2,…,p Among them, φ k (x) is the Lagrange interpolation function on the kth preset sub-template in the x-axis direction, φ k (x j+1 / 2 ) is the Lagrange interpolation function φ on the kth preset sub-template in the x-axis direction k (x) in x j+1 / 2 The general interpolation formula on the kth preset sub-template in the x-axis direction is: In the above formula are the Lagrange interpolation coefficients on the 2nd, 3rd, ..., p preset sub-templates in the x-axis direction, are the variable values at the solution points contained in the 1st, 2nd, 3rd, ..., pth preset sub-templates in the x-axis direction respectively; The Lagrange interpolation on the p preset sub-templates is weightedly combined using the linear weight and nonlinear weight of each preset sub-template in the x-axis direction to obtain a nonlinear interpolation format in the x-axis direction. The nonlinear interpolation format in the x-axis direction is: in, is the interpolation value at the j-th flux point in the x-axis direction, ω k is the nonlinear weight of the kth preset sub-template in the x-axis direction, γ k is the linear weight of the k-th preset sub-template in the x-axis direction, k = 1, 2, ..., p; In the y-axis and z-axis directions, nonlinear interpolation formats in the y-axis and z-axis directions are constructed based on the Lagrange interpolation functions, linear weights and nonlinear weights of multiple nested preset sub-templates in the y-axis and z-axis directions.

4. The method according to claim 3, characterized in that The Lagrangian interpolation on the p preset sub-templates is weightedly combined using the linear weight and nonlinear weight of each preset sub-template in the x-axis direction to obtain a nonlinear interpolation format in the x-axis direction. The step includes: a method for determining the nonlinear weight of the preset sub-template in the x-, y-, and z-axis directions; wherein the method for determining the nonlinear weight of the preset sub-template in the x-axis direction specifically includes: According to the Lagrangian interpolation function of the second to the p-th preset sub-template in the x-axis direction, a smooth indicator on the corresponding preset sub-template is obtained; According to the solution point in the first preset sub-template in the x-axis direction, the variable values at the two adjacent solution points, and the spatial step length in the x-axis direction, the smooth indicator on the first preset sub-template in the x-axis direction is determined; the smooth indicator on the first preset sub-template in the x-axis direction is: in, Obtaining an unnormalized nonlinear weight value corresponding to the preset sub-template in the x-axis direction according to the smoothness indicator, linear weight, and global smoothness metric factor on the preset sub-template in the x-axis direction; Normalize the unnormalized nonlinear weight values corresponding to all preset sub-templates to obtain the nonlinear weight corresponding to each preset sub-template in the x-axis direction; the expression of the nonlinear weight in the x-axis direction is: in, is the unnormalized nonlinear weight value, τ is the global smoothness factor, and ε is a constant; The preset sub-template in the x-axis direction is replaced by the preset sub-templates in the y and z-axis directions, and the nonlinear weights of the preset sub-templates in the y and z-axis directions are obtained by using the smooth indicator determination method on the preset sub-template in the x-axis direction.

5. The method according to claim 4, characterized in that According to the Lagrange interpolation function of the second to the p-th preset sub-template in the x-axis direction, the smooth indicator on the corresponding preset sub-template in the x-axis direction is obtained. In the step, the smooth indicator β on the k-th preset sub-template in the x-axis direction is k for: Where λ is a positive constant, The solution point x of the jth unit in the x-axis direction of the Lagrange interpolation function on the kth preset sub-template j The nth derivative at x j is the solution point of the jth element in the x-axis direction.

6. The method according to claim 4, characterized in that According to the smoothness indicator, linear weight and global smoothness metric factor on the preset sub-template in the x-axis direction, the unnormalized nonlinear weight value corresponding to the preset sub-template in the x-axis direction is obtained. In the step, the expression of the global smoothness metric factor in the x-axis direction is: Among them, β a is the smooth indicator on the a-th preset sub-template in the x-axis direction, a=1,2,…,p.

7. The method according to claim 1, characterized in that Interpolating each flux point according to the nonlinear interpolation format and a plurality of preset sub-templates to obtain a left upwind value and a right upwind value of the flux point includes: The flux points with computational unit number j satisfying p≤j≤N-p+1 are regarded as interior points, and the rest of the flux points are regarded as boundary points; When p≤j≤N-p+1, set the first preset sub-template to the single-point sub-template S I ={u j }, set the second preset sub-template to the three-point sub-template S II ={u j-1 ,u j ,u j+1 }, set the third preset sub-template to five-point sub-template S III ={u j-2 ,u j-1 ,u j ,u j+1 ,u j+2 }, and so on, set the p-th preset submodule to 2p-1 point sub-template S p ={u j-p+1 ,u j-p+2 ,…,u j+p-2 ,u j+p-1 };According to multiple preset sub-templates, a nonlinear interpolation format is used to obtain the left upwind value of the interior point; For the boundary point, when j = 0, the flux point on the leftmost end The value is set to the given boundary condition value, or the value obtained by interpolation from the adjacent area on the left in the case of multiple blocks of grid; When j=1,2,…,p-1, set the first preset sub-template to be the single-point sub-template S I ={u j }, the rest of the preset sub-templates S k ,k=2,3,…,p, when x j When there are enough solution points on the left, S k ={u j-k+1 ,u j-k+2 ,…,u j+k-2 ,u j+k-1 }, when the preset sub-template uses the leftmost u 1 / 2 When S k ={u 1 / 2 ,u1,…,u 2k-2 }, using the nonlinear interpolation format according to the preset sub-template to obtain the left upwind value at the j-th flux point; When j = N-p+2, N-p+3, ..., N, set the first preset sub-template to be the single-point sub-template S I ={u j }, the rest of the preset sub-templates S k ,k=2,3,…,p, when x j When there are enough solution points on the right side, S k ={u j-k+1 ,u j-k+2 ,…,u j+k-2 ,u j+k-1 }, when the preset sub-template has used the rightmost u N When S k ={u N-2p+2 ,u N-2p+3 ,…,u N-1 ,u N }, using the nonlinear interpolation format according to the preset sub-template to obtain the left upwind value at the j-th flux point; According to the left upwind value of the flux point, the right upwind value of the flux point is obtained by flipping the interpolation template.

8. The method according to claim 1, characterized in that According to the numerical flux at the flux point, a preset difference format is used to obtain the approximate derivative of the flux at the solution point. The preset difference format in the step is: Among them, θ j,k is the differential coefficient; f′ j is the approximate derivative of the flux at the jth solution point in the x-axis direction, Δx is the spatial step size in the x-axis direction, N is the total number of units in the calculation space, and f k-1 / 2 is the flux at the flux point with subscript k-1 / 2, f j+k-p-1 / 2 is the flux at the flux point with the subscript j+kp-1 / 2, f k+N-2p+3 / 2 is the flux at the flux point with subscript k+N-2p+3 / 2.

9. The method according to any one of claims 1 to 8, characterized in that Constructing a hyperbolic conservation law equation of the aircraft flow field to be simulated and setting a calculation interval, wherein the aircraft flow field to be simulated in the step includes at least one quadrilateral calculation grid.

Images