Abstract
A FEM for unsaturated transient seepage is established by using a quadrilateral isoparametric element, considering the fact that the main permeability does not coincide with the axis situation. It creates a function by using the element’s node hydraulic head and shape function instead of the real head in the Richard seepage control equation. With the help of the Galerkin weighted residual method, a FEM equation is given for analyzing 2-dimensional transient seepage problem. Further, based on the Jacobi matrix and Gauss numerical integral, it determines the elements of stiffness and capacitance matrices. This FEM equation considers not only the anisotropic of soil but also the uncoincidence between permeability and the axis. It is a common form of transient seepage. In the end, two examples illustrate the node accuracy of the quadrilateral element and the correctness of this FEM equation.
1. Introduction
Transient seepage of unsaturated soils is a hot issue of research currently [1–3], For instance, dam seepage [4, 5], subgrade moisture field, the steadiness of rainfall side slope [6], and other unsaturated seepage problems can be considered according to two-dimensional seepage problems. The finite element method (FEM) is a highly valued analysis method.
Lam et al. [7] have used the triangular element to obtain the finite element equation and stiffness and capacitance matrices of transient seepage of unsaturated soils. This research literature is widely referenced [8–10], and the establishment of a great deal of computational and analytical work is based on the triangular element. In the analysis of some symmetry problems, the triangular element has relatively large errors, which is caused by the asymmetry of the triangular element itself.
Under the condition of shape functions sharing the same power, approximate functions constructed by values on the nodes and shape functions are defective, while the quadrilateral element has no similar problem. In this regard, some scholars begin to use the rectangular element to calculate some complex engineering problems [11, 12], which demonstrates an excellent trend; however, the commercial software constituted of opaque principles is used for calculation, no determination method of stiffness and capacitance matrices is provided, and literature studies related to stiffness and capacitance matrices of the rectangular element are still scarce.
Further research [13, 14] gives the calculating idea and method of the finite element of homogeneous isotropic soils using the three-dimensional isoparametric finite element. Due to the anisotropy and heterogeneity of soils, the phenomenon of the main seepage coefficient not coinciding with the coordinate axis often appears. In addition to using the triangular element method, existing studies have not been involved in how to establish the finite element equation when the direction of the main seepage coefficient does not coincide with the coordinate axis. Its geometric shapes of the boundary are mostly oblique lines, so the one-dimensional isoparametric element is able to meet the precision requirements of simulated geometric conditions.
At present, triangular element is still the simplest element form in the FEM of unsaturated soil’s transient seepage problem. The stiffness matrix and capacitance matrix of the triangular element can be used to code a calculation program. But the triangular element has relatively large errors when analyzing some symmetry problems. The quadrilateral element is symmetric, so it is more accurate than the triangular element in the condition of symmetrical flow. Quadrilateral element is also widely used in various seepage analysis software, but its stiffness matrix and capacitance matrix have not been reported yet. This makes it difficult for geotechnical engineers who need to code a calculation program to analyze seepage problems. In addition, soil anisotropy also needs to be reflected.
Based on this, this paper attempts to consider the phenomenon of the anisotropy of soils and main seepage coefficient not coinciding with the coordinate axis and adopt one-dimensional quadrilateral isoparametric element to research on establishing the calculating method of the finite element for transient seepage of unsaturated soils.
On this basis, by using the Galerkin weighted residual method, this paper establishes a transient seepage finite element calculation method of unsaturated soils. This method uses the one-dimensional quadrilateral isoparametric element and gives its stiffness matrix and capacitance matrix. This method also considers the anisotropy coefficient of soil and seepage is not consistent with the coordinate axis. Two examples are given: one is to illustrate the difference between triangular element and quadrilateral element and the other one gives a calculation of certain dam seepage to clarify the accuracy of the proposed method in seepage flow.
2. Transient Seepage of Anisotropic Soils
Assume that, in the plane problem, the seepage coefficients of soils are, respectively, known to be k1 and k2 in two mutually perpendicular directions s1 and s2. They are functions of matric suction of unsaturated soils (ua−) and have an angle α to the coordinate axis, as shown in Figure 1.
According to Darcy’s law,where h is the total head, i.e., the sum of the gravity head and pore water pressure head. The continuous derivative rule and the simple relation of trigonometric functions are used to change the hydraulic gradient in the s1 direction
Formula 2 is introduced into formula (1); then, the flow velocity in the s1 direction can be written as
Similarly, the flow velocity in the s2 direction can be obtained as
and are orthogonally decomposed in the direction of the coordinate axis, and their components on the x-axis and y-axis make up the flow velocity and along the direction of the axis:
At present, introducing formula (3) and (4) into (5), the following forms are obtained: where
On the other hand, the continuity condition of seepage of unsaturated soils containswhere is called the change coefficient of the water yield associated with matric suction (ua − ) change. Formula (6) is introduced into (8), and the relation between pore water pressure and the total head is as follows:
Then, the control equation of transient seepage for anisotropic can be written as follows:where . For unsaturated soils, the seepage coefficient is not constant but synergistically varies with the suction head. Besides, the total head on the Dirichlet boundary (Г1) and the normal flow velocity on the Neumann boundary (Г2) are known aswhere n represents the normal direction.
3. Discreteness of Spatial Domain
When the quadrilateral one-dimensional isoparametric element is adopted, as shown in Figure 2, the shape function on each node and the displacement mode at any point in the element, respectively, arewhere (xi, yi) and (ζi, ηi) are global and local coordinate values of point i.
In the element, an approximate head function is constructed by means of the product of the node head and the shape function to approach the real head of any point on the nodes, such as
To replace formula (10) and the head h in formula (11) representing the boundary condition, for similar reasons, formula (10) will produce margin in the integral surface domain Ω, denoted as RΩ; formula (11) will produce margin on the boundary of Г2, denoted as RГ2. When shape functions of the quadrilateral element are used to multiply these margins, respectively, the integral sum of RΩ and RГ2 in each integral domain should be zero, and this method is called Galerkin weighted residual method. According to the Galerkin method, the following mathematical expression can be obtained:
The method of subsection integration is adopted to the first item of formula (14), and the surface integral is converted into the curve integral combined with the Green formula to obtain the “weak” form of its equivalent integral. After finishing formula (14) by the above method,where the subscripts m = x, y and n = x, y indicate that m and n should be taken as x and y, respectively, and [N] is the matrix of shape functions of the quadrilateral isoparametric element: .
In formula (15), stiffness matrix [D], capacitance matrix [E], and flux boundary condition [F] are stipulated as follows:
Then, formula (15) can be simplified as
3.1. Stiffness Matrix [D]
For the isoparametric quadrilateral element, there are a total of 4 × 4 = 16 elements in the stiffness matrix, among which
To desire for the expression of stiffness matrix, first of all, it is necessary to convert the integral into the local coordinate system, which requires the assistance of the transition of Jacobi matrix.
The partial derivative of the global coordinate to the local coordinate can be written as
The coefficients in the formula are
All of them are only related to the node coordinates, which are constant for the local coordinates. Then, the determinant of the Jacobi matrix is given by
The coefficients among them are
Similarly, these coefficients are constant.
The relation between the shape function and the partial derivative of global and local coordinates is obtained by using Jacobi matrix:where [J]−1 represents the inverse matrix of Jacobi matrix. Formula (23) is expansively written, in which the shape function on the partial derivative of x is
The coefficients in the formula are
Similarly, the shape function on the partial derivative of y can be obtained as
The coefficients in the formula are
The integral domain is transformed by Jacobi matrix
The results obtained by formulas (24)∼(28) are introduced into formula (18), and the integral form of Dij can be obtained as follows:
The integrand fDij (ζ, η) in formula (29) is expanded as
In formula (30), only ζ and η are variable, and the remaining coefficients are constant, in which the constants are only related to global and local coordinates of 4 nodes in the element, according to the foregoing calculation discussion of these coefficients, when the Gauss numerical method is used to evaluate the integral in formula (29).where hs and ht are weighting coefficients, ζs and ηt are integral points, and l is the number of integral points being adopted. Generally, for the 4-node isoparametric element, the number of integral points is l × l = 2 × 2, which is sufficient to meet the precision requirement. At this point, respectively, values of the weighting coefficient and the integral point should be
After being introduced into formula (31), the final form of the stiffness matrix is
Among them, the expression of the integrand fDij is shown in formula (30).
3.2. Capacitance Matrix [E]
The capacitance matrix is the same as the stiffness matrix, which also has 16 items. Among them, the calculation formula of the element at any position is
By using the Gauss integral method, each element in the capacitance matrix can be obtained. Considering that the form of the integrand is relatively simple in the capacitance matrix, the expression of its integral form can be presented by simplification.
For any element in the capacitance matrix, after the expression of formula (12) and the determinant of formula (21) of Jacobi matrix are introduced, the integrand is written as
If formula (35) is fully expanded, a high-order polynomial related to ζ and η can be obtained. Considering the symmetry of range of integration in the capacitance matrix, the odd polynomials containing ζ or η are equal to zero after the integral, then these items with the integral of zero are deleted directly after the expansion, and formula (35) can be simplified as follows:
The coefficients in the formula are
Formula (36) is introduced into formula (34) to conduct the integral, and after being finished, the expression for each element in the capacitance matrix is as follows:
After local coordinates on nodes are introduced, the capacitance matrix can be simplified as
Among which matrices [E5], [E6], and [E7] are all constant matrices, they, respectively, are
3.3. Flux Boundary Conditions [F]
Matrix [F] is the matrix reflecting Neumann boundary conditions, which represents equivalent nodal flow in the normal direction.
4. Analysis of Calculation Cases
The above method can be used to calculate unsaturated seepage under various conditions. Because all of the existing calculating methods are not involved in the problem of the main seepage coefficient not coinciding with the coordinate axis, for the convenience of comparative analysis, so does validate calculation examples, which can be achieved by demanding Corner α in formula (7) to be zero. When analyzing the isotropic problem, k1 = k2, and at this time, no matter how Corner α changes, it will not affect the seepage coefficient value. For the one-dimensional seepage problem, the seepage coefficient (k1 or k2) in a certain direction can be further demanded to be zero. The coefficients simplified by formula (7) under different operating conditions can also obtain different forms of stiffness matrix expressions after being introduced into formula (30). In the case of steady seepage, as the nodal head does not vary from the time, the partial derivative of the head to the time in formula (17) is zero, which means the item of capacitance matrix can be ignored.
At present, based on the formula derivation in this paper, Fortran language self-compiled program is applied to calculate two different problems.
4.1. Node Errors of the Triangular Element
Due to the asymmetry of the triangle element, node errors will appear in the analysis of some symmetric problems, while the quadrilateral element can overcome these errors. Considering the one-dimensional unsaturated seepage problem on a horizontal soil column, the linear triangular and quadrilateral elements are used to conduct discretization, as shown in Figure 3.
The length of the soil column is 1 m, and the boundary condition of the head is given on both sides: the suction head on the left side is maintained at zero, and the suction head on the right side is maintained at −10 m. Side lengths of the element are 0.1 m without exception. The quadrilateral element is composed of 10 elements and 22 nodes. The triangular element needs 20 elements and 22 nodes. In the process of seepage calculation, because the seepage coefficient k is the function of the suction head h, the value of the nodal head cannot be obtained when seepage is steady unless many iteration computations have been conducted. Among which, the relation between the seepage coefficient k and the suction head h can be obtained by referring to literature [15], as shown in Figure 4.
The two values of upper and lower nodal heads in two different elements away from the same position on the left side are shown in Figure 5. The tri-1 curve in the figure represents the variation conditions of the head of the triangular element beside the next row of nodes in Figure 3, while the tri-2 curve represents the head beside the top row of nodes at the same position. The tri-aver curve represents the average value of the two nodal heads. It can be seen that, for the triangular element, the two nodal heads at the same position are not the same, which does not accord with the law of symmetric seepage. The rec-1 and rec-2 in the figure are the values of nodal heads calculated by the quadrilateral element. At the same location, they are exactly the same, and the sizes of them are basically equal to the average value of the two nodal heads. It is shown that, for the triangular element, the two nodal water heads at the same position are not the same, which does not accord with the law of symmetric seepage. The rec-1 and rec-2 in the figure are the node head values calculated by the quadrilateral element. At the same location, they are exactly the same. And the sizes of them are basically equal to the average value of two nodal heads of the triangle, which indicates node errors of the triangular element objectively exist, and the reason of errors is mainly for the asymmetry of the triangle in the calculation of symmetry problem.
4.2. Seepage Caused by the Water Flowing through the Earth-Rock Dam
In order to verify the correctness of the formula deriving and program compiling in this paper, the condition of two-dimensional seepage of earth-rock dam in literature [15] is calculated particularly, whose calculation results are compared with the original paper. The quadrilateral isoparametric element is used to discretize the dam body, of which the body contains a total of 312 units and 351 nodes. Discrete graphics can be seen in Figure 6.
In the calculation, in order to be consistent with the conditions in the literature, it is assumed that the seepage coefficient of the soil is isotropic, and the seepage coefficient function used in the analysis is still given out by Figure 4. Select the dam bottom as a reference surface; in addition, the value of λ used in the analysis is adopted for 0.01 m–1. The water level of the initial reservoir is above the reference surface 4 m and the earth dam is in steady state. From the beginning of the zero moment, the water level of the water reservoir suddenly rises to the position above the reference surface 10 m. Water seeps into the dam under the new raised water pressure. The amount of evaporation and infiltration on the boundary is not considered in the process of dam seepage. Besides, in the subsequent transient process, the water level remains unchanged. The rise condition of the phreatic line within the dam from the initial steady state to the delaying 511 days is shown in Figure 7. The variation trend of the phreatic line reflected by the figure is basically the same as that in literature [15]. The velocity field and equipotential line can be calculated after the nodal head is obtained, and the same details will not be explained any more.
5. Conclusion
In view of the fact that the main seepage coefficient not coinciding with the coordinate axis is not considered in the existing numerical calculation of the moisture field of unsaturated soils, this paper adopts the calculating method of the finite element for transient seepage of unsaturated soils. By using the shape function on the element node and the head value, an approximate function can be constructed to replace the head function in the Richard seepage control equation. Galerkin weighted residual method is used to construct finite element equations by making the integral of the error generated by the approximation in the element equal to zero. The “weak” form of weighted residual method is obtained by the further application of multivariate function integral method and Green formula transforming the equation form, based on which, the finite element form is presented to analyze the seepage problem by using the quadrilateral isoparametric element, and stiffness matrix and capacitance matrix are determined by the Gauss numerical integration method. This finite element method takes the anisotropy of soils, the main seepage coefficient not coinciding with the coordinate axis, and other cases into account, which is the general form of transient seepage of unsaturated soil analysis.
Data Availability
All data and models generated or used during the study are available in the submitted article. All data included in this study are available upon request by contact with the corresponding author.
Conflicts of Interest
The authors have no conflicts of interest to declare.
Acknowledgments
This study was supported by the Fund Program for the National Natural Science Foundation of China (no. 51078309), Key Research and Development Program of Shaanxi (Program no. 2018ZDCXL-SF-30-9), and Shaanxi Key Laboratory of Geotechnical and Underground Space Engineering, XAUAT. The authors would like to thank the team members of the key discipline of geotechnical engineering and the workers of the Geotechnical Engineering Laboratory.