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 an improved multi-scale finite element method for simulating two-dimensional water flow motion in porous media. This method first converts the problem to be solved into a variational form; determines the boundary conditions, sets the grid unit scale h, subdivides the research area, and obtains coarse grid units; finely subdivides each coarse grid unit; according to The permeability coefficient K and the boundary conditions of the basis functions are used to solve the degenerated elliptic problem, and the basis functions are determined; the element stiffness matrix is obtained according to the basis functions, and the total stiffness matrix is obtained by adding them; the right-hand term is obtained according to the boundary conditions of the study area and the source and sink items; An effective calculation method solves the simultaneous equations of the total stiffness matrix and the right-hand term; obtains the hydraulic head of each node in the study area. Through various simulation experiments, the obtained results are in good agreement with the analytical solutions. Compared with the existing technology, the accuracy of this method is similar to it, but the calculation time is less than 10%; when solving large-scale, long-term or complex problems, the efficiency is greatly improved.

Description

An improved multiscale finite element method for simulating two-dimensional water flow in porous media

technical field

The invention relates to the field of hydraulics, in particular to an improved multi-scale finite element method for simulating two-dimensional water flow motion in porous media.

Background technique

The issue of water resources is an important issue closely related to human survival. Many cities around the world draw most of their water from groundwater. In addition, in geological engineering activities, the distribution of groundwater is also a factor that must be considered. Therefore, it is of great significance to study the calculation method and simulation of groundwater level for measuring and forecasting the distribution of groundwater.

The general equation of groundwater flow describes the distribution of steady flow by elliptic equation, and its two-dimensional form is:

$<span\; class="notranslate">\; -</span>\frac{<span\; class="notranslate">\; \&PartialD;</span>}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; xxx</span>}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\frac{<span\; class="notranslate">\; \&PartialD;</span>}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; yx</span>}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; yy</span>}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; W</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span>$

The general equation of groundwater flow is described by the parabolic equation for the distribution of steady flow, and its two-dimensional form is:

$\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; -</span>\frac{<span\; class="notranslate">\; \&PartialD;</span>}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; xxx</span>}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\frac{<span\; class="notranslate">\; \&PartialD;</span>}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; yx</span>}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; yy</span>}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; W</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span>$

Here H is the water head, K _xx , K _xy , K _yx , and K _yy are the permeability coefficients in the xx direction, xy direction, yx direction, and yy direction respectively, and W is the source-sink term.

The groundwater flow equation can be solved by conventional finite element method or finite difference method. However, these methods require that the medium in the cell grid is homogeneous. When solving heterogeneity problems, it must be finely divided to ensure that the rate of change of the permeability coefficient inside the cell is not large and can be approximated as a constant. When performing water flow simulation in a large-area research area, fine subdivision will generate a lot of nodes, which requires a large amount of computer storage and requires a lot of computing time. Therefore, hydrological researchers proposed a multi-scale finite element method [Hou, T.Y., and X.H.Wu (1997)] to solve this problem.

When the multi-scale finite element method divides the research area, it does not require that the permeability coefficient inside the unit must be approximately constant. This method constructs basis functions by solving simplified elliptic equations on the units by subdividing the units, which can well capture the heterogeneity of the medium through the basis functions. Especially for porous media, the heterogeneity of the media generally includes many scales, which is often reflected in the multi-scale fluctuation of the permeability coefficient of the media. If the finite element method is used to solve it, in order to ensure the accuracy, it needs to be solved on all small scales. This process requires a very large amount of calculation and calculation time. Since the basis function of the multiscale finite element method can capture the macroscopic information of the medium scale, the multiscale method does not need to be solved on a small scale, and is very suitable for engineering calculations. In addition, the multi-scale finite element method has been proved from mathematical derivation and numerical simulation that it can solve elliptic and parabolic problems very well, and it is convergent, accurate and 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 super-sampling technique in the multi-scale finite element method can avoid the resonance effect produced by the grid scale and the scale of the physical medium, improve the convergence speed, and reduce the error [Hou, T.Y., and X.H.Wu(1997)]. Compared with the finite element method and the finite difference method, the multi-scale finite element method has lower requirements on the storage capacity and calculation time of the computer, and can guarantee a certain accuracy. However, because the triangulation method of the traditional multi-scale finite element method will generate more points in the unit, it will require a large amount of calculation when solving the basis function. If the research area of the water flow problem is too large and the period is too long, Using the traditional multi-scale method requires a lot of time and calculation to solve the basis function, and the efficiency needs to be improved.

Contents of the invention

Aiming at the deficiencies in the prior art, the invention provides an improved multi-scale finite element method for simulating two-dimensional water flow motion in porous media, which can be applied to solve the water flow problems of steady flow and unsteady flow. The simulated results obtained by this method are in good agreement with the analytical solution. Given the number of subdivisions and the boundary conditions of basis functions, it is close to the accuracy of the traditional multi-scale finite element method, but the calculation time required is only 10% of the traditional method.

The improved multi-scale finite element method for simulating two-dimensional water flow motion in porous media according to the present invention comprises the following steps:

(1) Determine the boundary conditions according to the research area to be simulated, set the grid unit scale h, subdivide the research area, and obtain coarse grid units;

(2) In each coarse grid unit, with one or more internal points as the center, use radial triangular units for subdivision to obtain the fine grid unit of the coarse grid unit;

(3) According to the permeability coefficient K and the boundary conditions of the basis function, the degenerate ellipse problem is solved, the basis function is determined, and the finite element space is formed;

(4) Calculate the stiffness matrix of each coarse grid unit, and add the total stiffness matrix; according to the boundary conditions and source-sink items of the study area, calculate the right-hand end item to form a finite element equation;

(5) Provide effective solutions to finite element equations.

In the above step (1), triangular unit division is used for the division of the coarse grid unit.

In the above step (3), the permeability coefficient K and the source-sink item value on the fine grid unit are approximated by the average value of the permeability coefficient and source-sink item on all interior points of this unit

The present invention proposes a brand-new subdivision method of the multi-scale finite element method. By adopting the radial subdivision centered on the inner point, compared with the traditional method, when the unit is divided into the same number, the The invention requires fewer internal nodes (that is, unknowns), thereby reducing the calculation amount and calculation time required for calculating the basis functions. This method can effectively deal with the problem of groundwater flow in heterogeneous porous media, and is simple, efficient and accurate. This method is much more efficient when solving large-scale, long-term or complex problems.

Taking the area of the same triangle _Δijk as an example, Figure 1(a) divides it into 27 parts with one interior point, Figure 1(b) divides it into 25 parts with three interior points, Figure 1(c ) is a subdivision method of the traditional multi-scale finite element method, which requires 6 interior points, and divides _Δijk into 25 parts; the present invention is far greater than the traditional method in terms of subdivision efficiency. When calculating the multi-scale basis function of the triangular unit _Δijk , as shown in Figure 1(a), since there is only one interior point, the expression of a basis function can be directly obtained; using Figure 1(b), since there are 3 interior points point, it is necessary to solve a 3×3 matrix to obtain a basis function expression; using Figure 1(c), since there are 6 interior points, it is necessary to solve a 6×6 matrix to obtain a basis function expression; The present invention is chosen to be less than the traditional method in terms of calculation amount. Through the continuum model of two-dimensional steady flow and unsteady flow in porous media, the gradual change medium model of two-dimensional steady flow and unsteady flow, the sudden change medium model of two-dimensional steady flow and unsteady flow, two-dimensional steady flow diving Numerical simulation of the medium model and two-dimensional highly heterogeneous medium model (six different scales), it is found that the present invention is in good agreement with the analytical solution, the accuracy is close to the traditional method, and the calculation time is only 10% of the traditional method. The results show that: for the same water flow problem, the accuracy of the result obtained by subdividing Fig. 1(b) is higher than that obtained by using Fig. 1(a), but the calculation time is also slightly longer; The accuracy of the results obtained by subdivision is slightly higher than that obtained by using Figure 1(b), but the required calculation time is much longer than that required by using (b).

Description of drawings

Figure 1: (a): Improved fine subdivision method, using one interior point to divide a coarse grid unit into 27 parts; (b) Improved fine subdivision method, using three interior points to divide a coarse grid unit into 27 parts; The unit is divided into 25 parts; (c) The traditional subdivision method divides the coarse grid into 25 parts.

Figure 2: The continuum model of two-dimensional steady flow, the absolute error of the water head on the y=100m profile of each method.

Figure 3: Gradient medium model of two-dimensional unsteady flow (pumping model of flood plain), water head simulated by each method on y=5200m profile.

Figure 4: Two-dimensional unsteady flow abrupt change medium model, the water head simulated by each method on the y=5000m profile.

detailed description

The present invention is further explained below in conjunction with specific embodiment, wherein involves some abbreviation symbols, and the following are notes:

LFEM: Classical Finite Element Method.

LFEM-F: classical finite element method (fine subdivision).

MSFEM-L: Traditional multiscale finite element method, using linear boundary conditions.

MSFEM-O: Traditional multiscale finite element method, using oscillatory boundary conditions.

MSFEM-os-O: Traditional multiscale finite element method, using oscillatory boundary conditions, using supersampling techniques.

MMSFEM-p-L: Modified multiscale finite element method with linear boundary conditions and subdivision using p interior points.

MMSFEM-p-O: Modified multiscale finite element method with oscillatory boundary conditions and subdivision using p interior points.

MMSFEM-p-os--O: Modified multiscale finite element method, using oscillatory boundary conditions, subdivision with p interior points, using supersampling technique.

Embodiment 1: the continuum model of two-dimensional steady flow

The research area is a square unit: Ω=[50m, 150m]×[50m, 150m], permeability coefficient K(x,y)=x ² m/d. The water flow equation is formula (1), and the boundary condition is constant head boundary condition This model has an analytical solution: H=3x ² +y ² .

The model is solved by LFEM, LFEM-F, MSFEM-L, MSFEM-O, MMSFEM-1-L, MMSFEM-1-O, MMSFEM-1-os-O. Among them, LFEM-F divides the study area into 1800 parts, and other methods divide the study area into 200 coarse grid units. MSFEM uses the traditional triangulation method to subdivide each coarse grid into 9 units, and MMSFEM adopts the improved fine subdivision method to subdivide each coarse grid into 9 units. The super sample unit used by the super sample technology is 1.01 times of the coarse grid unit.

Figure 2 shows the absolute error of the water head calculated by the above method at the section y=100m. It can be seen from Figure 2 that the error of the FEM method is the largest; the accuracy of MSFEM-L and MMSFEM-L is very close; MMSFEM-1-O, LFEM-F, and MSFEM-O all obtain very accurate solutions, and their errors are very close ; MMSFEM-1-os-o obtained the best results. The result of MMSFEM-1-os-o is better than that of LFEM-F, which is consistent with the results obtained by T.Y.Hou and S.J.Ye when using MSFEM-os-O to simulate the water level model [Hou, T.Y., and X.H.Wu(1997), Ye , S., Y, Xue, and C. Xie (2004)].

Example 2: Gradient medium model of two-dimensional unsteady flow (pumping model of flood plain)

The research area is a square unit: Ω=[0,10km]×[0m,10km], the permeability coefficient increases from 1m/d to 250m/d from the left side to the right side of the boundary of the research area, namely K(x,y)= 1+x/40m/d. The water flow equation is formula (2), the left and right boundaries define the water head boundary, the left boundary is 10m, the right boundary is 0m, and the upper and lower water barrier boundaries. The thickness of the aquifer is 10m, and the water storage coefficient S=0.00001-0.000009x/1000/m. There is a pumping well at the coordinates (5200m, 5200m), the flow rate is 1000m ³ /d, the pumping time is 5 days, and the time step is 1 day. The water head at the initial moment H ₀ (x,y)=10-x/1000m. There is no analytical solution for this model, therefore, the solution of LFEM-F is used as the standard reference.

The model is solved by LFEM, LFEM-F, MSFEM-O, MMSFEM-3-O methods. Among them, LFEM-F divides the study area into 125,000 triangular units, and other methods divide the study area into 1250 units. In this model, MMSFEM-3-O adopts radial subdivision, and MSFEM-O adopts traditional subdivision to subdivide the coarse grid unit into 100 triangular units. Figure 3 shows the water head of each method at y=5200m section. The accuracy of MMSFEM is very close to that of MSFEM. In this model, the CPU time required by MMSFEM is 3 seconds, and that of MSFEM is 308 seconds. MMSFEM has higher efficiency.

Example 3: Two-dimensional unsteady flow abrupt change medium model

The research area is a square unit: Ω=[0,10km]×[0m,10km], the permeability coefficient changes abruptly at x=2480m, that is, when x<2480m, K=2m/d; x≥2480m, K =1000m/d. The water flow equation is formula (2), the left and right boundaries define the water head boundary, the left boundary is 10m, the right boundary is 0m, and the upper and lower water barrier boundaries. There is a pumping well at the coordinates (5000m, 5000m), the flow rate is 6000m ³ /d, the pumping time is 3 days, and the time step is 1 day. The thickness of the aquifer is 10m, the water storage coefficient x<2480m, S=0.000002/m; x≥2480m, S=0.0005/m. The water head at the initial moment H ₀ (x,y)=10-x/1000m. There is no analytical solution for this model, therefore, the solution of LFEM-F is used as the standard reference.

Use LFEM, LFEM-F, MSFEM-O, MMSFEM-3-O to solve this model. Among them, LFEM-F divides the study area into 80,000 triangular units, and other methods divide the study area into 400 units. In this model, MMSFEM adopts radial subdivision, and MSFEM adopts traditional subdivision to subdivide the coarse grid unit into 100 triangular units. Figure 4 shows the water head of each method at y=5000m section. The accuracy of MMSFEM is very close to that of MSFEM. In this model, the CPU time required by MMSFEM is 1 second, and that of MSFEM is 73 seconds.

Example 4: Two-dimensional Steady Flow Underwater Flow Model (Nonlinear Model)

The current equation is the submerged flow equation:

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

$\mathrm{<span\; class="notranslate">\; K</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left(\begin{array}{cc}\mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; b</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; b</span>}<span\; class="notranslate">\; )</span>\end{array}\right)$

Due to the nonlinear model, this method requires iteration:

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

And set the iteration error as η, that is, iterate until |H ⁽ⁿ⁾ -H ^(n-1) |<η.

The research area of this model is: Ω=[0,1m]×[0,1m], the boundary water head is the boundary of constant water head and both are 0.b=-4m, T=(1+x)(1+y).η =0.00001, H ₀ =0. This model has an analytical solution: H=xy(1+x)(1+y).

The model is solved using MMSFEM-3-L. The study area was divided into 800, 1800, and 3200 units for three simulations respectively. In each simulation, the coarse grid unit is subdivided into 25 units. The number of iterations required for the three simulations is 4, and the CPU time spent is 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.