Simulation of Groundwater Flow in Saturated Soil

Applicability in construction engineering

Ground water flow rates predefine largely predetermine the construction methods and materials used for footings, basement walls, underground constructions and many other in-ground works. The stability of beds and banks of water reservoirs and channels also depends on the filtration in coastal terrain. Factoring the flow of ground water improves the modeling accuracy of othe rphysical processes in the ground. When distribution of thermal fields takes place in the ground, the convective heat transfer is caused by groundwater flow.

Water is able to flow through the ground because of the presence of pores, which are voids of various diameters and shapes that appear due to the fact that structural elements produced during terrain formation don’t fit flush with each other. The approach for modeling water flow processes differ according to the degree of water porosity and the velocity of water in pores: if the pores are in a saturated condition, the groundwater flow process is simulated based on the Darcy differential equation; water flow in unsaturated ground is described by Richard’s or Brinkman’s equation.

The Darcy equation

Let us consider the Darcy equation in detail. This equation was named after French fluid power engineer, Henry Philibert Gaspard Darcy, and is a second order partial differential equation.

(1)   \begin{equation*}\bigtriangledown\ * (-K\bigtriangledown H)=Q_s\end{equation*}

where Q_s - water source, [Q_s]=c^{-1}, K – coefficient of groundwater flow, [K] = m/s; H - hydraulic pressure, [H] = m; \bigtriangledown=\frac{\partial} {\partial x}\vec{i}+\frac{\partial}{\partial y}\vec{j}+\frac{\partial}{\partial z}\vec{k} - nabla operator.

Water flow velocity vector v, according to the Darcy’s law is:

(2)   \begin{equation*}v= -K\bigtriangledown H\end{equation*}

The following boundary conditions are set for equation (1):

1) Boundary condition of the first order H=H_{0} , here H_{0} – hydraulic pressure given;

2) Boundary condition of the second order \vec{n}K\bigtriangledown H =q_{0}, here q_{0} – water flow through the boundary surface, [q_{0}] =m/s, and \vec{n} – unit outward normal to the surface.

3) If the boundary condition of the first order is set on the whole boundary, we have the Dirichlet problem for the Darcy equation; if the boundary condition of the second order is set on the whole boundary, we have Neumann’s problem for the Darcy equation. In the general case, we assume that boundary conditions of the first and second order are set at the boundary, and we must therefore consider the problem with mixed boundary conditions.

Due to the arbitrariness of the computational domain, where equation (1) is set, there is no analytical solution for this equation. Therefore, to solve the Darcy equation, numerical methods such as finite-difference approximation of differential operators and the finite element method are used. Here, the first of these methods is considered.

Numerical solution of the Darcy equation

Let the computational domain be a parallelepiped, and the computational mesh of nodes set \omega{{(x_{i}, y_{j}, z_{k})}}:i= 1,\ldots, N_{x}; j= 1,\ldots, N_{y}; k= 1,\ldots, N_{z}. The boundary is the boundary between a material (ground) and the environment as well as the boundary of the computational domain. Each node p={(x_{i}, y_{j}, z_{k})} of the computational mesh is assigned vector d^{p}={(d_{x}^{p}, d_{y}^{p}, d_{z}^{p})}, which is defined the following way:

1) p=(0,0,0), if the node p doesn’t pertain to the boundary;

2) if node p pertains to the boundary, d^{p}=(sign(n_{x}), sign(n_{y}), sign(n_{z}))}, where n_{x}, n_{y}, n_{z} – components of the normal of the vector \vec{n} in the node p.

Let us build the system of equations for numerical solution of problem (1). H_{i,j,k} is the approximate value of the pressure H (x,y,z) in the node p=(x_{i},y_{j},z_{k}). To resolve H_{i,j,k} we build a system of linear algebraic equations. Consider 2 variants according to the geometric location of node p. If node p=(x_{i},y_{j},z_{k}) doesn’t pertain to the boundary, equation (1) is approximated thus:

If node p is located on the boundary, there may be 2 variants depending on the order of the boundary condition set on node. Let the boundary condition of the first order be set on node ; in this case, it is necessary to use H_{i,j,k}= H_{0}(x_{i},y_{j},z_{k}) equation. If the boundary condition of the second order is set on node p - the following equation is used:

Combining the equations, written for every node of the computational mesh, we obtain a system of linear algebraic equations:

                                                          Ah=b,                                                                (4)

Where the matrix of A coefficients has dimensions N\times N,N=N_{x}\times N_{y}\times N_{z}.

Given that the order of system (4) can be several million for the general case, the system (4) is solved by iterative methods such as the method of simple iterations, the Seidel method, the conjugate gradient method and the biconjugate gradient stabilized method. To accelerate the solution of the system of equations, we use parallel computation on a multicore CPU and also on a GPU (in particular, CUDA technology is used).

After solving the system (4), the approximate value of water flow velocity vector is obtained through the finite-difference approximation of equation (2) by means of numerical differentiation of the calculated hydraulic pressure.

During simultaneous modeling of thermal fields and groundwater flow, it is necessary to consider the effect of ground freezing. When the ground is frozen, the pores are filled with ice, thereby reducing the rate of water flow. Assume, that water flow coefficient K(x,y,z) for the frozen ground is K(x,y,z)= (1-(w(x,y,z))K_{0}(x,y,z), where w(x,y,z) – ice content in the ground (x,y,z), and K_{0}(x,y,z) – water flow coefficient of the ground thaw at the same point.

The example of the computation of water flow velocities in Frost 3D Universal software

Here is an example of the calculation of water flow velocities, which demonstrates that the consistent account of water flow is essential for modeling thermal fields in the ground.

Let the discussed section of the ground have a parallelepiped shape with the dimensions 15m\times10m\times10m, we specify it as P (see figure 1). Initially, the ground in the parallelepiped is in the thawed state. There is a hydraulic pressure H_{x=0}= 0.5m on the left face of parallelepiped P (face x=0m); on the right face (x=15m) hydraulic pressure is H_{x=15}=0m. Under the pressure gradient, water flow with the following velocity takes place in the parallelepiped:

v_{x}= K \frac {\Delta H}{\Delta x} = 0.001 \frac {0.5-0}{15-0} = 3.33 \times 10^{-4} m/s

Here K is water flow coefficient.

Assume that in parallelepiped P there is a cylindrical section of ground С with the bottom radius (parallel to plane O_{xy}) r=4m and height h=8m, the center of the cylinder has coordinates (5m,5m,5m). Initially, the ground in this cylinder is in the frozen state. As the ground in section C is frozen, the water flow velocity in this section equals 0.

We determine the distribution of temperature fields in the discussed section of the ground in accordance with the set field of water flow velocities with the temperature value on the left face of the parallelepiped P (face x=0 m) being T_{x=0} =+2^\circ C, and the right face being (x=15m) T_{x=15} =-2^\circ C. The thermophysical properties of the ground for sections P and С set in the computation, are listed in the tables 1 and 2. As a result of the computation, we obtain the temperature distribution shown in figures 2 and 3.

Assessment model of the influence of water flow on heat transfer

Figure 1 – Computational domain

Table 1. – Thermophysical properties

Parallelepiped P Cylinder C
Volumetric heat capacity of the thawed ground (\frac{J}{m^3\times K})
2140000 2140000
Volumetric heat capacity of the frozen ground (\frac{J}{m^3\times K})
2140000 2140000
Thermal conductivity of the thawed ground (\frac{W}{m\times K})
1.8 1.8
Thermal conductivity of the frozen ground (\frac{W}{m\times K})
1.8 1.8
Total gravimetric moisture of soil 1 1
Density of dry soil (\frac{kg}{m^3})
1000 1000
The amount of unfrozen water 01temperaturephasetransfer.png
Phase transition temperature ^\circ C
0 0
Water flow coefficient (\frac{m}{c})
0.001 0.001
Initial temperature ^\circ C
+4 -4

Table 2. – Boundary conditions

Boundary condition of a heat problem Boundary condition of a water flow problem
face x=0 m Heat exchange by Newton:
temperature – +2^\circ C
heat exchange coefficient – 10(\frac{W}{m^2\times K})
Hydraulic pressure – 0.5 m
face x=15 m Heat exchange by Newton:
temperature – -2^\circ C
heat exchange coefficient – 10(\frac{W}{m^2\times K})
Hydraulic pressure – 0.0 m
Other faces Heat flow –0(\frac{W}{m^2\timesK}) Water flow velocity – 0 m/s

 

Temperature distribution with statically set water flow values

Figure 2 – Distribution of thermal fields in degrees Celsius at a depth of 5 meters after 10 days for the given water flow velocity.

Thermal fields with statically set water flow values

Figure 3 – Temperature isolines in degrees Celsius at a depth of 5 meters after 10 days for the given water flow velocity.

In reality, water flow velocity in parallelepiped P due to the existence of section C isn’t only equal to 3.33x10-4 m/s, but also has a complex structure. The distribution of water velocity can be obtained by solving the Darcy equation. We obtain the distribution of water flow velocity fields in Frost 3D Universal software (4 and 5) and analyze how the thermal field changes with more precise values for groundwater flow velocities. The obtained computational results are given on figures 6 and 7.

Groundwater flows according to the Darcy equation

Figure 4 – Projection of the water flow velocity on X axis (micrometers per second), computed on the basis of the Darcy equation

Groundwater flow velocities during heat transfer

Figure 5 – Projection of the water flow velocity on Y axis (micrometers per second), computed on the basis of the Darcy equation

Temperature distribution if the groundwater flow values were computed correctly

Figure 6 – Distribution of thermal fields in Celsius degrees at a depth of 5 meters after 10 days for the given water flow velocity based on the Darcy equation

Temperature fields when the groundwater flow values were computed by means of the Darcy equation

Figure 7 – Temperature isolines in degrees Celsius degrees at a depth of 5 meters after 10 days for the given water flow velocity based on the Darcy equation

Comparing the results of the thermal field computation with the computation of water flow velocity on the basis of the Darcy equation (figures 6 and 7), we see that, without any computation (figures 2 and 3), there is a significant difference in temperature distribution. Particularly critical is the fact that the temperature difference in the same nodes for these computations is more than 10 0C in absolute terms.

Conclusion

In the above article we analyzed the approach of the computation of the groundwater flow, based on the numerical solution of the Darcy equation. This equation is justified for applications involving groundwater seepage.

We have given an intuitive demonstration of how essential it is to consistently factor in water flow when modeling thermal fields in the ground.

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>