- Email: [email protected]
Modeling unconfined seepage flow in soil-rock mixtures using the numerical manifold method
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Modeling unconfined seepage flow in soil-rock mixtures using the numerical manifold method Yongtao Yang a, Guanhua Sun a,∗, Hong Zheng a,b a b
State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China Key Laboratory of Urban Security and Disaster Engineering, Ministry of Education, Beijing University of Technology, Beijing 100124, China
a r t i c l e
i n f o
Keywords: NMM Free surface SRMs (Soil-rock mixtures) Unconfined seepage problems
a b s t r a c t Recently, fruitful results of geotechnical problems solved by the NMM (abbreviated form for the numerical manifold method) have been gained. In the present work, a free surface updating strategy is implemented into the NMM to simulate unconfined seepage flow in the soil-rock mixtures (SRMs). Several examples are employed to validate the effectiveness and correctness of the proposed numerical model. Results from the numerical examples indicate that the proposed model can predict the location of free surface with high accuracy, as well as capture the main feature of the unconfined seepage flow. In addition, the influence of rock blocks on the seepage flow in the SRMs can be captured.
1. Introduction In regarding to physical and mechanical properties, there are significant differences between the soil and rock materials. The soil-rock mixtures (SRMs) [1–3] composed by these two materials hence show extremely heterogeneous character. In Chinese southwestern mountainous area, a lot of SRM slopes (see Fig. 1 for an example) are found and the failure of these slopes many result in great economic losses and casualties. Hence, evaluating the stability and proposing reinforcing measures for the weak areas of a SRM slope are of great importance. Before taking limit equilibrium methods or numerical methods for further investigation, the unconfined seepage analysis is needed since seepage flow may have adverse influence to the stability of SRM slopes. For an unconfined seepage analysis, the free surface location should be found. However, the free surface location is usually not known beforehand. Therefore, unconfined seepage analysis exhibits strong nonlinearity where an iterative process should be conducted to locate the free surface. Due to SRMs’ complex structural characteristics, it is almost impossible to locate the free surface in the unconfined seepage analysis by using the analytical methods. Therefore, numerical methods become alternative choices. In the past fifty years, a lot of numerical methods have been used to carry out seepage analysis, such as the BEM (abbreviated form for boundary element method) [4], the FEM (abbreviated form for finite element method) [5] and the MMs (meshfree methods) [6]. Although the BEM is able to reduce the problem dimension and simplify the solution, it has difficulties in solving problems with inhomogeneous domain and
∗
nonlinear characteristic [7]. The FEM [8] has been the most powerful and widely used numerical method in many types of engineering problems. The FEM has also been employed to calculate the equivalent permeability of SRMs [9]. However, the mesh used in the FEM has to match to the problem domain boundaries and material interfaces [10]. In regarding to the MMs [11,12], they are very suitable to solve problems with moving boundaries since the problems related to meshes have been overcame. However, computational cost for most of the MMs is usually very high. Apart from the MMs, the XFEM (abbreviated form for the extended finite element method) [13] and GFEM (abbreviated form for the generalized finite element method) [14] can also be employed to conduct seepage analysis. However, when dealing with hydro-mechanical coupling problems with rock blocks, the XFEM and the GFEM still have difficulties in simulating the movements of discrete rock block system. Compared to the methods mentioned above, the NMM (abbreviated form for the numerical manifold method) [15] has its unique characteristic. In the NMM, deployment of mathematical meshes is very easy, since there is no need to deploy the mathematical mesh that should match the boundaries of problem domain or interfaces of difference materials. This advantage of NMM means regular mesh, which usually results in satisfactory calculation accuracy, can always be used. In addition, the NMM has matured contact approach [16–18] to deal with contact problems which widely exist in the geotechnical engineering problems [19–22]. More importantly, the NMM can easily capture discontinuity without using the step function (or jump function). Because of the advantages mentioned above, many problems have been successfully treated using the NMM [23–50].
Corresponding author. E-mail address: [email protected] (G. Sun).
https://doi.org/10.1016/j.enganabound.2019.08.023 Received 30 May 2019; Received in revised form 11 August 2019; Accepted 13 August 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
where 𝛾 w denotes the unit weight of water, y denotes the coordinate in the vertical direction; p denotes the pore water pressure. The seepage flow velocity(v) in Ωw is assumed to satisfy Darcy’s law: 𝒗 = 𝑫𝒊
(2)
in which D denotes permeability matrix in the principal material directions, i.e., [ ] 𝑘 0 𝑫= 𝑥 (3) 0 𝑘𝑦 while i denotes hydraulic gradient, i.e., 𝒊 = −∇𝜙
(4)
in which ∇ denotes gradient operator. The continuity equation defined in the Ωw is ( ) ∇ ⋅ 𝒗 = 0, in Ω𝑤
(5)
The boundary conditions are: (1) boundary condition of water head: ( ) ̄ on Γ𝜙 = 𝐴𝐵 + 𝐶𝐷 𝜙 = 𝜙,
(6)
in which 𝜙̄ is the known water head. Fig. 1. Cross section of a SRM slope [3].
(2) boundary condition of flux: ( ) 𝑞𝑛 = −𝒏𝑇 𝒗 = 𝑞̄ = 0, on Γ𝑞 = 𝐵𝐶
A review of the literature has shown that the NMM can deal with geotechnical engineering problems with complex geometry. However, application of the NMM for unconfined seepage problems about SRMs has not been reported. In this study, a numerical model which is capable of simulating unconfined seepage flow in SRMs, is developed. A free surface updating strategy is implemented into the NMM code to obtain the free surface of SRMs. Note that there are also other methods [51–56,66] that have been implemented into the context of FEM to obtain the free surface in the unconfined seepage analysis.
Here, n and 𝑞̄ separately denote the unit normal outward vector and known flow rate.
For the convenience of description, the unconfined seepage flow in a SRM dam shown in Fig. 2 is taken as an example. The water in the problem domain (Ω = Ωw + Ωd ) is assumed to only flow within the saturated domain (Ωw ). The problem domain is divided by the free surface (Γf = AE) into two parts including the Ωw and the dry domain(Ωd ). Once the location of Γf is determined, the Ωw can also be found immediately. Therefore, the most important task in the analysis of unconfined seepage flow is the determination of free surface location. The water head 𝜙 at an arbitrary point of Ωw for a SRM dam can be expressed as: ) 𝑝 ( in Ω𝑤 𝛾𝑤 ,
(3) boundary condition corresponding to free surface: ( ) 𝑞𝑛 = 0,𝜙 = 𝑦, on Γ𝑓 = 𝐴𝐸
(8)
(4) boundary condition corresponding to seepage surface: ( ) 𝑞𝑛 ≤ 0,𝜙 = 𝑦, on Γ𝑠 = 𝐷𝐸
(9)
(5) continuity condition of material interface: { + ( ) 𝜙 = 𝜙− , on Γ𝑚 𝒏𝑇 𝒗+ = 𝒏𝑇 𝒗−
2. PDE form for the unconfined seepage problem
𝜙=𝑦+
(7)
(10)
where ‘‘+’’ and ‘‘−’’ respectively represent quantity on Γm belonging to two different material domains. 3. NMM for unconfined seepage analysis 3.1. Brief the basic concepts of NMM Since the NMM has been introduced in detail in [57], only basic concepts about NMM are introduced in the following. We take Fig. 3 as an example to explain these concepts. One salient feature of NMM is the adoption of two independent cover systems, i.e., the MC (mathematical cover) and PC (physical cover).
(1)
Fig. 2. Schematic for unconfined seepage flow in a SRM dam.
61
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 3. Physical mesh, mathematical patches, physical patches and manifold elements in numerical manifold method.
Fig. 4. Updating paths for the free surface.
In the present work, rectangular mathematical meshes are used to construct MC, although other forms such as triangular and polygonal mathematical meshes can also be used. A MP (mathematical patch) such as MP1 is the sum of several rectangles having a node in common. One MP can overlap another MP partially. And the union of all the MPs forms the mathematical cover. From Fig. 3, we can see that there is no need for the MC coinciding with the physical meshes. This merit of NMM implies that NMM could always employ regular mesh to discretize the problem domain. Physical meshes mentioned above are the sum of all
components including the problem boundaries, holes and fractures, which are capable of defining the integration areas. The PC is the sum of all the PPs(physical patches). Cutting MPs using the physical meshes will obtain the PPs. Cutting MP2 using the physical meshes will obtain PP2 ; Cutting MP3 using the physical meshes will obtain PP3 and PP4 . Note that each PP corresponding to a generalized node (GN) or NMM node. The “NMM nodes” plays similar role as the FEM nodes, since the DOFs (degrees of freedom) are attached to the NMM nodes. 62
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 5. Updating a point on the assumed free surface.
A manifold element (ME) which is the basic unit to conduct integration is the common part of four partially overlapping PPs. For instance, the manifold element E1 is the common areas of PP5 , PP6 , PP7 and PP8 . The sum of all the manifold element is exactly coincide with the domain of problem. The global displacement function uh (x) in a manifold element is expressed as 𝒖 𝒉 (𝒙 ) =
4 ∑ 𝑘=1
𝑤 𝑘 (𝒙 )𝒖 𝑘 (𝒙 )
(11)
in which wk (x) denote the PU functions (or weight functions), and uk (x) denote the cover functions (or local approximation functions). The following three conditions should be satisfied for the PU functions: 𝑤𝑖 (𝒙) = 0, 𝑖𝑓 𝒙 ∉ 𝑃 𝑃𝑖 ,
(12)
0 ≤ 𝑤𝑖 (𝒙) ≤ 1, 𝑖𝑓 𝒙 ∈ 𝑃 𝑃𝑖 ,
(13)
𝑝𝑛
𝑛 ∑ 𝑖=1
𝑤𝑖 (𝒙) = 1, 𝑖𝑓 𝒙 ∈ 𝑃 𝑃𝑖 .
(14)
where npn is the total number of PPs in the problem domain. In the present work, shape functions of the four-node isoparametric element are adopted to construct the PU functions. And constant is used to for the cover functions. The NMM global displacement function uh (x) can also be expressed as: 𝒖 𝒉 (𝒙 ) = 𝑵 𝒅
(15) Fig. 7. Several possible cases about the relative location between a free surface segment and an intersection manifold element.
where N is the matrix of shape function, while d is the vector corresponding to the DOFs of element.
Fig. 6. Three types of manifold elements.
63
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 10. Convergence of locations of the free surface from the NMM.
In the present work, the water head boundary conditions and continuity condition on the material interface are imposed using the penalty function method. The functional form presented in Eq. (16) can therefore be modified as
Fig. 8. A sketch for the homogeneous rectangular dam.
𝐼 ∗ (𝜙) = 𝐼 (𝜙) + 3.2. Weak form
1 ∇𝜙 ⋅ (𝑫 ∇(𝜙))dΩ − 𝑞̄𝜙dΓ ∬Ω 2 ∫Γ𝑞
(17)
in which k𝜙 denotes the penalty parameter related to the conditions of water head boundary, and km denotes the penalty parameter related to the continuity condition on the material interface. The following equations in matrix form can be obtained by letting 𝛿I∗ (𝜙) = 0,
According to Eq. (5) and boundary conditions for unconfined seepage problems, we can obtain the following functional form: 𝐼 (𝜙) =
1 ‖ 1 ‖ + 2 − ‖2 𝑘 𝜙 − 𝜙̄ ‖ ‖ d𝑆 + ∫ 2 𝑘𝑚 ‖𝜙 − 𝜙 ‖ d𝑆 ∫Γ𝜙 2 𝜙 ‖ Γ𝑚
(16)
𝑲𝜙 = 𝒒
(18) Fig. 9. Discretized NMM models for the homogeneous rectangular dam.
64
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 11. Water head contours for the homogeneous rectangular dam.
Fig. 12. A sketch for the rectangular dam with material interface.
in which K, q and 𝜙 denote the generalized global conductivity matrix, generalized flux term, and all degrees of freedom, respectively. The generalized global conductivity matrix K and generalized flux term q are separately generated by assembling Ke and qe , which can be obtained through: 𝑲𝑒 =
∫Ω𝑒
𝑩 𝑇 𝑫 𝑩 𝒅 Ω + 𝑘𝜙
+ 𝑘𝑚
∫Γ𝒆
Fig. 13. Discretized NMM model for the rectangular dam with two vertical layers (Model III 480 MEs, 558 PPs).
𝑵 𝑇 𝑵 d𝑆
𝜙
∫Γ𝒆
( + )𝑇 ( + ) 𝑵 − 𝑵− 𝑵 − 𝑵 − d𝑆
𝒒𝑒 =
(19)
∫Γ𝒆
𝒒
𝒎
65
𝑵 𝑇 𝑞̄𝒅 𝑆 + 𝑘𝜙
∫Γ𝒆
𝜙
̄ 𝑆 𝑵 𝑇 𝜙𝒅
(20)
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 14. Comparison of locations of the free surface assessed from different methods.
where: 𝑩 = 𝑳𝑑 𝑵 [ 𝑳𝒅 =
𝜕 ] 𝜕𝑥 𝜕 𝜕𝑦
(21)
(22)
3.3. The updating strategy for free surface Since natural boundary conditions qn = 0 corresponding to the boundaries BC and AE are imposed by not prescribing any normal flow to the surfaces in Eq. (18), the solution to the seepage problem could be obtained if the free surface AE could be determined. To obtain the free surface, the following strategy [58,59] is employed: (1) Determine the line segment that represents the potential seepage face, such as DF shown in Fig. 4. Note that the exit point of the free surface should be on DF, while the entrance point A is locating on BG. (2) Choose a point E that locates on DF as the initial guess of the exit point of the free surface. (3) Link point A and point E, and yield a line segment AE. Here, AE will be considered as the initial guess of the free surface. (4) Divide the initial free surface AE into a series of points Q1 , Q2 ,…, Qn . Determine the line segment that represents the base line, such as BC shown in Fig. 4. Then, divide the base line BC into a series of points B1 , B2 ,…, Bn . (5) Link point B1 and point Q1 , point B2 and point Q2 ,…, and point Bn and point Qn , so as to obtain a series of rays B1 Q1 , B2 Q2 ,…, Bn Qn . During the iteration steps, points Q1 , Q2 ,…, Qn which form the free surface, will move along these rays.
Fig. 15. Water head contours for the rectangular dam with two vertical layers.
is hold. Here, 𝜀 stands for a user-specified parameter, which is very small. And ( ) ( ) 𝛿𝑖 = 𝜙 𝑄𝑖 − 𝑦 𝑄𝑖 , 𝑖 = 1, 2, … , 𝑛. (24) where 𝜑(Qi ) stands for the water head at Qi , while y(Qi ) stands for the vertical coordinate of Qi . If Eq. (23) is satisfied, then the iteration process is terminated. Otherwise, point Qi will move along Bi Qi to point 𝑄̄ 𝑖 . The vertical coordinate of 𝑄̄ 𝑖 is expressed as follows (Fig. 5): ( ) ( ) 𝑦 𝑄̄ 𝑖 = 𝑦 𝑄𝑖 + 𝜔𝛿𝑖 , 𝑖 = 1, 2, … , 𝑛. (25)
Once the preparation listed above is finished, Eq. (18) is solved with the assumed free surface AE, and the water head of each NMM node is obtained. Then, we check if 𝛿 = max ||𝛿𝑖 || < 𝜀 𝑖
(23) 66
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 17. Location of the free surfaces for the rectangular dam with rock blocks.
Two discretized NMM models with different number of physical patches for this dam are adopted (Fig. 9). Free surface locations assessed from the NMM, as well as the analytical solutions are drawn together in Fig. 10. As can be seen from Fig. 10, the result predicted with the NMM based on Model II agrees well with the analytical solution, and better than that based on Model I. In addition, the water head contours predicted with the NMM based on the Model I and the Model II are drawn in Fig. 11.
Fig. 16. Discretized NMM model for the rectangular dam with rock blocks (542 MEs, 680 PPs).
Here, 𝜔 is a scaling factor which is set as 0.5 in this study. More details about 𝜔 can be found in [58]. Note that the above strategy is only applicable for updating the inner points of the free surface, such as Q1 , Q2 ,…, Qn , but cannot be used to update the location of the exit point. This is because the exit point E has been given a water head equaling to its vertical ordinate when solving Eq. (18). Therefore, the exit point 𝐸̄ is determined by extrapolating the last two points, 𝑄̄ 𝑛−1 and 𝑄̄ 𝑛 . More ingenious approaches can be found in [58]. It is noticed that during the iteration, the contribution of domain areas locating above the free surface to the generalized global conductivity matrix will not be considered. As shown in Fig. 6, there are three types of manifold elements, namely, the internal manifold elements, the external manifold elements and the intersection manifold elements. The internal manifold elements are totally under the free surface, while the external manifold elements are totally above the free surface. The intersection manifold elements are intersected by the free surface. During the calculation, the contribution of external manifold elements to the generalized global conductivity matrix will not be considered. For the intersection manifold elements, only the parts bellow the free surface will be considered. Fig. 7 shows several possible (but not limited) cases about the relative location between a free surface and an intersection manifold element.
4.2. Verification of NMM using a rectangular dam with material interface To verify the NMM for treating continuity condition for unconfined seepage problem, a dam with material interface shown in Fig. 12 is considered. In the computation, a water head with H1 = 10 is imposed on the left boundary, and that with H2 = 2 is applied to the right boundary. A discretized NMM model with 480 manifold elements and 558 physical patches is plotted in Fig. 13. Note that this example has already been used by many researchers to test their numerical methods to treat drastic singularity in the free surface [59,62–65]. Fig. 14 shows the positions of the free surface assessed from different methods. It can be seen in Fig. 14 that the result of NMM is close to the solutions of KazemzadehParsi and Daneshmand [59], and Bardet and Tobita [65]. This means the improved NMM can treat the continuity condition on the material interface for unconfined seepage problem. The water head contours predicted with the NMM is drawn in Fig. 15. 4.3. A SRM dam
4. Numerical tests
In soil-rock mixtures (SRMs), rock blocks may influence the location of free surface. Therefore, it is meaningful to conduct a test to analyze this influence. As the last example, we consider a SRM dam, as shown Fig. 16. Rock blocks with different shapes are randomly distributed in the soil matrix. The boundary conditions for this SRM dam is the same to the rectangular dam presented in Section 4.1. In the computation, the permeability coefficient of soil is set as ks = 1.0 × 10−4 , while different values of permeability coefficients for the rock blocks are considered, which are 𝑘𝑅 = 1.0 × 10−7 ∼ 1.0 × 10−4 . Fig. 17 shows the locations of free surface corresponding to different values of rock block permeability coefficients. It can be seen clearly that when the permeability coefficient of rock blocks is the same to that of soil, the free surface of the SRM dam is very smooth, which means the rock blocks have no influence on the free
In the following, the proposed numerical model will be validated by three numerical examples about unconfined seepage analysis. 4.1. Verification of NMM using a rectangular dam As the first example, a homogeneous rectangular dam is considered (Fig. 8). The height and width of the dam are 1 and 0.5, respectively. The permeability coefficient of the dam is set as 1. In the computation, a water head with H1 = 1 is imposed on the left boundary, and that with H2 = 0.5 is applied to the right boundary. The bottom boundary of the dam is impervious. The analytical solution of the location of free surface for this dam can be found in [60,61]. 67
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 18. Water head contours for the rectangular dam with rock blocks.
surface. As the permeability coefficient of rock blocks becomes smaller, the influence of rock blocks can no longer be neglected. Fig. 18 plots the water head contours of the SRM dam corresponding to different values of rock block permeability coefficients.
The proposed numerical model is validated by three numerical examples about unconfined seepage analysis. The first example is about a homogeneous rectangular dam, which is employed to test the accuracy of the proposed numerical model in predicting the free surface location. The second example is about a rectangular dam with two vertical layers, which is employed to test the ability of the proposed numerical model in treating continuity condition. The third example is about a rectangular SRM dam, which is employed to evaluate the rock-block influence on the seepage path of seepage flow.
5. Conclusions In this paper, a free surface updating strategy is implemented into the NMM to model the unconfined seepage flow in soil-rock mixtures (SRMs). 68
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
The simulation results of the first two examples are very close to the theoretical or reference solution, which show that the main features of the unconfined seepage flow problems can be captured using the proposed numerical model. The simulation results of the last example show the ability of the proposed model in predicting the free surface location of SRMs. Note that the developed model can only consider unconfined seepage flow in SRMs. In our future work, the model will be further improved to conduct unsaturated seepage analysis in SRMs.
[28] Yang YT, Zheng H. A three-node triangular element fitted to numerical manifold method with continuous nodal stress for crack analysis. Eng Fract Mech 2016;162:51–75. [29] Yang YT, Xu DD, Sun GH, Zheng H. Modeling complex crack problems using the three-node triangular element fitted to numerical manifold method with continuous nodal stress. Sci ChinaTechnol Sci 2017;60(10):1537–47. [30] Yang YT, Sun GH, Zheng H, Fu XD. A four-node quadrilateral element fitted to numerical manifold method with continuous nodal stress for crack analysis. Comput Struct 2016;177:69–82. [31] Yang YT, Tang XH, Zheng H. A three-node triangular element with continuous nodal stress. Comput Struct 2014;141:46–58. [32] Yang YT, Sun GH, Cai KJ, Zheng H. A high order numerical manifold method and its application to linear elastic continuous and fracture problems. Sci China Technol Sci 2018;61(3):346–58. [33] Yang YT, Xu DD, Zheng H. Explicit discontinuous deformation analysis method with lumped mass matrix for highly discrete block system. Int J Geomech 2018;18(9):04018098. [34] Zhang GH, Yang YT, Wang H. Improved numerical manifold method (iNMM)-An extra-DOF free and interpolating NMM with continuous nodal stress. Eng Anal Bound Elem 2017;84:117–28. [35] Liu F, Yu CY, Yang YT. An edge-based smoothed numerical manifold method and its application to static, free and forced vibration analyses. Eng Anal Bound Elem 2018;86:19–30. [36] Chen T, Yang YT, Zheng H, Wu ZJ. Numerical determination of the effective permeability coefficient of soil-rock mixtures using the numerical manifold method. Int J Numer Anal Methods Geomech 2019;43(1):381–414. [37] Zhang GH, Yang YT, Zheng H. A mass lumping scheme for the second-order numerical manifold method. Comput Struct 2019;213:23–39. [38] Yang YT, Sun GH, Zheng H. Investigation of the sequential excavation of a soil-rock-mixture slope using the numerical manifold method. Eng Geol 2019;256:93–109. [39] Zheng H, Yang YT. On generation of lumped mass matrices in partition of unity based methods. Int J Numer Methods Eng 2017;112(8):1040–69. [40] Wu ZJ, Xu XY, Liu QS, Yang YT. A zero-thickness cohesive element-based numerical manifold method for rock mechanical behavior with micro-Voronoi grains. Eng Anal Bound Elem 2018;96:94–108. [41] Zheng H, Yang YT, Shi GH. Reformulation of dynamic crack propagation using the numerical manifold method. Eng Anal Bound Elem 2019;105:279–95. [42] Jiang QH, Deng SS, Zhou CB, Lu WB. Modeling unconfined seepage flow using three-dimensional numerical manifold method. J Hydrodyn 2010;22:554–61. [43] Zheng H, Liu F, Du XL. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method. Comput Methods Appl Mech Eng 2015;295:150–71. [44] Wu Z, Wong LNY. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput Geotech 2012;39(1):38–53. [45] Hu M, Wang Y, Rutqvist J. Fully coupled hydro-mechanical numerical manifold modeling of porous rock with dominant fractures. Acta Geotech 2017;12(2):231–52. [46] Hu M, Wang Y, Rutqvist J. Development of a discontinuous approach for modeling fluid flow in heterogeneous media using the numerical manifold method. Int J Numer Anal Methods Geomech 2015;39:1932–52. [47] Hu M, Wang Y, Rutqvist J. An effective approach for modeling fluid flow in heterogeneous media using numerical manifold method. Int J Numer Methods Fluids 2015;77:459–76. [48] Wang Y, Hu MS, Zhou QL, Rutqvist J. A new second-order numerical manifold method model with an efficient scheme for analysing free surface flow with inner drains. Appl Math Model 2016;40:1427–45. [49] Wang Y, Hu MS, Zhou QL, Rutqvist J. Energy-work-based numerical manifold seepage analysis with an efficient scheme to locate the phreatic surface. Int J Numer Anal Methods Geomech 2014;38(15):1633–50. [50] Zhang GX, Li X, Li HF. Simulation of hydraulic fracture utilizing numerical manifold method. Sci China Technol Sci 2015;58(9):1542–57. [51] Bathe KJ, Khoshgoftaar MR. Finite element free surface seepage analysis without mesh iteration. Int J Numer Anal Methods Geomech 1979;3(1):13–22. [52] Desai CS, Li GC. A residual flow procedure and application for free surface in porous media. Adv Water Resour 1983;6(1):27–35. [53] Wang Y. The modified initial flow method for 3-D unconfined seepage computation. J Hydraul Eng 1998;3(3):68–73. [54] Zhou CB, Xiong WL, Liang YG. A new method for the solution unconfined seepage field. J Hydrodyn Ser A 1996;11(5):528–34. [55] Chen SH, Wang WM, She CX, et al. Unconfined seepage analysis of discontinuous rock slope. J Hydrodyn Ser B 2000;12(3):75–86. [56] Zheng H, Liu DF, Lee CF, et al. A new formulation of Signorini’s type for seepage problems with free surfaces. Int J Numer Methods Eng 2005;64(1):1–16. [57] Zheng H, Xu DD. New strategies for some issues of numerical manifold method in simulation of crack propagation. Int J Numer Methods Eng 2014;97(13):986–1010. [58] Zheng H, Liu F, Li CG. Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method. Appl Math Model 2015;39:794–808. [59] Kazemzadeh-Parsi MJ, Daneshmand F. Unconfined seepage analysis in earth dams using smoothed fixed grid finite element method. Int J Numer Anal Methods Geomech 2012;36:780–97. [60] Polubarinova-Kochina P. Theory of ground water movement. Princeton: Princeton University Press; 1962. [61] Hornung U, Krueger T. Evaluation of the Polubarinova-Kochina formula for the dam problem. Water Resour Res 1985;21:395–8. [62] Lacy SJ, Prevost JH. Flow through porous media: a procedure for locating the free surface. Int J Numer Anal Methods Geomech 1987;11:585–601.
Acknowledgements This study is supported by the National Natural Science Foundation of China, under the grant nos. 51609240, 11572009 and 51538001. References [1] Zhou Z, Fu HL, Liu BC, et al. Orthogonal tests on permeability of soil-rock-mixture. Chin J Geotech Eng 2006;28(9):1134–8. [2] Chen L, Yang YT, Zheng H. Numerical study of soil-rock mixture: generation of random aggregate structure. Sci China Technol Sci 2018;61(3):359–69. [3] Xu WJ, Wang YJ, Chen ZY, Hu RL. Stability analysis of soil-rock mixed slope based on digital image technology. Rock Soil Mech 2008;28:341–6. [4] Chugh AK, Falvey HT. Seepage analysis in a zoned anisotropic medium by the boundary element method. Int J Numer Anal Methods Geomech 1984;8:399–407. [5] Gioda G, Gentile C. A nonlinear programming analysis of unconfined steady-state seepage. Int J Numer Anal Methods Geomech 1987;11:283–305. [6] Li GX, Ge JH, Jie YX. Free surface seepage analysis based on the element-free method. Mech Res Commun 2003;30:9–19. [7] Kazemzadeh-Parsi MJ, Daneshmand F. Three dimensional smoothed fixed grid finite element method for the solution of unconfined seepage problems. Finite Elem Anal Des 2013;64:24–35. [8] Zienkiewicz OC, Taylor RL. The finite element method. 5th ed. Oxford: Butterworth-Heinemann; 2000. [9] Xu WJ, Wang YG. Meso-structural permeability of S-RM based on numerical tests. Chin J Geotech Eng 2010;32(4):542–50. [10] Lee NS, Bathe KJ. Effects of element distortions on the performance of isoparametric elements. Int J Numer Methods Eng 1993;36(20):3553–76. [11] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Methods Eng 2004;61(13):2316–43. [12] Zhuang X, Augarde CE, Mathisen KM. Fracture modeling using meshless methods and level sets in 3D: framework and modeling. Int J Numer Methods Eng 2012;92:969–98. [13] Mohammadnejad T, Khoei AR. An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model. Finite Elem Anal Des 2013;73:77–95. [14] Meschke G, Leonhart D. A generalized finite element method for hydro-mechanically coupled analysis of hydraulic fracturing problems using space-time variant enrichment functions. Comput Methods Appl Mech Eng 2015;290:438–65. [15] Shi GH. Manifold method of material analysis. In: Proceedings of the transactions of the ninth army conference on applied mathematics and computing; 1991. p. 57–76. [16] Choo LQ, Zhao Z, Chen H, Tian Q. Hydraulic fracturing modeling using the discontinuous deformation analysis (DDA) method. Comput Geotech 2016;76:12–22. [17] Jiao YY, Zhang HQ, Zhang XL, Li HB, Jiang QH. A two-dimensional coupled hydromechanical discontinuum model for simulating rock hydraulic fracturing. Int J Numer Anal Methods Geomech 2015;39(5):457–81. [18] Shi GH. Discontinuous deformation analysis—a new numerical model for the statics and dynamics of block systems. Berkeley: Department of civil engineering, University of California; 1988. [19] Zheng Y, Chen CX, Liu TT, Song DR, Meng F. Stability analysis of anti-dip bedding rock slopes locally reinforced by rock bolts. Eng Geol 2019;251:228–40. [20] Zheng Y, Chen CX, Liu TT, Zhang HN, Xia KZ, Liu F. Study on the mechanisms of flexural toppling failure in anti-inclined rock slopes using numerical and limit equilibrium models. Eng Geol 2018;237:116–28. [21] Sun GH, Yang YT, Cheng SG, et al. Phreatic line calculation and stability analysis of slopes under the combined effect of reservoir water level fluctuations and rainfall. Can Geotech J 2017;54(5):631–45. [22] Sun GH, Lin S, Jiang W, Yang YT. A simplified solution for calculating the phreatic line and slope stability during a sudden drawdown of the reservoir water level. Geofluids 2018 UNSP 1859285. doi:10.1155/2018/1859285. [23] Yang YT, Tang XH, Zheng H, Liu QS, Liu ZJ. Hydraulic fracturing modeling using the enriched numerical manifold method. Appl Math Model 2018;53:462–86. [24] Yang YT, Zheng H, Sivaselvan MV. A rigorous and unified mass lumping scheme for higher-order elements. Comput Methods Appl Mech Eng 2017;319:491–514. [25] Yang YT, Guo HW, Fu XD, Zheng H. Boundary settings for the seismic dynamic response analysis of rock masses using the numerical manifold method. Int J Numer Anal Methods Geomech 2018;42(9):1095–122. [26] Yang YT, Tang XH, Zheng H, Liu QS, He L. Three-dimensional fracture propagation with numerical manifold method. Eng Anal Bound Elem 2016;72:65–77. [27] Yang YT, Zheng H. Direct approach to treatment of contact in numerical manifold method. Int J Geomech 2017;17(5):1–14 UNSP E4016012. 69
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
[63] Oden JT, Kikuchi N. Recent advances: theory of variational inequalities with applications to problems of flow through porous media. Int J Eng Sci 1980;18:1173–284. [64] Borja RI, Kishnani SS. On the solution of elliptic free boundary problems via Newton’s method. Comput Methods Appl Mech Eng 1991;88:341–61.
[65] Bardet JP, Tobita T. A practical method for solving free-surface seepage problems. Comput Geotech 2002;29:451–75. [66] Zheng H, Yang YT, Shi GH. Reformulation of dynamic crack propagation using the numerical manifold method. Eng Anal Bound Elem 2019;105:279–95.
70
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Modeling unconfined seepage flow in soil-rock mixtures using the numerical manifold method Yongtao Yang a, Guanhua Sun a,∗, Hong Zheng a,b a b
State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China Key Laboratory of Urban Security and Disaster Engineering, Ministry of Education, Beijing University of Technology, Beijing 100124, China
a r t i c l e
i n f o
Keywords: NMM Free surface SRMs (Soil-rock mixtures) Unconfined seepage problems
a b s t r a c t Recently, fruitful results of geotechnical problems solved by the NMM (abbreviated form for the numerical manifold method) have been gained. In the present work, a free surface updating strategy is implemented into the NMM to simulate unconfined seepage flow in the soil-rock mixtures (SRMs). Several examples are employed to validate the effectiveness and correctness of the proposed numerical model. Results from the numerical examples indicate that the proposed model can predict the location of free surface with high accuracy, as well as capture the main feature of the unconfined seepage flow. In addition, the influence of rock blocks on the seepage flow in the SRMs can be captured.
1. Introduction In regarding to physical and mechanical properties, there are significant differences between the soil and rock materials. The soil-rock mixtures (SRMs) [1–3] composed by these two materials hence show extremely heterogeneous character. In Chinese southwestern mountainous area, a lot of SRM slopes (see Fig. 1 for an example) are found and the failure of these slopes many result in great economic losses and casualties. Hence, evaluating the stability and proposing reinforcing measures for the weak areas of a SRM slope are of great importance. Before taking limit equilibrium methods or numerical methods for further investigation, the unconfined seepage analysis is needed since seepage flow may have adverse influence to the stability of SRM slopes. For an unconfined seepage analysis, the free surface location should be found. However, the free surface location is usually not known beforehand. Therefore, unconfined seepage analysis exhibits strong nonlinearity where an iterative process should be conducted to locate the free surface. Due to SRMs’ complex structural characteristics, it is almost impossible to locate the free surface in the unconfined seepage analysis by using the analytical methods. Therefore, numerical methods become alternative choices. In the past fifty years, a lot of numerical methods have been used to carry out seepage analysis, such as the BEM (abbreviated form for boundary element method) [4], the FEM (abbreviated form for finite element method) [5] and the MMs (meshfree methods) [6]. Although the BEM is able to reduce the problem dimension and simplify the solution, it has difficulties in solving problems with inhomogeneous domain and
∗
nonlinear characteristic [7]. The FEM [8] has been the most powerful and widely used numerical method in many types of engineering problems. The FEM has also been employed to calculate the equivalent permeability of SRMs [9]. However, the mesh used in the FEM has to match to the problem domain boundaries and material interfaces [10]. In regarding to the MMs [11,12], they are very suitable to solve problems with moving boundaries since the problems related to meshes have been overcame. However, computational cost for most of the MMs is usually very high. Apart from the MMs, the XFEM (abbreviated form for the extended finite element method) [13] and GFEM (abbreviated form for the generalized finite element method) [14] can also be employed to conduct seepage analysis. However, when dealing with hydro-mechanical coupling problems with rock blocks, the XFEM and the GFEM still have difficulties in simulating the movements of discrete rock block system. Compared to the methods mentioned above, the NMM (abbreviated form for the numerical manifold method) [15] has its unique characteristic. In the NMM, deployment of mathematical meshes is very easy, since there is no need to deploy the mathematical mesh that should match the boundaries of problem domain or interfaces of difference materials. This advantage of NMM means regular mesh, which usually results in satisfactory calculation accuracy, can always be used. In addition, the NMM has matured contact approach [16–18] to deal with contact problems which widely exist in the geotechnical engineering problems [19–22]. More importantly, the NMM can easily capture discontinuity without using the step function (or jump function). Because of the advantages mentioned above, many problems have been successfully treated using the NMM [23–50].
Corresponding author. E-mail address: [email protected] (G. Sun).
https://doi.org/10.1016/j.enganabound.2019.08.023 Received 30 May 2019; Received in revised form 11 August 2019; Accepted 13 August 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
where 𝛾 w denotes the unit weight of water, y denotes the coordinate in the vertical direction; p denotes the pore water pressure. The seepage flow velocity(v) in Ωw is assumed to satisfy Darcy’s law: 𝒗 = 𝑫𝒊
(2)
in which D denotes permeability matrix in the principal material directions, i.e., [ ] 𝑘 0 𝑫= 𝑥 (3) 0 𝑘𝑦 while i denotes hydraulic gradient, i.e., 𝒊 = −∇𝜙
(4)
in which ∇ denotes gradient operator. The continuity equation defined in the Ωw is ( ) ∇ ⋅ 𝒗 = 0, in Ω𝑤
(5)
The boundary conditions are: (1) boundary condition of water head: ( ) ̄ on Γ𝜙 = 𝐴𝐵 + 𝐶𝐷 𝜙 = 𝜙,
(6)
in which 𝜙̄ is the known water head. Fig. 1. Cross section of a SRM slope [3].
(2) boundary condition of flux: ( ) 𝑞𝑛 = −𝒏𝑇 𝒗 = 𝑞̄ = 0, on Γ𝑞 = 𝐵𝐶
A review of the literature has shown that the NMM can deal with geotechnical engineering problems with complex geometry. However, application of the NMM for unconfined seepage problems about SRMs has not been reported. In this study, a numerical model which is capable of simulating unconfined seepage flow in SRMs, is developed. A free surface updating strategy is implemented into the NMM code to obtain the free surface of SRMs. Note that there are also other methods [51–56,66] that have been implemented into the context of FEM to obtain the free surface in the unconfined seepage analysis.
Here, n and 𝑞̄ separately denote the unit normal outward vector and known flow rate.
For the convenience of description, the unconfined seepage flow in a SRM dam shown in Fig. 2 is taken as an example. The water in the problem domain (Ω = Ωw + Ωd ) is assumed to only flow within the saturated domain (Ωw ). The problem domain is divided by the free surface (Γf = AE) into two parts including the Ωw and the dry domain(Ωd ). Once the location of Γf is determined, the Ωw can also be found immediately. Therefore, the most important task in the analysis of unconfined seepage flow is the determination of free surface location. The water head 𝜙 at an arbitrary point of Ωw for a SRM dam can be expressed as: ) 𝑝 ( in Ω𝑤 𝛾𝑤 ,
(3) boundary condition corresponding to free surface: ( ) 𝑞𝑛 = 0,𝜙 = 𝑦, on Γ𝑓 = 𝐴𝐸
(8)
(4) boundary condition corresponding to seepage surface: ( ) 𝑞𝑛 ≤ 0,𝜙 = 𝑦, on Γ𝑠 = 𝐷𝐸
(9)
(5) continuity condition of material interface: { + ( ) 𝜙 = 𝜙− , on Γ𝑚 𝒏𝑇 𝒗+ = 𝒏𝑇 𝒗−
2. PDE form for the unconfined seepage problem
𝜙=𝑦+
(7)
(10)
where ‘‘+’’ and ‘‘−’’ respectively represent quantity on Γm belonging to two different material domains. 3. NMM for unconfined seepage analysis 3.1. Brief the basic concepts of NMM Since the NMM has been introduced in detail in [57], only basic concepts about NMM are introduced in the following. We take Fig. 3 as an example to explain these concepts. One salient feature of NMM is the adoption of two independent cover systems, i.e., the MC (mathematical cover) and PC (physical cover).
(1)
Fig. 2. Schematic for unconfined seepage flow in a SRM dam.
61
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 3. Physical mesh, mathematical patches, physical patches and manifold elements in numerical manifold method.
Fig. 4. Updating paths for the free surface.
In the present work, rectangular mathematical meshes are used to construct MC, although other forms such as triangular and polygonal mathematical meshes can also be used. A MP (mathematical patch) such as MP1 is the sum of several rectangles having a node in common. One MP can overlap another MP partially. And the union of all the MPs forms the mathematical cover. From Fig. 3, we can see that there is no need for the MC coinciding with the physical meshes. This merit of NMM implies that NMM could always employ regular mesh to discretize the problem domain. Physical meshes mentioned above are the sum of all
components including the problem boundaries, holes and fractures, which are capable of defining the integration areas. The PC is the sum of all the PPs(physical patches). Cutting MPs using the physical meshes will obtain the PPs. Cutting MP2 using the physical meshes will obtain PP2 ; Cutting MP3 using the physical meshes will obtain PP3 and PP4 . Note that each PP corresponding to a generalized node (GN) or NMM node. The “NMM nodes” plays similar role as the FEM nodes, since the DOFs (degrees of freedom) are attached to the NMM nodes. 62
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 5. Updating a point on the assumed free surface.
A manifold element (ME) which is the basic unit to conduct integration is the common part of four partially overlapping PPs. For instance, the manifold element E1 is the common areas of PP5 , PP6 , PP7 and PP8 . The sum of all the manifold element is exactly coincide with the domain of problem. The global displacement function uh (x) in a manifold element is expressed as 𝒖 𝒉 (𝒙 ) =
4 ∑ 𝑘=1
𝑤 𝑘 (𝒙 )𝒖 𝑘 (𝒙 )
(11)
in which wk (x) denote the PU functions (or weight functions), and uk (x) denote the cover functions (or local approximation functions). The following three conditions should be satisfied for the PU functions: 𝑤𝑖 (𝒙) = 0, 𝑖𝑓 𝒙 ∉ 𝑃 𝑃𝑖 ,
(12)
0 ≤ 𝑤𝑖 (𝒙) ≤ 1, 𝑖𝑓 𝒙 ∈ 𝑃 𝑃𝑖 ,
(13)
𝑝𝑛
𝑛 ∑ 𝑖=1
𝑤𝑖 (𝒙) = 1, 𝑖𝑓 𝒙 ∈ 𝑃 𝑃𝑖 .
(14)
where npn is the total number of PPs in the problem domain. In the present work, shape functions of the four-node isoparametric element are adopted to construct the PU functions. And constant is used to for the cover functions. The NMM global displacement function uh (x) can also be expressed as: 𝒖 𝒉 (𝒙 ) = 𝑵 𝒅
(15) Fig. 7. Several possible cases about the relative location between a free surface segment and an intersection manifold element.
where N is the matrix of shape function, while d is the vector corresponding to the DOFs of element.
Fig. 6. Three types of manifold elements.
63
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 10. Convergence of locations of the free surface from the NMM.
In the present work, the water head boundary conditions and continuity condition on the material interface are imposed using the penalty function method. The functional form presented in Eq. (16) can therefore be modified as
Fig. 8. A sketch for the homogeneous rectangular dam.
𝐼 ∗ (𝜙) = 𝐼 (𝜙) + 3.2. Weak form
1 ∇𝜙 ⋅ (𝑫 ∇(𝜙))dΩ − 𝑞̄𝜙dΓ ∬Ω 2 ∫Γ𝑞
(17)
in which k𝜙 denotes the penalty parameter related to the conditions of water head boundary, and km denotes the penalty parameter related to the continuity condition on the material interface. The following equations in matrix form can be obtained by letting 𝛿I∗ (𝜙) = 0,
According to Eq. (5) and boundary conditions for unconfined seepage problems, we can obtain the following functional form: 𝐼 (𝜙) =
1 ‖ 1 ‖ + 2 − ‖2 𝑘 𝜙 − 𝜙̄ ‖ ‖ d𝑆 + ∫ 2 𝑘𝑚 ‖𝜙 − 𝜙 ‖ d𝑆 ∫Γ𝜙 2 𝜙 ‖ Γ𝑚
(16)
𝑲𝜙 = 𝒒
(18) Fig. 9. Discretized NMM models for the homogeneous rectangular dam.
64
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 11. Water head contours for the homogeneous rectangular dam.
Fig. 12. A sketch for the rectangular dam with material interface.
in which K, q and 𝜙 denote the generalized global conductivity matrix, generalized flux term, and all degrees of freedom, respectively. The generalized global conductivity matrix K and generalized flux term q are separately generated by assembling Ke and qe , which can be obtained through: 𝑲𝑒 =
∫Ω𝑒
𝑩 𝑇 𝑫 𝑩 𝒅 Ω + 𝑘𝜙
+ 𝑘𝑚
∫Γ𝒆
Fig. 13. Discretized NMM model for the rectangular dam with two vertical layers (Model III 480 MEs, 558 PPs).
𝑵 𝑇 𝑵 d𝑆
𝜙
∫Γ𝒆
( + )𝑇 ( + ) 𝑵 − 𝑵− 𝑵 − 𝑵 − d𝑆
𝒒𝑒 =
(19)
∫Γ𝒆
𝒒
𝒎
65
𝑵 𝑇 𝑞̄𝒅 𝑆 + 𝑘𝜙
∫Γ𝒆
𝜙
̄ 𝑆 𝑵 𝑇 𝜙𝒅
(20)
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 14. Comparison of locations of the free surface assessed from different methods.
where: 𝑩 = 𝑳𝑑 𝑵 [ 𝑳𝒅 =
𝜕 ] 𝜕𝑥 𝜕 𝜕𝑦
(21)
(22)
3.3. The updating strategy for free surface Since natural boundary conditions qn = 0 corresponding to the boundaries BC and AE are imposed by not prescribing any normal flow to the surfaces in Eq. (18), the solution to the seepage problem could be obtained if the free surface AE could be determined. To obtain the free surface, the following strategy [58,59] is employed: (1) Determine the line segment that represents the potential seepage face, such as DF shown in Fig. 4. Note that the exit point of the free surface should be on DF, while the entrance point A is locating on BG. (2) Choose a point E that locates on DF as the initial guess of the exit point of the free surface. (3) Link point A and point E, and yield a line segment AE. Here, AE will be considered as the initial guess of the free surface. (4) Divide the initial free surface AE into a series of points Q1 , Q2 ,…, Qn . Determine the line segment that represents the base line, such as BC shown in Fig. 4. Then, divide the base line BC into a series of points B1 , B2 ,…, Bn . (5) Link point B1 and point Q1 , point B2 and point Q2 ,…, and point Bn and point Qn , so as to obtain a series of rays B1 Q1 , B2 Q2 ,…, Bn Qn . During the iteration steps, points Q1 , Q2 ,…, Qn which form the free surface, will move along these rays.
Fig. 15. Water head contours for the rectangular dam with two vertical layers.
is hold. Here, 𝜀 stands for a user-specified parameter, which is very small. And ( ) ( ) 𝛿𝑖 = 𝜙 𝑄𝑖 − 𝑦 𝑄𝑖 , 𝑖 = 1, 2, … , 𝑛. (24) where 𝜑(Qi ) stands for the water head at Qi , while y(Qi ) stands for the vertical coordinate of Qi . If Eq. (23) is satisfied, then the iteration process is terminated. Otherwise, point Qi will move along Bi Qi to point 𝑄̄ 𝑖 . The vertical coordinate of 𝑄̄ 𝑖 is expressed as follows (Fig. 5): ( ) ( ) 𝑦 𝑄̄ 𝑖 = 𝑦 𝑄𝑖 + 𝜔𝛿𝑖 , 𝑖 = 1, 2, … , 𝑛. (25)
Once the preparation listed above is finished, Eq. (18) is solved with the assumed free surface AE, and the water head of each NMM node is obtained. Then, we check if 𝛿 = max ||𝛿𝑖 || < 𝜀 𝑖
(23) 66
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 17. Location of the free surfaces for the rectangular dam with rock blocks.
Two discretized NMM models with different number of physical patches for this dam are adopted (Fig. 9). Free surface locations assessed from the NMM, as well as the analytical solutions are drawn together in Fig. 10. As can be seen from Fig. 10, the result predicted with the NMM based on Model II agrees well with the analytical solution, and better than that based on Model I. In addition, the water head contours predicted with the NMM based on the Model I and the Model II are drawn in Fig. 11.
Fig. 16. Discretized NMM model for the rectangular dam with rock blocks (542 MEs, 680 PPs).
Here, 𝜔 is a scaling factor which is set as 0.5 in this study. More details about 𝜔 can be found in [58]. Note that the above strategy is only applicable for updating the inner points of the free surface, such as Q1 , Q2 ,…, Qn , but cannot be used to update the location of the exit point. This is because the exit point E has been given a water head equaling to its vertical ordinate when solving Eq. (18). Therefore, the exit point 𝐸̄ is determined by extrapolating the last two points, 𝑄̄ 𝑛−1 and 𝑄̄ 𝑛 . More ingenious approaches can be found in [58]. It is noticed that during the iteration, the contribution of domain areas locating above the free surface to the generalized global conductivity matrix will not be considered. As shown in Fig. 6, there are three types of manifold elements, namely, the internal manifold elements, the external manifold elements and the intersection manifold elements. The internal manifold elements are totally under the free surface, while the external manifold elements are totally above the free surface. The intersection manifold elements are intersected by the free surface. During the calculation, the contribution of external manifold elements to the generalized global conductivity matrix will not be considered. For the intersection manifold elements, only the parts bellow the free surface will be considered. Fig. 7 shows several possible (but not limited) cases about the relative location between a free surface and an intersection manifold element.
4.2. Verification of NMM using a rectangular dam with material interface To verify the NMM for treating continuity condition for unconfined seepage problem, a dam with material interface shown in Fig. 12 is considered. In the computation, a water head with H1 = 10 is imposed on the left boundary, and that with H2 = 2 is applied to the right boundary. A discretized NMM model with 480 manifold elements and 558 physical patches is plotted in Fig. 13. Note that this example has already been used by many researchers to test their numerical methods to treat drastic singularity in the free surface [59,62–65]. Fig. 14 shows the positions of the free surface assessed from different methods. It can be seen in Fig. 14 that the result of NMM is close to the solutions of KazemzadehParsi and Daneshmand [59], and Bardet and Tobita [65]. This means the improved NMM can treat the continuity condition on the material interface for unconfined seepage problem. The water head contours predicted with the NMM is drawn in Fig. 15. 4.3. A SRM dam
4. Numerical tests
In soil-rock mixtures (SRMs), rock blocks may influence the location of free surface. Therefore, it is meaningful to conduct a test to analyze this influence. As the last example, we consider a SRM dam, as shown Fig. 16. Rock blocks with different shapes are randomly distributed in the soil matrix. The boundary conditions for this SRM dam is the same to the rectangular dam presented in Section 4.1. In the computation, the permeability coefficient of soil is set as ks = 1.0 × 10−4 , while different values of permeability coefficients for the rock blocks are considered, which are 𝑘𝑅 = 1.0 × 10−7 ∼ 1.0 × 10−4 . Fig. 17 shows the locations of free surface corresponding to different values of rock block permeability coefficients. It can be seen clearly that when the permeability coefficient of rock blocks is the same to that of soil, the free surface of the SRM dam is very smooth, which means the rock blocks have no influence on the free
In the following, the proposed numerical model will be validated by three numerical examples about unconfined seepage analysis. 4.1. Verification of NMM using a rectangular dam As the first example, a homogeneous rectangular dam is considered (Fig. 8). The height and width of the dam are 1 and 0.5, respectively. The permeability coefficient of the dam is set as 1. In the computation, a water head with H1 = 1 is imposed on the left boundary, and that with H2 = 0.5 is applied to the right boundary. The bottom boundary of the dam is impervious. The analytical solution of the location of free surface for this dam can be found in [60,61]. 67
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
Fig. 18. Water head contours for the rectangular dam with rock blocks.
surface. As the permeability coefficient of rock blocks becomes smaller, the influence of rock blocks can no longer be neglected. Fig. 18 plots the water head contours of the SRM dam corresponding to different values of rock block permeability coefficients.
The proposed numerical model is validated by three numerical examples about unconfined seepage analysis. The first example is about a homogeneous rectangular dam, which is employed to test the accuracy of the proposed numerical model in predicting the free surface location. The second example is about a rectangular dam with two vertical layers, which is employed to test the ability of the proposed numerical model in treating continuity condition. The third example is about a rectangular SRM dam, which is employed to evaluate the rock-block influence on the seepage path of seepage flow.
5. Conclusions In this paper, a free surface updating strategy is implemented into the NMM to model the unconfined seepage flow in soil-rock mixtures (SRMs). 68
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
The simulation results of the first two examples are very close to the theoretical or reference solution, which show that the main features of the unconfined seepage flow problems can be captured using the proposed numerical model. The simulation results of the last example show the ability of the proposed model in predicting the free surface location of SRMs. Note that the developed model can only consider unconfined seepage flow in SRMs. In our future work, the model will be further improved to conduct unsaturated seepage analysis in SRMs.
[28] Yang YT, Zheng H. A three-node triangular element fitted to numerical manifold method with continuous nodal stress for crack analysis. Eng Fract Mech 2016;162:51–75. [29] Yang YT, Xu DD, Sun GH, Zheng H. Modeling complex crack problems using the three-node triangular element fitted to numerical manifold method with continuous nodal stress. Sci ChinaTechnol Sci 2017;60(10):1537–47. [30] Yang YT, Sun GH, Zheng H, Fu XD. A four-node quadrilateral element fitted to numerical manifold method with continuous nodal stress for crack analysis. Comput Struct 2016;177:69–82. [31] Yang YT, Tang XH, Zheng H. A three-node triangular element with continuous nodal stress. Comput Struct 2014;141:46–58. [32] Yang YT, Sun GH, Cai KJ, Zheng H. A high order numerical manifold method and its application to linear elastic continuous and fracture problems. Sci China Technol Sci 2018;61(3):346–58. [33] Yang YT, Xu DD, Zheng H. Explicit discontinuous deformation analysis method with lumped mass matrix for highly discrete block system. Int J Geomech 2018;18(9):04018098. [34] Zhang GH, Yang YT, Wang H. Improved numerical manifold method (iNMM)-An extra-DOF free and interpolating NMM with continuous nodal stress. Eng Anal Bound Elem 2017;84:117–28. [35] Liu F, Yu CY, Yang YT. An edge-based smoothed numerical manifold method and its application to static, free and forced vibration analyses. Eng Anal Bound Elem 2018;86:19–30. [36] Chen T, Yang YT, Zheng H, Wu ZJ. Numerical determination of the effective permeability coefficient of soil-rock mixtures using the numerical manifold method. Int J Numer Anal Methods Geomech 2019;43(1):381–414. [37] Zhang GH, Yang YT, Zheng H. A mass lumping scheme for the second-order numerical manifold method. Comput Struct 2019;213:23–39. [38] Yang YT, Sun GH, Zheng H. Investigation of the sequential excavation of a soil-rock-mixture slope using the numerical manifold method. Eng Geol 2019;256:93–109. [39] Zheng H, Yang YT. On generation of lumped mass matrices in partition of unity based methods. Int J Numer Methods Eng 2017;112(8):1040–69. [40] Wu ZJ, Xu XY, Liu QS, Yang YT. A zero-thickness cohesive element-based numerical manifold method for rock mechanical behavior with micro-Voronoi grains. Eng Anal Bound Elem 2018;96:94–108. [41] Zheng H, Yang YT, Shi GH. Reformulation of dynamic crack propagation using the numerical manifold method. Eng Anal Bound Elem 2019;105:279–95. [42] Jiang QH, Deng SS, Zhou CB, Lu WB. Modeling unconfined seepage flow using three-dimensional numerical manifold method. J Hydrodyn 2010;22:554–61. [43] Zheng H, Liu F, Du XL. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method. Comput Methods Appl Mech Eng 2015;295:150–71. [44] Wu Z, Wong LNY. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput Geotech 2012;39(1):38–53. [45] Hu M, Wang Y, Rutqvist J. Fully coupled hydro-mechanical numerical manifold modeling of porous rock with dominant fractures. Acta Geotech 2017;12(2):231–52. [46] Hu M, Wang Y, Rutqvist J. Development of a discontinuous approach for modeling fluid flow in heterogeneous media using the numerical manifold method. Int J Numer Anal Methods Geomech 2015;39:1932–52. [47] Hu M, Wang Y, Rutqvist J. An effective approach for modeling fluid flow in heterogeneous media using numerical manifold method. Int J Numer Methods Fluids 2015;77:459–76. [48] Wang Y, Hu MS, Zhou QL, Rutqvist J. A new second-order numerical manifold method model with an efficient scheme for analysing free surface flow with inner drains. Appl Math Model 2016;40:1427–45. [49] Wang Y, Hu MS, Zhou QL, Rutqvist J. Energy-work-based numerical manifold seepage analysis with an efficient scheme to locate the phreatic surface. Int J Numer Anal Methods Geomech 2014;38(15):1633–50. [50] Zhang GX, Li X, Li HF. Simulation of hydraulic fracture utilizing numerical manifold method. Sci China Technol Sci 2015;58(9):1542–57. [51] Bathe KJ, Khoshgoftaar MR. Finite element free surface seepage analysis without mesh iteration. Int J Numer Anal Methods Geomech 1979;3(1):13–22. [52] Desai CS, Li GC. A residual flow procedure and application for free surface in porous media. Adv Water Resour 1983;6(1):27–35. [53] Wang Y. The modified initial flow method for 3-D unconfined seepage computation. J Hydraul Eng 1998;3(3):68–73. [54] Zhou CB, Xiong WL, Liang YG. A new method for the solution unconfined seepage field. J Hydrodyn Ser A 1996;11(5):528–34. [55] Chen SH, Wang WM, She CX, et al. Unconfined seepage analysis of discontinuous rock slope. J Hydrodyn Ser B 2000;12(3):75–86. [56] Zheng H, Liu DF, Lee CF, et al. A new formulation of Signorini’s type for seepage problems with free surfaces. Int J Numer Methods Eng 2005;64(1):1–16. [57] Zheng H, Xu DD. New strategies for some issues of numerical manifold method in simulation of crack propagation. Int J Numer Methods Eng 2014;97(13):986–1010. [58] Zheng H, Liu F, Li CG. Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method. Appl Math Model 2015;39:794–808. [59] Kazemzadeh-Parsi MJ, Daneshmand F. Unconfined seepage analysis in earth dams using smoothed fixed grid finite element method. Int J Numer Anal Methods Geomech 2012;36:780–97. [60] Polubarinova-Kochina P. Theory of ground water movement. Princeton: Princeton University Press; 1962. [61] Hornung U, Krueger T. Evaluation of the Polubarinova-Kochina formula for the dam problem. Water Resour Res 1985;21:395–8. [62] Lacy SJ, Prevost JH. Flow through porous media: a procedure for locating the free surface. Int J Numer Anal Methods Geomech 1987;11:585–601.
Acknowledgements This study is supported by the National Natural Science Foundation of China, under the grant nos. 51609240, 11572009 and 51538001. References [1] Zhou Z, Fu HL, Liu BC, et al. Orthogonal tests on permeability of soil-rock-mixture. Chin J Geotech Eng 2006;28(9):1134–8. [2] Chen L, Yang YT, Zheng H. Numerical study of soil-rock mixture: generation of random aggregate structure. Sci China Technol Sci 2018;61(3):359–69. [3] Xu WJ, Wang YJ, Chen ZY, Hu RL. Stability analysis of soil-rock mixed slope based on digital image technology. Rock Soil Mech 2008;28:341–6. [4] Chugh AK, Falvey HT. Seepage analysis in a zoned anisotropic medium by the boundary element method. Int J Numer Anal Methods Geomech 1984;8:399–407. [5] Gioda G, Gentile C. A nonlinear programming analysis of unconfined steady-state seepage. Int J Numer Anal Methods Geomech 1987;11:283–305. [6] Li GX, Ge JH, Jie YX. Free surface seepage analysis based on the element-free method. Mech Res Commun 2003;30:9–19. [7] Kazemzadeh-Parsi MJ, Daneshmand F. Three dimensional smoothed fixed grid finite element method for the solution of unconfined seepage problems. Finite Elem Anal Des 2013;64:24–35. [8] Zienkiewicz OC, Taylor RL. The finite element method. 5th ed. Oxford: Butterworth-Heinemann; 2000. [9] Xu WJ, Wang YG. Meso-structural permeability of S-RM based on numerical tests. Chin J Geotech Eng 2010;32(4):542–50. [10] Lee NS, Bathe KJ. Effects of element distortions on the performance of isoparametric elements. Int J Numer Methods Eng 1993;36(20):3553–76. [11] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Methods Eng 2004;61(13):2316–43. [12] Zhuang X, Augarde CE, Mathisen KM. Fracture modeling using meshless methods and level sets in 3D: framework and modeling. Int J Numer Methods Eng 2012;92:969–98. [13] Mohammadnejad T, Khoei AR. An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model. Finite Elem Anal Des 2013;73:77–95. [14] Meschke G, Leonhart D. A generalized finite element method for hydro-mechanically coupled analysis of hydraulic fracturing problems using space-time variant enrichment functions. Comput Methods Appl Mech Eng 2015;290:438–65. [15] Shi GH. Manifold method of material analysis. In: Proceedings of the transactions of the ninth army conference on applied mathematics and computing; 1991. p. 57–76. [16] Choo LQ, Zhao Z, Chen H, Tian Q. Hydraulic fracturing modeling using the discontinuous deformation analysis (DDA) method. Comput Geotech 2016;76:12–22. [17] Jiao YY, Zhang HQ, Zhang XL, Li HB, Jiang QH. A two-dimensional coupled hydromechanical discontinuum model for simulating rock hydraulic fracturing. Int J Numer Anal Methods Geomech 2015;39(5):457–81. [18] Shi GH. Discontinuous deformation analysis—a new numerical model for the statics and dynamics of block systems. Berkeley: Department of civil engineering, University of California; 1988. [19] Zheng Y, Chen CX, Liu TT, Song DR, Meng F. Stability analysis of anti-dip bedding rock slopes locally reinforced by rock bolts. Eng Geol 2019;251:228–40. [20] Zheng Y, Chen CX, Liu TT, Zhang HN, Xia KZ, Liu F. Study on the mechanisms of flexural toppling failure in anti-inclined rock slopes using numerical and limit equilibrium models. Eng Geol 2018;237:116–28. [21] Sun GH, Yang YT, Cheng SG, et al. Phreatic line calculation and stability analysis of slopes under the combined effect of reservoir water level fluctuations and rainfall. Can Geotech J 2017;54(5):631–45. [22] Sun GH, Lin S, Jiang W, Yang YT. A simplified solution for calculating the phreatic line and slope stability during a sudden drawdown of the reservoir water level. Geofluids 2018 UNSP 1859285. doi:10.1155/2018/1859285. [23] Yang YT, Tang XH, Zheng H, Liu QS, Liu ZJ. Hydraulic fracturing modeling using the enriched numerical manifold method. Appl Math Model 2018;53:462–86. [24] Yang YT, Zheng H, Sivaselvan MV. A rigorous and unified mass lumping scheme for higher-order elements. Comput Methods Appl Mech Eng 2017;319:491–514. [25] Yang YT, Guo HW, Fu XD, Zheng H. Boundary settings for the seismic dynamic response analysis of rock masses using the numerical manifold method. Int J Numer Anal Methods Geomech 2018;42(9):1095–122. [26] Yang YT, Tang XH, Zheng H, Liu QS, He L. Three-dimensional fracture propagation with numerical manifold method. Eng Anal Bound Elem 2016;72:65–77. [27] Yang YT, Zheng H. Direct approach to treatment of contact in numerical manifold method. Int J Geomech 2017;17(5):1–14 UNSP E4016012. 69
Y. Yang, G. Sun and H. Zheng
Engineering Analysis with Boundary Elements 108 (2019) 60–70
[63] Oden JT, Kikuchi N. Recent advances: theory of variational inequalities with applications to problems of flow through porous media. Int J Eng Sci 1980;18:1173–284. [64] Borja RI, Kishnani SS. On the solution of elliptic free boundary problems via Newton’s method. Comput Methods Appl Mech Eng 1991;88:341–61.
[65] Bardet JP, Tobita T. A practical method for solving free-surface seepage problems. Comput Geotech 2002;29:451–75. [66] Zheng H, Yang YT, Shi GH. Reformulation of dynamic crack propagation using the numerical manifold method. Eng Anal Bound Elem 2019;105:279–95.
70