A new WENO format construction method for numerical simulation based on trigonometric function framework

Abstract

The present invention relates to a method for constructing a new WENO format under a trigonometric function framework for numerical simulation. The weighted substantially non-oscillating format constructed with trigonometric function polynomials is easier to simulate wave-type or high-frequency oscillation-type problems than the classical substantially non-oscillating format constructed with algebraic polynomials. High-order numerical accuracy can be obtained in smooth areas, and substantially non-oscillating properties can be maintained at shock waves and contact discontinuities. Although the new TWENO format uses the same five-point information as the classical fifth-order WENO format, it can obtain a lower global L 1 and L ∞ Norm truncation error. The linear weights used in the new TWENO format no longer need to be calculated tediously to obtain the optimal solution, and can be set to any positive number that satisfies the sum of 1. Compared with the classic WENO format, it is simpler, more robust, and easier to generalize to high-dimensional space. The new TWENO format effectively simulates several classic Euler problems and fully verifies its effectiveness.

Description

[0001] This application is a divisional application of a patent application entitled "A method for constructing a new WENO format under a trigonometric function framework". The application date of the original application is May 16, 2018, and the application number is 201810472192.4. Technical Field

[0002] The invention belongs to the technical field of computational fluid dynamics engineering, and in particular relates to a new WENO format construction method under a trigonometric function framework for numerical simulation. Background Art

[0003] In engineering applications, flow field problems often occur, such as gas dynamic systems and shallow water modeling. Therefore, it is crucial to develop robust, accurate, and efficient numerical simulation methods to solve such problems, which has also attracted the interest of many researchers. In 1959, Godunov proposed a first-order accuracy numerical simulation format for solving flow field problems. The first-order accuracy numerical simulation method will not produce non-physical numerical oscillations when capturing shock waves, but will over-smooth strong discontinuities, which are often of great significance to the subsequent research of the problem. Therefore, it is necessary to introduce a high-precision numerical calculation format to simulate strong discontinuity problems.

[0004] In order to improve the accuracy of the format, simulate the structure of the solution and accurately capture the shock wave position, Harten first proposed the TVD (Total Variation Diminishing) format in 1983, and on this basis, he and Osher proposed the ENO (Essentially Non-Oscillatory) high-precision format in 1987. The main idea of ​​the ENO format is to select the smoothest template in the successively expanded templates to construct a polynomial to obtain the value at the unit boundary, thereby achieving high-order accuracy and high resolution in the smooth area, and achieving a basically non-oscillation effect near the discontinuity. However, in the implementation of the method, the ENO format ultimately only selects the best template among all candidate templates, resulting in a waste of calculation results, and the higher the accuracy of the construction, the more waste there is, resulting in low calculation efficiency. Therefore, Liu, Osher and Chan et al. proposed the WENO (Weighted Essentially Non-oscillatory) format in 1994, which improved the utilization rate of the calculation results and improved the r-order accuracy of the ENO format to r+1-order accuracy. In 1996, Jiang and Shu further improved the WENO format, so that the numerical accuracy can be improved to 2r-1 order, and designed a new smoothing factor and nonlinear weight construction framework. The main idea of ​​the WENO format is to obtain a high-order approximation through a linear convex combination of low-order reconstruction fluxes. However, in the implementation of the classic WENO format, the linear weights depend on the mother template, and the solution process is quite complicated. Therefore, in 2016, Zhu and Qiu improved the WENO format by randomly selecting linear weights greater than zero and summing to one while maintaining the accuracy. These formats have been successfully used in many application fields, especially problems involving shock waves and complex solution structures, such as simulating compressible turbulent systems and aeroacoustic systems.

[0005] Wave and high-frequency oscillation problems often occur in engineering applications. However, there are few studies on the more suitable trigonometric polynomial interpolation WENO format for simulating such problems. Although Baron presented the Neville-type method in 1976 and Muhlbach proposed Newton's trigonometric interpolation, these results cannot be directly applied to the ENO-type interpolation format. For this reason, Christofi proposed a trigonometric reconstruction method that can be directly used in the ENO format in 1996. Zhu and Qiu proposed a method to reconstruct the WENO format using trigonometric polynomials in 2010, but the calculation is complex and difficult to implement. Summary of the invention

[0006] In view of the deficiencies in the prior art, the present invention provides a new WENO format construction method under a trigonometric function framework for numerical simulation, which can perform high-precision numerical simulation on various compressible flow field problems.

[0007] To achieve the above object, the present invention provides the following solutions:

[0008] A new WENO format construction method for numerical simulation under the framework of trigonometric functions is provided. In a Cartesian coordinate system, the TWENO format is used to perform numerical simulation on a compressible flow field problem, and the method comprises the following steps:

[0009] Step 1: Discretize the hyperbolic conservation law equation into a spatial semi-discrete finite difference format and reconstruct the approximate value of the flux using the TWENO format;

[0010] Step 2: Use the third-order TVD Runge-Kutta discretization formula to discretize the spatial semi-discrete finite difference format into a space-time fully discrete finite difference format for the time derivative in the control equation;

[0011] Step 3: Obtain the approximate value on the next time layer according to the fully discrete finite difference format of time and space, iterate in sequence, and obtain the numerical simulation value of the flow field in the calculation area at the end time;

[0012] Among them, the approximate values ​​of the flux reconstructed using the TWENO format include:

[0013] Discretize space into grid cells of uniform length Unit length Unit Center Where i is the coordinate number, there are:

[0014]

[0015] Where L(U) represents -f(U) x The spatial discretization form of , U = (ρ, ρu, E) T represents the conserved variable, f(U)=(ρu,ρu 2 +p,u(E+p)) T represents the flux, f(U) x represents the derivative of f(U) with respect to x, x represents the spatial variable, ρ, u, p, E represent the fluid density, velocity, pressure, and energy respectively, T represents the transposition, and U i (t) represents the position of U in the target grid unit I i Interior point x i The value U(x i ,t), and They represent the flux f(U) in the target grid cell I i The border and The numerical flux of the fifth-order approximation at ;

[0016] Find the flux f(U) in the target grid cell Ii The border and The numerical flux of the fifth-order approximation at and The specific steps are as follows:

[0017] Step 1: Use Lax-Friedrichs splitting to split the flux into in,

[0018] Step 2: Set the target grid cell I i And the five grid cells around it form a large template T1 = [I i-2 ,I i-1 ,I i ,I i+1 ,I i+2 ], select two small templates containing two grid cells from the large template T2 = [I i-1 ,I i ] and T3=[I i ,I i+1 ];

[0019] Step 3: Reconstruct the trigonometric polynomials p1(x), p2(x) and p3(x) on each template of T1, T2 and T3 respectively, so that:

[0020]

[0021] p2(x),p3(x)∈span{1,sin(x i )};

[0022] Step 4: Randomly select three groups of linear weights, the sum of each group of linear weights is 1;

[0023] Step 5: Calculate the smoothness indicator based on the trigonometric polynomial to measure the trigonometric polynomial p. l (x) Smoothness on the target grid cell; l = 1, 2, 3 represents the corresponding template number;

[0024] Step 6, calculating the nonlinear weight through the linear weight and the smooth indicator;

[0025] Step 7: Calculate the numerical flux splitting f based on the nonlinear weight and trigonometric polynomial + (U) at point TWENO reconstruction value at ;

[0026] Similarly, find the numerical flux splitting f - (U) at point TWENO reconstruction value at, numerical flux splitting f + (U) at point TWENO reconstruction value at, numerical flux splitting f - (U) at point TWENO reconstruction value at ;

[0027] Substituting the calculated results into the spatial semi-discrete finite difference format containing the time derivative term, the ordinary differential equation about the time derivative is obtained.

[0028] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0029] Compared with the WENO format, the TWENO format simulates wave-like and high-frequency oscillation-like compressible flow field problems by using trigonometric polynomials instead of algebraic polynomials as the building blocks of the finite difference TWENO format, and can achieve high-order accuracy in smooth areas; compared with the existing trigonometric polynomial reconstruction format, the global L obtained by the TWENO format is 1 Truncation error and L ∞ The truncation error is smaller, and non-physical oscillations at strong shock waves and contact discontinuities are avoided. The relevant linear weights in the new fifth-order TWENO format no longer need to be obtained through complex calculations but are set to arbitrary positive numbers that sum to one. Therefore, the new TWENO format has the advantage of being simpler and easier to expand to high-dimensional space. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0031] Figure 1a-1c This is the step problem in the first embodiment. The density contour map obtained by the finite difference TWENO format of the present invention adopts linear weights ①, ②, and ③ respectively.

[0032] Figure 2a-2c This is the double Mach problem in the second embodiment. The density contour map obtained by the finite difference TWENO format of the present invention adopts linear weights ①, ②, and ③ respectively.

[0033] Figure 3a-3c This is the mutual interference problem between shock waves and eddies in Example 3. The pressure contour map at t=0.35 obtained by the finite difference TWENO format of the present invention uses linear weights ①, ②, and ③ respectively.

[0034] Figure 4a-4cThis is the mutual interference problem between shock waves and eddies in Example 3. The pressure contour map at t=0.6 obtained by the finite difference TWENO format of the present invention uses linear weights ①, ②, and ③ respectively.

[0035] Figure 5a-5c This is the mutual interference problem between shock waves and eddies in Example 3. The pressure contour map at t=0.8 obtained by the finite difference TWENO format of the present invention uses linear weights ①, ②, and ③ respectively.

[0036] Figure 6a-6c It is the two-dimensional Euler Riemann problem with the initial value condition (15) in the fourth embodiment. The density contour map at t=0.25 is obtained by using the finite difference TWENO format of the present invention, and linear weights ①, ②, and ③ are used respectively.

[0037] Figure 7a-7c It is the two-dimensional Euler Riemann problem with the initial value condition (16) in the fourth embodiment. The density contour map at t=0.25 is obtained by using the finite difference TWENO format of the present invention, and linear weights ①, ②, and ③ are used respectively.

[0038] Figure 8a-8c It is the two-dimensional Euler Riemann problem with the initial value condition (17) in the fourth embodiment. The density contour map at t=0.3 is obtained by using the finite difference TWENO format of the present invention, and linear weights ①, ②, and ③ are used respectively.

[0039] Figure 9a-9c It is the two-dimensional Euler Riemann problem with the initial value condition (18) in the fourth embodiment. The density contour map at t=0.2 is obtained by using the finite difference TWENO format of the present invention, and linear weights ①, ②, and ③ are used respectively.

[0040] Figure 10a-10c It is the two-dimensional Euler Riemann problem with the initial value condition (19) in the fourth embodiment. The density contour map at t=0.3 is obtained by using the finite difference TWENO format of the present invention, and linear weights ①, ②, and ③ are used respectively. DETAILED DESCRIPTION

[0041] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0042] The purpose of the present invention is to provide a new WENO format construction method under the framework of trigonometric functions for numerical simulation, which can perform high-precision numerical simulation for various compressible flow field problems.

[0043] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0044] The present invention provides a construction process of a new fifth-order finite difference TWENO high-precision numerical calculation format for solving compressible flow field problems under Cartesian grids. Compared with the classic WENO format, the TWENO format solves the numerical simulation of compressible flow field problems of wave type and high-frequency oscillation type by using reconstructed trigonometric function polynomials instead of reconstructed algebraic polynomials as the building blocks of the finite difference WENO format, and can achieve high-order precision approximation in smooth areas, capturing the transition of sharp and non-oscillating shock waves. The linear weights used in the new TWENO format no longer need to be obtained through tedious numerical calculations, and can be set to any positive number that satisfies the sum of unity. This format method is simple, highly accurate, and can be easily extended to multi-dimensional space. The method uses the TWENO format to perform numerical simulation of compressible flow field problems in a Cartesian coordinate system, and the specific steps are as follows:

[0045] Step 1: Discretize the hyperbolic conservation law equations into a spatial semi-discrete finite difference format and reconstruct the approximate value of the flux using the TWENO format.

[0046] Consider the one-dimensional hyperbolic conservation law equation:

[0047]

[0048] Its semi-discrete format is:

[0049]

[0050] Where U = (ρ, ρu, E) T represents the conserved variable, f(U)=(ρu,ρu 2 +p,u(E+p)) T represents flux, U t represents the derivative of U with respect to t, f(U) x represents the derivative of f(U) with respect to x, t represents the time variable, x represents the space variable, ρ, u, p, E represent the fluid density, velocity, pressure, and energy respectively, T represents the transposition, U0 represents the initial state value, and L(U) represents -f(U) x The spatial discretization form of .

[0051] Discretize space into grid cells of uniform length Unit length The unit center is Where i is the coordinate number, there are:

[0052]

[0053] in, and They represent the flux f(U) in the target grid cell I i The border and The numerical flux of the fifth-order approximation at i (t) represents the position of U in the target grid unit I i Interior point x i The value U(x i ,t).

[0054] Find the flux f(U) in the target grid cell I i The border and The fifth-order approximation at and The specific steps are as follows:

[0055] Step 1: Use the simplest Lax-Friedrichs splitting to split the flux into in For the sake of simplicity, the present invention only describes + (U) at point The reconstruction process at and defines it as

[0056] Step 2: Set the target grid cell I i And the 5 grid cells around it form a large template T1 = [I i-2 ,I i-1 ,I i ,I i+1 ,I i+2 ], select two small templates containing two grid cells from the large template T2 = [I i-1 ,I i ] and T3=[I i ,I i+1 ], where I i is the grid unit with the corresponding serial number.

[0057] Step 3: Reconstruct the trigonometric polynomials p1(x), p2(x) and p3(x) on each template respectively, so that:

[0058]

[0059] p2(x),p3(x)∈span{1,sin(x i )}.

[0060] The specific process is as follows:

[0061] Step 3.1, construct trigonometric polynomials p1(x), p2(x) and p3(x) on three templates T1, T2 and T3 respectively, so that they satisfy:

[0062]

[0063]

[0064]

[0065] Step 3.2, obtain the trigonometric interpolation polynomials p1(x), p2(x) and p3(x) on each template as follows:

[0066]

[0067]

[0068]

[0069] Among them, I i-2 ,I i-1 ,I i ,I i+1 ,I i+2 Respectively represent the i-2, i-1, i, i+1, and i+2 grid cells, Respectively represent f + (U) at point x i-2 、x i-1 、x i 、x i+1 、x i+2 The value of , h is the grid step size.

[0070] Step 4: Take three sets of linear weights at random:

[0071] ①γ1=0.98, γ2=0.01, γ3=0.01;

[0072] ②γ1=1 / 3, γ2=1 / 3, γ3=1 / 3;

[0073] ③γ1=0.01, γ2=0.495, γ3=0.495.

[0074] Step 5: Calculate the smoothness indicator β l , used to measure the reconstruction polynomial p l (x) The smoothness on the target grid cell is calculated as:

[0075]

[0076] Among them, l = 1, 2, 3 represents the corresponding template number, Represents the polynomial p l (x) The α-order derivative of x, r1=4, r2=1, r3=1.

[0077] Step 6: Use linear weight γ l and smooth indicator β l Calculate the nonlinear weight ω l , and its calculation formula is:

[0078]

[0079] Among them, l = 1, 2, 3 represents the corresponding template number, τ is the transition value in the calculation process, β l is a smooth indicator, ε=10 -6 Prevents the denominator from being zero.

[0080] Step 7: Find the numerical flux splitting f + (U) at point TWENO reconstruction value at:

[0081]

[0082] Secondly, the calculation results are substituted into the semi-discrete finite difference format containing the time derivative term to obtain the ordinary differential equation about the time derivative.

[0083] Step 2: Use the third-order TVD Runge-Kutta discretization formula to discretize the semi-discrete finite difference format into a time-space fully discrete finite difference format for the time derivative in the control equation.

[0084] Using the third-order TVD Runge-Kutta discretization formula:

[0085]

[0086] The space-time fully discrete finite difference format is obtained, where U (1) ,U (2) is the intermediate transition value, Δt is the time step, the superscript n represents the nth time layer, L(U n ),L(U (1) ),L(U (2) ) is -f(U) x An approximation of the higher-order spatial discretization form of .

[0087] Step 3: Obtain the approximate value on the next time layer according to the fully discrete finite difference format in time and space, iterate in sequence, and obtain the numerical simulation value of the flow field in the calculation area at the termination time.

[0088] The space-time fully discrete finite difference format is an iterative formula for the time layer. The initial state value is known. The approximate value of the next time layer is obtained through the iterative formula, and the numerical simulation value of the flow field in the calculation area at the termination time is obtained in turn. For two-dimensional problems, the above reconstruction process is used dimension by dimension.

[0089] Several examples are given below as specific embodiments of the method disclosed in the present invention.

[0090] Embodiment 1, step problem. This problem is a classic example proposed by Emery in 1968 for verifying the nonlinear hyperbolic conservation law format. The initial data is that the horizontal inflow Mach number is 3, the density is 1.4, the horizontal velocity is 3, the vertical velocity is 0, the pressure is 1, the pipe area is [0,3]×[0,1], and there is a step with a height of 0.2 at a distance of 0.6 from the left boundary, and the step extends to the end of the pipe. The upper and lower boundaries are reflection boundaries, the left boundary is the inflow boundary, and the right boundary is the outflow boundary. Figure 1a-1c The density contour map at t=4 is given.

[0091] Example 2: Double Mach reflection problem. This problem describes the changes that occur when a strong shock wave with an angle of 60° to the x-axis hits a reflection wall. The incoming flow is a strong shock wave with a Mach number of 10. The calculation area is [0,4]×[0,1], and the bottom of the area is The reflection boundary condition starts at y=0, and the other bottom boundaries (from x=0 to That part) is the wavefront condition. Figure 2a-2c The density contour map in the [0,3]×[0,1] region at t=0.2 is given.

[0092] Example 3: Mutual interference between shock wave and vortex. A shock wave with a Mach number of 1.1 is located at x=0.5 and is perpendicular to the x-axis. The initial state of the shock wave is The small vortex is located to the left of the shock wave and its center is located at (x c ,y c )=(0.25,0.5), the eddy current can be regarded as the disturbance of velocity, temperature and entropy of the mean flow, expressed as:

[0093]

[0094] in, ε=0.3,r c =0.05, α=0.204, γ=1.4, and the calculation area is [0,2]×[0,1]. Figure 3a-3c The pressure contour map in the [0,1]×[0,1] region at t=0.35 is given. Figure 4a-4c The pressure contour map in the [0.4,1.45]×[0,1] region at t=0.6 is given. Figure 5a-5cThe pressure contour map in the [0,2]×[0,1] region at t=0.8 is given.

[0095] Example 4: Two-dimensional Euler Riemann problem. The computational domain is [0,1]×[0,1], and the initial value conditions are set as:

[0096]

[0097]

[0098]

[0099]

[0100]

[0101] Figure 6a-6c The density contour diagram of the two-dimensional Euler Riemann problem at t=0.25 is given when the initial value condition is (15). Figure 7a-7c The density contour diagram of the two-dimensional Euler Riemann problem at t=0.25 is given when the initial value condition is (16). Figure 8a-8c The density contour diagram of the two-dimensional Euler Riemann problem at time t=0.3 is given when the initial value condition is (17). Figure 9a-9c The density contour diagram of the two-dimensional Euler Riemann problem at time t=0.2 is given when the initial value condition is (18). Figure 10a-10c The density contour map of the two-dimensional Euler Riemann problem at t = 0.3 when the initial value condition is (19) is given. It can be seen from the figure that the finite difference TWENO format based on the trigonometric function polynomial space can well capture most of the flow characteristics of the Riemann problem for any linear weight obtained by the present invention.

[0102] Each embodiment in this specification focuses on the differences from other embodiments, and the same or similar parts between the embodiments can be referenced to each other.

[0103] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A new WENO format construction method for numerical simulation under the framework of trigonometric functions, using the TWENO format to perform numerical simulation on compressible flow field problems in a Cartesian coordinate system, characterized in that: The following steps are involved: Step 1: Discretize the hyperbolic conservation law equation into a spatial semi-discrete finite difference format and reconstruct the approximate value of the flux using the TWENO format; Step 2: Use the third-order TVD Runge-Kutta discretization formula to discretize the spatial semi-discrete finite difference format into a space-time fully discrete finite difference format for the time derivative in the control equation; Step 3: Obtain the approximate value on the next time layer according to the fully discrete finite difference format of time and space, iterate in sequence, and obtain the numerical simulation value of the flow field in the calculation area at the end time; Among them, the approximate values ​​of the flux reconstructed using the TWENO format include: Discretize space into grid cells of uniform length Unit length Unit Center Where i is the coordinate number, there are: Where L(U) represents -f(U) x The spatial discretization form of , U=(ρ,ρu,E) T represents the conserved variable, f(U)=(ρu,ρu 2 +p,u(E+p)) T represents the flux, f(U) x represents the derivative of f(U) with respect to x, x represents the spatial variable, ρ, u, p, E represent the fluid density, velocity, pressure, and energy respectively, T represents the transposition, and U i (t) represents the position of U in the target grid unit I i Interior point x i The value U(x i ,t), and They represent the flux f(U) in the target grid cell I i The border and The numerical flux of the fifth-order approximation at ; Find the flux f(U) in the target grid cell I i The border and The numerical flux of the fifth-order approximation at and The specific steps are as follows: Step 1: Use Lax-Friedrichs splitting to split the flux into in, Step 2: Set the target grid cell I i And the five grid cells around it form a large template T1 = [I i-2 ,I i-1 ,I i ,I i+1 ,I i+2 ], select two small templates containing two grid cells from the large template T2 = [I i-1 ,I i ] and T3=[I i ,I i+1 ]; Step 3: Reconstruct the trigonometric polynomials p1(x), p2(x) and p3(x) on each template of T1, T2 and T3 respectively, so that: p2(x),p3(x)∈span{1,sin(x-x i )}; Step 4: Randomly select three groups of linear weights, the sum of each group of linear weights is 1; Step 5: Calculate the smoothness indicator based on the trigonometric polynomial to measure the trigonometric polynomial p. l (x) Smoothness on the target grid cell; l = 1, 2, 3 represents the corresponding template number; Step 6, calculating the nonlinear weight through the linear weight and the smooth indicator; Step 7: Calculate the numerical flux splitting f based on the nonlinear weight and trigonometric polynomial + (U) at point TWENO reconstruction value at ; Similarly, find the numerical flux splitting f - (U) at point TWENO reconstruction value at, numerical flux splitting f + (U) at point TWENO reconstruction value at, numerical flux splitting f - (U) at point TWENO reconstruction value at ; Substituting the calculated results into the spatial semi-discrete finite difference format containing the time derivative term, the ordinary differential equation about the time derivative is obtained.

2. The method according to claim 1, characterized in that: In step 1, the hyperbolic conservation law equation is: Among them, U t It represents the derivative of U with respect to t, t represents the time variable, and U0 represents the initial state value.

3. The method according to claim 2, characterized in that: In the step 1, the spatial semi-discrete finite difference format is:

4. The method according to claim 1, characterized in that: The specific steps of step 3 are as follows: Step 3.1: Construct trigonometric polynomials p1(x) and p2(x) on three templates T1, T2 and T3 respectively. and p3(x) such that: Step 3.2, obtain the trigonometric interpolation polynomials p1(x), p2(x) and p3(x) on each template as follows: Among them, I i-2 ,I i-1 ,I i ,I i+1 ,I i+2 Respectively represent the i-2, i-1, i, i+1, and i+2 grid cells, f i + , Respectively represent f + (U) at point x i-2 、x i-1 、x i 、x i+1 、x i+2 The value of .

5. The method according to claim 1, characterized in that: In step 4, the three sets of linear weights are: γ1=0.98, γ2=0.01, γ3=0.01; γ1=1 / 3, γ2=1 / 3, γ3=1 / 3; γ1=0.01, γ2=0.495, γ3=0.

495.

6. The method according to claim 1, characterized in that: In step 5, the calculation formula for calculating the smoothness indicator is: Among them, β l For smooth indicators, Represents the trigonometric polynomial p l (x) The α-order derivative of x, r1=4, r2=1, r3=1.

7. The method according to claim 5, characterized in that: In step 6, the calculation formula for calculating the nonlinear weight is: Among them, ω l is the nonlinear weight, τ is the transition value in the calculation process, γ l is the linear weight, β l is a smooth indicator, ε=10 -6 .

8. The method according to claim 7, characterized in that: In step 7, the numerical flux splitting f is obtained. + (U) at point The calculation formula of the TWENO reconstruction value at is:

9. The method according to claim 1, characterized in that: The third-order TVD Runge-Kutta discrete formula in step 2 is: Among them, U (1) ,U (2) is the intermediate transition value, Δt is the time step, the superscript n represents the nth time layer, L(U n ),L(U (1) ),L(U (2) ) is -f(U) x An approximation of the higher-order spatial discretization form of .

10. The method according to claim 1, characterized in that: In the step three, the space-time fully discrete finite difference format is an iterative formula for the time layer. The initial state value is known, and the approximate value of the next time layer is calculated through the iterative formula to obtain the numerical simulation value of the flow field in the calculation area at the termination time.

Images