CN103778298A - Improved multi-scale finite element method for stimulating two-dimensional water flow movement in porous media - Google Patents

Abstract

The invention discloses an improved multi-scale finite element method for stimulating two-dimensional water flow movement in porous media. The method comprises the steps as follows: firstly, converting a problem to be solved into a variational form; determining a boundary condition, setting a grid cell size h, subdividing a research area, and obtaining coarse grid cells; subdividing each coarse grid cell; according to a permeability coefficient K and the boundary condition of a basis function, solving the problem of degenerating elliptic type, and determining the basis function; according to the basis function, obtaining cell stiffness matrixes, and adding the same to obtain a total stiffness matrix; according to the boundary condition of the research area and a source sink term, obtaining a right-hand side; adopting an effective calculating method to solve the simultaneous equations of the total stiffness matrix and the right-hand side; and obtaining the hydraulic heads of all nodes in the research area. Through various simulation tests, the obtained result coincides with the analytical solution. Compared with the prior art, the method disclosed by the invention is similar to the same in precision, but the calculating time of the method is less than 10% of the calculating time in the prior art. The efficiency is greatly improved when the method is used for solving wide-range, long-time or complicated problems.

Description

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

Technical field

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

Background technology

Water resources problems is current and a most close major issue of human survival relation.Countries in the world have most of water in a lot of cities all to take from underground water.In addition,, in Geological Engineering activity, the distribution of underground water is also the factor that must consider.Therefore, the computing method of Study of The Underground water level and simulation, have very important significance for distribution situation and the forecast of measuring underground water.

The general equation of groundwater flow is described the distribution of steady flow by elliptic equation, 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 groundwater flow is described the distribution of steady flow by parabolic equation, 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 _yybe respectively xx direction, xy direction, yx direction, yy direction infiltration coefficient, W is source sink term.

Groundwater flow equation can solve with conventional Finite Element Method or finite difference method.But it is homogeneous that these methods require the medium in unit grid, in the time solving heterogeneous body problem, necessary fine dissection is little to guarantee the infiltration coefficient rate of change of inside, unit, can be approximated to be constant.In the time of the water flow simulation of carrying out in large-area survey region, fine dissection can produce very many nodes, and the memory space of computing machine is required very large and needs a large amount of computing times.Therefore, hydrologic research worker has proposed multi-scale finite elements method [Hou, T.Y., and X.H.Wu (1997)] and has solved this problem.

Multi-scale finite elements method, in the time of subdivision survey region, does not require that the infiltration coefficient of inside, unit must be approximately constant.The method, by the thin subdivision to unit, is constructed basis function by the elliptic equation that solves simplification on unit, can well catch by basis function the heterogeneous body character of medium.Particularly, for porous medium, the nonuniformity of this medium has generally comprised a lot of yardsticks, is usually reflected in the multiple dimensioned fluctuation of infiltration coefficient of medium.If utilize Finite Element to solve, in order to guarantee degree of accuracy, need in all small scales, solve the very large calculated amount of this process need and computing time.Because the basis function of multi-scale finite elements method can be caught the macroscopic information of the yardstick of medium, multi-scale method need to not solve in small scale, is highly suitable for engineering calculation.In addition, multi-scale finite elements method has proved that from mathematical derivation and numerical simulation it can be good at solving ellipse and Solving Parabolic Problems, and convergence, 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)].In addition, the resonance effects that the super sample technology in multi-scale finite elements method can avoid the yardstick of mesh scale and physical medium to produce, improves speed of convergence, reduces error [Hou, T.Y., and X.H.Wu (1997)].Compare with finite difference method with Finite Element, memory space and the requirement of computing time of multi-scale finite elements method to computing machine is lower, and can guarantee certain precision.But, because the thin subdivision method of triangle of traditional multi-scale finite elements method can produce point in more unit, cause in the time solving basis function the calculated amount of needs larger, if the survey region of current problem is too huge, excessive cycle, adopt traditional multi-scale method to need plenty of time and calculated amount to remove to solve basis function, efficiency needs to improve.

Summary of the invention

The present invention is directed to the deficiencies in the prior art, a kind of improved multi-level finite element modeling method of simulating two dimension flow motion in porous medium is provided, it can be applied to the current problem that solves steady flow and unsteady fluid flow.It is very identical that the method is simulated the result and the analytic solution that obtain.The in the situation that of given thin subdivision umber, basis function boundary condition, the precision of it and traditional multi-scale finite elements method approaches, but only has 10% of classic method the computing time needing.

The multi-level finite element modeling method of two dimension flow motion in improved simulation porous medium of the present invention, comprises the following steps:

(1) determine boundary condition according to the survey region that will simulate, set grid cell yardstick h, this survey region of subdivision, obtains coarse grid unit;

(2) in each coarse grid unit, centered by one or more interior points, adopt radial triangular element to carry out thin subdivision, obtain the refined net unit of this coarse grid unit;

(3) according to the boundary condition of coefficient of permeability K and basis function, solve the elliptic problem of degeneration, determine basis function, form Finite Element Space;

(4) calculate the stiffness matrix of each coarse grid unit, be added to obtain global stiffness matrix; According to the boundary condition of survey region, source sink term, calculate right-hand vector, form finite element equation;

(5) provide the efficient solution method of finite element equation.

In above-mentioned steps (1), what the subdivision of described formation coarse grid unit adopted is triangular element subdivision.

In above-mentioned steps (3), the infiltration coefficient on the coefficient of permeability K on refined net unit, the approximate all interior point of getting this unit of source sink term value, the mean value of source sink term

The present invention proposes a kind of thin subdivision method of brand-new multi-scale finite elements method, by adopting the radial thin subdivision centered by interior point, compared with classic method, in the time that element subdivision is same umber, internal node (the being unknown number) number that the present invention needs still less, is calculated basis function required calculated amount and computing time thereby reduced.The method can effectively be processed non-homogeneous porous medium Groundwater Flow Problems, and simple, efficiently accurate.Solving on a large scale, when long-time or challenge, it is high a lot of that the efficiency of the method is wanted.

With same triangle Δ _ijkregion is example, and Fig. 1 (a) is for being 27 parts with its subdivision of naming a person for a particular job in, and Fig. 1 (b) is for being 25 parts with its subdivision of naming a person for a particular job in three, and Fig. 1 (c) is the thin subdivision method of traditional multi-scale finite elements method, and 6 interior points of needs, by Δ _ijksubdivision is 25 parts; Subdivision efficiency, the present invention will be much larger than classic method.Calculating triangular unit Δ _ijkmultiple-Scale basis function time, as adopt Fig. 1 (a), owing to only having an interior point, can be directly to obtain the expression formula of a basis function; Adopt Fig. 1 (b), owing to there being 3 interior points, need to solve the matrix of 3 × 3 to obtain the expression formula of a basis function; Adopt Fig. 1 (c), owing to there being 6 interior points, need to solve the matrix of 6 × 6 to obtain the expression formula of a basis function; Calculated amount, the present invention will select and be less than classic method.By the gradual change dielectric model of the continuum Model of the two-dimentional steady flow under porous medium and unsteady fluid flow, two-dimentional steady flow and unsteady fluid flow, the sudden change dielectric model of two dimension steady flow and unsteady fluid flow, two dimension steady flow diving dielectric model, the numerical simulation of two dimension height nonisotropic medium model (six kinds of different scales), find that the present invention and analytic solution are coincide fine, precision and classic method approach, and only have 10% of classic method computing time.Result shows: for same current problem, and the precision of the result that the result precision that adopts the thin subdivision of Fig. 1 (b) to obtain obtains higher than employing Fig. 1 (a), but computing time is also bigger; Adopt result precision that the thin subdivision of Fig. 1 (c) obtains a little more than with the precision that adopts Fig. 1 (b) to obtain, but required computing time is much larger than (b) the required time of employing.

Accompanying drawing explanation

Fig. 1: (a): improved thin subdivision method, with the coarse grid element subdivision of naming a person for a particular job in be 27 parts; (b) improved thin subdivision method, with the coarse grid element subdivision of naming a person for a particular job in three be 25 parts; (c) traditional subdivision method is 25 parts by coarse grid subdivision.

Fig. 2: the continuum Model of two-dimentional steady flow, the absolute error of each method head on y=100m section.

Fig. 3: the gradual change dielectric model of two-dimentional unsteady fluid flow (impact Plain draw water model), the head that each method is simulated on y=5200m section.

Fig. 4: two-dimentional unsteady fluid flow sudden change dielectric model, each method on y=5000m section, simulate head.

Embodiment

Below in conjunction with specific embodiment, the present invention will be further explained, wherein relates to some shorthand notations, below for explaining:

LFEM: classic Finite Element.

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

MSFEM-L: traditional multi-scale finite elements method, uses linear barrier's condition.

MSFEM-O: traditional multi-scale finite elements method, uses oscillating edge movement condition.

MSFEM-os-O: traditional multi-scale finite elements method, use oscillating edge movement condition, use super sample technology.

MMSFEM-p-L: improved multi-scale finite elements method, use linear barrier's condition, use point in p to carry out thin subdivision.

MMSFEM-p-O: improved multi-scale finite elements method, use oscillating edge movement condition, use point in p to carry out thin subdivision.

MMSFEM-p-os--O: improved multi-scale finite elements method, use oscillating edge movement condition, use point in p to carry out thin subdivision, use super sample technology.

Embodiment 1: the continuum Model of two-dimentional steady flow

Study area is a square shaped cells: Ω=[50m, 150m] × [50m, 150m], coefficient of permeability K (x, y)=x ²m/d.Current equation is formula (1), and boundary condition is for determining head boundary condition

this model has analytic solution: H=3x ²+ y ².

Adopt LFEM, LFEM-F, MSFEM-L, MSFEM-O, MMSFEM-1-L, MMSFEM-1-O, MMSFEM-1-os-O solves this model.Wherein, LFEM-F is 1800 parts by study area subdivision, and additive method is 200 coarse grid unit by study area subdivision.It is 9 unit by thin each coarse grid subdivision that MSFEM adopts traditional triangle subdivision method, and MMSFEM adopts and improves thin subdivision method is 9 unit by thin each coarse grid subdivision.The super sample unit that super sample technology adopts is 1.01 times of coarse grid unit.

Fig. 2 is the head that said method calculates, in the absolute error of this section of y=100m.As can be seen from Figure 2, the error of FEM method is maximum; The precision of MSFEM-L and MMSFEM-L is very approaching; MMSFEM-1-O, LFEM-F, MSFEM-O obtains point-device solution, and their error is very approaching; MMSFEM-1-os-o has obtained best result.The result of MMSFEM-1-os-o is better than LFEM-F, result consistent [Hou, T.Y. that this obtains with S.J.Ye with T.Y.Hou in the time adopting MSFEM-os-O mimic water-depth model, and X.H.Wu (1997), Ye, S., Y, Xue, and C.Xie (2004)].

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

Study area is a square shaped cells: Ω=[0,10km] × [0m, 10km], and it is K (x, y)=1+x/40m/d that infiltration coefficient is increased to 250m/d to right side from 1m/d from the left side, border of study area.Current equation is formula (2), and head boundary is determined on border, left and right, and left margin is 10m, and right margin is 0m, upper the next water proof border.Water-bearing zone thickness is 10m, water storage coefficient S=0.00001-0.000009x/1000/m.Located a pumped well at coordinate (5200m, 5200m), flow is 1000m ³/ d, time of pumping is 5 days, time step is 1 day.The head H of initial time ₀(x, y)=10-x/1000m.This model does not have analytic solution, therefore, adopts the solution of LFEM-F as standard reference.

Adopt LFEM, LFEM-F, MSFEM-O, MMSFEM-3-O method solves this model.Wherein, LFEM-F is 125000 triangular elements by study area subdivision, and additive method is 1250 unit by study area subdivision.In this model, MMSFEM-3-O adopts radial thin subdivision, and it is all 100 triangular elements by thin coarse grid unit subdivision that MSFEM-O adopts the thin subdivision of tradition.Fig. 3 is the head of each method y=5200m section.The precision of MMSFEM and MSFEM is very approaching.In this model, the required CPU time of MMSFEM is 3 seconds, and MSFEM is 308 seconds, and MMSFEM has higher efficiency.

Embodiment 3: two-dimentional unsteady fluid flow sudden change dielectric model

Study area is a square shaped cells: Ω=[0,10km] × [0m, 10km], and infiltration coefficient is undergone mutation on this section of x=2480m, as x < 2480m, and K=2m/d; X >=2480m, K=1000m/d.Current equation is formula (2), and head boundary is determined on border, left and right, and left margin is 10m, and right margin is 0m, upper the next water proof border.Located a pumped well at coordinate (5000m, 5000m), flow is 6000m ³/ d, time of pumping is 3 days, time step is 1 day.Water-bearing zone thickness is 10m, water storage coefficient x < 2480m, S=0.000002/m; X >=2480m, S=0.0005/m.The head H of initial time ₀(x, y)=10-x/1000m.This model does not have analytic solution, therefore, adopts the solution of LFEM-F as standard reference.

Adopt LFEM, LFEM-F, MSFEM-O, MMSFEM-3-O solves this model.Wherein, LFEM-F is 80000 triangular elements by study area subdivision, and additive method is 400 unit by study area subdivision.In this model, MMSFEM adopts radial thin subdivision, and it is all 100 triangular elements by thin coarse grid unit subdivision that MSFEM adopts the thin subdivision of tradition.Fig. 4 is the head of each method y=5000m section.The precision of MMSFEM and MSFEM is very approaching.In this model, the required CPU time of MMSFEM is 1 second, and MSFEM is 73 seconds.

Embodiment 4: two-dimentional 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 set iteration error be η, iteration until | H ⁽ⁿ⁾-H ^(n-1)| < η.

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

Adopt MMSFEM-3-L to solve this model.Be that 800,1800,3200 unit carry out respectively three simulations by study area subdivision.In each simulation, be all 25 unit by thin coarse grid unit subdivision.The iterations that three simulations need is 4 times, and the CPU time spending is all less than 1s.

Claims (3)

1. a multi-level finite element modeling method for two dimension flow motion in improved simulation porous medium, is characterized in that, comprises the following steps:

(1) determine boundary condition according to the survey region that will simulate, set grid cell yardstick h, this survey region of subdivision, obtains coarse grid unit;

(2) in each coarse grid unit, centered by one or more interior points, adopt radial triangular element to carry out thin subdivision, obtain the refined net unit of this coarse grid unit;

(3) according to the boundary condition of coefficient of permeability K and basis function, solve the elliptic problem of degeneration, determine basis function, form Finite Element Space;

(4) calculate the stiffness matrix of each coarse grid unit, be added to obtain global stiffness matrix; According to the boundary condition of survey region, source sink term, calculate right-hand vector, form finite element equation;

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

2. the multi-level finite element modeling method of two dimension flow motion in improved simulation porous medium according to claim 1, is characterized in that in step (1), what the subdivision of described formation coarse grid unit adopted is triangular element subdivision.

3. according to the multi-level finite element modeling method of two dimension flow motion in improved simulation porous medium described in claim 1 or 2, it is characterized in that: in step (3), the infiltration coefficient on the coefficient of permeability K on refined net unit, the approximate all interior point of getting this unit of source sink term value, the mean value of source sink term.

Images