CN103778298B - The multi-level finite element modeling method of two dimension flow motion in the simulation porous media improved - Google Patents

Abstract

The invention discloses the multi-level finite element modeling method of two dimension flow motion in the simulation porous media of a kind of improvement.First the problem that needs solve is converted into variational form by the method；Determine boundary condition, set grid cell yardstick h, subdivision survey region, obtain coarse grid cell；Each coarse grid cell is carried out thin subdivision；According to coefficient of permeability K and the boundary condition of basic function, solve the elliptic problem of degeneration, determine basic function；Obtain element stiffness matrix according to basic function, be added and obtain global stiffness matrix；Right-hand vector is obtained according to survey region boundary condition and source sink term；Effective computational methods are used to solve global stiffness matrix and right-hand vector Simultaneous Equations；Try to achieve the head of each node in survey region.Being tested by multiple analog, the result obtained matches with analytic solutions.Compared with prior art, the precision of the method is close with it, but the time of calculating is less than its 10%；Solving on a large scale, during long-time or challenge, efficiency increases substantially.

Description

The multi-level finite element modeling method of two dimension flow motion in the simulation porous media improved

Technical field

The present invention relates to field of hydraulics, be specifically related to a kind of improvement simulating the motion of two dimension flow in porous media Multi-level finite element modeling method.

Background technology

Water resources problems is current and that human survival relation the is the closest major issue.Countries in the world have very Most of water in many cities all takes from subsoil water.Additionally, in Geological Engineering activity, the distribution of subsoil water is also It it is the factor that must take into.Therefore, the computational methods of Study of The Underground water level and simulation, for measuring subsoil water Distribution situation has very important significance with forecast.

The general equation of subsurface flow is described the distribution of stationary flow by elliptic equation, and its two dimensional form is:

$-\frac{\∂}{\∂x}(K_{\mathrm{xx}}\frac{\∂H}{\∂x}+K_{\mathrm{xy}}\frac{\∂H}{\∂y})-\frac{\∂}{\∂y}(K_{\mathrm{yx}}\frac{\∂H}{\∂x}+K_{\mathrm{yy}}\frac{\∂H}{\∂y})=W,---\left(1\right),$

The general equation of subsurface flow is described the distribution of stationary flow by parabolic equation, and its two dimensional form is:

$\frac{\∂H}{\∂t}-\frac{\∂}{\∂x}(K_{\mathrm{xx}}\frac{\∂H}{\∂x}+K_{\mathrm{xy}}\frac{\∂H}{\∂y})-\frac{\∂}{\∂y}(K_{\mathrm{yx}}\frac{\∂H}{\∂x}+K_{\mathrm{yy}}\frac{\∂H}{\∂y})=W,---\left(2\right),$

Here H is head, K_xx、K_xy、K_yx、K_yyIt is respectively xx direction, xy direction, yx direction, yy Direction infiltration coefficient, W is source sink term.

Subsurface flow equation can solve with conventional Finite Element Method or finite difference method.But these Method requires that the medium in unit grid is homogenizing, when solving heterogeneous body problem, then must fine dissection with Ensure that the infiltration coefficient rate of change within unit is little, can be approximated to be constant.Carrying out large-area study area During water flow simulation in territory, fine dissection can produce very many nodes, requires the biggest to the amount of storage of computer And need the substantial amounts of calculating time.Therefore, hydrologic research worker propose Multi-scale remote sensing [Hou, T.Y., and X.H.Wu (1997)] solve this problem.

Multi-scale remote sensing is when subdivision survey region, it is not required that the infiltration coefficient within unit must approximate For constant.The method, by the thin subdivision to unit, constructs by solving the elliptic equation of simplification on unit Basic function, can well catch the heterogeneous body character of medium by basic function.Especially for porous media, The anisotropism of this medium has generally comprised a lot of yardstick, is usually reflected in the multiple dimensioned ripple of the infiltration coefficient of medium On Dong.Solve if, with Finite Element, in order to ensure degree of accuracy, need to solve on all little yardsticks, This process needs the biggest amount of calculation and the time of calculating.Owing to the basic function of Multi-scale remote sensing is permissible Catching the macroscopic information of the yardstick of medium, multi-scale method need not solve on little yardstick, is highly suitable for work Journey calculates.Additionally, Multi-scale remote sensing from mathematical derivation and numerical simulation, demonstrate it can be very Good ellipse and the Solving Parabolic Problems of solving, and restrain, accurately, efficient [Hou, T.Y., and X.H. Wu(1997),X.H.Wu,and Z.Cai(1999),W.Deng et al.(2008),Ye,S.,Y,Xue,and C. Xie(2004)].Additionally, super sample technology in Multi-scale remote sensing is it can be avoided that mesh scale and physics The resonance effect that the yardstick of medium produces, improves convergence rate, reduces error [Hou, T.Y., and X.H. Wu(1997)].Compared with Finite Element and finite difference method, Multi-scale remote sensing is to computer The requirement of amount of storage and the time of calculating is relatively low, and ensure that certain precision.But, due to traditional many chis The triangle thin subdivision method of degree Finite Element can produce point in more unit, causes the needs when solving basic function Amount of calculation relatively big, if the survey region of water flow problems is the hugest, excessive cycle, use traditional many chis Degree method needs plenty of time and amount of calculation to remove to solve basic function, and efficiency needs to improve.

Summary of the invention

The present invention is directed to the deficiencies in the prior art, it is provided that two dimension flow fortune in a kind of simulation porous media The multi-level finite element modeling method of dynamic improvement, its current that can apply to solve stationary flow and unsteady fluid flow are asked Topic.The result that the method simulation obtains is the most identical with analytic solutions.On given thin subdivision number, basic function border In the case of condition, it is close with the precision of traditional multiscale transform Finite Element, but the calculating time needed only has The 10% of traditional method.

The multi-level finite element modeling method of two dimension flow motion, bag in the simulation porous media of improvement of the present invention Include following steps:

(1) determine boundary condition according to survey region to be simulated, set grid cell yardstick h, subdivision This survey region, obtains coarse grid cell；

(2) in each coarse grid cell, centered by one or more interior points, radial triangle is used Shape unit carries out thin subdivision, obtains the refined net unit of this coarse grid cell；

(3) according to coefficient of permeability K and the boundary condition of basic function, the elliptic problem of degeneration is solved, really Determine basic function, form Finite Element Space；

(4) calculate the stiffness matrix of each coarse grid cell, be added to obtain global stiffness matrix；According to survey region Boundary condition, source sink term, calculate right-hand vector, forms finite element equation；

(5) efficient solution method of finite element equation is provided.

In above-mentioned steps (1), the subdivision of described formation coarse grid cell uses triangular element subdivision.

In above-mentioned steps (3), coefficient of permeability K, source sink term value approximation on refined net unit take this unit All interior point on infiltration coefficient, the meansigma methods of source sink term

The present invention proposes the thin subdivision method of a kind of brand-new Multi-scale remote sensing, by point within using Centered by radial thin subdivision, compared with traditional method, when being same number by element subdivision, this Internal node (i.e. unknown number) number of bright needs is less, thus decreases the amount of calculation calculated needed for basic function With the time of calculating.The method can effectively process non-homogeneous porous medium Groundwater Flow Problems, and the easiest OK, the most accurately.Solving on a large scale, during long-time or challenge, the efficiency of the method is much higher.

With same triangle Δ_ijkAs a example by region, Fig. 1 (a) for being 27 parts with putting its subdivision in one, Fig. 1 B (), for being 25 parts with putting its subdivision in three, Fig. 1 (c) is carefully cuing open of traditional multiscale transform Finite Element Divide method, need 6 interior points, by Δ_ijkSubdivision is 25 parts；From subdivision efficiency, the present invention is significantly larger than and passes System method.Calculating triangular unit Δ_ijkMultiple-Scale basic function time, as used Fig. 1 (a), due to only One interior point, can be directly to obtain the expression formula of a basic function；Use Fig. 1 (b), in having 3 Point, needs the matrix solving 3 × 3 to obtain the expression formula of a basic function；Use Fig. 1 (c), due to There are 6 interior points, need the matrix solving 6 × 6 to obtain the expression formula of a basic function；From amount of calculation The present invention to select less than traditional method.Continuous by the Two-dimecnsional steady flow under porous media and unsteady fluid flow The gradual change dielectric model of dielectric model, Two-dimecnsional steady flow and unsteady fluid flow, Two-dimecnsional steady flow and unsteady fluid flow Sudden change dielectric model, Two-dimecnsional steady flow diving dielectric model, (six kinds are not or not two dimensional height nonisotropic medium model Same yardstick) numerical simulation, find identical fine of the present invention and analytic solutions, precision is close with traditional method, The calculating time only has the 10% of traditional method.Result shows: for same water flow problems, uses Fig. 1 (b) The result precision that thin subdivision obtains is higher than the precision of the result using Fig. 1 (a) to obtain, but the time of calculating is the most bigger； The precision that the result precision using the thin subdivision of Fig. 1 (c) to obtain slightly above obtains with employing Fig. 1 (b), but required The calculating time is much larger than the time used needed for (b).

Accompanying drawing explanation

Fig. 1: (a): the thin subdivision method of improvement, was 27 parts with in one o'clock by a coarse grid cell subdivision;(b) The thin subdivision method improved, was 25 parts with in three o'clock by a coarse grid cell subdivision;(c) tradition subdivision side Method, is 25 parts by coarse grid subdivision.

The continuum Model of Fig. 2: Two-dimecnsional steady flow, each method the absolute of head on y=100m section misses Difference.

The gradual change dielectric model (draw water model in impact Plain) of Fig. 3: two dimension unsteady fluid flow, each method is at y=5200m The head of simulation on section.

Fig. 4: two dimension unsteady fluid flow sudden change dielectric model, what each method was simulated on y=5000m section head.

Detailed description of the invention

Below in conjunction with specific embodiment, the present invention will be further explained, is directed to some shorthand notations, with Lower for explaining:

LFEM: classic Finite Element.

LFEM-F: classic Finite Element (fine dissection).

MSFEM-L: traditional multiscale transform Finite Element, uses linear barrier's condition.

MSFEM-O: traditional multiscale transform Finite Element, uses oscillating edge movement condition.

MSFEM-os-O: traditional multiscale transform Finite Element, uses oscillating edge movement condition, uses super sample technology.

The Multi-scale remote sensing of MMSFEM-p-L: improvement, uses linear barrier's condition, clicks in using p The thin subdivision of row.

The Multi-scale remote sensing of MMSFEM-p-O: improvement, uses oscillating edge movement condition, clicks in using p The thin subdivision of row.

The Multi-scale remote sensing of MMSFEM-p-os--O: improvement, uses oscillating edge movement condition, in using p Point carries out thin subdivision, uses super sample technology.

Embodiment 1: the continuum Model of Two-dimecnsional steady flow

Study area is a square shaped cells: Ω=[50m, 150m] × [50m, 150m], infiltration coefficient K (x, y)=x²m/d.Current equation is formula (1), and boundary condition is for determine head boundary conditionThis model has analytic solutions: H=3x²+y²。

Use LFEM, LFEM-F, MSFEM-L, MSFEM-O, MMSFEM-1-L, MMSFEM-1-O, MMSFEM-1-os-O solves this model.Wherein, study area subdivision is 1800 parts by LFEM-F, its other party Study area subdivision is 200 coarse grid cell by method.MSFEM use traditional triangle subdivision method by each slightly The thin subdivision of grid is 9 unit, and MMSFEM employing improves thin subdivision method and by thin for each coarse grid subdivision is 9 unit.Super sample unit is coarse grid cell 1.01 times of super sample technology employing.

Fig. 2 is the head that said method calculates, in the absolute error of this section of y=100m.Can from Fig. 2 Knowing, the error of FEM method is maximum；The precision of MSFEM-L Yu MMSFEM-L is closely； MMSFEM-1-O, LFEM-F, MSFEM-O obtain point-device solution, and their error connects very much Closely；MMSFEM-1-os-o obtains best result.The result of MMSFEM-1-os-o is better than LFEM-F, this and T.Y.Hou and S.J.Ye obtain when using MSFEM-os-O mimic water-depth model Result is consistent [Hou, T.Y., and X.H.Wu (1997), Ye, S., Y, Xue, and C.Xie (2004)].

Embodiment 2: the gradual change dielectric model (draw water model in impact Plain) of two dimension unsteady fluid flow

Study area is a square shaped cells: Ω=[0,10km] × [0m, 10km], and infiltration coefficient is from the limit of study area 250m/d i.e. K (x, y)=1+x/40m/d is increased to right side from 1m/d on the left of boundary.Current equation is public Formula (2), right boundary determines head boundary, and left margin is 10m, and right margin is 0m, upper the next water proof border. Water-bearing layer thickness is 10m, water storage coefficient S=0.00001-0.000009x/1000/m.At coordinate There is a pumped well at (5200m, 5200m) place, and flow is 1000m³/ d, time of pumping is 5 days, time step A length of 1 day.The head H of initial time₀(x, y)=10-x/1000m.This model does not has analytic solutions, because of This, use the solution of LFEM-F as standard reference.

Using LFEM, LFEM-F, MSFEM-O, MMSFEM-3-O method solves this model.Wherein, Study area subdivision is 125000 triangular elements by LFEM-F, and study area subdivision is 1250 by additive method Individual unit.In this model, MMSFEM-3-O uses radial thin subdivision, and MSFEM-O uses tradition thin Thin for coarse grid cell subdivision is all 100 triangular elements by subdivision.Fig. 3 is each method y=5200m section Head.The precision of MMSFEM Yu MSFEM is closely.In this model, needed for MMSFEM CPU time is 3 seconds, and MSFEM is 308 seconds, and MMSFEM has higher efficiency.

Embodiment 3: two dimension unsteady fluid flow sudden change dielectric model

Study area is a square shaped cells: Ω=[0,10km] × [0m, 10km], and infiltration coefficient is at x=2480m Undergo mutation on this section, i.e. as x ＜ 2480m, K=2m/d;X >=2480m, K=1000m/d.Current Equation is formula (2), and right boundary determines head boundary, and left margin is 10m, and right margin is 0m, upper bottom every Water boundaries.Having a pumped well at coordinate (5000m, 5000m) place, flow is 6000m³/ d, time of pumping Being 3 days, time step is 1 day.Water-bearing layer thickness is 10m, water storage coefficient X ＜ 2480m, S=0.000002/m;X >=2480m, S=0.0005/m.The head of initial time H₀(x, y)=10-x/1000m.This model does not has analytic solutions, therefore, uses the solution of LFEM-F as mark Quasi-reference.

Using LFEM, LFEM-F, MSFEM-O, MMSFEM-3-O solve this model.Wherein, Study area subdivision is 80000 triangular elements by LFEM-F, and study area subdivision is 400 by additive method Unit.In this model, MMSFEM uses radial thin subdivision, MSFEM use the thin subdivision of tradition all by The thin subdivision of coarse grid cell is 100 triangular elements.Fig. 4 is the head of each method y=5000m section. The precision of MMSFEM Yu MSFEM is closely.CPU time in this model, needed for MMSFEM Being 1 second, MSFEM is 73 seconds.

Embodiment 4: Two-dimecnsional steady flow underground flow model (nonlinear model)

Current equation is underground flow equation:

-K (x, y, H) H=W,

$K(x,y,H)=\left(\begin{array}{cc}T(H-b)& 0\\ 0& T(H-b)\end{array}\right)$

Owing to being nonlinear model, the method needs iteration:

-▽·K(x,y,H^(n-1))▽H⁽ⁿ⁾=W,

And the iteration error set as η, i.e. iteration until | H⁽ⁿ⁾-H^(n-1)| ＜ η.

This scale-model investigation district is: Ω=[0,1m] × [0,1m], border head is for determine head boundary and to be 0.b=-4m, T=(1+x) (1+y). η=0.00001, H₀=0. this model has analytic solutions: H=xy (1+x) (1+y).

MMSFEM-3-L is used to solve this model.It is 800,1800,3200 unit by study area subdivision Carry out three simulations respectively.Every time in simulation, it is all 25 unit by thin for coarse grid cell subdivision.Three moulds The iterations intending needing is 4 times, and the CPU time spent is respectively less than 1s.

Claims (3)

1. one kind improve simulation porous media in two dimension flow motion multi-level finite element modeling method, it is characterised in that include with Lower step:

(1) determine boundary condition according to survey region to be simulated, set grid cell yardstick h, this survey region of subdivision, Obtain coarse grid cell；

(2) in each coarse grid cell, centered by one or more interior points, radial triangular element is used to carry out carefully Subdivision, obtains the refined net unit of this coarse grid cell；

(3) according to coefficient of permeability K and the boundary condition of basic function, solve the elliptic problem of degeneration, determine basic function, shape Become Finite Element Space；

(4) calculate the stiffness matrix of each coarse grid cell, be added to obtain global stiffness matrix；Boundary condition according to survey region, source Remittance item, calculates right-hand vector, forms finite element equation；

(5) efficient solution method of finite element equation is provided, tries to achieve the head of each node in survey region.

The most according to claim 1, the multi-level finite element modeling method of two dimension flow motion in the simulation porous media improved, it is special Levying and be in step (1), the subdivision of described formation coarse grid cell uses triangular element subdivision.

The multi-level finite element modeling method of two dimension flow motion in the simulation porous media of improvement the most according to claim 1 or claim 2, It is characterized in that: in step (3), the coefficient of permeability K on refined net unit, source sink term value take on all interior point of this unit Infiltration coefficient, the meansigma methods of source sink term.