High-order nonlinear weighted difference numerical simulation method based on accurate dissipative control

Abstract

The present invention discloses a high-order nonlinear weighted difference numerical simulation method based on precise dissipative control, which relates to the field of compressible flow numerical simulation. The method comprises: step 1: constructing a computational grid based on aircraft shape parameters and computational conditions; step 2: assigning initial values ​​to nodes of the computational grid based on aircraft numerical simulation conditions; step 3: setting a time step based on aircraft numerical simulation conditions and CFL conditions; step 4: cyclically executing the following time advancement process until the termination time or convergence condition required by the aircraft numerical simulation is reached, so as to more accurately control the numerical dissipation of the WENO format during the aircraft numerical simulation, and improve the flow field simulation capability in small-scale areas without substantially reducing the shock wave capture stability.

Description

Technical Field

[0001] The present invention relates to the field of compressible flow numerical simulation, and in particular to a high-order nonlinear weighted difference numerical simulation method based on precise dissipative control. Background Art

[0002] When an aircraft is flying at high speed, shock waves, contact discontinuities, and turbulence may exist in the flow field around it. In order to obtain a more detailed flow field, a numerical method is needed to simulate both discontinuities and continuous fluctuations accurately when simulating compressible flows. The high-order precision WENO (Weighted essentially non-oscillatory) format can capture discontinuities and achieve high-order accuracy in smooth areas through nonlinear weighting technology, and is often used for detailed simulation of such problems.

[0003] For the conservative difference method solution of compressible flow, the typical numerical solution steps are as follows:

[0004] Step 1: Mesh generation. Generate mesh according to the geometric shape and calculation conditions of the simulated problem;

[0005] Step 2: Assign initial values. According to the given conditions of the simulated problem, give the initial values ​​of the nodes at the grid points;

[0006] Step 3: Set the time step. Set the time step according to the problem and the method used in combination with the CFL condition;

[0007] Step 4: Start to cycle the following time advancement process until the termination time or convergence condition required by the problem is reached:

[0008] (1) For interior points, the spatial format is used for numerical calculation;

[0009] (2) Adopt appropriate numerical boundary treatment according to boundary conditions;

[0010] (3) After the numerical calculation of the interior points and boundaries is completed in spatial discretization, it is advanced using a suitable time format.

[0011] The shock wave capture numerical simulation method mentioned here involves the spatial format in step (1) of step 4. The process of numerical simulation of the spatial format is explained by taking the one-dimensional Euler equation as an example. The governing equation of the one-dimensional Euler equation is:

[0012] ;

[0013] ;

[0014] in, U is a conserved variable, Fis the flux, and the subscripts t and x represent the partial derivatives of the corresponding variables. F ( U )express F is a function of U, T is the transposition operator, ρ,μ,p,E are density, velocity, pressure and total energy respectively, and the state equation is:

[0015] ;

[0016] Where γ is the specific heat ratio of the gas. The semi-discrete conservation difference form of the Euler equation is:

[0017] ;

[0018] in, is the numerical flux, j is the grid point subscript, h is the grid spacing. The spatial format in step (1) of step 4 is the right-hand side of the above equation. for The same method is used to obtain the value after translating one grid unit. The calculation steps are:

[0019] (a) Get the calculated mesh half-node j Variable values ​​and flux values ​​at grid points near +1 / 2;

[0020] (b) Calculation of grid points j and j +1Roe average or algebraic average of the flow field variable values, and then obtain the left characteristic matrix L and the right characteristic matrix R;

[0021] (c) Combining the left feature matrix and the grid point variable value, the flux value of the grid point is subjected to flux splitting based on the feature projection;

[0022] (d) The positive and negative parts of the flux split are calculated using the WENO format according to the components, and then the positive and negative parts are added together, and finally multiplied by the right characteristic matrix R to obtain .

[0023] The high-order WENO format has been widely used in the fine numerical simulation of compressible flows. However, in the continuous region, especially in the small-scale fluctuation region, the dissipation caused by the WENO format's linear basis format and nonlinear weighting mechanism often makes the calculation results less than ideal. When the common technology of reducing dissipation is used, the stability of shock wave capture will be affected, such as obvious non-physical oscillations near the shock wave. Summary of the invention

[0024] The purpose of the present invention is to more accurately control the numerical dissipation of the WENO format during numerical simulation of aircraft, and to improve the flow field simulation capability in small-scale areas without substantially reducing the shock wave capture stability.

[0025] To achieve the above object, the present invention provides a high-order nonlinear weighted difference numerical simulation method based on precise dissipative control, the method comprising:

[0026] Step 1: Construct the computational grid based on the aircraft shape parameters and computational conditions;

[0027] Step 2: assigning initial values ​​to the nodes of the computational grid based on the numerical simulation conditions of the aircraft;

[0028] Step 3: Set the time step based on the aircraft numerical simulation conditions and CFL conditions;

[0029] Step 4: The following time advancement process is executed cyclically. The loop ends when the termination time or convergence condition required by the numerical simulation of the aircraft is reached. The time advancement process includes:

[0030] (1) For the interior points in the computational grid, the conservative derivative approximation is used: Perform discrete calculation of the inviscid term to obtain the discrete calculation result of the inviscid term of the interior point, where: h is the grid spacing, j is the grid point subscript, Half Node j +1 / 2 corresponds to the numerical flux, Half Node j -1 / 2 corresponds to the numerical flux;

[0031] (2) Discretize the viscous term to obtain the discrete result of the viscous term;

[0032] (3) Combining the discrete calculation results of the inviscid term, the discrete results of the viscous term, the turbulence model and the numerical boundary to complete the spatial discretization and obtain the spatial discretization results;

[0033] (4) Based on the spatial discrete results, advance using the preset time format.

[0034] One or more technical solutions provided by the present invention have at least the following technical effects or advantages:

[0035] This method constructs a nonlinear weighted format by carefully designing the basis weight coefficients d k , which makes the dissipation of the basis format very low in areas with good continuity, and maintains sufficient dissipation in areas with poor continuity. Combined with the S-type smoothing factor, the final nonlinear weighted format retains the properties of the basis format to a large extent and has good shock wave capturing stability. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] The drawings described herein are used to provide a further understanding of the embodiments of the present invention, constitute a part of the present invention, and do not constitute a limitation on the embodiments of the present invention;

[0037] Figure 1 It is a flow chart of the high-order nonlinear weighted difference numerical simulation method based on precise dissipative control;

[0038] Figure 2 Schematic diagram of the real part of the approximate dispersion relation;

[0039] Figure 3 Schematic diagram of the imaginary part of the approximate dispersion relation;

[0040] Figure 4 This is a locally enlarged schematic diagram of the imaginary part of the approximate dispersion relation;

[0041] Figure 5 This is a schematic diagram of the calculation results for the modulated Gaussian wave propagation problem;

[0042] Figure 6 A schematic diagram of density comparison of the calculation results of the shock wave entropy wave interaction problem;

[0043] Figure 7 Schematic diagram of density comparison in the high-frequency vibration zone of the calculation results of the shock wave entropy wave interaction problem. DETAILED DESCRIPTION

[0044] In order to more clearly understand the above-mentioned purpose, features and advantages of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be noted that the embodiments of the present invention and the features in the embodiments can be combined with each other without conflict.

[0045] In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention may also be implemented in other ways different from those within the scope of this description. Therefore, the protection scope of the present invention is not limited to the specific embodiments disclosed below.

[0046] Embodiment 1;

[0047] Please refer to Figure 1 , Figure 1 The present invention provides a high-order nonlinear weighted difference numerical simulation method based on precise dissipative control, and the method includes:

[0048] Step 1: construct a computational grid based on the aircraft shape parameters and the computational conditions, such as: construct a corresponding computational domain geometric model based on the aircraft shape parameters and the computational conditions, define the boundary conditions, and then use the structured grid type to generate the grid required for the computation. For problems that require parallel computation, the grid needs to be partitioned. The existing computational grid construction method can be referred to, and the embodiments of the present invention will not be described in detail. The aircraft shape parameters and the computational conditions can be set or adjusted according to actual needs. The embodiments of the present invention will not be specifically limited. Specifically, the grid generation software can be used. The grid generation software can use some commonly used software on the market, such as ANSYS ICEM CFD, PointWise, etc. The embodiments of the present invention will not be specifically limited. The computational grid construction method can refer to the prior art such as CN109918808B. The embodiments of the present invention will not be correspondingly limited and can be adjusted according to actual needs.

[0049] Step 2: Based on the numerical simulation conditions of the aircraft, assign initial values ​​to the nodes of the computational grid, such as selecting the control equations and state equations used for numerical simulation based on the numerical simulation conditions of the aircraft, and considering the turbulence model for turbulence simulation, and then assign initial values ​​to the nodes of the computational grid according to the numerical simulation conditions of the aircraft, such as performing full-field assignment based on the physical state value of the flow; Step 2 is to obtain corresponding parameter values ​​using the numerical simulation conditions of the aircraft, and then assign initial values ​​to the nodes of the computational grid in the grid software, and the numerical simulation conditions of the aircraft can be adjusted accordingly according to the simulation requirements, and the embodiment of the present invention does not specifically limit this;

[0050] Step 3: Set the time step based on the numerical simulation conditions of the aircraft, the control equations used (including the turbulence model adopted), the grid used for calculation, and the CFL condition; the CFL condition is the convergence condition for judging the calculation in computational fluid dynamics;

[0051] Step 4: The following time advancement process is executed cyclically, and the loop ends when the termination time or convergence condition required by the numerical simulation of the aircraft is reached.

[0052] Among them, the present invention aims to solve the problem that it is difficult to simultaneously ensure the stability of shock wave capture and the high-fidelity simulation of small-scale fluctuations in the numerical simulation of compressible flow. Based on this, a nonlinear weighted format with more precise control of dissipation is proposed to improve the flow field simulation capability in small-scale areas without reducing the stability of shock wave capture.

[0053] The time advancement process in step 4 includes:

[0054] (1) For the interior points in the computational grid, the conservative derivative approximation is used: Perform discrete calculation of the inviscid term to obtain the discrete calculation result of the inviscid term of the interior point, where: his the grid spacing, j is the grid point subscript, Half Node j +1 / 2 corresponds to the numerical flux, Half Node j -1 / 2 corresponds to the numerical flux;

[0055] (2) Using the sixth-order central difference method or other methods to discretize the viscosity term and obtain the discrete result of the viscosity term;

[0056] (3) Combining the discrete calculation results of the inviscid term, the discrete results of the viscous term, the turbulence model and the numerical boundary to complete the spatial discretization and obtain the spatial discretization result; the turbulence model is such as SA-IDDES or PANS, etc., which is not specifically limited in the embodiment of the present invention;

[0057] (4) Based on the spatial discrete results, advance using the preset time format.

[0058] Among them, in the embodiment of the present invention, The calculation method is:

[0059] (a) Get half node j +1 / 2 the variable value of the neighboring grid point U l And the flux value F l , l = j -3,…, j +4, l are the numbers of the 8 selected nodes; the adjacent grid points can be: j +1 / 2 is the center, and four nodes are taken on each side. j +1 / 2 is located j and j +1 between;

[0060] (b) Calculation of grid points j and grid points j +1 The Roe average or algebraic average of the flow field variable values ​​is used to obtain the average value, and the left characteristic matrix is ​​obtained based on the average value. L and the right characteristic matrix R ;

[0061] (c) Based on the left feature matrix L , the variable value of the grid point U l And the flux value F l Perform flux splitting based on feature projection to obtain the forward flux splitting result and negative flux splitting results ;

[0062] (d) Taking the feature direction as the positive direction Each component of is reconstructed by WENO to obtain the first reconstruction value ;

[0063] (e) Taking the feature direction as the negative pair Each component of is reconstructed by WENO to obtain the second reconstruction value ;

[0064] (f) The first reconstruction value and the second reconstruction value Add to get the added value, multiply the added value by the left and right characteristic matrix R get ;

[0065] in, The calculation method of Roughly the same, just turn the whole j The negative direction of is translated by one grid unit. The calculation method of includes the above steps (a) to (f), and the steps (a) to (f) are j Replace with j -1, and l Replace with l = j -4,…, j +3.

[0066] Among them, in the embodiment of the present invention, is the numerical flux of the nonlinear weighted format. The nonlinear weighted format in the present invention corresponds to the WENO reconstruction method in step (d). To make the formula more concise, f l and replace and . , where the feature direction is Each component of is reconstructed by WENO to obtain the first reconstruction value ,use f l express For any component, the reconstructed value is The corresponding component of ,but The calculation method is:

[0067] ; (1)

[0068] in, ω 0 to ω 4 are the first to fifth nonlinear weight coefficients, to They are the numerical fluxes of the first to fifth sub-templates respectively.

[0069] In the embodiment of the present invention, the calculation methods of the numerical fluxes of the first to fifth sub-templates are respectively:

[0070] The calculation methods of the numerical flux of the first to fifth sub-templates are:

[0071] ; (2)

[0072] ; (3)

[0073] ; (4)

[0074] ; (5)

[0075] ; (6)

[0076] in, f j-3 to f j+4 for The corresponding amount.

[0077] In this embodiment of the present invention, the first to fifth nonlinear weight coefficients are calculated as follows:

[0078] ; (7)

[0079] in, m To sum the variables, is the unnormalized weight coefficient, d k is the basis weight coefficient, γ k is the intermediate variable, τ is the intermediate quantity, β k is the smoothing factor, ε To prevent the denominator from being a positive number of 0. Among them, ε is 10 -40 . Smoothness factor β k A four-point S-shaped smoothing factor is used.

[0080] Among them, in the embodiment of the present invention, β k The calculation method is:

[0081] ; (8)

[0082] ;

[0083] ;

[0084] ;

[0085] in, f j+k-3 to f j+k for The corresponding amount, n , y and z All are intermediate amounts.

[0086] In the embodiment of the present invention, in order to ensure the stability of the format, β 4 The following correction is made, and the corresponding nonlinear weight coefficient is calculated using the corrected smoothing factor. The correction method is:

[0087] ; (9)

[0088] in, β 4r is the corrected smoothing factor.

[0089] In the embodiment of the present invention, the parameter τ adopts an eight-point S-shaped smoothing factor, and the calculation method of τ is:

[0090] ; (10)

[0091] ;

[0092] ;

[0093] ;

[0094] in, f j-3 to f j+4 for The corresponding amount, v , w and e All are intermediate amounts.

[0095] In this embodiment of the present invention, the base weight coefficient d k is the coefficient of the seventh-order upwind format and the coefficients of the eighth-order central format A nonlinear combination of d k The calculation method is:

[0096] ; (11)

[0097] Among them, λ is the weight, are the coefficients of the seventh-order upwind format, are the coefficients of the eighth-order central format.

[0098] In this embodiment of the present invention, the specific calculation method is:

[0099] ; (12)

[0100] ; (13)

[0101] in, All are coefficients;

[0102] λ is calculated as:

[0103] ; (14)

[0104] in, γ m is an intermediate variable.

[0105] The calculation formula in step (e) is symmetric with the formula in step (d) about j+1 / 2, with the characteristic direction as the negative pair. Each component of is reconstructed by WENO to obtain the second reconstruction value , still use f l express For any component, the reconstructed value is The corresponding component of ,but The calculation method of is mirror symmetric with formula (1)-formula (14).

[0106] The present invention constructs a nonlinear weighted format by carefully designing the basis weight coefficients d k , which makes the base format dissipate very little in areas with good continuity, and maintains sufficient dissipation in areas with poor continuity. Combined with the S-type smoothing factor, the final nonlinear weighted format retains the properties of the base format to a large extent and has good shock wave capture stability. The comparative analysis with several comparable formats is as follows:

[0107] This method has excellent resolution properties. The relationship between the corrected wave number and the dimensionless wave number can be obtained through approximate dispersion relation analysis. Figure 2-Figure 4 As shown, Figure 2 is a schematic diagram of the real part of the approximate dispersion relation, Figure 3 is a schematic diagram of the imaginary part of the approximate dispersion relation, Figure 4 This is a locally enlarged schematic diagram of the imaginary part of the approximate dispersion relation. Figure 2The horizontal axis is the dimensionless wave number, and the vertical axis is the real part of the corrected wave number. Figure 3 The horizontal axis is the dimensionless wave number, and the vertical axis is the imaginary part of the corrected wave number. Figure 4 The horizontal axis is the dimensionless wave number, and the vertical axis is the imaginary part of the corrected wave number. Figure 2-Figure 4 WENO8-SU is used to represent the format of the present invention. Its real part curve is consistent with the WENO7-S format, which is better than the WENO7-JS format and the WENO7-Z format, indicating its good dispersion properties. The imaginary part curve is similar to WENO7-S, and its absolute value is smaller than the other two comparison formats. In addition, the imaginary part of the ADR value of this format at the dimensionless wave number π is 0, indicating that there is no dissipation at this time. Only when the waveform completely conforms to this situation will there be no dissipation, which almost never happens in actual calculations. Even if it occurs occasionally, as the calculation proceeds, the slight change in the waveform causes its dissipation to increase sharply, so this situation does not cause adverse effects. In the enlarged figure, you can see that , its absolute value is significantly smaller than other formats, while , and its curve coincides with that of WENO7-S. This result shows that the dissipation of this format is smaller in the low wavenumber region, while the dissipation in the high wavenumber region is similar to that of WENO7-S.

[0108] Good shock wave capture stability. According to the construction method of this format, the weight λ It takes 1 or very close to 1 near the discontinuity, that is, its basis format is basically consistent with WENO7-S. In addition, its weighting method is similar to WENO7-S, but due to the use of a larger template, when the template contains a discontinuity, τ is usually slightly larger than the corresponding value in WENO7-S. This makes the dissipation of this format near the discontinuity often slightly larger than that of the WENO7-S format, ensuring its good shock wave capture stability.

[0109] The following example introduces this method:

[0110] Example 1: Propagation problem of modulated Gaussian wave;

[0111] The governing equation for this problem is the linear convection equation with a velocity of 1. The computational domain is [0,1], a periodic boundary is used, and the initial value condition is:

[0112] ;

[0113] in, u 0( x ) is the waveform function, x is the spatial coordinate;

[0114] The number of grid points N is 120, and the number of CFL is 0.5. The calculation end time is T=1. The calculation results of several WENO formats are shown in Figure 5 , Figure 5The solid line Exact in the middle represents the exact solution. Obviously, the WENO8-SU format proposed in the present invention has the smallest dissipation and obtains the best result.

[0115] Example 2: The problem of shock wave-entropy wave interaction;

[0116] The governing equation for this problem is the one-dimensional Euler equation, where a shock wave with a Mach number of 3 interferes with the entropy wave, inducing high-frequency vibration. Here the number of grid points is N = 200, and the CFL number is 0.5. Figure 6 A comparison chart of the calculation results of several WENO formats when T=1.8 is given. Since there is no analytical solution to this problem, the calculation results of the classic fifth-order WENO format when N=8000 are used as the reference exact solution, which is represented by the solid line Exact in the figure. Figure 6 You can see that several formats capture shock waves very well. Figure 7 This is a local enlarged view of the high-frequency vibration area. It can be seen from the figure that the WENO8-SU format proposed in the present invention achieves the best result.

[0117] Although the preferred embodiments of the present invention have been described, those skilled in the art may make other changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present invention.

[0118] Obviously, those skilled in the art can make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include these modifications and variations.

Claims

1. A high-order nonlinear weighted difference numerical simulation method based on precise dissipative control, characterized in that: The method comprises: Step 1: Construct the computational grid based on the aircraft shape parameters and computational conditions; Step 2: assigning initial values ​​to the nodes of the computational grid based on the numerical simulation conditions of the aircraft; Step 3: Set the time step based on the aircraft numerical simulation conditions and CFL conditions; Step 4: The following time advancement process is executed cyclically. The loop ends when the termination time or convergence condition required by the numerical simulation of the aircraft is reached. The time advancement process includes: (1) For the interior points in the computational grid, the conservative derivative approximation is used: Perform discrete calculation of the inviscid term to obtain the discrete calculation result of the inviscid term of the interior point, where: h is the grid spacing, j is the grid point subscript, Half Node j +1 / 2 corresponds to the numerical flux, Half Node j -1 / 2 corresponds to the numerical flux; (2) Discretize the viscous term to obtain the discrete result of the viscous term; (3) Combining the discrete calculation results of the inviscid term, the discrete results of the viscous term, the turbulence model and the numerical boundary to complete the spatial discretization and obtain the spatial discretization results; (4) Based on the spatial discretization results, advance using a preset time format; The calculation method is: (a) Get half node j +1 / 2 the variable value of the neighboring grid point U l And the flux value F l , l = j -3,…, j +4, l is the number of the 8 selected nodes; (b) Calculation of grid points j and grid points j +1 The Roe average or algebraic average of the flow field variable values ​​is used to obtain the average value, and the left characteristic matrix is ​​obtained based on the average value. L and the right characteristic matrix R ; (c) Based on the left feature matrix L , the variable value of the grid point U l And the flux value F l Perform flux splitting based on feature projection to obtain the forward flux splitting result and negative flux splitting results ; (d) Taking the feature direction as the positive direction Each component of is reconstructed by WENO to obtain the first reconstruction value ; (e) Taking the feature direction as the negative pair Each component of is reconstructed by WENO to obtain the second reconstruction value ; (f) The first reconstruction value and the second reconstruction value Add to get the added value, multiply the added value by the left and right characteristic matrix R get ; The calculation method of includes the above steps (a) to (f), and the steps (a) to (f) are j Replace with j -1, and l Replace with l = j -4,…, j +3; Taking the feature direction as positive Each component of is reconstructed by WENO to obtain the first reconstruction value ,use f l express For any component, the reconstructed value is The corresponding component of , The calculation method is: ； in, ω 0 to ω 4 are the first to fifth nonlinear weight coefficients, to They are the numerical fluxes of the first to fifth sub-templates respectively; The calculation method of the first to fifth nonlinear weight coefficients is: ； in, m To sum the variables, is the unnormalized weight coefficient, d k is the basis weight coefficient, γ k is the intermediate variable, τ is the intermediate quantity, β k is the smoothing factor, ε is a positive number; β k The calculation method is: ； ； ； ； in, f j+k-3 to f j+k for The corresponding amount, n , y and z All are intermediate amounts.

2. The high-order nonlinear weighted difference numerical simulation method based on precise dissipative control according to claim 1, characterized in that: The calculation methods of the numerical flux of the first to fifth sub-templates are: ； ； ； ； ； in, f j-3 to f j+4 for The corresponding amount.

3. The high-order nonlinear weighted difference numerical simulation method based on precise dissipative control according to claim 1, characterized in that: right β 4 The following correction is made, and the corresponding nonlinear weight coefficient is calculated using the corrected smoothing factor. The correction method is: ； in, β 4r is the corrected smoothing factor.

4. The high-order nonlinear weighted difference numerical simulation method based on precise dissipative control according to claim 1, characterized in that: τ is calculated as: ； ； ； ； in, f j-3 to f j+4 for The corresponding amount, v , w and e All are intermediate amounts.

5. The high-order nonlinear weighted difference numerical simulation method based on precise dissipative control according to claim 1, characterized in that: d k The calculation method is: ； Among them, λ is the weight, are the coefficients of the seventh-order upwind format, are the coefficients of the eighth-order central format.

6. The high-order nonlinear weighted difference numerical simulation method based on precise dissipative control according to claim 5, characterized in that: ； ； in, All are coefficients; λ is calculated as: ； in, γ m is an intermediate variable.

Images