Method for solving physical quantities in fluid problems

Abstract

The application discloses a kind of physical quantity solving method in fluid problem, in actual flow problem calculation region, rectangular coordinate system is generated and grid is divided;Based on the calculation region after division, the spatial discretization equation of Euler or Navier-Stokes equation set describing fluid motion is obtained;The grid points in the calculation region are classified and marked, to obtain smooth grid points and discontinuous grid points;Linear format is used at smooth grid point respectively and non-equidistant finite difference WENO format is used at discontinuous grid point to carry out spatial discretization to its spatial derivative term;The classical Runge-Kutta time-discrete method is used to the time direction of fluid mechanics control equation set to advance, to obtain the physical quantity value at grid discrete point at any time, and the state of each time in actual flow problem is numerically simulated.The application can realize accurate, efficient and robust solution of compressible flow problem, and greatly improve its calculation efficiency.

Description

Technical Field

[0001] This invention belongs to the field of fluid mechanics calculation technology, specifically relating to a method for solving physical quantities in fluid problems. Background Technology

[0002] High-precision numerical algorithms possess low numerical dissipation and, under the same conditions, can yield more refined flow structures compared to lower-order schemes. Therefore, they have always been a key research topic in computational fluid dynamics, computational mathematics, computational electromagnetics, and computational acoustics. Among numerous high-precision algorithms, the WENO scheme plays a crucial role. The WENO scheme was first proposed by Liu et al. in 1994 and refined by Jiang and Shu in 1996, who provided a unified framework for arbitrary high-order WENO schemes. Based on this framework, many scholars have conducted subsequent research. For example, Balsara et al. constructed monotone-preserving WENO schemes with 7th-13th order accuracy; Henrick et al. designed the WENO-M scheme to recover the numerical accuracy of traditional WENO schemes near extreme points; and Dumbser et al. developed a finite-volume WENO method for unstructured mesh systems and implemented the engineering application of unstructured WENO schemes using UCNS3D code developed by Tsoutsanis et al. The construction process of the aforementioned WENO schemes is relatively complex, involving tedious linear weight calculations. Therefore, in 2016, Zhu et al. designed a scheme called the non-equidistant WENO. The linear weights of this scheme can be any set of positive numbers that sum to one, which greatly simplifies the implementation process (especially in unstructured grid systems) and has better steady-state convergence for steady problems.

[0003] Undoubtedly, the WENO scheme is a popular high-precision shock-capturing scheme with excellent engineering application prospects. However, its computational efficiency is low due to the eigenvalue decomposition and complex nonlinear weight calculations involved in its construction. Therefore, how to inherit the excellent characteristics of the WENO scheme while reducing its computational cost has always been a research hotspot. The hybrid WENO scheme is a feasible alternative, which classifies mesh cells through a discontinuity identifier, then uses a low-cost and efficient linear method in smooth regions, while using the expensive but highly effective WENO scheme in discontinuous regions. This combination allows numerical algorithms to fully leverage the advantages of both linear methods and nonlinear WENO schemes. However, as Pirozzoli stated, the role of the discontinuity identifier in classifying mesh cells is particularly important in the hybrid WENO scheme process. How to quickly, accurately, efficiently, and without human intervention distinguish between smooth and discontinuous flow field cells remains an open research topic. Summary of the Invention

[0004] The technical problem to be solved by this invention is to address the shortcomings of the prior art by providing a method for solving physical quantities in fluid problems. This method involves designing an accurate and efficient discontinuity identifier that requires no human intervention, classifying the grid cells within the computational domain, and discretizing the fluid dynamics control equations using both inexpensive and efficient linear methods and the non-equidistant WENO scheme with excellent shock wave capture capabilities. Finally, a third-order Runge-Kutta time discretization method with TVD properties is used to advance the control equations in the time direction, thereby obtaining the physical quantities at all grid points within the computational domain at any given time. The effectiveness of the method is verified through inviscid and viscous compressible flow problems.

[0005] To achieve the above-mentioned technical objectives, the technical solution adopted by the present invention is as follows:

[0006] A method for solving physical quantities in fluid problems includes the following steps:

[0007] Step 1: Within the computational domain of the actual flow problem, generate a Cartesian coordinate system and perform mesh generation;

[0008] Step 2: Based on the partitioned computational domain, obtain the spatial discretization equations of the Euler or Navier-Stokes equations describing fluid motion.

[0009] Step 3: Based on the spatial discretization equation, classify and label the grid points in the computational region to obtain smooth grid points and discontinuous grid points;

[0010] Step 4: Discretize the spatial derivative terms at smooth grid points using the linear scheme and at discontinuous grid points using the non-equidistant finite difference WENO scheme, respectively.

[0011] Step 5: Based on Step 4, the classic Runge-Kutta time discretization method is used to advance the fluid dynamics control equations in the time direction, thereby obtaining the physical quantity values ​​at the grid discrete points at any time, and then numerically simulating the state of the actual flow problem at each time.

[0012] To optimize the above technical solution, the specific measures also include:

[0013] The aforementioned practical flow problems include the front step reflection problem, the Rayleigh-Taylor instability problem, the viscous shock tube problem, the three-dimensional shock-shear layer interaction problem, and the Taylor-Green vortex evolution problem.

[0014] In step one above, a Cartesian coordinate system is constructed for the computational domain, and a mesh is generated. The corresponding mesh point coordinates are denoted as I. i,j,k =(x i ,yj ,z k ), where x, y, z represent spatial variables, and subscripts i, j, k represent grid point numbers.

[0015] In step two above, the Euler or Navier-Stokes equations describing fluid motion are:

[0016]

[0017] Where t represents the time variable, x = (x, y, z) represents the spatial variable vector, and U(x, t) = (ρ, ρu, ρv, ρw, E) T Let T denote the vector of conserved variables;

[0018] The total energy is represented by ρ, u, v, w, and e, which are functions of the time variable t and the spatial variables x, y, and z, respectively, representing the fluid density, velocity in the x-direction, velocity in the y-direction, velocity in the z-direction, and internal energy.

[0019] p = (γ-1)ρe represents pressure, where γ represents the specific heat constant;

[0020] ▽ denotes the divergence operator, σ is the conversion symbol between the Euler equation and the Navier-Stokes equation, σ=0 represents the Euler equation, and σ=1 represents the Navier-Stokes equation;

[0021] Re is the Reynolds number. Represents the convection flow vector. This represents the viscous flux vector, with superscripts x, y, and z indicating the corresponding spatial directions, and subscripts c and v indicating convection and viscosity, respectively.

[0022] F c (U) and F v Each component in (U,▽U) and The specific expression is

[0023]

[0024]

[0025]

[0026]

[0027]

[0028]

[0029] Where, τxx τ yy τ zz Let τ represent the normal stresses in the x, y, and z directions, respectively. xy τ xz τ yz The table represents the shear stress in the xy, xz, and yz planes, q x q y q z Represents the heat flux in the x, y, and z directions;

[0030] Based on the mesh in a spatial rectangular coordinate system, the following spatial discretization equations can be obtained.

[0031]

[0032] Among them, U i,j,k This indicates that the conserved variable U(x,t) is at grid point (x i ,y j ,z k Discrete values ​​at ) Represents the flow rate in the x-direction. At the grid boundary (x i±1 / 2 ,y j ,z k Discrete values ​​at ) F represents the flux in the y-direction. c y At the grid boundary (x i ,y j±1 / 2 ,z k Discrete values ​​at ) Represents the flow rate in the z-direction At the grid boundary (x) i ,y j ,z k±1 / 2 Discrete values ​​at ) Represents viscous flux in the x-direction At the grid boundary (x) i±1 / 2 ,y j ,z k Discrete values ​​at ) Represents viscous flux in the y direction At the grid boundary (x) i ,y j±1 / 2 ,z k Discrete values ​​at ) Represents the viscous flux in the z-direction At the grid boundary (x) i ,y j ,z k±1 / 2 The discrete value at ).

[0033] In step three above, the specific process of classifying and marking grid points within the computational region is as follows:

[0034] Step 31: For the target cell, i.e., target grid point I i,j,k First, fix the y and z directions, that is, keep the subscripts j and k unchanged, and select five adjacent elements in the x direction to form the template. We obtain a quartic polynomial p1(x) that satisfies the following conditions:

[0035]

[0036] Where Δx represents the spatial grid step size in the x-direction;

[0037]

[0038] in, Let a represent the basis function. n n = 0, 1, 2, 3, 4 represents the polynomial coefficients and the expression is as follows:

[0039]

[0040]

[0041]

[0042]

[0043]

[0044] Step 32: Find the first derivative of the polynomial p1(x) to obtain the derivative polynomial p1′(x):

[0045]

[0046] The superscript ' indicates the derivative symbol;

[0047] Step 33: Find the weak forms of the maximum and minimum values ​​of the derivative polynomial p1′(x);

[0048] Step 34: Classify all mesh subdivision points according to the weak form of the maximum and minimum values;

[0049] If (MAX1±Δx) and (MIN1±Δx) have opposite signs or (MAX2±Δx) and (MIN2±Δx) have opposite signs, then the grid point is marked as a discontinuous grid point containing discontinuities, denoted as Flag=1; otherwise, if they have the same sign, it is marked as a smooth grid point, denoted as Flag=0.

[0050] Step 35: Fix the x and z directions and the x and y directions respectively, and repeat steps 31-34 to classify and label all units.

[0051] In step 33 above, the derivative polynomial in In the interval [x i-5 / 2 ,x i ] and [x i ,x i+5 / 2 Each interval has monotonicity, so its maximum and minimum values ​​are obtained at the endpoints of its interval. Specifically:

[0052] In the interval [x i-5 / 2 ,x i [Contains]

[0053]

[0054]

[0055] In the interval [x i ,x i+5 / 2 ]have

[0056]

[0057]

[0058] Step four above is as follows:

[0059] Step 41: Fix the y and z directions, i.e., keep the subscripts j and k unchanged. For smooth grid points, use a linear format to obtain the convective flux at the grid boundary. The fifth-order approximation is:

[0060]

[0061]

[0062] Step 42: For discontinuous grid points, obtain the flux using a nonlinear, non-equidistant WENO scheme. The fifth-order approximation;

[0063] Step 43: Fix the x and z directions and the x and y directions respectively, and repeat steps 41-42 to achieve spatial discretization of the convection terms of all grid cells.

[0064] Step 42 above includes:

[0065] a) First, given the polynomial p1(x), for the target unit I... i,j,k And its two grid points to the left and right, and two more templates S are selected. (2) ={I i-1,I i} and S (3) ={I i ,I i+1 We can obtain two linear polynomials p2(x) and p3(x) respectively, which satisfy the following conditions.

[0066]

[0067]

[0068] The specific expressions are as follows:

[0069]

[0070]

[0071] b) Randomly select three positive numbers whose sum is 1 as linear weights, denoted as γ1, γ2, γ3;

[0072] c) Based on the three existing polynomials p1(x), p2(x), and p3(x), calculate their corresponding smoothness factors using the following formulas:

[0073]

[0074] Where n represents the index of the corresponding polynomial, α represents the summation index, and r represents the degree of the corresponding polynomial. Represents the polynomial p n Find the α-th derivative of (x) with respect to the independent variable x;

[0075] d) Calculate the nonlinear weights based on the linear weights and the smoothing factor, as follows:

[0076]

[0077] Where τ and For the transition value in the calculation, ε = 10 -6 To prevent small quantities with a denominator of zero;

[0078] e) Obtain the convective flux at the grid boundary using the non-equidistant WENO format. The fifth-order nonlinear approximation:

[0079]

[0080]

[0081] In step five above, assuming the time step is Δt, t n Indicates the nth time level. Represents the grid point (x) at time n. i ,y j,z k U(x) at point i ,y j ,z k ,t n The value of ) in t n The spatial derivatives of the fluid dynamics governing equations are discretized in step four at the time level, and denoted as... Then, the semi-discrete finite difference scheme is discretized using the Runge-Kutta time discretization method that satisfies the TVD property, thereby obtaining the spatiotemporally fully discrete finite difference scheme, as follows:

[0082]

[0083] in, and To calculate intermediate transition values, the spatiotemporal fully discrete finite difference scheme obtained after discretization in the spatial and temporal directions is an iterative formula for the time layer. Given the initial flow density ρ, x-direction velocity u, y-direction velocity v, z-direction velocity w, and internal energy e, the numerical approximation of the physical quantities of the flow field at a certain moment or in a steady state can be obtained by continuously iterating according to the iterative formula.

[0084] The present invention has the following beneficial effects:

[0085] To address the issue that existing non-equidistant WENO schemes can accurately simulate compressible flow problems, but their computational cost is relatively high, the hybrid non-equidistant WENO scheme of this invention can inherit all the excellent characteristics of the non-equidistant WENO scheme itself, while significantly improving its computational efficiency (roughly estimated to be about 30%).

[0086] Furthermore, this invention classifies computational grid points through a novel hybridization strategy, enabling the use of an efficient linear scheme to simulate the refined structure of the flow field in smooth regions, and a non-equidistant WENO scheme with essentially non-oscillatory characteristics to simulate discontinuous structures in discontinuous regions. This achieves accurate, efficient, and robust solutions to compressible flow problems. Finally, this invention demonstrates its superiority in simulating fluid mechanics problems through various inviscid and viscous compressible flow problems.

[0087] The discontinuity identifier developed in this invention has the following advantages compared to existing research results:

[0088] 1) It reduces mesh dependency; its construction process is only related to the reconstruction polynomial of the target element, without needing to consider the mesh type.

[0089] 2) It does not contain any manual parameters, and there is no need to adjust parameters according to specific problems, which greatly improves its practicality;

[0090] 3) It can automatically, accurately, and efficiently distinguish between smooth and discontinuous elements, thus effectively combining linear methods and nonlinear schemes. Attached Figure Description

[0091] Figure 1 Density contour map of the front step reflection problem;

[0092] Figure 2 A diagram showing the grid cells for the front step problem;

[0093] Figure 3 Density contour map of the inviscid Rayleigh-Taylor instability problem;

[0094] Figure 4 A diagram showing the grid cell markings for the inviscid Rayleigh-Taylor instability problem;

[0095] Figure 5 Density contour map for the viscous Rayleigh-Taylor instability problem;

[0096] Figure 6 A diagram showing the grid cell markings for the viscous Rayleigh-Taylor instability problem;

[0097] Figure 7 Density contour plot for a viscous shock tube problem;

[0098] Figure 8 A diagram showing the grid cell markings for the viscous shock tube problem;

[0099] Figure 9 The flow field structure diagram for the three-dimensional viscous shock wave-shear layer interaction problem;

[0100] Figure 10 The curves showing the kinetic energy change over time in the three-dimensional inviscid / viscous Taylor-Green eddy evolution problem;

[0101] Figure 11 This is a flowchart of a method for solving physical quantities in a fluid problem. Detailed Implementation

[0102] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0103] Although the steps in this invention are arranged by reference numerals, this is not intended to limit the order of the steps. Unless the order of the steps is explicitly stated or the execution of a step requires other steps as a basis, the relative order of the steps can be adjusted. It is understood that the term "and / or" as used herein refers to and covers any and all possible combinations of one or more of the associated listed items.

[0104] like Figure 11 As shown, this invention provides a method for solving physical quantities in fluid problems. It establishes a spatial rectangular coordinate system, divides the computational domain using a uniform grid, and constructs a fifth-order precision hybrid fundamental weighted non-oscillating scheme within the formed computational grid to perform fluid dynamics simulation within the computational domain. Specifically, it includes the following steps:

[0105] Step 1: Within the computational domain of the actual flow problem, generate a Cartesian coordinate system and perform mesh generation;

[0106] Step 2: Based on the partitioned computational domain, obtain the spatial discretization equations of the Euler or Navier-Stokes equations describing fluid motion.

[0107] Step 3: Based on the spatial discretization equation, classify and label the grid points in the computational region to obtain smooth grid points and discontinuous grid points;

[0108] Step 4: Discretize the spatial derivative terms at smooth grid points using the linear scheme and at discontinuous grid points using the non-equidistant finite difference WENO scheme, respectively.

[0109] Step 5: Based on Step 4, the classic Runge-Kutta time discretization method is used to advance the fluid dynamics control equations in the time direction, thereby obtaining the physical quantities such as flow field density, velocity in each direction, pressure and energy at any time at the grid discrete point, and then numerically simulating the state of the actual flow problem at each time.

[0110] In the embodiments, the actual flow problems include the front step reflection problem, the Rayleigh-Taylor instability problem, the viscous shock tube problem, the three-dimensional shock-shear layer interaction problem, and the Taylor-Green vortex evolution problem.

[0111] In this embodiment, in step one, a Cartesian coordinate system is constructed for the computational domain, and a mesh is generated, with the corresponding mesh point coordinates denoted as I. i,j,k =(x i ,y j ,z k ), where x, y, z represent spatial variables, and subscripts i, j, k represent grid point numbers.

[0112] In step two, the Euler or Navier-Stokes equations describing fluid motion are:

[0113]

[0114] Where t represents the time variable, x = (x, y, z) represents the spatial variable vector, and U(x, t) = (ρ, ρu, ρv, ρw, E) T Let T denote the vector of conserved variables;

[0115] The total energy is represented by ρ, u, v, w, and e, which are functions of the time variable t and the spatial variables x, y, and z, respectively, representing the fluid density, velocity in the x-direction, velocity in the y-direction, velocity in the z-direction, and internal energy.

[0116] p = (γ-1)ρe represents pressure, where γ represents the specific heat constant (γ = 1.4 in gases);

[0117] ▽ denotes the divergence operator, σ is the conversion symbol between the Euler equation and the Navier-Stokes equation, σ=0 represents the Euler equation, and σ=1 represents the Navier-Stokes equation;

[0118] Re is the Reynolds number. Represents the convection flow vector. This represents the viscous flux vector, with superscripts x, y, and z indicating the corresponding spatial directions, and subscripts c and v indicating convection and viscosity, respectively.

[0119] F c (U) and F v Each component in (U,▽U) and The specific expression is

[0120]

[0121]

[0122]

[0123]

[0124]

[0125]

[0126] Where, τ xx τ yy τ zz Let τ represent the normal stresses in the x, y, and z directions, respectively. xyτ xz τ yz The table represents the shear stress in the xy, xz, and yz planes, q x q y q z Represents the heat flux in the x, y, and z directions;

[0127] It is worth noting that all of the above physical quantities are functions of time and space variables, that is, they change with time and spatial location.

[0128] Based on the mesh in a spatial rectangular coordinate system, the following spatial discretization equations can be obtained.

[0129]

[0130] Among them, U i,j,k This indicates that the conserved variable U(x,t) is at grid point (x i ,y j ,z k Discrete values ​​at ) Represents the flow rate in the x-direction. At the grid boundary (x) i±1 / 2 ,y j ,z k Discrete values ​​at ) Indicates the flow rate in the y-direction At the grid boundary (x) i ,y j±1 / 2 ,z k Discrete values ​​at ) Represents the flow rate in the z-direction At the grid boundary (x) i ,y j ,z k±1 / 2 Discrete values ​​at ) Represents viscous flux in the x-direction At the grid boundary (x) i±1 / 2 ,y j ,z k Discrete values ​​at ) Represents viscous flux in the y direction At the grid boundary (x) i ,y j±1 / 2 ,z k Discrete values ​​at ) Represents the viscous flux in the z-direction At the grid boundary (x) i ,y j ,z k±1 / 2 The discrete value at ).

[0131] In step three, the specific process of classifying and marking the grid points within the computational area is as follows:

[0132] Step 31: For the target cell, i.e., target grid point I i,j,k First, fix the y and z directions, i.e., keep the subscripts j and k unchanged (which can be ignored for now), and select five adjacent elements in the x direction to form the template. We obtain a quartic polynomial p1(x) that satisfies the following conditions:

[0133]

[0134] Where Δx represents the spatial grid step size in the x-direction; its specific expression is as follows:

[0135]

[0136] in, Let a represent the basis function. n n = 0, 1, 2, 3, 4 represents the polynomial coefficients and the expression is as follows:

[0137]

[0138]

[0139]

[0140]

[0141]

[0142] Step 32: Find the first derivative of the polynomial p1(x) to obtain the derivative polynomial p1′(x):

[0143]

[0144] The superscript ' indicates the derivative symbol;

[0145] Step 33: Find the weak forms of the maximum and minimum values ​​of the derivative polynomial p1′(x);

[0146] It is easy to see that in the derivative polynomial In the interval [x i-5 / 2 ,x i ] and [x i ,x i+5 / 2 Each of the intervals has monotonicity, so its maximum and minimum values ​​must be obtained at the endpoints of the interval.

[0147] First, in the interval [x i-5 / 2 ,x i [Contains]

[0148]

[0149]

[0150] Secondly, in the interval [x i ,x i+5 / 2 ]have

[0151]

[0152]

[0153] Step 34: Classify all mesh subdivision points according to the weak form of the maximum and minimum values;

[0154] If (MAX1±Δx) and (MIN1±Δx) have opposite signs or (MAX2±Δx) and (MIN2±Δx) have opposite signs, then the grid point is marked as a discontinuous grid point containing discontinuities, denoted as Flag=1; otherwise, if they have the same sign, it is marked as a smooth grid point, denoted as Flag=0.

[0155] Step 35: Fix the x and z directions (i.e., keep the subscripts i and k unchanged) and the x and y directions (i.e., keep the subscripts i and j unchanged), and repeat steps 31-34 to classify and label all units.

[0156] Step four is described in detail below:

[0157] Step 41: Fix the y and z directions, i.e., keep the subscripts j and k unchanged. For the points identified as smooth grid points in Step 3, use a linear format to obtain the convective flux at the grid boundary. The fifth-order approximation; temporarily ignoring the subscripts j and k, we have

[0158]

[0159]

[0160] Step 42: For points identified as discontinuous grid points, the flux is obtained using a nonlinear, non-equidistant WENO scheme. The fifth-order approximation; temporarily ignoring the subscripts j and k, we have

[0161] a) First, given the polynomial p1(x), for the target unit I... i,j,k And its two grid points to the left and right, and two more templates S are selected. (2) ={I i-1 ,I i} and S (3) ={I i ,I i+1 We can obtain two linear polynomials p2(x) and p3(x) respectively, which satisfy the following conditions.

[0162]

[0163]

[0164] The specific expressions are as follows:

[0165]

[0166]

[0167] b) Randomly select three positive numbers whose sum is 1 as linear weights, denoted as γ1, γ2, γ3;

[0168] c) Based on the existing three polynomials p1(x), p2(x), and p3(x), calculate their corresponding smoothness factors, denoted as β1, β2, and β3, respectively, using the following formulas:

[0169]

[0170] Where n represents the index of the corresponding polynomial, α represents the summation index, and r represents the degree of the corresponding polynomial. Represents the polynomial p n Find the α-th derivative of (x) with respect to the independent variable x;

[0171] d) Calculate the nonlinear weights based on the linear weights and the smoothing factor, as follows:

[0172]

[0173] Where τ and For the transition value in the calculation, ε = 10 -6 To prevent small quantities with a denominator of zero;

[0174] e) Obtain the convective flux at the grid boundary using the non-equidistant WENO format. The fifth-order nonlinear approximation:

[0175]

[0176]

[0177] Step 43: Fix the x and z directions (i.e., keep the indices i and k unchanged) and the x and y directions (i.e., keep the indices i and j unchanged), and repeat steps 41-42 to achieve spatial discretization of the convection terms of all grid cells.

[0178] In step five, assuming the time step is Δt, t n Indicates the nth time level. Represents the grid point (x) at time n. i ,yj ,z k U(x) at point i ,y j ,z k ,t n The value of ) in t n The spatial derivatives of the fluid dynamics governing equations are discretized in step four at the time level, and denoted as... Then, the semi-discrete finite difference scheme is discretized using the Runge-Kutta time discretization method that satisfies the TVD property, thereby obtaining the spatiotemporally fully discrete finite difference scheme, as follows:

[0179]

[0180] in, and To calculate intermediate transition values, the spatiotemporal fully discrete finite difference scheme obtained after discretization in the spatial and temporal directions is an iterative formula for the time layer. Given the initial flow state density ρ, x-direction velocity u, y-direction velocity v, z-direction velocity w, internal energy e, and other physical quantities, the numerical approximations of the physical quantities of the flow field at a certain moment or in a steady state can be obtained by continuously iterating according to the iterative formula.

[0181] Five examples are given below as specific embodiments of the method disclosed in this invention.

[0182] Example 1: Front step reflection problem.

[0183] Consider a two-dimensional Euler equation within the region [0,3]×[0,1], with a step of height 0.2 at x=0.6. The initial condition is a free flow with a Mach number of 3. The upper boundary and the step are reflection boundaries, the left boundary is the incoming flow boundary, and the right boundary is the outgoing flow boundary. The computational mesh is 1200×400, and the computation terminates at t=4. Figure 1 The calculation results of the non-equidistant WENO scheme and the hybrid non-equidistant WENO scheme under the same conditions are presented. It can be seen that the method of the present invention can well simulate the complex flow phenomena involved in this problem, including incident shock waves, reflected shock waves, secondary shock waves, and their interactions. Figure 2 The discontinuous grid points identified by the method of this invention are presented. For this problem, the original non-equidistant WENO format takes 16447 seconds to compute, while the method of this invention only takes 11562 seconds, improving efficiency by 30%.

[0184] Example 2: Rayleigh-Taylor instability problem.

[0185] This problem involves the complex mixing phenomenon of light and heavy fluids under the influence of gravity. The Euler or Navier-Stokes equations are considered in the region [0,0.25]×[0,1], with the following initial conditions:

[0186]

[0187] in It is the speed of sound. Figure 3 and Figure 5 Density contour lines of the original non-equidistant WENO scheme and the hybrid non-equidistant WENO scheme of this invention are given for both non-viscous and viscous cases. Figure 4 and Figure 6 The grid points requiring a nonlinear format identified by the method of this invention are presented. The results show that the method of this invention can effectively capture flow instability interfaces and complex flow phenomena, and can accurately identify discontinuous elements. In terms of computational efficiency, the original method requires 7898 seconds and 9925 seconds in inviscid and viscous cases, respectively, while the method of this invention requires only 5268 seconds and 7062 seconds, respectively, improving computational efficiency by approximately 30%.

[0188] Example 3: Problems with viscous shock tubes.

[0189] This problem involves the interaction between shock waves and the viscous boundary layer. Due to the interaction between the shock wave and the lower and right walls, complex multi-scale structures such as shock waves, secondary shock waves, flow separation, and contact discontinuities are generated. The Navier-Stokes equations are considered within the computational domain [0,1]×[0,0.5], with a Reynolds number of 200 and a computation termination time of 1. The initial state is as follows:

[0190]

[0191] Figure 7 Density contour results obtained in the original non-equidistant WENO and hybrid non-equidistant WENO schemes are presented in an 800×400 grid. Figure 8 The discontinuous units identified by the method of this invention are presented. Results show that the method of this invention can accurately simulate this type of viscous problem. In terms of computational efficiency, the original non-equidistant WENO scheme requires 12784 seconds, while the hybrid non-equidistant WENO scheme requires only 9138 seconds, representing an efficiency improvement of 28.5%.

[0192] Example 4: Three-dimensional shock-shear layer interaction problem.

[0193] This problem describes the flow phenomenon generated by the interaction between an obliquely upward incident shock wave and a free shear layer. The computational domain is [0,200]×[-20,20]×[-20,20], the computational mesh is 600×120×120, the computation ends at t=120, the Reynolds number is 500, and the governing equations are the three-dimensional Navier-Stokes equations. Figure 9 Density contour plots calculated using two methods are presented. The results show that the method of this invention can effectively simulate the incident shock wave, reflected shock wave, and downstream vortex structure generated after the shear layer action involved in this problem. Under the same conditions, the original non-equidistant WENO scheme requires 262,314 seconds, while the method of this invention requires 180,295 seconds, representing a 31% improvement in efficiency.

[0194] Example 5: Taylor-Green vortex evolution problem.

[0195] The computational domain is [0,2π]×[0,2π]×[0,2π], the computational grid is 128×128×128, and the initial expression for the vortex is as follows:

[0196] ρ0 = 1,

[0197]

[0198]

[0199] w0(x,y,z)0,

[0200]

[0201] Where L = 1, V0 = 1, p0 = 100. Figure 10 The kinetic energy versus time curves for inviscid and viscous cases are presented. As can be seen from the figures, the method of this invention is closer to the reference solution than the original non-equidistant WENO method. Furthermore, in terms of computational efficiency, the original method takes 137,669 seconds and 327,132 seconds for inviscid and viscous cases, respectively, while the method of this invention requires only 74,925 seconds and 219,358 seconds, representing efficiency improvements of 45% and 33%, respectively.

[0202] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

[0203] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A method for solving physical quantities in fluid problems, characterized in that, Includes the following steps: Step 1: Within the computational domain of the actual flow problem, generate a Cartesian coordinate system and perform mesh generation; Step 2: Based on the partitioned computational domain, obtain the spatially discretized equations of the Euler or Navier-Stokes equations describing fluid motion; wherein, the Euler or Navier-Stokes equations describing fluid motion are: in, Represents a time variable. Represents a vector of spatial variables. Represents a vector of conserved variables. Indicates transpose; It means always, , , , , They are all about time variables. and spatial variables , , The functions represent fluid density, Directional velocity, Directional velocity, Directional velocity and internal energy; Indicates pressure, among which Represents the specific heat ratio constant; Denotes the divergence operator, The conversion symbol between the Euler equation and the Navier-Stokes equation. Representing the Euler equation, Representing the Navier-Stokes equations; For Reynolds numbers, Represents the convection flow vector. Represents the viscous flux vector, with superscript , , Indicates the corresponding spatial direction, subscript and Indicates convection and viscosity; and Each component in and The specific expression is in, , , They represent , , Normal stress in the direction, , , Tables are displayed , , Shear stress in a plane , , express , , Heat flux in the direction; Based on the mesh in a spatial rectangular coordinate system, the following spatial discretization equations can be obtained. in, Representing conserved variables At grid points discrete values ​​at that point express Direction to Flow At grid boundaries discrete values ​​at that point express Direction to Flow At grid boundaries discrete values ​​at that point express Direction to Flow At grid boundaries discrete values ​​at that point express Directional viscous flux At grid boundaries discrete values ​​at that point express Directional viscous flux At grid boundaries discrete values ​​at that point express Directional viscous flux At grid boundaries Discrete values ​​at; Step 3: Based on the spatial discretization equation, classify and label the grid points within the computational region to obtain smooth grid points and discontinuous grid points; the specific process of classifying and labeling the grid points within the computational region is as follows: Step 31: For the target cell, i.e. the target grid point First fix direction and Direction, i.e., subscript and Unchanged, select The template consists of five adjacent units on the left and right sides. This yields a quartic polynomial. And it satisfies the following conditions: ， in, express Spatial grid step size in the direction; ， in, Denotes basis functions. Denotes the polynomial coefficients and their expression is as follows: Step 32: Find the polynomial The first derivative yields the derivative polynomial. : ， Among them, superscript Indicates the derivative symbol; Step 33: Find the derivative polynomial The weak forms of the maximum and minimum values; Step 34: Classify all mesh subdivision points according to the weak form of the maximum and minimum values; If satisfied and opposite signs or and If the signs are opposite, the grid point is marked as a discontinuous grid point containing discontinuities, denoted as . Conversely, if they have the same sign, they are marked as smooth grid points, denoted as... ; in, , For interval The weak forms of the maximum and minimum values; , for The weak forms of the maximum and minimum values; Step 35: Fix them separately direction and direction and direction and In this direction, repeat steps 31-34 to classify and label all units; Step 4: Discretize the spatial derivative terms at smooth grid points using the linear scheme and at discontinuous grid points using the non-equidistant finite difference WENO scheme, as follows: Step 41, Fix direction and Direction, i.e., subscript and For smooth grid points, the convective flux at the grid boundaries is obtained using a linear format, remaining unchanged. The fifth-order approximation is: Step 42: For discontinuous grid points, obtain the flux using a nonlinear, non-equidistant WENO scheme. The fifth-order approximation; Step 43, fix them respectively direction and direction and direction and By repeating steps 41-42, spatial discretization of convection terms in all grid cells can be achieved. Step 5: Based on Step 4, the classic Runge-Kutta time discretization method is used to advance the fluid dynamics control equations in the time direction, thereby obtaining the physical quantity values ​​at the grid discrete points at any time, and then numerically simulating the state of the actual flow problem at each time.

2. The method for solving physical quantities in fluid problems according to claim 1, characterized in that, The actual flow problems include the front step reflection problem, the Rayleigh-Taylor instability problem, the viscous shock tube problem, the three-dimensional shock-shear layer interaction problem, and the Taylor-Green vortex evolution problem.

3. The method for solving physical quantities in fluid problems according to claim 1, characterized in that, In step one, a Cartesian coordinate system is constructed for the computational domain, and a mesh is generated. The corresponding mesh point coordinates are as follows: ,in , , Indicates spatial variables, subscript , , Indicates the grid point number.

4. A method for solving physical quantities in fluid problems according to claim 1, characterized in that, In step 33, the derivative polynomial in In the interval and Each interval possesses monotonicity, therefore its maximum and minimum values ​​are obtained at the endpoints of its interval, specifically: In the interval Inside ， In the interval have 。 5. A method for solving physical quantities in fluid problems according to claim 1, characterized in that, Step 42 includes: a) First, the polynomial is known. For the target unit And its two grid points to the left and right, and two more templates are selected. and Two linear polynomials can be obtained respectively. and The following conditions must be met. ， ； The specific expressions are as follows: b) Randomly select three positive numbers whose sum is 1 as linear weights, denoted as... , , ; c) Based on the three existing polynomials , and The smoothness factor is calculated using the following formulas: in, Indicates the index of the corresponding polynomial. Indicates the summation index. Indicates the degree of the corresponding polynomial. Representing a polynomial For independent variable beg Second derivative; d) Calculate the nonlinear weights based on the linear weights and the smoothing factor, as follows: in, and These are transitional values ​​in the calculation. To prevent small quantities with a denominator of zero; e) Obtain the convective flux at the grid boundary using the non-equidistant WENO format. The fifth-order nonlinear approximation: 。 6. The method for solving physical quantities in a fluid problem according to claim 1, characterized in that, In step five, it is assumed that the time step is... , Indicates the first Time layer, Indicates the first Time layer grid points Place The value, in The spatial derivatives of the fluid dynamics governing equations are discretized in step four at the time level, and denoted as... Then, the semi-discrete finite difference scheme is discretized using the Runge-Kutta time discretization method that satisfies the TVD property, thereby obtaining the spatiotemporally fully discrete finite difference scheme, as follows: in, and To calculate the intermediate transition values, the spatiotemporal fully discrete finite difference scheme obtained after discretization in both spatial and temporal directions is an iterative formula with respect to the time layer, given the initial flow density. , Directional velocity , Directional velocity , Directional velocity , internal energy Then, by continuously iterating according to the iterative formula, the numerical values ​​of the physical quantities of the flow field at a certain moment or in a steady state can be obtained.

Images