CN104091065A - Intermittent flow numerical simulation method for solving shallow water problem - Google Patents

Abstract

The invention discloses an intermittent flow numerical simulation method for solving the shallow water problem and relates to the field of computational fluid mechanics, in particular to a numerical method for solving the two-dimensional hyperbolic system Riemannian problem. The method is applied to the solution of a shallow-water wave equation for numerical simulation of intermittent flow. According to the method, in a two-dimensional structureless triangular net, the Riemannian problem of a hyperbolic type shallow-water equation is solved, numerical precision and resolution are improved, and hydrodynamic force features of intermittent flow are discussed; a high-resolution numerical model applicable to a complex terrain and the phenomenon that a large gradient or intermittent flow is contained in a complex calculation area is established based on the structureless triangular net, errors can be reduced to a large extent, and the resolution of simulation of the hydrodynamic force features is improved.

Description

A kind of method that solves shallow water problems simulation interruption current numerical value

Technical field

The invention discloses a kind of method that shallow water problems simulation is interrupted current numerical value that solves, be specifically related to Fluid Mechanics Computation technical field.

Background technology

Numerical simulation is one of main path of understanding HYDRODYNAMIC CHARACTERISTICS.In coastal waters, the large area such as river mouth, lake, large reservoir opens wide waters, can adopt Wave of Two-Dimension Shallow Water Equation to describe the feature of hydrodynamic parameter.Because shallow water equation has nonlinear characteristic, generally cannot obtain analytic solution, mainly solve by numerical simulation.

Numerical method is the core of numerical simulation.River mouth, seashore and offshore waters, be subject to the impact of the factors such as off-lying sea tidal wave, water front and landform, hydrodynamic environment is very complicated, except these physical factors, also has the mathematical theory characteristic of shallow water equation itself, if numerical method selects improper meeting to cause false concussion, and the not all numerical method that is applicable to solving shallow water current can be used for simulation and is interrupted current.In general, even if given initial condition is very simple, the solution of shallow water equation also may develop into discontinuity wave or shock wave, and in fact, this is the general characteristic without the non-linear hyperbolic conservation law of dissipativeness.

The mobile numerical simulations of strong discontinuity such as dam break ripple, tidal bore be all the time important in hydrodynamics research, be also difficult problem.A lot of traditional numerical methods that solve shallow water equation are all the governing equations based on differential form, require physical descriptor everywhere continuous, can be micro-, therefore while there is be interrupted mobile, just there will be some problems at this class water level of simulation, flow velocity.Some general hydrodynamics numerical simulation softwares, are applied to general river course all very successful both at home and abroad, flow but can't well simulate discontinuous flow, particularly strong discontinuity.The main difficulty of discontinuous flow numerical simulation is to set up the computation scheme of " harmony ", shallow water equation flux term will with source item " harmony "; Before and after being interrupted, water level, current gradient are very large, require computation scheme can simulate large gradient and flow, and form is stable again.Thereby the theoretical research of discontinuous solution and numerical simulation just have the academic significance of particular importance.

Theoretical analysis shows, utilizes high-order, high precision numeric format not only can reduce the requirement to grid, can also accurately catch flow phenomenon complicated in flow field, and therefore developing high order accurate numerical method is one of important channel of improving numerical simulation calculation precision.In addition, have the numerical simulation of the physical phenomenon of multiple dimensional variation for room and time, require multiple yardstick information to be caught in numerical evaluation, this also needs high-accuracy high-resolution form.Therefore in shallow water current, the high resolution numerical simulation of discontinuous flow has more important learning value.

Correct simulation is interrupted current and first needs correct simulation Rankine-Hugoniot.In the time not considering viscous effect, discontinuity wave is a discontinuity surface.Current is discontinuity surface to be included in zoning and to be calculated by unified method, be i.e. so-called shock capturing method compared with effective method.Its basic point of departure is to use the conservativeness difference scheme compatible with conservation law differential equation group, the Difference Solution of gained meets discontinuity condition automatically in interruption both sides, no matter thereby whether there is interruption in solution, can not unify to calculate with making any distinction between, but require computing grid closeer, and require computation scheme to there is the ability of the large gradient of simulation.Single order precision format in common trapping, what can make that shock wave smears is excessively flat, and second order or the classical form numerical solution higher than second order accuracy easily produce non-physics concussion near shock wave.In recent years for the resolution and the false concussion that improve near numerical solution shock wave have made great progress, some are produced and have solved the high resolving power numeric format of shallow water discontinuity wave.

Sharp keen and the verisimilitude of so-called " high resolving power " exponential quantity solution figure, shows that the numerical dissipation of numerical method is faint and moderate.To ignore the conservation form shallow water equation U of source item _t+ F (U) _x=0 is example, and wherein U is conservativeness physical quantity, and F (U) is flux term.If the discontinuous point x of true solution _sin interval in, wherein subscript i is grid node numbering, as x < x _stime U=U _l, as x > x _stime U=U _r, subscripting L, the symbol U of R _l, U _rrepresent that respectively U is at x _sthe value on the some left side and the right.

The Difference Solution of conservation scheme is:

$U_j=\left\{\begin{array}{cc}{U}_L& (j<i)\\ U_M& (j=i)\\ U_R& (j>i)\end{array}\right.;$

Wherein, subscript j is grid node numbering, U _mfor x _japproximate U value on point.

In order to meet conservativeness, need following formula to set up:

$U_M(x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}})=U_L(x_s-x_{i-\frac{1}{2}})+U_R(x_{i+\frac{1}{2}}-x_s).$

This formula and difference scheme have been determined x uniquely _sand U _m.Therefore, numerical value is interrupted and must at least relates to three lattice points (i.e. two intervals).The interruption resolution of some good forms has reached two lattice to three lattice at present.Realizing the design of high resolution scheme, is exactly under the prerequisite of the compatibility that ensures form, stability, in the time of design and structure form, can weaken as far as possible its numerical dispersion, keeps numerical dissipation sexual clorminance.

It is the basis of numerical simulation that grid generates, and its essence is that physics solves territory and calculates the conversion that solves territory.In the time solving the flow field with complex geometric shapes, it is a very important problem that suitable grid generates, and the quality of mesh quality directly has influence on convergence situation and the precision of result of calculation.Non-structure grid has that complex region adaptability is good, local refinement flexibly and be convenient to adaptive advantage, simulating nature border and complicated underwater topography, improve boundary condition simulation precision well.

Summary of the invention

Technical matters to be solved by this invention is: for the defect of prior art, a kind of method that shallow water problems simulation is interrupted current numerical value that solves is provided, specifically a kind of numerical method that solves Wave of Two-Dimension Shallow Water Equation Riemann Solution, can be applied to simulation dam break ripple, tidal bore etc. and can not press interruption current, and then understand the HYDRODYNAMIC CHARACTERISTICS of being interrupted current.

The present invention is for solving the problems of the technologies described above by the following technical solutions:

Solve the method that shallow water problems simulation is interrupted current numerical value, set up the high resolving power numerical model that simulation is interrupted current problem, utilize in described modeling landform and computational fields with the current phenomenon of large Gradient Solution or interruption, concrete steps are as follows:

Step 1, set up computing grid, non-structure triangular mesh is carried out in physical computing region and cut apart, obtain a plurality of triangular mesh unit, using single grid cell as controlling unit, wherein i triangle control unit is designated as Ω _i;

Step 2, reflect by conservation form physical quantity and under Descartes's rectangular coordinate system, by the Representation Equation of planar motion of shallow conservation form be the HYDRODYNAMIC CHARACTERISTICS of each computing grid:

U _t+F(U) _x+G(U) _y＝S (1)；

Wherein, U is conserved quantity vector; F, G is flux vector; S is source item vector, and t is the time, and x and y are cartesian coordinate system coordinate, U _t, F (U) _x, G (U) _ybe expressed as time t and space x, the partial derivative of y direction, is specifically expressed as:

U＝[h,uh,vh] ^T；

F＝[uh,u ²h+gh ²/2,uvh] ^T；

G＝[vh,uvh,v ²h+gh ²/2] ^T；

S＝[0,-ghb _x,-ghb _y] ^T；

Wherein, T represents vectorial transposition; U, v be respectively x, y direction along the average velocity component of the depth of water, and be about x, y, the function of t; H is the depth of water; G is acceleration of gravity; b _x, b _ybe respectively the component of bed sloped in x, y direction, b (x, y) is bed level of the river function;

Solve the Riemann Solution of above-mentioned equation, to system of equations (1), make vector function E=(F, G), for gradient operator, formula (1) is expressed as:

U _t+▽·E＝S (2)；

Adopt the triangular mesh unit of dividing in step 1 to carry out zoning discrete, at Ω _iabove system of equations (2) is carried out to integration, and by utilizing counterclockwise Green formula that area is divided into line integral:

$\frac{d}{\mathrm{dt}}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Ud\&Omega;}=-\underset{m=1}{\overset{3}{\mathrm{\&Sigma;}}}{\&Integral;}_{L_m}E\&CenterDot;n_m\mathrm{dl}+{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Sd\&Omega;}---\left(3\right);$

Wherein, for the ordinary differential operator of time orientation, d Ω is that area divides infinitesimal, and dl is line integral infinitesimal, L _mrepresent i triangular element Ω _im article of limit, m=1,2,3, n _m=(n _mx, n _my)=(cos θ, sin θ) be L _munit outside normal vector, n _mx, n _myfor n _mat the component of x, y direction, θ is the angle of outer normal direction and horizontal direction;

Definition unit mean value: ${\stackrel{\&OverBar;}{U}}_i=\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Ud\&Omega;},{\stackrel{\&OverBar;}{S}}_i=\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Sd\&Omega;},$ Wherein, | Ω _i| be grid cell area;

Utilize 2 Gauss quadrature formulas to solve the line integral in formula (3):

Wherein, | L _m| be triangle edges boundary line L _mlength, G _qfor triangular element limit L _mon Gauss point, for corresponding weight coefficient, q is Gauss point numbering,

Construct a numerical flux form approximate substitution E (U (G _q, t)), variate-value inside and outside on triangular mesh elementary boundary is designated as respectively to U _l, U _r, will be by Gauss point G _qflux be designated as wherein, U _l(G _q, t), U _r(G _q, t) be respectively U in the inside and outside value in Gauss point triangular mesh unit, draw thus the limited bulk semidiscrete formulation of governing equation (3):

$\frac{d}{\mathrm{dt}}{\stackrel{\&OverBar;}{U}}_i=-\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}\underset{m=1}{\overset{3}{\mathrm{\&Sigma;}}}\left|L_m\right|\underset{q=1}{\overset{2}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_q\hat{E}\left(U_R(G_q,t)U_L(G_q,t)\right)\&CenterDot;n_m+{\stackrel{\&OverBar;}{S}}_i---\left(5\right);$

Step 3, solve the HLL numerical flux on triangular mesh elementary boundary

Step 4, employing three rank Runge-Kutta time discrete methods, to discrete processing the on time orientation, draw the time discrete form of third-order;

Step 5, employing three rank WENO form reconstructing methods solve U _l(G _q, t), U _r(G _q, t);

Step 6, to each grid cell, utilize the conserved quantity vector integral mean of each unit on current time layer, solve the conserved quantity vector integral mean in next moment:

Through on step 5 convection cell equation time orientation discrete after, then convection cell equation carries out spatial spreading, expression formula is:

$\frac{\mathrm{dU}}{\mathrm{dt}}=L\left(U\right);$

Wherein, L (U) is spatial discretization operator;

Set n the U value on time horizon for U ⁿ, after a time step Δ t, the approximate value U of next time horizon ⁿ⁺¹be expressed as:

U ⁽¹⁾＝U ⁿ+ΔtL(U ⁿ)；

$U^{\left(2\right)}=\frac{3}{4}{U}^n+\frac{1}{4}{U}^{\left(1\right)}+\frac{1}{4}\mathrm{\&Delta;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{\&Delta;tL}\left(U^{\left(2\right)}\right);$

Wherein, U ⁽¹⁾, U ⁽²⁾for middle transition vector.

As present invention further optimization scheme, in described step 3, the HLL flux expression formula of equation is specific as follows:

${\hat{E}}^{\mathrm{HLL}}(U_L,U_R)=\left\{\begin{array}{cc}E\left(U_L\right)& s_L\&GreaterEqual;0\\ \frac{s_{R}E\left(U_L\right)-s_{L}E\left(U_R\right)+s_R{s}_L(U_R-U_L)}{s_R-s_L}& s_L\&le;0\&le;s_R\\ E\left(U_R\right)& s_R\&le;0\end{array}\right.---\left(6\right);$

Wherein, s _l, s _rrespectively the upper infimum of wave velocity: s _l=u _l-a _lq _l, s _r=u _r+ a _rq _r;

Wherein, the velocity of wave of the right and left is expressed as:

Set q _k, h ^*for middle transition amount, K=L, R, is specially:

$q_K=\left\{\begin{array}{cc}\sqrt{\frac{1}{2}[\frac{(h_{*}+h_K)h_{*}}{h_K^2}}& h_{*}>h_K\\ 1& h_{*}\&le;h_K\end{array}\right.;$

$h^{*}=\frac{1}{g}{[\frac{1}{2}(a_L+a_R)+\frac{1}{4}(u_L-u_R)]}^2.$

As present invention further optimization scheme, described step 4 comprises: for triangular mesh unit structure triangular element Ω ₀polynomial expression, concrete steps comprise:

(501) structure of linear polynomial:

Based on triangular element Ω ₀and nine adjacent triangular elements, set nine templates, be specially:

B ₁＝{Ω ₀,Ω _j,Ω _k}；B ₂＝{Ω ₀,Ω _k,Ω _i}；B ₃＝{Ω ₀,Ω _i,Ω _j}；B ₄＝{Ω ₀,Ω _i,Ω _ia}；B ₅＝{Ω ₀,Ω _i,Ω _ib}；B ₆＝{Ω ₀,Ω _j,Ω _ja}；B ₇＝{Ω ₀,Ω _j,Ω _jb}；B ₈＝{Ω ₀,Ω _k,Ω _ka}；B ₉＝{Ω ₀,Ω _k,Ω _kb}；

Wherein, subscript i, j, k, ia, ib, ja, jb, ka, the numbering that kb is different adjacent cells;

Definition of T=B ₁∪ B ₂∪ ... ∪ B ₉for stencil sets, wherein ∪ is union symbol.In each little template, construct linear polynomial P (x, y), make mean value and known units Ω on three grid cells that P (x, y) covers in little template ₀mean value equate;

(502) the linear power of structure:

Be constructed as follows quadratic polynomial: Q (x, y)=b ₁+ b ₂x+b ₃y+b ₄xy+b ₅x ²+ b ₆y ²;

Wherein, b ₁, b ₂, b ₃, b ₄, b ₅, b ₆for polynomial undetermined coefficient, at { Ω ₀, Ω _i, Ω _j, Ω _k, Ω _ia, Ω _ib, Ω _ja, Ω _jb, Ω _ka, Ω _kbmeet on these ten unit: s=0, i, j, k, ia, ib, ja, jb, ka, kb; Wherein represent unit Ω _sthe mean value of upper U.

Diabolo unit Ω ₀each Gauss point (x ^g, y ^g), set linear power for γ _l(l=1 ..., 9), γ _lbe the constant by the unique decision of geometric properties of triangular mesh, weight coefficient is non-negative;

By linear power γ for linear polynomial P (x, y) _lmerge:

$R(x,y)=\underset{l=1}{\overset{9}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_l{P}_l(x,y)---\left(7\right)$

γ _lmeet following equation: R (x ^g, y ^g)=Q (x ^g, y ^g) (8)

(503) structure nonlinear weight and the smooth factor:

On the basis of linearity power, construct nonlinear weight, based on consistance, weight coefficient and be 1,

Nine linear polynomials are divided into four groups: $\underset{l=1}{\overset{9}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_l{P}_l(x,y)=\underset{\tilde{l}=1}{\overset{4}{\mathrm{\&Sigma;}}}{\widetilde{\mathrm{\&gamma;}}}_{\tilde{l}}{\tilde{P}}_{\tilde{l}}(x,y);$

First Gauss point on i limit:

First group contains equation P ₂(0, k, i), P ₄(0, i, ia), P ₅(0, i, ib),

${\tilde{P}}_1=({\mathrm{\&gamma;}}_2{P}_2+{\mathrm{\&gamma;}}_4{P}_4+{\mathrm{\&gamma;}}_5{P}_5)/({\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_4+{\mathrm{\&gamma;}}_5),{\widetilde{\mathrm{\&gamma;}}}_1={\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_4+{\mathrm{\&gamma;}}_5;$

Second group contains equation P ₃(0, i, j), P ₆(0, j, ja), P ₇(0, j, jb),

${\tilde{P}}_2=({\mathrm{\&gamma;}}_3{P}_3+{\mathrm{\&gamma;}}_6{P}_6+{\mathrm{\&gamma;}}_7{P}_7)/({\mathrm{\&gamma;}}_3+{\mathrm{\&gamma;}}_6+{\mathrm{\&gamma;}}_7),{\widetilde{\mathrm{\&gamma;}}}_2={\mathrm{\&gamma;}}_3+{\mathrm{\&gamma;}}_6+{\mathrm{\&gamma;}}_7;$

The 3rd group contains equation P ₁(0, j, k), P ₈(0, k, ka), P ₉(0, k, kb),

${\tilde{P}}_3=({\mathrm{\&gamma;}}_1{P}_1+{\mathrm{\&gamma;}}_8{P}_8+{\mathrm{\&gamma;}}_9{P}_9)/({\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_8+{\mathrm{\&gamma;}}_9),{\widetilde{\mathrm{\&gamma;}}}_3={\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_8+{\mathrm{\&gamma;}}_9;$

The 4th group contains equation P ₁(0, j, k), P ₂(0, k, i), P ₃(0, i, j),

${\tilde{P}}_4=({\mathrm{\&gamma;}}_1{P}_1+{\mathrm{\&gamma;}}_1{P}_2+{\mathrm{\&gamma;}}_3{P}_3)/({\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_3),{\widetilde{\mathrm{\&gamma;}}}_3={\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_3;$

Draw linear polynomial: for polynomial non-negative solution;

Order: ${\tilde{P}}_1(x^{G_1},y^{G_1})=a_1\stackrel{\&OverBar;}{u_0}+a_2\stackrel{\&OverBar;}{u_i}+a_3\stackrel{\&OverBar;}{u_k}+a_4\stackrel{\&OverBar;}{u_{\mathrm{ia}}}+a_5\stackrel{\&OverBar;}{u_{\mathrm{ib}}},$ To a polynomial of degree n p (x, y), define smooth measurement: $\mathrm{\&Psi;}=\underset{1\&le;\left|\mathrm{\&alpha;}\right|\&le;n}{\mathrm{\&Sigma;}}{\&Integral;}_{\mathrm{\&Omega;}}{\left|\mathrm{\&Omega;}\right|}^{\left|\mathrm{\&alpha;}\right|-1}{\left(D^{\mathrm{\&alpha;}}p(x,y)\right)}^2\mathrm{dxdy};$

Wherein, α is a complex exponential, and D is partial derivative, and nonlinear weight is:

${\mathrm{\&omega;}}_{\tilde{l}}=\frac{{\widetilde{\mathrm{\&gamma;}}}_{\tilde{l}}}{{(\mathrm{\&epsiv;}+{\mathrm{\&Psi;}}_{\tilde{l}})}^2}$

Wherein, individual polynomial expression smoothness, ε is little positive number.

As present invention further optimization scheme, the span of described ε is (10 ^-6, 10 ^-2).

As present invention further optimization scheme, the value of described ε is ε=10 ^-6.

As present invention further optimization scheme, in described step 6, the concrete account form of time step Δ t is:

$\mathrm{\&Delta;t}=\mathrm{CFL}\&CenterDot;\underset{i}{\min }\left(\frac{\underset{m}{\min }\left|L_m\right|}{\sqrt{g{h}_i}+\sqrt{u_i^2+v_i^2}}\right);$

Wherein, CFL is CFL (Courant-Friedrichs-Lewy) conditional number.

The present invention adopts above technical scheme compared with prior art, there is following technique effect: the present invention is based on non-structure triangular mesh and set up and can adapt in complex-terrain and complicated calculations region with large gradient or be interrupted the high resolving power numerical model of current phenomenon, the method can farthest reduce error, improve the resolution of Simulated Water dynamic characteristic, capturing discontinuous ripple more accurately.

Brief description of the drawings

Fig. 1 is triangular mesh control unit and outer normal direction schematic diagram thereof.

Fig. 2 is WENO form typical case network of triangle grid template.

Fig. 3 is the HLL numerical flux Approximate Riemann Solution of shallow water equation.

Fig. 4 is the geometrical plane figure of local dam break in one embodiment of the invention.

Fig. 5 calculates water level three-dimensional distribution map while one embodiment of the invention intermediate cam shape mesh generation and time being 7.2 seconds.

Fig. 6 calculates the depth of water elevation isoline schematic diagram when time is 7.2 seconds in one embodiment of the invention.

Fig. 7 calculates water flow field distribution plan when the time is 7.2 seconds in one embodiment of the invention.

Embodiment

Below in conjunction with accompanying drawing, technical scheme of the present invention is described in further detail:

Under Descartes's rectangular coordinate system, the equation of the conservation form of planar motion of shallow is typically expressed as:

U _t+F(U) _x+G(U) _y＝S (1)

Wherein, U is conserved quantity vector; F, G is flux vector; S is source item vector, is specifically expressed as:

U＝[h,uh,vh] ^T，

F＝[uh,u ²h+gh ²/2,uvh] ^T，

G＝[vh,uvh,v ²h+gh ²/2] ^T，

S＝[0,-ghb _x,-ghb _y] ^T，

Wherein, subscript T represents vectorial transposition, and t is the time; U (x, y, t), v (x, y, t) be respectively x, y direction along the average velocity component of the depth of water, and be about x, y, the function of t; H (x, y) is the depth of water; G is acceleration of gravity; b _x, b _ybe respectively the component of bed sloped in x, y direction, b (x, y) is bed level of the river function, and subscript x, y are illustrated respectively in the partial derivative of x, y direction.

Wish simulation shallow water is interrupted current, will solve the Riemann Solution of above-mentioned shallow water equation.In the structure of numerical method, to hyperbolic conservation equation group (1), make vector function E=(F, G), for gradient operator, formula (1) can be expressed as:

U _t+▽·E＝S (2)

For the convenience on realizing and calculating, it is discrete that the present invention adopts unified triangular mesh unit to carry out zoning, and according to counterclockwise record cell information, using single grid cell as controlling unit, physical descriptor is configured in the center of each unit.I triangle control unit is designated as to Ω _i, at Ω _ithe upper system of equations (2) to vector mode is carried out integration, and by utilizing counterclockwise Green formula that area is divided into line integral:

$\frac{d}{\mathrm{dt}}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Ud\&Omega;}=-\underset{m=1}{\overset{3}{\mathrm{\&Sigma;}}}{\&Integral;}_{L_m}E\&CenterDot;n_m\mathrm{dl}+{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Sd\&Omega;}---\left(3\right)$

Wherein, for the ordinary differential operator of time orientation, d Ω is that area divides infinitesimal, and dl is line integral infinitesimal, L _mrepresent triangular element Ω _im article of limit, m=1,2,3, vector n _m=(n _mx, n _my)=(cos θ, sin θ) be L _munit outside normal vector, n _mx, n _myfor n _mat the component of x, y direction, θ is the angle of outer normal direction and horizontal direction, as shown in Figure 1.

Utilize the thought definition unit mean value of following integral mean:

${\stackrel{\&OverBar;}{U}}_i=\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Ud\&Omega;},{\stackrel{\&OverBar;}{S}}_i=\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Sd\&Omega;}$

Wherein, | Ω _i| be grid cell area.

To the line integral in formula (3), the present invention utilizes 2 Gauss quadrature formulas, can obtain:

Wherein, | L _m| be triangle edges boundary line L _mlength, G _qfor triangular element limit L _mon Gauss point, corresponding weight coefficient is

What finite volume method obtained is the mean value on grid cell by can reconstruct high-order moment R (x, y), be similar to this value that U is ordered at Gauss.Because can not ensure that R (x, y) is continuous in net boundary, so can not, directly to formula (4) integration, therefore need to construct a numerical flux form here and be similar to E (U (G _q, t)).

Suppose that considered a problem time step is enough little, the information of remaining element not yet propagates into investigated boundary surface, according to the thought of being interrupted, physical descriptor value is considered to different on elementary boundary both sides, and variate-value inside and outside on triangle border is designated as respectively U _l, U _r, the outside of triangular element is also just equivalent to adjacent triangular element inside.In this case, the calculating of numerical flux also must be considered this interruption.Note for by Gauss point G _qflux, as shown in Figure 2, U _l(G _q, t), U _r(G _q, t) be respectively U in the inside and outside value of Gauss point triangle.

Limited bulk semidiscrete formulation that so just can controlled equation (3):

$\frac{d}{\mathrm{dt}}{\stackrel{\&OverBar;}{U}}_i=-\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}\underset{m=1}{\overset{3}{\mathrm{\&Sigma;}}}\left|L_m\right|\underset{q=1}{\overset{2}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_q\hat{E}\left(U_R(G_q,q)U_L(G_q,q)\right)\&CenterDot;n_m+{\stackrel{\&OverBar;}{S}}_i---\left(5\right)$

After so discrete, the numerical solution problem of partial differential equations is just converted into the Solve problems of an ordinary differential equation group.So there are three major issues to need to solve in ensuing departure process: the one, numerical flux form; The 2nd, U _r(G _q, t), U _l(G _q, computing method (structure of weighing again) t); The 3rd, time discrete mode.

What model of the present invention adopted respectively is HLL numerical flux, three rank WENO form reconstructing methods and three step Runge-kutta time discrete forms, next will introduce respectively.

For the form of numerical flux, the present invention adopts HLL flux, and HLL flux is the numerical flux based on Approximate Riemann Solution, and it is divided into three states with two kinds of ripples by solution space, as shown in Figure 3.The HLL flux expression formula of shallow water equation is as follows:

${\hat{E}}^{\mathrm{HLL}}(U_L,U_R)=\left\{\begin{array}{cc}E\left(U_L\right)& s_L\&GreaterEqual;0\\ \frac{s_{R}E\left(U_L\right)-s_{L}E\left(U_R\right)+s_R{s}_L(U_R-U_L)}{s_R-s_L}& s_L\&le;0\&le;s_R\\ E\left(U_R\right)& s_R\&le;0\end{array}\right.---\left(6\right)$

Wherein, s _l, s _rbe respectively the upper infimum of wave velocity, expression is:

s _L＝u _L-a _Lq _L，s _R＝u _R+a _Rq _R

Wherein, the velocity of wave of the right and left is expressed as q _k, (K=L, R), h ^*for middle transition amount, be specially:

$q_K=\left\{\begin{array}{cc}\sqrt{\frac{1}{2}[\frac{(h_{*}+h_K)h_{*}}{h_K^2}}& h_{*}>h_K\\ 1& h_{*}\&le;h_K\end{array}\right.$

$h^{*}=\frac{1}{g}{[\frac{1}{2}(a_L+a_R)+\frac{1}{4}(u_L-u_R)]}^2$

The Runge-Kutta discrete scheme on discrete employing three rank of time orientation, after considering that convection cell equation carries out spatial spreading, become shape as: ordinary differential equation group, wherein L (U) is spatial discretization operator.The object of time discrete form is the equation discrete by time horizon solution room, supposes n the U value U on time horizon ⁿknown, after a time step Δ t, the approximate value U of next time horizon ⁿ⁺¹can try to achieve as follows:

U ⁽¹⁾＝U ⁿ+ΔtL(U ⁿ)

$U^{\left(2\right)}=\frac{3}{4}{U}^n+\frac{1}{4}{U}^{\left(1\right)}+\frac{1}{4}\mathrm{\&Delta;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{\&Delta;tL}\left(U^{\left(2\right)}\right)$

Wherein, U ⁽¹⁾, U ⁽²⁾for middle transition vector.

Three step Runge-Kutta time discrete methods are to ensure TVD (Total Variation Diminishing) character, therefore the problem of (as shock wave or contact discontinuity) suppresses numerical value concussion to be suitable for having strong discontinuity, and combining for discontinuous solution and stable spatial spreading is nonlinear stability.The present invention selects the Runge-Kutta time discrete method of three steps to obtain the time discrete form of third-order.

Next will be for following net template structure triangular element Ω ₀polynomial expression, to reaching space third-order.

Introduce the polynomial method of structure below:

(1) structure of linear polynomial

First, based on triangular element Ω ₀and nine adjacent triangular elements, as shown in Figure 2, wherein subscript is just described different adjacent cells numberings, gives the coldest days of the year end individual little template, as follows:

B ₁＝{Ω ₀,Ω _j,Ω _k}，B ₂＝{Ω ₀,Ω _k,Ω _i}，B ₃＝{Ω ₀,Ω _i,Ω _j}，B ₄＝{Ω ₀,Ω _i,Ω _ia}，B ₅＝{Ω ₀,Ω _i,Ω _ib}，B ₆＝{Ω ₀,Ω _j,Ω _ja}，B ₇＝{Ω ₀,Ω _j,Ω _jb}，B ₈＝{Ω ₀,Ω _k,Ω _ka}，B ₉＝{Ω ₀,Ω _k,Ω _kb}。Definition of T=B ₁∪ B ₂∪ ... ∪ B ₉for stencil sets, wherein ∪ is union symbol.

Next, in each little template, construct linear polynomial, below with template B ₁for example is introduced construction process, the linear polynomial P of structure ₁(x, y) meets:

$\stackrel{\&OverBar;}{P_1{(x,y)}_{{\mathrm{\&Omega;}}_r}}={\stackrel{\&OverBar;}{U}}_r,r=0,j,k$

Be P ₁mean value on three grid cells that (x, y) covers in little template equate with known cell-average value.

Linear equation in two unknowns contains three unknown quantitys, above three conditions can uniquely determine a linear polynomial.In like manner, in each little template, can construct a linear polynomial.

(2) the linear power of structure

Use linear polynomial reconstruct three rank WENO forms below.Have third-order in order to reach result of calculation, the present invention must construct the quadratic polynomial that U is satisfied and be similar to three rank that reach variable on each triangular element.The present invention is constructed as follows quadratic polynomial:

Q(x,y)＝b ₁+b ₂x+b ₃y+b ₄xy+b ₅x ²+b ₆y ²

Wherein, b ₁, b ₂, b ₃, b ₄, b ₅, b ₆for polynomial undetermined coefficient, at { Ω ₀, Ω _i, Ω _j, Ω _k, Ω _ia, Ω _ib, Ω _ja, Ω _jb, Ω _ka, Ω _kbmeet on these ten unit:

${\stackrel{\&OverBar;}{Q(x,y)}}_{\mathrm{\&Omega;s}}={\stackrel{\&OverBar;}{U}}_s,$ s＝0,i,j,k,ia,ib,ja,jb,ka,kb

Because there are six unknowm coefficients, ten equations, so system of equations is overdetermination, can solve by least square method.Because for same grid cell, some adjacent cells status should be identical, as Ω _iaand Ω _ibequal; Ω _jaand Ω _jbequal; Ω _kaand Ω _kbbe equal, this is not on affecting with least square method structure quadratic equation.

Next determine linear power γ _l(l=1 ..., 9), reconstruct linear polynomial.In order to meet stability and compatibility, claim coefficient is non-negative, and and is 1 equally.

Diabolo unit Ω ₀each Gauss point (x ^g, y ^g), the present invention will find one group of linear power γ _l, they are some constants by the unique decision of geometric properties of grid.The linear polynomial R (x, y) of reconstruct merges linear polynomial to get by these power.

$R(x,y)=\underset{l=1}{\overset{9}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_l{P}_l(x,y)---\left(7\right)$

Weight function γ _l(l=1 ..., 9) selection to make following equation meet:

$R(x^G,y^G)=Q(x^G,y^G)---\left(8\right)$

Each Gauss point difference, the weight coefficient obtaining is also different.In the time that Gauss point is fixing, formula (8) left end and right-hand member are all $\{\stackrel{\&OverBar;}{U_0},{\stackrel{\&OverBar;}{U}}_i,\stackrel{\&OverBar;}{U_j},\stackrel{\&OverBar;}{U_k},\stackrel{\&OverBar;}{U_{\mathrm{ia}}},\stackrel{\&OverBar;}{U_{\mathrm{ja}}},\stackrel{\&OverBar;}{U_{\mathrm{ib}}},\stackrel{\&OverBar;}{U_{\mathrm{jb}}},\stackrel{\&OverBar;}{U_{\mathrm{ka}}},\stackrel{\&OverBar;}{U_{\mathrm{kb}}}\}$ Linear combination because two ends, left and right equate, the coefficient of so each grid mean value equates, the present invention has just obtained 10 equations like this.9 unknown number γ _l(l=1 ..., 9), ten equations, system of equations is overdetermination on the surface, and in fact the order of the matrix of coefficients of system of equations is 8, and solution of equations space is one dimension like this, has Basic Solutions system, and contains a solution vector in a Basic Solutions system.

Obviously, the polynomial expression R (x, y) of reconstruct is linear like this, and can reach the precision of quadratic polynomial, therefore at Gauss point $G_q=(x^{G_q},y^{G_q})$ On $U_{L,R}(G_q,t)=R(x^{G_q},y^{G_q})$ Can be to exact solution U (G _q, the approximate third-order that reaches t).

(3) structure nonlinear weight and the smooth factor

On the basis of linearity power, construct nonlinear weight below.Upper surface construction be linear power, but when the problem that contains shock wave class in calculating, must employing nonlinear weight.Select the suitable nonlinear weight factor, make the interpolation polynomial of separating smooth domain structure contribute large to combination polynomial expression, be that weighting factor is larger, the interpolation polynomial of constructing on the region of containing interruption is got to minimum weighting factor, while making to approach the circulation on interface with this weighted interpolation polynomial expression, numerical solution is not shaken substantially.In order to ensure to stablize near shock wave, need to construct non-negative power, therefore need to be from constructing non-negative nonlinear weight.Because when structure linearity is weighed above, system of equations contains a free variable, solution space is one dimension, from going out to send the non-negative nonlinear weight of structure here.

In order to keep consistency, to need weight coefficient and be 1, next, nine linear polynomials are divided into four groups:

$\underset{l=1}{\overset{9}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_l{P}_l(x,y)=\underset{\tilde{l}=1}{\overset{4}{\mathrm{\&Sigma;}}}{\widetilde{\mathrm{\&gamma;}}}_{\tilde{l}}{\tilde{P}}_{\tilde{l}}(x,y)$

Each polynomial expression remain linear polynomial, time still second approximation in U.The present invention still needs rational allocation template to carry out corresponding four new polynomial expressions in normal situation, with in ensureing that working as concussion occurs, not every template is all shaken.

For first Gauss point on i limit, see Fig. 2:

First group contains equation P ₂(0, k, i), P ₄(0, i, ia), P ₅(0, i, ib):

${\tilde{P}}_1=({\mathrm{\&gamma;}}_2{P}_2+{\mathrm{\&gamma;}}_4{P}_4+{\mathrm{\&gamma;}}_5{P}_5)/({\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_4+{\mathrm{\&gamma;}}_5),{\widetilde{\mathrm{\&gamma;}}}_1={\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_4+{\mathrm{\&gamma;}}_5;$

Second group contains equation P ₃(0, i, j), P ₆(0, j, ja), P ₇(0, j, jb):

${\tilde{P}}_2=({\mathrm{\&gamma;}}_3{P}_3+{\mathrm{\&gamma;}}_6{P}_6+{\mathrm{\&gamma;}}_7{P}_7)/({\mathrm{\&gamma;}}_3+{\mathrm{\&gamma;}}_6+{\mathrm{\&gamma;}}_7),{\widetilde{\mathrm{\&gamma;}}}_2={\mathrm{\&gamma;}}_3+{\mathrm{\&gamma;}}_6+{\mathrm{\&gamma;}}_7;$

The 3rd group contains equation P ₁(0, j, k), P ₈(0, k, ka), P ₉(0, k, kb):

${\tilde{P}}_3=({\mathrm{\&gamma;}}_1{P}_1+{\mathrm{\&gamma;}}_8{P}_8+{\mathrm{\&gamma;}}_9{P}_9)/({\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_8+{\mathrm{\&gamma;}}_9),{\widetilde{\mathrm{\&gamma;}}}_3={\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_8+{\mathrm{\&gamma;}}_9;$

The 4th group contains equation P ₁(0, j, k), P ₂(0, k, i), P ₃(0, i, j):

${\tilde{P}}_4=({\mathrm{\&gamma;}}_1{P}_1+{\mathrm{\&gamma;}}_1{P}_2+{\mathrm{\&gamma;}}_3{P}_3)/({\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_3),{\widetilde{\mathrm{\&gamma;}}}_3={\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_3;$

The present invention obtains linear polynomial:

$\tilde{R}(x,y)=\underset{\tilde{l}=1}{\overset{4}{\mathrm{\&Sigma;}}}{\widetilde{\mathrm{\&gamma;}}}_{\tilde{l}}{\tilde{P}}_{\tilde{l}}(x,y);$

There is non-negative solution at above formula situation under, order:

${\tilde{P}}_1(x^{G_1},y^{G_1})=a_1\stackrel{\&OverBar;}{u_0}+a_2\stackrel{\&OverBar;}{u_i}+a_3\stackrel{\&OverBar;}{u_k}+a_4\stackrel{\&OverBar;}{u_{\mathrm{ia}}}+a_5\stackrel{\&OverBar;}{u_{\mathrm{ib}}};$

The smooth factor and nonlinear weight are with reference to the document of Jiang, and to a polynomial of degree n p (x, y), the present invention defines smooth measurement below:

$\mathrm{\&Psi;}=\underset{1\&le;\left|\mathrm{\&alpha;}\right|\&le;n}{\mathrm{\&Sigma;}}{\&Integral;}_{\mathrm{\&Omega;}}{\left|\mathrm{\&Omega;}\right|}^{\left|\mathrm{\&alpha;}\right|-1}{\left(D^{\mathrm{\&alpha;}}p(x,y)\right)}^2\mathrm{dxdy}$

Wherein, α is a complex exponential, and D is partial derivative, and nonlinear weight is defined as:

${\widetilde{\mathrm{\&omega;}}}_{\tilde{l}}=\frac{{\mathrm{\&omega;}}_{\tilde{l}}}{\underset{\tilde{l}}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_{\tilde{l}}},{\mathrm{\&omega;}}_{\tilde{l}}=\frac{{\widetilde{\mathrm{\&gamma;}}}_{\tilde{l}}}{{(\mathrm{\&epsiv;}+{\mathrm{\&Psi;}}_{\tilde{l}})}^2}$

Wherein, be $\tilde{R}(x,y)=\underset{\tilde{l}=0}{\overset{4}{\mathrm{\&Sigma;}}}{\widetilde{\mathrm{\&gamma;}}}_{\tilde{l}}{\tilde{P}}_{\tilde{l}}(x,y)$ In weight coefficient, individual polynomial expression smoothness; ε is little positive number, and object is to prevent that denominator from being zero, gets ε=10 in the present invention ^-6.ε value is (10 ^-6, 10 ^-2) logarithm value result do not have too large impact.

Use nonlinear weight replace linear power complete WENO reconstruct.

The method of introducing according to the present invention, provide below one complete, simulate local dam break and produce the example of dam break ripple by Wave of Two-Dimension Shallow Water Equation.

This example adopts the model of rectangle computational fields, and as shown in Figure 4, length and width is respectively 200 meters, and reservoir dam is positioned at x=100 rice, 10 meters of upstream, the initial time dam depth of waters, and 5 meters of dam downstream water depths, a certain moment dam bursts suddenly between y=95～170 meter.The dam break ripple forming is afterwards interrupted current exactly, still meets shallow water equation.

The present invention adopt the numerical method that solves shallow water equation Riemann Solution that the present invention provides simulate this part dam body collapse after the situation of current.

Calculating parameter is as follows, CFL=0.6; Non-structure triangle gridding subdivision is carried out in zoning, be split into altogether 17,658 triangular elements; All borders are adopted to wall boundary condition; Calculating and stopping the moment is 7.2 seconds.When Fig. 5, Fig. 6, Fig. 7 are followed successively by 7.2 seconds, triangular mesh subdivision and calculating water level three-dimensional distribution map, calculating depth of water elevation isoline schematic diagram, calculating water flow field distribution plan, particularly, calculate and adopt following steps:

(1) mesh generation.Non-structure triangular mesh subdivision is carried out in square zoning, be split into altogether M grid, and record the relevant information of each grid cell; Calculate the geometric properties of each grid cell;

(2) initialization.Set and calculate termination time T _end, set starting condition the present invention adopts cold start starting condition, and the depth of water is according to specified criteria, and in the time of x≤100, the depth of water is h=10, depth of water h=5 in the time of x≤100.Flow velocity is zero at initial time, i.e. u=v=0.

(4) computing time step delta t.For meeting stable condition, time step choose CFL (Courant-Friedrichs-Lewy) condition that is controlled by.For finite volume method or finite difference arbitrarily, CFL condition is a necessary condition that ensures that algorithm is stable and converge to the solution of partial differential equation.The present invention often adopts CFL condition number to be one to be greater than zero definite value that is less than 1, and the present invention gets CFL condition and counts CFL=0.6.The computing formula of time step can be expressed as so:

$\mathrm{\&Delta;t}=\mathrm{CFL}\&CenterDot;\underset{i}{\min }\left(\frac{\underset{m}{\min }\left|L_m\right|}{\sqrt{g{h}_i}+\sqrt{u_i^2+v_i^2}}\right),$ m＝1,2,3，i＝1,2,...,M；

Wherein, L _m(m=1,2,3) are three articles of length of sides of i triangular element.

(5) utilize the reconstruct of WENO form.Utilize and the information of nine adjacent triangular elements, the Gauss point value on three limits of each triangular element is carried out to the reconstruct of WENO form, calculate the value of each Gauss point left-right dots.

(6) calculate flux, renewal function value.Utilize HLL numerical flux expression formula to calculate flux term, calculate the variate-value of next time point according to three step Runge-Kutta time discrete formula

(7) make t=t+ Δ t, if t < is T _end, will value be assigned to turn back to (4), continue executive routine; If t < is T _end, export last result of calculation, program is carried out and is finished.

Fig. 6 has provided and has calculated the water level isoline that stops the moment.From all results, longitudinal and transverse propagation downstream when indentation, there dam break earial drainage, and have rarefaction wave upstream to propagate, dam break ripple to form harmony to bank high, because dam break breach is asymmetric, be subject to the cause that bank Yong Shui affects and flow velocity is larger, form two asymmetrical vortexs on both sides, dam site place.Wavefront is propagated with the form of discontinuity wave, and discontinuities has Yong Shui, realistic physical phenomenon.Result of calculation of the present invention will be with to have been reported result all more identical.Only, from figure, the result of the form based on various numerical fluxs is more or less the same.

By reference to the accompanying drawings embodiments of the present invention are explained in detail above, but the present invention is not limited to above-mentioned embodiment, in the ken possessing those of ordinary skill in the art, can also under the prerequisite that does not depart from aim of the present invention, makes a variety of changes.

Claims (6)

1. one kind solves the method that shallow water problems simulation is interrupted current numerical value, it is characterized in that, set up the high resolving power numerical model that simulation is interrupted current problems, utilize in described modeling landform and computational fields with the current phenomenon of large Gradient Solution or interruption, concrete steps are as follows:

Step 1, set up computing grid, non-structure triangular mesh is carried out in physical computing region and cut apart, obtain a plurality of triangular mesh unit, using single grid cell as controlling unit, wherein i triangle control unit is designated as Ω _i;

Step 2, reflect by conservation form physical quantity and under Descartes's rectangular coordinate system, by the Representation Equation of planar motion of shallow conservation form be the HYDRODYNAMIC CHARACTERISTICS of each computing grid:

U _t+F(U) _x+G(U) _y＝S (1)；

Wherein, U is conserved quantity vector; F, G is flux vector; S is source item vector, and t is the time, and x and y are cartesian coordinate system coordinate, U _t, F (U) _x, G (U) _ybe expressed as time t and space x, the partial derivative of y direction, is specifically expressed as:

U＝[h,uh,vh] ^T；

F＝[uh,u ²h+gh ²/2,uvh] ^T；

G＝[vh,uvh,v ²h+gh ²/2] ^T；

S＝[0,-ghb _x,-ghb _y] ^T；

Wherein, T represents vectorial transposition; U, v be respectively x, y direction along the average velocity component of the depth of water, and be about x, y, the function of t; H is the depth of water; G is acceleration of gravity; b _x, b _ybe respectively the component of bed sloped in x, y direction;

Solve the Riemann Solution of above-mentioned equation, to formula (1), make vector function E=(F, G), for gradient operator, formula (1) is expressed as:

U _t+▽·E＝S (2)；

Adopt the triangular mesh unit of dividing in step 1 to carry out zoning discrete, at Ω _iabove system of equations (2) is carried out to integration, and by utilizing counterclockwise Green formula that area is divided into line integral:

$\frac{d}{\mathrm{dt}}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Ud\&Omega;}=-\underset{m=1}{\overset{3}{\mathrm{\&Sigma;}}}{\&Integral;}_{L_m}E\&CenterDot;n_m\mathrm{dl}+{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Sd\&Omega;}---\left(3\right);$

Wherein, for the ordinary differential operator of time orientation, d Ω is that area divides infinitesimal, and dl is line integral infinitesimal, L _mrepresent i triangular element Ω _im article of limit, m=1,2,3, n _m=(n _mx, n _my)=(cos θ, sin θ) be L _munit outside normal vector, n _mx, n _myfor n _mat the component of x, y direction, θ is the angle of outer normal direction and horizontal direction;

Definition unit mean value: ${\stackrel{\&OverBar;}{U}}_i=\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Ud\&Omega;},{\stackrel{\&OverBar;}{S}}_i=\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}{\&Integral;}_{{\mathrm{\&Omega;}}_i}\mathrm{Sd\&Omega;},$ Wherein, | Ω _i| be grid cell area;

Utilize 2 Gauss quadrature formulas to solve the line integral in formula (3):

Wherein, | L _m| be triangle edges boundary line L _mlength, G _qfor triangular element limit L _mon Gauss point, for corresponding weight coefficient, q is Gauss point numbering,

Construct a numerical flux form approximate substitution E (U (G _q, t)), variate-value inside and outside on triangular mesh elementary boundary is designated as respectively to U _l, U _r, will be by Gauss point G _qflux be designated as wherein, U _l(G _q, t), U _r(G _q, t) be respectively U in the inside and outside value in Gauss point triangular mesh unit, draw thus the limited bulk semidiscrete formulation of governing equation (3):

$\frac{d}{\mathrm{dt}}{\stackrel{\&OverBar;}{U}}_i=-\frac{1}{\left|{\mathrm{\&Omega;}}_i\right|}\underset{m=1}{\overset{3}{\mathrm{\&Sigma;}}}\left|L_m\right|\underset{q=1}{\overset{2}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_q\hat{E}\left(U_R(G_q,t)U_L(G_q,t)\right)\&CenterDot;n_m+{\stackrel{\&OverBar;}{S}}_i---\left(5\right);$

Step 3, solve the HLL numerical flux on triangular mesh elementary boundary

Step 4, employing three rank Runge-Kutta time discrete methods, to discrete processing the on time orientation, draw the time discrete form of third-order;

Step 5, employing three rank WENO form reconstructing methods solve U _l(G _q, t), U _r(G _q, t);

Step 6, to each grid cell, utilize the conserved quantity vector integral mean of each unit on current time layer, solve the conserved quantity vector integral mean in next moment:

Through on step 5 convection cell equation time orientation discrete after, then convection cell equation carries out spatial spreading, expression formula is:

$\frac{\mathrm{dU}}{\mathrm{dt}}=L\left(U\right);$

Wherein, L (U) is spatial discretization operator;

Set n the U value on time horizon for U ⁿ, after a time step Δ t, the approximate value U of next time horizon ⁿ⁺¹be expressed as:

U ⁽¹⁾＝U ⁿ+ΔtL(U ⁿ)；

$U^{\left(2\right)}=\frac{3}{4}{U}^n+\frac{1}{4}{U}^{\left(1\right)}+\frac{1}{4}\mathrm{\&Delta;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{\&Delta;tL}\left(U^{\left(2\right)}\right);$

Wherein, U ⁽¹⁾, U ⁽²⁾for middle transition vector.

2. a kind of method that shallow water problems simulation is interrupted current numerical value that solves as claimed in claim 1, is characterized in that, in described step 3, the HLL flux expression formula of equation is specific as follows:

${\hat{E}}^{\mathrm{HLL}}(U_L,U_R)=\left\{\begin{array}{cc}E\left(U_L\right)& s_L\&GreaterEqual;0\\ \frac{s_{R}E\left(U_L\right)-s_{L}E\left(U_R\right)+s_R{s}_L(U_R-U_L)}{s_R-s_L}& s_L\&le;0\&le;s_R\\ E\left(U_R\right)& s_R\&le;0\end{array}\right.---\left(6\right);$

Wherein, s _l, s _rrespectively the upper infimum of wave velocity: s _l=u _l-a _lq _l, s _r=u _r+ a _rq _r;

Wherein, the velocity of wave of the right and left is expressed as:

Set q _k, h ^*for middle transition amount, K=L, R, is specially:

$q_K=\left\{\begin{array}{cc}\sqrt{\frac{1}{2}[\frac{(h_{*}+h_K)h_{*}}{h_K^2}}& h_{*}>h_K\\ 1& h_{*}\&le;h_K\end{array}\right.;$

$h^{*}=\frac{1}{g}{[\frac{1}{2}(a_L+a_R)+\frac{1}{4}(u_L-u_R)]}^2.$

3. a kind of method that shallow water problems simulation is interrupted current numerical value that solves as claimed in claim 1, is characterized in that, described step 4 comprises: for triangular mesh unit structure triangular element Ω ₀polynomial expression, concrete steps comprise:

(501) structure of linear polynomial:

Based on triangular element Ω ₀and nine adjacent triangular elements, set nine templates, be specially:

B ₁＝{Ω ₀,Ω _j,Ω _k}；B ₂＝{Ω ₀,Ω _k,Ω _i}；B ₃＝{Ω ₀,Ω _i,Ω _j}；B ₄＝{Ω ₀,Ω _i,Ω _ia}；B ₅＝{Ω ₀,Ω _i,Ω _ib}；B ₆＝{Ω ₀,Ω _j,Ω _ja}；B ₇＝{Ω ₀,Ω _j,Ω _jb}；B ₈＝{Ω ₀,Ω _k,Ω _ka}；B ₉＝{Ω ₀,Ω _k,Ω _kb}；

Wherein, subscript i, j, k, ia, ib, ja, jb, ka, the numbering that kb is different adjacent cells;

Definition of T=B ₁∪ B ₂∪ ... ∪ B ₉for stencil sets, wherein ∪ is union symbol, in each little template, constructs linear polynomial P (x, y), makes mean value and known units Ω on three grid cells that P (x, y) covers in little template ₀mean value equate;

(502) the linear power of structure:

Be constructed as follows quadratic polynomial: Q (x, y)=b ₁+ b ₂x+b ₃y+b ₄xy+b ₅x ²+ b ₆y ²;

Wherein, b ₁, b ₂, b ₃, b ₄, b ₅, b ₆for polynomial undetermined coefficient, at { Ω ₀, Ω _i, Ω _j, Ω _k, Ω _ia, Ω _ib, Ω _ja, Ω _jb, Ω _ka, Ω _kbmeet on these ten unit: s=0, i, j, k, ia, ib, ja, jb, ka, kb; Wherein represent unit Ω _sthe mean value of upper U,

Diabolo unit Ω ₀each Gauss point (x ^g, y ^g), set linear power for γ _l(l=1 ..., 9), γ _lbe the constant by the unique decision of geometric properties of triangular mesh, weight coefficient is non-negative;

By linear power γ for linear polynomial P (x, y) _lmerge:

$R(x,y)=\underset{l=1}{\overset{9}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_l{P}_l(x,y)---\left(7\right)$

γ _lmeet following equation: R (x ^g, y ^g)=Q (x ^g, y ^g) (8)

(503) structure nonlinear weight and the smooth factor:

On the basis of linearity power, construct nonlinear weight, based on consistance, weight coefficient and be 1,

Nine linear polynomials are divided into four groups: $\underset{l=1}{\overset{9}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_l{P}_l(x,y)=\underset{\tilde{l}=1}{\overset{4}{\mathrm{\&Sigma;}}}{\widetilde{\mathrm{\&gamma;}}}_{\tilde{l}}{\tilde{P}}_{\tilde{l}}(x,y);$

First Gauss point on i limit:

First group contains equation P ₂(0, k, i), P ₄(0, i, ia), P ₅(0, i, ib),

${\tilde{P}}_1=({\mathrm{\&gamma;}}_2{P}_2+{\mathrm{\&gamma;}}_4{P}_4+{\mathrm{\&gamma;}}_5{P}_5)/({\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_4+{\mathrm{\&gamma;}}_5),{\widetilde{\mathrm{\&gamma;}}}_1={\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_4+{\mathrm{\&gamma;}}_5;$

Second group contains equation P ₃(0, i, j), P ₆(0, j, ja), P ₇(0, j, jb),

${\tilde{P}}_2=({\mathrm{\&gamma;}}_3{P}_3+{\mathrm{\&gamma;}}_6{P}_6+{\mathrm{\&gamma;}}_7{P}_7)/({\mathrm{\&gamma;}}_3+{\mathrm{\&gamma;}}_6+{\mathrm{\&gamma;}}_7),{\widetilde{\mathrm{\&gamma;}}}_2={\mathrm{\&gamma;}}_3+{\mathrm{\&gamma;}}_6+{\mathrm{\&gamma;}}_7;$

The 3rd group contains equation P ₁(0, j, k), P ₈(0, k, ka), P ₉(0, k, kb),

${\tilde{P}}_3=({\mathrm{\&gamma;}}_1{P}_1+{\mathrm{\&gamma;}}_8{P}_8+{\mathrm{\&gamma;}}_9{P}_9)/({\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_8+{\mathrm{\&gamma;}}_9),{\widetilde{\mathrm{\&gamma;}}}_3={\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_8+{\mathrm{\&gamma;}}_9;$

The 4th group contains equation P ₁(0, j, k), P ₂(0, k, i), P ₃(0, i, j),

${\tilde{P}}_4=({\mathrm{\&gamma;}}_1{P}_1+{\mathrm{\&gamma;}}_1{P}_2+{\mathrm{\&gamma;}}_3{P}_3)/({\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_3),{\widetilde{\mathrm{\&gamma;}}}_3={\mathrm{\&gamma;}}_1+{\mathrm{\&gamma;}}_2+{\mathrm{\&gamma;}}_3;$

Draw linear polynomial: for polynomial non-negative solution;

Order: ${\tilde{P}}_1(x^{G_1},y^{G_1})=a_1\stackrel{\&OverBar;}{u_0}+a_2\stackrel{\&OverBar;}{u_i}+a_3\stackrel{\&OverBar;}{u_k}+a_4\stackrel{\&OverBar;}{u_{\mathrm{ia}}}+a_5\stackrel{\&OverBar;}{u_{\mathrm{ib}}},$ To a polynomial of degree n p (x, y), define smooth measurement: $\mathrm{\&Psi;}=\underset{1\&le;\left|\mathrm{\&alpha;}\right|\&le;n}{\mathrm{\&Sigma;}}{\&Integral;}_{\mathrm{\&Omega;}}{\left|\mathrm{\&Omega;}\right|}^{\left|\mathrm{\&alpha;}\right|-1}{\left(D^{\mathrm{\&alpha;}}p(x,y)\right)}^2\mathrm{dxdy};$

Wherein, α is a complex exponential, and D is partial derivative, and nonlinear weight is:

${\mathrm{\&omega;}}_{\tilde{l}}=\frac{{\widetilde{\mathrm{\&gamma;}}}_{\tilde{l}}}{{(\mathrm{\&epsiv;}+{\mathrm{\&Psi;}}_{\tilde{l}})}^2}$

Wherein, individual polynomial expression smoothness, ε is little positive number.

4. a kind of method that shallow water problems simulation is interrupted current numerical value that solves as claimed in claim 3, is characterized in that: the span of described ε is (10 ^-6, 10 ^-2).

5. a kind of method that shallow water problems simulation is interrupted current numerical value that solves as claimed in claim 4, is characterized in that: the value of described ε is ε=10 ^-6.

6. a kind of method that shallow water problems simulation is interrupted current numerical value that solves as claimed in claim 1, is characterized in that: in described step 6, the concrete account form of time step Δ t is:

$\mathrm{\&Delta;t}=\mathrm{CFL}\&CenterDot;\underset{i}{\min }\left(\frac{\underset{m}{\min }\left|L_m\right|}{\sqrt{g{h}_i}+\sqrt{u_i^2+v_i^2}}\right);$

Wherein, CFL represents CFL condition number.