Discretization of (1) seepage region (taking two-dimensional flow as an example)
Using numerical simulation technology to study groundwater movement, the aquifer in the hydrogeological model to be studied is discretized first. The so-called discretization means that the heterogeneous anisotropic aquifer in the seepage area to be studied is divided into many interconnected small equilibrium areas in a certain way, and each small equilibrium area is homogeneous and isotropic. In each small equilibrium region, aquifer parameters are regarded as constants; The central head value or the average head value under the condition is regarded as the representative head value of the small equilibrium area. There are usually two forms of subdivision (rectangle and polygon).
1. Rectangular equilibrium domain
It divides the balance into many small rectangular balance domains with two orthogonal parallel lines, as shown in Figure 7-3. When cutting, it is agreed that: ① the boundary of constant head or known head (a kind of boundary) should pass through the center of small equilibrium domain; ② The water-resisting boundary (the second boundary) coincides with the boundary of the small equilibrium region. This subdivision method is similar to rectangular coordinate system, and small areas and nodes (nodes) (center points of small balanced areas) are calibrated with appropriate numbers. Commonly used terms are:
Figure 7-3 Seepage area is divided into small rectangular equilibrium areas.
(According to Li Junting et al., 1987)
1) points, rows and columns, where the point (node) is the center point of a small area, and the grid is called a row horizontally and a column vertically.
2) Step size, which is divided into space step size (δ x, δ y, δ z) (Figure 7-3) and time step size (δ t).
3) The numbers of small areas and nodes are uniformly marked as (i, j), which means that small areas and nodes are located in the I-th row and the J-th column. ..
2. Polygon equilibrium domain
Because the polygon equilibrium domain is close to the complex boundary geometry, it is used more. It first divides the seepage basin into triangles, and then forms a polygonal equilibrium domain based on the triangles, as shown in Figure 7-4. Commonly used terms and precautions:
1) point element, bin, line element, the edge of triangle is called line element, the vertex of triangle is called point element (node or node), and the area of triangle is called bin;
2) The single internal angle of the triangle is required to be 30 ~ 90 when splitting;
3) The area behind the leakage section shall be consistent with the original area, and shall not be duplicated or cracked.
(2) Basic equilibrium discrete equation (taking regular grid finite difference method as an example)
The water exchange relationship between the equilibrium zone in (i, j) in Figure 7-3 and the adjacent equilibrium zone is shown in Figure 7-5.
Figure 7-4 Triangle of Seepage Area
Fig. 7-5(i, j) Schematic diagram of flow relationship in equilibrium area.
(According to Li Junting et al., 1987)
1) The equilibrium period is δTN+ 1:
δTN+ 1 = TN+ 1-tn0
Represents the head of point (i, j) at time tn.
2) If there is no vertical water alternation in the (i, j) equilibrium area, according to the principle of water balance, there are:
Groundwater hydraulics
The different balance periods in the X-axis direction are:
Groundwater hydraulics
Among them:
Groundwater hydraulics
3) Considering the difference between formula (7- 14) and formula (7- 15), different calculation results will be produced. The calculation scheme (difference scheme) will be written as the following general formula:
Groundwater hydraulics
Where: 0≤θ≤ 1. θ usually takes three cases: ① when θ = 0, it is called the display difference format of finite difference method; ② Symmetric (central) difference scheme of finite difference method when θ =1/2; ③θ= 1 is an implicit difference scheme of finite difference method.
The finite difference equation is actually an approximate expression of the basic differential equation, and its approximate degree can be analyzed by Taylor series. Through the difference expression of differential equation, it can be seen that there are errors in replacing differential equation with difference scheme, that is, simulating groundwater flow system with finite difference equation will produce errors.
(3) boundary conditions and treatment of vertical water exchange
No matter the boundary of known water head or the boundary of known discharge, the calculation point falls on the boundary, so it does not need to be included in the equilibrium equation. The same is true for the treatment of vertical water exchange. If the point coincides with the pumping well and the point is included in the equilibrium discrete equation, the pumping capacity will directly participate in the water balance in the equilibrium area where the point is located.
(4) Solving the equilibrium discrete equation
Obviously, under the condition of given aquifer parameters and boundary conditions, as long as the head values of the flow field at a certain moment are known, the head values of all points in the next time step can be calculated. That is to say, on the basis of knowing the initial conditions, the head value and flow field at different times can be calculated. In order to solve this kind of problem, from the point of view of widely using microcomputer processing, the over-relaxation iteration method has been adopted by many researchers.
(5) Application
To sum up, when the initial conditions, boundary conditions, vertical water exchange and aquifer parameters are known, the head values at different times and nodes in the seepage area can be calculated. At present, numerical methods have been widely used in the calculation of water quantity in groundwater resources evaluation, the calculation of mine groundwater drainage or the prevention and control of geological disasters caused by large-scale groundwater level decline.
At present, there are many groundwater numerical calculation softwares with strong adaptability and high simulation, which are widely used, such as MOP-FLOW (three-dimensional finite difference pore water numerical calculation software) and GWMS-3D (two-dimensional or three-dimensional groundwater flow and pollutant migration numerical simulation software).
(6) Examples
Through the study of examples, students can understand the process of solving by numerical methods. This process includes: ① generalizing hydrogeological conditions and establishing a conceptual model; (2) establishing a numerical model according to the hydrogeological conceptual model; (3) dividing the calculation area and sorting out the calculation data; (4) modifying the numerical model; ⑤ Verify the numerical model; ⑥ Use the model to forecast.
This example is located at the junction of alluvial and diluvial fans at the eastern foot of Taihang Mountain. The aquifer is Quaternary unconsolidated layer, with fine sand and silty sand in the upper part and gravel, coarse gravel plus soil layer and clay gravel layer in the lower part. The groundwater in the upper aquifer has been drained. At present, the mining depth is 40 ~ 80m, the water level is mostly below 10m, and the central area of the funnel has reached 30m. The buried depth of water level in some marginal areas is 2 ~ 10m.
1. Hydrogeological conceptual model
(1) The aquifer floor is a waterproof clay layer; (2) Aquifers are mainly heterogeneous and isotropic phreatic aquifers; ③ Three sides of the boundary of the calculation area are the first-class boundary with known water head, and the other side is the weak permeable layer with different degrees, with a calculation area of nearly 600km2；; (4) There are mining mines in this area; ⑤ The groundwater flow is an unsteady plane flow, and the flow conforms to Darcy flow.
2. Numerical model
1) differential equation:
Groundwater hydraulics
2) initial condition: H(x, y, 0)=H0(x, y)
3) A boundary condition: H(x, y, t) | γ 1 = H 1 (x, y, t).
4) The second boundary condition:
Groundwater hydraulics
Among them, w is the return term, which is obtained by the algebraic sum of precipitation infiltration and well production; N is the inner normal; Other symbols are the same as before.
3. Divide the calculation area and arrange the calculation data.
The calculation area is divided into 506 cells and 230 nodes, including 40 boundary points of the first type and 65,438+04 boundary points of the second type. Taking the dry season as the model correction period, the partition parameters of 10 are given, and the initial values of the parameters are given through experiments.
4. Modify the numerical model
The correction results show that the differential equation is consistent with the boundary conditions.
5. Verify the numerical model
Take the rainy season water level data to verify the water level of the seventh phase. According to the verification data, the fitting diagram of high and low water levels and other required fitting diagrams are drawn, which proves that the fitting degree is good and meets the requirements of the specification.
6. Model use
By using the verified realistic simulation model, the exploitation amount can be predicted according to the designed water level, or the water level drop in different periods can be predicted according to the designed exploitation amount, especially in the center of the funnel.