CN105760602A - Total flow field numerical simulation method for finite volume weighted essentially non-oscillatory scheme - Google Patents

Abstract

The invention discloses a total flow field numerical simulation method for a finite volume weighted essentially non-oscillatory scheme.Specific steps of the structured grid FVWENO5 (Finite Volume Weighted Essentially Non-oscillatory) scheme are given; then aiming at the Cartesian grid non-body-fitted feature, an object plane boundary is treated through a virtual unit immersed boundary method; a compressible inviscid flow problem of complex objects contained in a flow field can be solved under structured grid coordinates by combining the former operation.A numerical experiment is performed aiming at multiple classic complex analysis examples, and the reliability and effectiveness of the method can be fully verified.

Description

A kind of whole flow field method for numerical simulation of limited bulk weighted essential non-oscillatory schemes

Technical field

The invention discloses the whole flow field method for numerical simulation of a kind of limited bulk weighted essential non-oscillatory schemes, relate to calculating fluid Mechanics, partial differential equation numerical technology field.

Background technology

One main challenge of Fluid Mechanics Computation is to design effective numerical method to simulate whole flow field accurately.And be interrupted Existence, the shock wave in such as high-speed aerodynamics and contact discontinuity, logarithm value method proposes the highest requirement.In recent years On the basis of TVD concept, Harten A. and Osher S. proposes the concept of ENO first and constructs high-order ENO form.Its After, Liu X.D. et al. proposes the first weighted ENO scheme approached based on point, and the ENO form of r rank precision can be entered by this form Row weighted average processes and obtains r+1 rank numerical precision.1996, Jiang G.S. and Shu C.W. provided a new structure weighting The basic framework of ENO form, can make the numerical precision of form reach 2r-1 rank.Along with the development of finite difference method, limited body Long-pending method also emerges in an endless stream.Harten A. proposes a two-dimensional finite volume ENO form in the literature.Thereafter, Hu C.Q. and Shu C.W., the ordinary construction rule of two-dimensional finite volume WENO form is proposed.Qiu J.X. and Shu C.W. proposes under one-dimensional case HWENO form, and use it as RKDG limiter and can press flow field problem in order to calculate.Method is promoted the use of two by them In the structure of dimension limited bulk WENO form.These classical big multipair mesh qualities of numerical computation method require higher, it is impossible to Directly apply to computational flow comprises the pressed circumferential motion problem on complicated object plane border.Complex geometry profile object is contained for these Flow Field Calculation, use cartesian grid method will necessarily produce difform cutter unit, this results in the firm of equation system Property and produce the pseudo-vibration solved at object plane boundary.To this end, there is following several solution, mainly there is Mixed grid method, Integrated unit method, and immerse boundary method or dummy unit method.

Single order numerical method does not haves the numerical oscillation of non-physical but understands excessive floating strong discontinuity when capturing shock.Need for this Imported high precision numerical computations form carries out correlation values simulation.When zoning occurs complicated object plane, many based on structure Mesh design high-order computational scheme out is not directly applicable correlation values and calculates, and needs first to generate the body fitted grids of complexity, Reconstructing high-order computational scheme based on this body fitted grids, universality is poor.

Summary of the invention

The technical problem to be solved is: for the defect of prior art, it is provided that a kind of limited bulk weights substantially without shaking Swing the whole flow field method for numerical simulation of form, structured grid and height can be used uniformly across for various complex object circumferential motion problem whole flow fields Accuracy computation form.First the present invention provides structured grid limited bulk and adds weighted essentially non-oscillatory highly-accurate nephelometric titrimetry form, and For the non-fit characteristic of cartesian grid, use suitable boundary processing method to process Solid boundary condition, both can be had Effect combines and is applied to the viscous numerical simulation that can press circumferential motion problem of nothing with complex geometry profile object plane.

The present invention solves above-mentioned technical problem by the following technical solutions:

By " dimension in the whole flow field method for numerical simulation of a kind of limited bulk weighted essential non-oscillatory schemes, with finite difference By dimension fashion " mode structuring one-dimensional Polynomial Reconstructing have certain difference, the especially determination of smoothing factor has bigger Different.This way more conforms to conservation variable actual change under two-dimensional case in flow field.

Concrete steps include:

Step 1. is at the boundary point of object elementWithFunction u is reconstructed by place, wherein, and x, y For cartesian coordinate, i, j are respectively the mesh point sequence number on x, y direction；

1.1 use least square method strategy obtain on each Gauss point four times reconstruct multinomial Q (x, y)；

Large form T is divided into a plurality of little template by 1.2, in each little template, and reconstruction of function u at each Gauss point respectively Value；

Step 2. calculates half Discrete Finite volume format；

Step 3. uses three rank TVD Runge-Kutta time discrete formula to obtain space-time full discrete scheme.

As present invention further optimization scheme, non-Body fit to due to cartesian grid, the shearing list produced at object plane Unit, uses two kinds of dummy unit methods of ST and GBCM to process boundary condition, and concrete steps include:

Step one, at object plane boundary, determine virtual grid unit I₁Central point A is about the symmetric points B of object plane；

Step 2, determine the virtual grid unit I that B point falls into₂, and virtual grid unit I around₃,I₄,I₅；

Step 3, utilize virtual grid unit I₂,I₃,I₄,I₅, by interpolation formula, determine B point flow field value；

Step 4, determine dummy unit I by Solid boundary condition and momentum relationship₁With the relation of B point, thus obtain virtual list Unit I₁Flow field value.

As present invention further optimization scheme, in step 1.1, (x is y) that two-Dimensional Algebraic is many to described four reconstruct multinomial Q Item formula, (x, y) and function u is at object element I for Q₁₃Mean value equal:

$\frac{1}{\left|I_{13}\right|}\underset{I_{13}}{\∫}Q(x,y)dxdy={\stackrel{\‾}{u}}_{13}$

Wherein,Represent object element I₁₃Cell-average value, | I₁₃| represent object element I₁₃Area；

Further, (x y) meets the mean value of other unit superior function u of large form T under least square meaning to Q.

As present invention further optimization scheme, in step 1.2, the multinomial that each little template is constructed is: p_l(x, y), (l=1 ..., n), wherein p_l(x,y)∈span{1,x,y,x²,xy,y²,x²y,xy²,x²y², i.e. belong to these substrates and open into Linear space, subscript l represent correspondence little template sequence number, n is little template number；

Multinomial p_l(x, y) must be with u in same little template S_lThe mean value of upper unit is equal；

At object element I₁₃Gaussian quadrature nodeOn function u is reconstructed, and solve at Gauss pointLinearly weigh γ_l；

Subsequently, smoothing factor β is calculated_l, smoothing factor β_lIt is used for weighing multinomial p_l(x, y), (l=1 ..., n) at object element I₁₃On Smooth degree, its computing formula is:

${\mathrm{\β}}_l=\underset{\left|\mathrm{\α}\right|=1}{\overset{r}{\Σ}}\underset{I_{13}}{\∫}|I_{13}{|}^{\left|\mathrm{\α}\right|-1}{\left(D^{\mathrm{\α}}{p}_l\right(x,y\left)\right)}^{2}dxdy$

Wherein, α is multiple indexes, and D is derivation operator, and r is multinomial p_l(x, number of times y)；

Calculate nonlinear weight on this basis:

${\mathrm{\ω}}_l=\frac{{\stackrel{\‾}{\mathrm{\ω}}}_l}{\underset{k=1}{\overset{9}{\Σ}}{\stackrel{\‾}{\mathrm{\ω}}}_k},{\stackrel{\‾}{\mathrm{\ω}}}_l=\frac{{\mathrm{\γ}}_l}{{({10}^{-7}+{\mathrm{\β}}_l)}^2},(l=1,...,n)$

Wherein, γ_lIt is linear power, β_lIt is smoothing factor,For calculating process median, ω_lFor nonlinear weight.Finally, approximation Result is:

$u_{i+\frac{1}{2},j+\frac{1}{2\sqrt{3}}}\≈\underset{l=1}{\overset{n}{\Σ}}{\mathrm{\ω}}_l{p}_l(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}}).$

Repeat above step, at object element I₁₃Other Gaussian quadrature nodes on function u is reconstructed.

As present invention further optimization scheme, in step 2, the calculating formula of half Discrete Finite volume format is:

$\frac{d}{dt}{\stackrel{\‾}{u}}_{i,j}\left(t\right)=L\left(u\right)=-\frac{1}{\mathrm{\Δ}y\mathrm{\Δ}y}{\∫}_{\∂I_{i,j}}(f,g)*\stackrel{\→}{n}ds$

Wherein, ${\stackrel{\‾}{u_i}}_{,j}\left(t\right)=\frac{1}{\mathrm{\Δ}x\mathrm{\Δ}y}{\∫}_{I_{i,j}}u(x,y,t)dxdy,$ L is right-hand vector,For outer normal unit vector；

Line integral formula is rewritten as 2 Gauss quadrature formulas:

${\∫}_{\∂I_{i,j}}(f,g)*\stackrel{\→}{n}ds\≈|\∂I_{i,j}|\underset{l=1}{\overset{2}{\Σ}}{\mathrm{\ω}}_l\left(f\right(u\left(x_{i\±\frac{1}{2}},y_{j+G_l},t\right)),g(u\left(x_{i+G_l},y_{j\±\frac{1}{2}},t\right)\left)\right)*\stackrel{\→}{n}$

Wherein,AndReplaced by numerical flux.

As present invention further optimization scheme, in step 3, calculating formula is:

$\left\{\begin{array}{c}{u}^{\left(1\right)}=u^n+\mathrm{\Δ}tL\left(u^n\right)\\ u^{\left(2\right)}=\frac{3}{4}{u}^n+\frac{1}{4}{u}^{\left(1\right)}+\frac{1}{4}\mathrm{\Δ}tL\left(u^{\left(1\right)}\right)\\ u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}{u}^{\left(2\right)}+\frac{2}{3}\mathrm{\Δ}tL\left(u^{\left(2\right)}\right)\end{array}\right..$

The present invention uses above technical scheme compared with prior art, has following technical effect that limited for structured grid high accuracy Volume weighting essentially non oscillatory scheme with immerse border dummy unit method and be combined, be applied to have complex geometry profile without viscosity flow In dynamic problem numerical simulation, can effectively reduce computation scheme structure and calculate the complexity of mess generation.Conventional method is in the face of tool It is that difficulty is the most insurmountable when having double Mach reflection problem and a peripheral flow problem on true inclined-plane, but the present invention can be sane, Accurately, capturing to basic dead-beat intense shock wave, contact is discontinuous and keeps consistent higher order accuracy at the smooth domain solved simultaneously.

Accompanying drawing explanation

Fig. 1 is large form T schematic diagram.

Fig. 2 is that in cartesian grid, grid cell intersects schematic diagram with object plane border.

Fig. 3 is in embodiment one, double Mach reflection problem density isogram.

Fig. 4 is in embodiment one, double Mach reflection problem pressure isogram.

Fig. 5 is in embodiment two, processes the peripheral flow problem density isogram of boundary condition by ST method.

Fig. 6 is in embodiment two, processes the peripheral flow problem Mach number isogram of boundary condition by ST method.

Fig. 7 is in embodiment two, processes the peripheral flow problem density isogram of boundary condition by GBCM method.

Fig. 8 is in embodiment two, processes the peripheral flow problem Mach number isogram of boundary condition by GBCM method.

Detailed description of the invention

Embodiments of the present invention are described below in detail, and the example of described embodiment is shown in the drawings, the most extremely Same or similar label represents same or similar element or has the element of same or like function eventually.Below by ginseng The embodiment examining accompanying drawing description is exemplary, is only used for explaining the present invention, and is not construed as limiting the claims.

Below in conjunction with the accompanying drawings technical scheme is described in further detail:

One, limited bulk weighted essential non-oscillatory schemes

Consideration two dimension Conservation Law Equations:

$\left\{\begin{array}{c}{u}_t+f{\left(u\right)}_x+g{\left(u\right)}_y=0\\ u(x,y,0)=u_0(x,y)\end{array}\right.---\left(1\right)$

Wherein, u_tRepresenting the conservation variable u partial derivative to time t, x, y are cartesian coordinate, and f, g are respectively x, y direction Numerical flux, u₀Represent the conservation variable u initial value in the t=0 moment.

It is consistent for setting grid cell: x_i+1/2-x_i-1/2=Δ x, y_j+1/2-y_j-1/2=Δ y, wherein, i, j are respectively x, y side Mesh point sequence number upwards；

Unit center is: (x_i,y_j)=((x_i-1/2+x_i+1/2)/2,(y_j-1/2+y_j+1/2)/2)。

Definition two-dimensional cell is: I_i,j=I_i×I_j=[x_i-1/2,x_i+1/2]×[y_j-1/2,y_j+1/2], wherein I is grid cell.

Set Δ t as time step, tⁿ⁺¹=tⁿ+Δt；

Represent tⁿThe value of point on time horizon,Represent tⁿCell-average value on time horizon.Net template Figure is as shown in Figure 1.

Step 1. is at object element I_i,j=I₁₃Boundary pointWithFunction u is reconstructed by place.

The least square method strategy that step 1.1. uses Barth T. and Frederickson P. to propose obtains four reconstruct multinomials, fixed Justice be Q (x, y).(x, y) is two-Dimensional Algebraic multinomial to Q, needs to meet following condition and determines it, is allowed to necessary and function u at mesh Mark unit I₁₃Mean value equal:

$\frac{1}{\left|I_{13}\right|}\underset{I_{13}}{\∫}Q(x,y)dxdy={\stackrel{\‾}{u}}_{13},---\left(2\right)$

Wherein,Represent object element I₁₃Cell-average value, | I₁₃| represent object element I₁₃Area.Same restriction is multinomial (x y) must meet large form T other unit I under least square meaning to formula Q_l(l=1 ..., 12,14 ..., 25) superior function u flat Average.The present invention is by expression formulaIt is abbreviated as u_l。

Step 1.1.1. is at object element I₁₃Gaussian quadrature nodeThe value of place's reconstruct u.

$\begin{array}{c}Q(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=\frac{1}{163598400}\left(\right(5421910+920782\sqrt{3})u1-2(110690+\\ 2290053\sqrt{3})u10-1052580u11-6277980u12+144923100u13+\\ 202093080u14+5569260u15-2168980u16+1565654\sqrt{3}u16-\\ 8871640u17+2424762\sqrt{3}u17-3865120u18+4501120\sqrt{3}u18+\\ 15278600u19+5863358\sqrt{3}u19-5712530u2+1228611\sqrt{3}u2-\\ 221380u20+4580106\sqrt{3}u20+5421910u21-920782\sqrt{3}u21-\\ 5712530u22-1228611\sqrt{3}u22-4520390u23+898060\sqrt{3}u23+\\ 11426350u24+1596491\sqrt{3}u24-6653210u25-2996058\sqrt{3}u25-\\ 4520390u3-898060\sqrt{3}u3+11426350u4-1596491\sqrt{3}u4-\\ 6653210u5+2996058\sqrt{3}u5-2168980u6-1565654\sqrt{3}u6-\\ 8871640u7-2424762\sqrt{3}u7-3865120u8-4501120\sqrt{3}u8+\\ 15278600u9-5863358\sqrt{3}u9)\end{array}---\left(3\right)$

Step 1.1.2.～1.1.8. are similar with step 1.1.1, at object element I₁₃Other Gaussian quadrature nodes at reconstruct u value.

Step 1.2. constructs corresponding 9 multinomial p in following 9 little templates_l(x, y), (l=1 ..., 9), wherein p_l(x,y)∈span{1,x,y,x²,xy,y²,x²y,xy²,x²y², i.e. belonging to the linear space that these substrates are opened, subscript l represents right The little template sequence number answered:

S₁={ I₁,I₂,I₃,I₆,I₇,I₈,I₁₁,I₁₂,I₁₃, S₂={ I₆,I₇,I₈,I₁₁,I₁₂,I₁₃,I₁₆,I₁₇,I₁₈,

S₃={ I₁₁,I₁₂,I₁₃,I₁₆,I₁₇,I₁₈,I₂₁,I₂₂,I₂₃, S₄={ I₂,I₃,I₄,I₇,I₈,I₉,I₁₂,I₁₃,I₁₄,

S₅={ I₇,I₈,I₉,I₁₂,I₁₃,I₁₄,I₁₇,I₁₈,I₁₉, S₆={ I₁₂,I₁₃,I₁₄,I₁₇,I₁₈,I₁₉,I₂₂,I₂₃,I₂₄,

S₇={ I₃,I₄,I₅,I₈,I₉,I₁₀,I₁₃,I₁₄,I₁₅, S₈={ I₈,I₉,I₁₀,I₁₃,I₁₄,I₁₅,I₁₈,I₁₉,I₂₀,

S₉={ I₁₃,I₁₄,I₁₅,I₁₈,I₁₉,I₂₀,I₂₃,I₂₄,I₂₅}。S_lRepresent l little template, and haveThe most multinomial Formula p_l(x, y) must be with u in same little template S_lThe mean value of upper unit is equal.

Such as:

$\frac{1}{\left|I_l\right|}\underset{I_l}{\∫}{p}_1(x,y)dxdy={\stackrel{\‾}{u}}_l(l=1,2,3,6,7,8,11,12,13).---\left(4\right)$

Can obtain at each Gaussian quadrature node p_l(x_*,y_*), (l=1 ..., 9) explicit expression.

Step 1.2.1. is at object element I₁₃Gaussian quadrature nodeFunction u is reconstructed by place.

$\begin{array}{c}{p}_1(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=(2\sqrt{3}{u}_1+6(4+\sqrt{3})u_{11}-21(4+\sqrt{3})u_{12}+33(4+\sqrt{3})u_{13}\\ -7\sqrt{3}{u}_2+11\sqrt{3}{u}_3-8\sqrt{3}{u}_6+28\sqrt{3}{u}_7-44\sqrt{3}{u}_8)/72,\end{array}---\left(5\right)$

$\begin{array}{c}{p}_2(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=(24u_{11}-84u_{12}+132u_{13}+2\sqrt{3}{u}_{16}-7\sqrt{3}{u}_{17}+11\sqrt{3}{u}_{18}\\ -2\sqrt{3}{u}_6+7\sqrt{3}{u}_7-11\sqrt{3}{u}_8)/72,\end{array}---\left(6\right)$

$\begin{array}{c}{p}_3(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=(-6(-4+\sqrt{3})u_{11}+21(-4+\sqrt{3})u_{12}+132u_{13}-33\sqrt{3}{u}_{13}+8\sqrt{3}{u}_{16}\\ -28\sqrt{3}{u}_{17}+44\sqrt{3}{u}_{18}-2\sqrt{3}{u}_{21}+7\sqrt{3}{u}_{22}-11\sqrt{3}{u}_{23})/72,\end{array}---\left(7\right)$

$\begin{array}{c}{p}_4(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=(-3(4+\sqrt{3})u_{12}+15(4+\sqrt{3})u_{13}+24u_{14}+6\sqrt{3}{u}_{14}-\sqrt{3}{u}_2\\ +5\sqrt{3}{u}_3+2\sqrt{3}{u}_4+4\sqrt{3}{u}_7-20\sqrt{3}{u}_8-8\sqrt{3}{u}_9)/72,\end{array}---\left(8\right)$

$\begin{array}{c}{p}_5(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=(-12u_{12}+60u_{13}+24u_{14}-\sqrt{3}{u}_{17}+5\sqrt{3}{u}_{18}+2\sqrt{3}{u}_{19}\\ +\sqrt{3}{u}_7-5\sqrt{3}{u}_8-2\sqrt{3}{u}_9)/72,\end{array}---\left(9\right)$

$\begin{array}{c}{p}_6(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=\left(3\right(-4+\sqrt{3})u_{12}-15(-4+\sqrt{3})u_{13}+6(4-\sqrt{3})u_{14}-4\sqrt{3}{u}_{17}\\ +20\sqrt{3}{u}_{18}+8\sqrt{3}{u}_{19}+\sqrt{3}{u}_{22}-5\sqrt{3}{u}_{23}-2\sqrt{3}{u}_{24})/72,\end{array}---\left(10\right)$

$\begin{array}{c}{p}_7(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=(4\sqrt{3}{u}_{10}+6(4+\sqrt{3})u_{13}+60u_{14}+15\sqrt{3}{u}_{14}-3(4+\sqrt{3})u_{15}\\ +2\sqrt{3}{u}_3+5\sqrt{3}{u}_4-\sqrt{3}{u}_5-8\sqrt{3}{u}_8-20\sqrt{3}{u}_9)/72,\end{array}---\left(11\right)$

$\begin{array}{c}{p}_8(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=(\sqrt{3}{u}_{10}+24u_{13}+60u_{14}-12u_{15}+2\sqrt{3}{u}_{18}+5\sqrt{3}{u}_{19}\\ -\sqrt{3}{u}_{20}-2\sqrt{3}{u}_8-5\sqrt{3}{u}_9)/72,\end{array}---\left(12\right)$

$\begin{array}{c}{p}_9(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}})=(-6(-4+\sqrt{3})u_{13}-15(-4+\sqrt{3})u_{14}-3(4-\sqrt{3})u_{15}+8\sqrt{3}{u}_{18}\\ +20\sqrt{3}{u}_{19}-4\sqrt{3}{u}_{20}-2\sqrt{3}{u}_{23}-5\sqrt{3}{u}_{24}+\sqrt{3}{u}_{25})/72.\end{array}---\left(13\right)$

It is then determined linearly weigh γ₁,...,γ₉.Obtain at Gaussian quadrature nodeLinear power:

$\begin{array}{c}{\mathrm{\γ}}_1=\frac{-1+70\sqrt{3}}{3600\sqrt{3}},{\mathrm{\γ}}_2=\frac{11}{180},{\mathrm{\γ}}_3=\frac{1+70\sqrt{3}}{3600\sqrt{3}},\\ {\mathrm{\γ}}_4=\frac{210-\sqrt{3}}{1800},{\mathrm{\γ}}_5=\frac{11}{30},{\mathrm{\γ}}_6=\frac{210+\sqrt{3}}{1800},\\ {\mathrm{\γ}}_7=\frac{-1+70\sqrt{3}}{1200\sqrt{3}},{\mathrm{\γ}}_8=\frac{11}{60},{\mathrm{\γ}}_9=\frac{1+70\sqrt{3}}{1200\sqrt{3}}.\end{array}---\left(14\right)$

Then, calculate smoothing factor, be designated as β₁,...,β₉, they are corresponding to different templates S₁,...,S₉.Smoothing factor is used for weighing Amount multinomial p_l(x, y), (l=1 ..., 9) at object element I₁₃On smooth degree.

The computing formula of smoothing factor is:

${\mathrm{\β}}_l=\underset{\left|\mathrm{\α}\right|=1}{\overset{r}{\Σ}}\underset{I_{13}}{\∫}|I_{13}{|}^{\left|\mathrm{\α}\right|-1}{\left(D^{\mathrm{\α}}{p}_l\right(x,y\left)\right)}^{2}dxdy,---\left(15\right)$

Wherein, α is multiple indexes, and D is derivation operator, and r is multinomial p_l(x, number of times y).

Calculate nonlinear weight on this basis:

${\mathrm{\ω}}_l=\frac{{\stackrel{\‾}{\mathrm{\ω}}}_l}{\underset{k=1}{\overset{9}{\Σ}}{\stackrel{\‾}{\mathrm{\ω}}}_k},{\stackrel{\‾}{\mathrm{\ω}}}_l=\frac{{\mathrm{\γ}}_l}{{({10}^{-7}+{\mathrm{\β}}_l)}^2},(l=1,...,9).---\left(16\right)$

Wherein, γ_lIt is linear power obtained above, β_lIt is smoothing factor obtained above,For calculating process median, ω_lFor Nonlinear weight.Final approximation is:

$u_{i+\frac{1}{2},j+\frac{1}{2\sqrt{3}}}\≈\underset{l=1}{\overset{9}{\Σ}}{\mathrm{\ω}}_l{p}_l(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}}).---\left(17\right)$

Step 1.2.2. repeats above step, at object element I₁₃Other Gaussian quadrature nodes at equally function u is reconstructed.

Step 2. obtains half Discrete Finite volume format:

$\frac{d}{dt}{\stackrel{\‾}{u}}_{i,j}\left(t\right)=L\left(u\right)=-\frac{1}{\mathrm{\Δ}y\mathrm{\Δ}y}{\∫}_{\∂I_{i,j}}(f,g)*\stackrel{\→}{n}ds,---\left(18\right)$

WhereinL is right-hand vector,For outer normal unit vector.By line integral formula It is rewritten as 2 Gauss quadrature formulas:

${\∫}_{\∂I_{i,j}}(f,g)*\stackrel{\→}{n}ds\≈|\∂I_{i,j}|\underset{l=1}{\overset{2}{\Σ}}{\mathrm{\ω}}_l\left(f\right(u\left(x_{i\±\frac{1}{2}},y_{j+G_l},t\right)),g(u\left(x_{i+G_l},y_{j\±\frac{1}{2}},t\right)\left)\right)*\stackrel{\→}{n},---\left(19\right)$

ThereinAndReplaced by numerical flux.

Step 3. uses three rank TVD Runge-Kutta time discrete formula:

$\left\{\begin{array}{c}{u}^{\left(1\right)}=u^n+\mathrm{\Δ}tL\left(u^n\right),\\ u^{\left(2\right)}=\frac{3}{4}{u}^n+\frac{1}{4}{u}^{\left(1\right)}+\frac{1}{4}\mathrm{\Δ}tL\left(u^{\left(1\right)}\right),\\ u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}{u}^{\left(2\right)}+\frac{2}{3}\mathrm{\Δ}tL\left(u^{\left(2\right)}\right)\·\end{array}\right.---\left(20\right)$

Obtain space-time full discrete scheme.

Two, dummy unit method

Non-Body fit due to cartesian grid, the place intersected with object plane at grid necessarily leads to cut cells, and this can cause The Nonphysical Oscillation solved is produced at object plane.In order to solve this problem, present invention employs immersion border dummy unit method and process Solid boundary condition.As in figure 2 it is shown, in cartesian grid, grid cell intersects with object boundary.In figure, circle represents thing The grid cell in internal portion, i.e. dummy unit；Filled squares represents the dummy unit central point symmetric points about object plane border. We are it needs to be determined that dummy unit I in fig. 2₁Flow field value, then we are firstly the need of determining unit I₁A point in center is about thing The symmetric points B on border, face；It is then determined that the flow field value of B point；Finally by Solid boundary condition and momentum relationship etc., determine Relation between A point flow field value and B point flow field value.Determine A point flow field value and B point flow field value the most in two different ways Between relation.

2.1, ST method

For determining the relation between A point flow field value and B point flow field value, the simplest most straightforward approach is exactly to use general symmetry Reflective boundary condition ST, relational expression is as follows:

In above formula, ρ is density, and pre is pressure,For velocity,Represent the unit outer normal vector at object plane.The most true Fixed boundary condition disclosure satisfy that wall is without penetrating condition.

2.2, GBCM method

Dadone et al. proposes a kind of new boundary condition treatment method, and named CCST method for body fitted grids.? In CCST method, the flow field value on dummy unit point is to have Local Symmetric by the normal direction assumed near object plane to be distributed The vortex field of entropy and total enthalpy obtain.This flow model meets the following normal direction equation of momentum:

$\frac{\∂pre}{\∂n}=-\mathrm{\ρ}\frac{{\tilde{u}}^2}{R},---\left(22\right)$

Wherein, R is signed local radius of curvature, if the center of curvature is just at interior of articles, otherwise is negative；It is to cut To speed, meet wall without penetrating conditionOwing to forcing antisymmetric entropy and total enthalpy along the surface normal of object Normal derivativeThis entropy makes when irrotational flow with total enthalpy distribution, and entropy is zero with the normal derivative of total enthalpy, Crocco relation even also can be met when whirlpool occurs in object plane.CCST method is generalized to non-fit flute card by Dadone et al. That grid, and by named for the method GBCM method, as in figure 2 it is shown, in order to seek dummy unit I₁The flow field value of center A point, Give following expression formula:

$\left\{\begin{array}{c}{\mathrm{pre}}_A={\mathrm{pre}}_B-{\mathrm{\ρ}}_B({\tilde{u}}_B^2/R)\mathrm{\Δ}n,\\ {\mathrm{\ρ}}_A={\mathrm{\ρ}}_B{\left(\frac{{\mathrm{pre}}_A}{{\mathrm{pre}}_B}\right)}^{1/\mathrm{\γ}},\\ {\tilde{u}}_A^2={\tilde{u}}_B^2+\frac{2\mathrm{\γ}}{\mathrm{\γ}-1}(\frac{{\mathrm{pre}}_B}{{\mathrm{\ρ}}_B}-\frac{{\mathrm{pre}}_A}{{\mathrm{\ρ}}_A}),\\ {\tilde{v}}_A=-{\tilde{v}}_B\·\end{array}\right.---\left(23\right)$

Three, limited bulk weighted essential non-oscillatory schemes and the combination of dummy unit method

From the constitution step of above form, limited bulk weighted essential non-oscillatory schemes is a kind of limited bulk high accuracy numerical value Computation scheme.Finite Volume Method from the flow equation of integration conservation form, by the equation on grid cell is carried out from Dissipate and construct integral form discrete equation.Compared to finite difference method, Finite Volume Method is more flexible on grid processes, and keeps Perseverance is good.But, on structured grid development limited bulk weighted essential non-oscillatory schemes, solve double Mach reflection, It is difficult during the complicated interface circumferential motion problems such as peripheral flow.Because limited bulk weighted essential non-oscillatory schemes is to grid cell Need to use (including self) 25 grid cells altogether around when flow field value carries out Polynomial Reconstructing, especially at grid and object plane The place that border is intersected, problem just seems especially prominent.In the face of these problems, in double Mach reflection numerical simulation, traditional Solution is to make incident shock tilt, and inclined-plane just overlaps with grid lines, and numerical simulation is just apparent from and simple.Though so So can preferably capture the reflection case of shock wave, but the most inherently have different with experimental simulation. And in peripheral flow problem, traditional solution is to use hybrid grid or unstrctured grid, this increases net the most undoubtedly The difficulty that lattice generate.Structured grid limited bulk weighted essential non-oscillatory schemes is tied by the present invention with immersing border dummy unit method Close, effectively simulate the double Mach reflection problem with true inclined-plane, accurately capture the feelings of shock wave and ramp effect back reflection Condition.Similarly, High Mach number peripheral flow have also been obtained preferable numerical simulation result.Highlight limited bulk below Weighted essential non-oscillatory schemes is combined how to process Solid boundary condition with dummy unit method.Using, limited bulk weighting is basic When non-oscillatory scheme carries out Polynomial Reconstructing to the flow field value of some object element, need to use 25 grid cells around (big Template T, is shown in Fig. 1), in the place near object plane border, the grid cell in large form is positioned at inside object plane, as single in Fig. 2 Unit I₁。

The present invention first determines grid cell I₁A point in center is about the symmetric points B point of object plane, it is then determined that the grid list that B point falls into Unit I₂.So formula (21), the flow field value of some A, B in (23) just can replace by corresponding grid cell flow field value respectively. Numerical experiment shows, when B point flow field value I₂,I₃,I₄,I₅When the weighted average of four unit replaces, numerical simulation result is optimal.

Two examples specific embodiment as presently disclosed method is given below.

Embodiment one, double Mach reflection problem.This problem is after one intense shock wave of description is incident on the slope becoming 30 ° of angles with plane The change occurred, carrys out the intense shock wave that stream is Mach 2 ship 10.Zoning is taken as [0,3] × [0,2], starts fromShock wave It is perpendicular to x-axis.Primary data is:

$u(x,y,0)=\left\{\begin{array}{cc}{u}_L,& x\≤\sqrt{3}/12,\\ u_R,& x>\sqrt{3}/12,\end{array}\right.,$

Left and right therein state value: u_L=(8,66.009,0,563.544)^T,u_R=(1.4,0,0,2.5)^T.Boundary condition treatment is: left Right margin is taken as left and right state value respectively.Lower boundary: whenTime be reflecting boundary；Time, it is set to left shape State.Coboundary: as x ＞ g (t), takes right state；During x≤g (t), take left state.Wherein: Owing to, in this example, object plane border is straight line, radius of curvature is infinitely great, and at this moment ST method is consistent with GBCM method. It is engraved in the density isopleth of region [0,3] × [0,2] when that Fig. 3 reflects is t=0.2, when that Fig. 4 reflects is t=0.2, is engraved in region [0,3] The pressure isogram of × [0,2].

Embodiment two, peripheral flow problem.Consider that supersonic speed inviscid fluid flows to unit cylinder.Primary data is: level flows Mach 2 ship 3, density is 1, and pressure is 1, and region is [-10,10] × [-10,10], and left margin is for flowing border, and right margin is super Velocity of sound border.Fig. 5 reflection is the t=100 moment to process, by ST method, the density isogram that boundary condition obtains, and Fig. 6 reflects To be the t=100 moment process, by ST method, the Mach number isogram that boundary condition obtains, Fig. 7 reflection is to use in the t=100 moment GBCM method processes the density isogram that boundary condition obtains, and Fig. 8 reflection is to process by GBCM method in the t=100 moment The Mach number isogram that boundary condition obtains.

Above in conjunction with accompanying drawing, embodiments of the present invention are explained in detail, but the present invention are not limited to above-mentioned embodiment, In the ken that those of ordinary skill in the art are possessed, it is also possible to make various on the premise of without departing from present inventive concept Change.The above, be only presently preferred embodiments of the present invention, and the present invention not makees any pro forma restriction, although The present invention is disclosed above with preferred embodiment, but is not limited to the present invention, any those skilled in the art, In the range of without departing from technical solution of the present invention, when the technology contents of available the disclosure above makes a little change or is modified to equivalent The Equivalent embodiments of change, as long as being without departing from technical solution of the present invention content, according to the technical spirit of the present invention, in the present invention Spirit and principle within, any simple amendment that above example is made, equivalent and improvement etc., all still fall within this Within the protection domain of inventive technique scheme.

Claims (6)

1. the whole flow field method for numerical simulation of a limited bulk weighted essential non-oscillatory schemes, it is characterised in that: at 2-d polynomial weight During structure, according to the actual change under two-dimensional case of the conservation variable in flow field and then determine smoothing factor, concrete steps bag Include:

Step 1. is at the boundary point of object elementWithFunction u is reconstructed by place, wherein, and x, y For cartesian coordinate, i, j are respectively the mesh point sequence number on x, y direction；

1.1 use least square method strategies obtain on each Gauss point four times reconstruct multinomial Q (x, y)；

Large form T is divided into a plurality of little template by 1.2, in each little template, and reconstruction of function u at each Gauss point respectively Value；

Step 2. calculates half Discrete Finite volume format；

Step 3. uses three rank TVD Runge-Kutta time discrete formula to obtain space-time full discrete scheme.

The whole flow field method for numerical simulation of a kind of limited bulk weighted essential non-oscillatory schemes the most as claimed in claim 1, its feature exists In: non-Body fit to due to cartesian grid, the cut cells produced at object plane, use ST and GBCM two kinds virtual Element method processes boundary condition, and concrete steps include:

Step one, at object plane boundary, determine virtual grid unit I₁Central point A is about the symmetric points B of object plane；

Step 2, determine the virtual grid unit I that B point falls into₂, and virtual grid unit I around₃,I₄,I₅；

Step 3, utilize virtual grid unit I₂,I₃,I₄,I₅, by interpolation formula, determine B point flow field value；

Step 4, determine dummy unit I by Solid boundary condition and momentum relationship₁With the relation of B point, thus obtain virtual list Unit I₁Flow field value.

The whole flow field method for numerical simulation of a kind of limited bulk weighted essential non-oscillatory schemes the most as claimed in claim 1, its feature exists In: in step 1.1, (x, y) is two-Dimensional Algebraic multinomial to described four reconstruct multinomial Q, and (x, y) and function u is at target list for Q Unit I₁₃Mean value equal:

$\frac{1}{\left|I_{13}\right|}\underset{I_{13}}{\∫}Q(x,y)dxdy={\stackrel{\‾}{u}}_{13}$

Wherein,Represent object element I₁₃Cell-average value, | I₁₃| represent object element I₁₃Area；

Further, (x y) meets the mean value of other unit superior function u of large form T under least square meaning to Q.

The whole flow field method for numerical simulation of a kind of limited bulk weighted essential non-oscillatory schemes the most as claimed in claim 1, its feature exists In: in step 1.2, the multinomial that each little template is constructed is: p_l(x, y), (l=1 ..., n), wherein p_l(x,y)∈span{1,x,y,x²,xy,y²,x²y,xy²,x²y², i.e. belonging to the linear space that these substrates are opened, subscript l represents right The little template sequence number answered, n is little template number；

Multinomial p_l(x, y) must be with u in same little template S_lThe mean value of upper unit is equal；

At object element I₁₃Gaussian quadrature nodeOn function u is reconstructed, and solve at Gauss pointLinearly weigh γ_l；

Subsequently, smoothing factor β is calculated_l, smoothing factor β_lIt is used for weighing multinomial p_l(x, y), (l=1 ..., n) at object element I₁₃On Smooth degree, its computing formula is:

${\mathrm{\β}}_l=\underset{\left|\mathrm{\α}\right|=1}{\overset{r}{\Σ}}\underset{I_{13}}{\∫}|I_{13}{|}^{\left|\mathrm{\α}\right|-1}{\left(D^{\mathrm{\α}}{p}_l\right(x,y\left)\right)}^{2}dxdy$

Wherein, α is multiple indexes, and D is derivation operator, and r is multinomial p_l(x, number of times y)；

Calculate nonlinear weight on this basis:

${\mathrm{\ω}}_l=\frac{{\stackrel{\‾}{\mathrm{\ω}}}_l}{\underset{k=1}{\overset{9}{\Σ}}{\stackrel{\‾}{\mathrm{\ω}}}_k},{\stackrel{\‾}{\mathrm{\ω}}}_l=\frac{{\mathrm{\γ}}_l}{{({10}^{-7}+{\mathrm{\β}}_l)}^2},(l=1,\mathrm{...},n)$

Wherein, γ_lIt is linear power, β_lIt is smoothing factor,For calculating process median, ω_lFor nonlinear weight.Finally, approximation Result is:

$u_{i+\frac{1}{2},j+\frac{1}{2\sqrt{3}}}\≈\underset{l=1}{\overset{n}{\Σ}}{\mathrm{\ω}}_l{p}_l(x_{i+\frac{1}{2}},y_{j+\frac{1}{2\sqrt{3}}}).$

Repeat above step, at object element I₁₃Other Gaussian quadrature nodes on function u is reconstructed.

The whole flow field method for numerical simulation of a kind of limited bulk weighted essential non-oscillatory schemes the most as claimed in claim 1, its feature exists In: in step 2, the calculating formula of half Discrete Finite volume format is:

$\frac{d}{dt}{\stackrel{\‾}{u}}_{i,j}\left(t\right)=L\left(u\right)=-\frac{1}{\mathrm{\Δ}x\mathrm{\Δ}y}{\∫}_{\∂I_{i,j}}(f,g)*\stackrel{\→}{n}ds$

Wherein, ${\stackrel{\‾}{u}}_{i,j}\left(t\right)=\frac{1}{\mathrm{\Δ}x\mathrm{\Δ}y}{\∫}_{I_{i,j}}u(x,y,t)dxdy,$ L is right-hand vector,For outer normal unit vector；

Line integral formula is rewritten as 2 Gauss quadrature formulas:

${\∫}_{\∂I_{i,j}}(f,g)*\stackrel{\→}{n}ds\≈|\∂I_{i,j}|\underset{l=1}{\overset{2}{\Σ}}{\mathrm{\ω}}_l\left(f\right(u\left(x_{i\±\frac{1}{2}},y_{j+G_l},t\right)),g(u\left(x_{i+G_l},y_{j\±\frac{1}{2}},t\right)\left)\right)*\stackrel{\→}{n}$

Wherein,AndReplaced by numerical flux.

The whole flow field method for numerical simulation of a kind of limited bulk weighted essential non-oscillatory schemes the most as claimed in claim 1, its feature exists In: in step 3, calculating formula is:

$\left\{\begin{array}{c}{u}^{\left(1\right)}=u^n+\mathrm{\Δ}tL\left(u^n\right)\\ u^{\left(2\right)}=\frac{3}{4}{u}^n+\frac{1}{4}{u}^{\left(1\right)}+\frac{1}{4}\mathrm{\Δ}tL\left(u^{\left(1\right)}\right)\\ u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}{u}^{\left(2\right)}+\frac{2}{3}\mathrm{\Δ}tL\left(u^{\left(2\right)}\right)\end{array}\right..$