Hydraulic fracture at the dam-foundation interface using the scaled boundary finite element method coupled with the cohesive crack model

Engineering Analysis with Boundary Elements 88 (2018) 41–53

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Hydraulic fracture at the dam-foundation interface using the scaled boundary finite element method coupled with the cohesive crack model Hong Zhong a,b,∗, Hongjun Li a, Ean Tat Ooi c, Chongmin Song d a

State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China c School of Science, Information Technology and Engineering, Federation University, Ballarat VIC 3353, Australia d School of Civil and Environmental Engineering, The University of New South Wales, Sydney NSW 2052, Australia b

a r t i c l e

i n f o

Keywords: Scaled boundary finite element method Polygon element Cohesive crack model Hydraulic fracture Interfacial fracture Crack propagation

a b s t r a c t The scaled boundary finite element method coupled with the cohesive crack model is extended to investigate the hydraulic fracture at the dam-foundation interface. The concrete and rock bulk are modeled by the scaled boundary polygons. Cohesive interface elements model the fracture process zone between the crack faces. The cohesive tractions are modeled as side-face tractions in the scaled boundary polygons. The solution of the stress field around the crack tip is expressed semi-analytically as a power series. Accurate displacement field, stress field and stress intensity factors can be obtained without asymptotic enrichment or local mesh refinement. The proposed procedure is verified by the hydraulic fracture of a rectangular embankment on rigid foundation and applied to the modeling of hydraulic fracture on the dam-foundation interface of a benchmark dam. Different distributions of water pressure inside the crack are investigated. It is found that the water pressure inside the crack decreases the peak overflow to less than 20% of the case without water in the crack. Considering the water lag or not is significant to the response, while different distribution of pressure following the water lag region in the fracture process zone has negligible influence. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction The interface between concrete dam and its underlying rock foundation is a weak link of the dam-foundation system. The stress concentration at the interface is a source for nucleation of micro-cracks. In practice, micro-cracks are inherently present at the dam heel. These cracks can potentially coalesce to develop a major crack that propagates along the concrete-rock interface, which increases the permeability of concrete, harms the drainage system and reduces the shear bearing capacity of the dam. Moreover, water inflow from the reservoir into the crack faces further promotes the crack propagation through hydraulic mechanisms. Rescher [1] reported that the penetration of water into the concrete-rock interface exerts an uplift pressure on the dam, which affects its structural stability. Therefore, hydraulic fracture at the damfoundation interface is an important phenomenon for the safety evaluation of concrete dams. The distribution of the water pressure within the crack is an important parameter in hydraulic fracture modeling. Generally, two approaches are usually used to determine the water pressure distribution within the crack during hydraulic fracture. The first approach employs ∗

Corresponding author. E-mail address: [email protected] (H. Zhong).

coupled hydro-mechanical simulations, in which, the water pressure distribution is determined from the equilibrium conditions of fluid flow in a control volume [2]. The second approach idealizes the water pressure distribution as a purely mechanical load [3,4]. A mathematical model for the water pressure distribution is usually derived analytically based on experimental observations. This approach is suitable for crack propagation in materials with low permeability. The latter approach is adopted in this paper. In treating the water pressure as a purely mechanical load, initial studies using this approach assumed a uniform distribution of the water pressure within the crack [5]. Independent numerical studies by Dewey et al. [6] and Plizzari [7], however, reported that the stress field in the vicinity of the crack is sensitive to the assumed water pressure distribution within the crack. Bruhwiler and Saouma [8] studied the distribution of water pressure within a crack by experimenting on wedge splitting specimens. They reported that the water pressure increased from zero to full hydrostatic pressure along the fracture process zone. An empirical model for the water pressure distribution expressed as a function of the crack opening displacement was proposed. Slowik and Saouma [9] studied the effect of water pressure in fast propagating cracks using wedge splitting tests subjected to dynamic loads. They reported that the rate of the COD influences the water pressure distribution inside the

https://doi.org/10.1016/j.enganabound.2017.11.009 Received 19 September 2017; Received in revised form 24 November 2017; Accepted 25 November 2017 0955-7997/© 2017 Elsevier Ltd. All rights reserved.

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

crack and that the trapped water has a wedging effect that increases the magnitude of the tensile stresses in the vicinity of the crack tip. A model for the fluid flow in concrete cracks was proposed based on experimental observations and simplifying assumptions based on, for example, the fluid flow velocity, the size of crack widths and the impermeability of the concrete at the fracture process zone. Hydraulic fracture at the dam-joint foundation is usually investigated using the cohesive zone model [10–12] in combination with the finite element method (FEM). Bolzon and Cocchetti [3] developed a friction-cohesive softening interface model with coupled degradation of normal and tangential strength to model the joint in low permeability materials. Two models of water pressure distribution within the crack was studied i.e. piecewise linear and exponential. Alfano et al. [4] developed a multi-scale interface model to simulate crack propagation in dams considering the effect of water pressure in the crack. The water pressure distribution was assumed as an exponential function of the crack opening displacement, consistent with the experimental observations of Bruhwiler and Saouma [8]. Barpi and Valente [13] modeled the water penetration at a dam-foundation joint using the cohesive zone model. The model for the water pressure distribution inside the crack was described by polynomial functions and was based on the models proposed in [8,14,15]. Barpi and Valente [16] also derived the asymptotic expansion of the stress field around the crack tip of a damfoundation joint. Consistent with the model developed by Desroches et al. [17], their analysis assumed water lag near the crack tip where the hydrostatic pressure penetrates only to the knee-point of the cohesive law. Using a similar approach, Alberto and Valente [18] derived the asymptotic fields around a cohesive frictional crack at the bi-material interface between a dam and the rock foundation. Efficient simulations of hydraulic fracture at the dam-foundation interface require accurate modeling of the stress field in the vicinity of the crack. Compared with homogeneous materials, the stress field at the crack on a concrete-rock interface is more complex and is highly influenced by the mismatch of the two constituent materials. Although the FEM can be used to model such behavior, its polynomial-based shape functions necessitate the use of a sufficiently dense mesh in the vicinity of the crack tip, as was demonstrated in the numerical examples reported in e.g. [4,13]. To improve the performance of the FEM in modeling quasi-brittle fracture, the finite element shape functions can be enriched with the asymptotic expansion of a cohesive crack [19,20]. The scaled boundary finite element method (SBFEM) is an alternative approach that can effectively model the fracture behavior of quasibrittle materials. The SBFEM is a semi-analytical method and was developed by Song and Wolf [21]. Its unique characteristics make it very efficient in modeling problems of unbounded media [22,23], thin plates [24], heat problem [25] and fracture problems [26–28]. Application of the SBFEM to model quasi-brittle fracture was reported in [29,30,33]. When modeling fracture, the SBFEM does not require a fine mesh at the crack tip such as the FEM e.g. [13] [34] or asymptotic enrichment functions such as the extended FEM [19,35]. The asymptotic fields in the vicinity of crack tips, notches, material junctions and cohesive cracks are analytically represented in the formulation. In this paper, the SBFEM is applied to model hydraulic fracture at the dam-foundation interface. The concrete and rock are assumed to be linear elastic and impermeable. The fracture process zone and the water pressure distribution in the vicinity of the crack tip is modeled using the cohesive zone model via interface elements. The water pressure distribution is applied to the crack surfaces as a purely mechanical load and is assumed to be an exponential function of the crack opening displacement, similar to [3,4]. The interface elements are coupled to the SBFEM by a shadow domain procedure [29]. The SBFEM solution using generalized coordinates [36] is adopted to model the cohesive tractions acting on the crack faces. To facilitate the modeling of crack propagation, the SBFEM is used together with polygon meshes having arbitrary number of sides [44]. Crack propagation is modeled by a simple modification of the polygons along the interface.

This paper is organized as follows. Section 2 describes the cohesive zone model and interface elements used to model the cohesive tractions and water pressure distribution acting in the fracture process zone of the interface crack. Section 3 details the fundamental theory of the SBFEM and its solution using generalized coordinates. Section 4 describes the shadow domain procedure used to couple the SBFEM and the interface elements, the stress field in the vicinity of the crack tip and the condition of stability of the cohesive crack. In Section 5, the developed approach is validated by standard numerical bench marks. The influence of the water pressure on the fracture behavior at the interface crack is investigated numerically. The major conclusions of this study are summarized in Section 6. 2. Modeling of cohesive crack and water pressure in the process 2.1. Cohesive traction distribution The cohesive zone model [10–12,49] is commonly used to model the fracture behavior of quasi-brittle materials such as concrete, rocks and ceramics. This model idealizes the fracture behavior in the process zone as a line that is characterized by nonlinear constitutive relations of the cohesive traction and the relative displacements of the crack surfaces. The constitutive relation at the fracture process zone is determined by the mode I and mode II fracture energy, GfI and GfII , and the tensile and shear strength, tn (u) and ts (u) . The exact relation between the normal and tangential cohesive tractions, tn and ts , with the relative crack opening and sliding displacements, wn and ws , are not known explicitly. For homogeneous materials, the experimentally measured GfI , GfII , tn and ts are fitted into various simplified models assuming linear, bi-linear, exponential and trapezoidal decaying relations have been proposed and studied e.g [37]. For bi-material interfacial fracture, there is a lack of consensus on the form of the constitutive relation at the interface and its related fracture parameters. For example, the determination of GfII , the coupling between the constitutive relation in the normal direction and that in the shear direction, still remain to be further investigated. Mahalingam presented a cohesive law for the NFU-SiN interface, in which the parameters for the normal component and shear component were obtained through independent experiments. The shape of the traction-displacement law for both normal and tangential components was assumed to be the same [38]. Toftegaard et al. studied the interface between two dissimilar composite materials by mixed mode cohesive laws. They considered linear and exponential softening cohesive laws for both the normal and shear components, and a parametric analysis was provided [39]. The authors previously investigated the quasi-brittle fracture behavior at the interface between concrete and rock using four-point-shear specimens [40]. By assuming a set of unique parameters, the constitutive relations relating the normal and tangential cohesive traction to the relative displacements of the crack faces at the interface were determined by calibrating the experimentally measured load-crack opening and sliding displacements with numerical simulations. In order to be able to use the interface constitutive relations previously determined in [40] to model interface fracture of practical dams, the experimentally determined relations between the cohesive traction and the crack opening and sliding displacements, and the fracture energies, GfI and GfII , have to be augmented to account for the size effect [41]. Due to the lack of experimental results on specimens with aggregate sizes typical in practical dams, the tails of the cohesive curves determined in [40] are elongated so that the energy dissipation in the process zone due to larger crack opening and shearing can be modeled. The augmented constitutive relations of the concrete-rock interface are shown in Fig. 1. Note that for the shear constitutive relation, only the positive part of the anti-symmetric curve is shown. In these curves, the knee points remain unchanged. For the case of shear case, since there is no knee point, the first half of the descending section remains unchanged. The size of the fracture process zone, wn and ws , is assumed 42

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

Fig. 1. Bi-linear cohesive constitutive relation.

Table 1 Parameters for cohesive traction law. t n −w n

GfI (N/m) 207.00

tn (u) (MPa) 2.00

tn1 (u) (MPa) 0.67

wn0 (mm) 0.001

wn1 (mm) 0.040

wnc (mm) 0.500

t s −w s

GfII (N/m) 425.00

ts (u) (MPa) 3.50

ts1 (u) (MPa) 1.17

ws0 (mm) 0.001

ws1 (mm) 0.076

wsc (mm) 0.500

to be 0.5 mm, which is commonly used in cohesive fracture modeling of concrete dams. It is to be noted that both the fracture energies GfI and GfII have been increased as a result of augmenting the constitutive relations at the interface. The fracture parameters defining the constitutive relations shown in Fig. 1 are listed in Table 1. In addition to this bi-linear cohesive law, in order to verify the proposed procedure, an exponential cohesive law is also assumed in the first numerical example following that in the literature. The cohesive crack model is included through the cohesive interface element (CIE) [29,30] which embodies the cohesive laws. Then the stiffness matrices of CIEs are assembled with those of the bulk material.

Fig. 2. Combination of normal traction and water pressure.

crack. As 𝜌 decreases, the hydrostatic pressure decreases as it approaches the crack tip. For 𝜌 = 0, the water pressure is zero, corresponding to a perfect drainage system. In this equation, the water pressure distribution is controlled by only a single parameter, i.e. 𝜌. This enables convenient study on the effect of water pressure on response of the structure.

2.2. Water pressure distribution When a crack propagates along the dam-foundation joint, water from the reservoir penetrates into the crack, creating a wedging effect which further promotes crack propagation. Therefore, the water pressure distribution within the crack has to be taken into account in order to accurately model the structural response. The interaction between the water pressure and the structural response at the interface is an inherently complex phenomenon involving hydro-mechanical coupling between the fluid and the surrounding solid. To simplify the modeling process, it is assumed that both concrete and rock are impermeable and that the loading is quasi-static. With these assumptions, the effect of fluid motion can be neglected and the water pressure distribution in the crack is a function of the crack opening displacement that increases from zero to full hydrostatic pressure along the fracture process zone [8]. In this study, the water pressure distribution p is assumed to be related to the crack opening displacement by an exponential function as ( ) 𝑝 = 𝑝0 1 − 𝑒−𝜌(𝑤𝑛 −𝑤𝑛1 )

2.3. Combined effect of water pressure and cohesive tractions on crack face In hydraulic fracture, the load on the crack face consists of two parts: (1) the cohesive traction due to aggregate interlocking and surface friction at the process zone and (2) the hydrostatic pressure. Two approaches can be used to model the combined effect of these loads on the crack faces. The first approach is based on experimental observations that both the cohesive tractions and hydrostatic pressure distribution in the crack face are functions of the crack opening displacements. A combined relation between the cohesive tractions, hydrostatic pressure and the crack opening displacement is obtained by superposition. The resulting constitutive relation between the cohesive traction, hydrostatic pressure and crack opening displacement is shown in Fig. 2 for a bilinear softening law. For a given crack opening wn , the secant normal stiffness of the interface element is decided by the corresponding value of (𝜎-p) divided by wn . When unloading takes place, only the normal cohesive traction 𝜎 is unloaded, the water pressure p is only dependent on the current wn . In the second approach, water pressure is applied on the crack face as an external distributed load. It varies with the opening of the crack, so during the nonlinear iteration process the water pressure is updated for each iteration.

(1)

where p0 denotes the full hydrostatic pressure, wn is the crack opening displacement, wn1 is the crack opening corresponding to the knee point on the normal traction curve as shown in Fig. 1. When wn is smaller than wn1 , p can be set to zero to consider the water lag near the crack tip. 𝜌 is a parameter that modulates the pressure distribution. When 𝜌 is sufficiently large, p approaches p0 , the pressure is independent of the crack opening and the water pressure is uniformly distributed inside the 43

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

The first approach is easy to implement since only the normal cohesive traction law needs to be updated with a combined one that included the influence of water pressure. Its disadvantage is that negative stiffness of the cohesive interface elements may be involved, which could influence the convergence of nonlinear iteration. The second approach is theoretically clearer, but the iteration between water pressure and crack opening displacement has to be explicitly included. Through practice it is found that both approaches lead to the same result and either one can be used. But in subsequent sections the first one is employed since it converges faster. 2.4. Interface elements The combined effect of water pressure and cohesive tractions on crack face can be implemented using interface elements e.g. [31,32]. These elements connect pairs of nodes on the crack faces. Both normal and shear tractions components can be incorporated. The element stiffness matrix KIE in local coordinate system is expressed as

Fig. 3. Scaled boundary finite-element modeling of a plate with a side crack.

𝑁

𝐊IE =

𝑔 𝐴∑ 𝑇 𝐌 𝐤𝐌𝑤 2 𝑖=1 𝑖 𝑖 𝑖 𝑖

In polygons modeling a crack, the scaling center O is selected at the crack tip. The crack faces are formed by scaling the nodes A and B on the boundary and are not discretized. The crack faces are denoted as “sidefaces” in the scaled boundary finite element convention. The material interface is also formed by scaling similarly. The global Cartesian coordinate of a point within a line sector covered by an element on the polygon boundary is expressed as { } { } 𝑥(𝜉, 𝜂) 𝐱 = 𝜉𝐍(𝜂) 𝑏 (4) 𝑦(𝜉, 𝜂) 𝐲𝑏

(2)

where A is the crack surface area, Ng is the number of Gaussian integration points, Mi is the shape function matrix, and wi is the weighting function. The matrix ki is computed from { } [ ]{ } 𝑡𝑠 𝑘 0 𝑤𝑠 = 𝑠 (3) 𝑡𝑛 0 𝑘𝑛 𝑤𝑛 where kn and ks are the secant stiffness of the softening functions shown in Figs. 1 and 2. The interface elements can be coupled with the finite element method or the SBFEM. When coupled with the SBFEM, a shadow domain procedure [29,30] has to be employed. This will be explained in Section 4.

where 𝐍(𝜂) =

[

𝑁1 (𝜂) 0 𝑁2 (𝜂) 0 ⋯ 0 𝑁𝑀 (𝜂) 0 0 𝑁1 (𝜂) 0 𝑁2 (𝜂) 0 ⋯ 0 𝑁𝑀 (𝜂)

] (5)

is the shape function matrix with M number of nodes, xb and yb are the vector of nodal coordinates of the element. 3.2. Displacement and stress interpolation

3. Fundamentals of the scaled boundary finite element method for fracture problems

The displacement field u(𝜉,𝜂) at any point in a sector covered by a line element on a polygon boundary is expressed as

The SBFEM is a semi-analytical technique that has unique advantages when applied to problems involving fracture. In the vicinity of the crack tip, the stress fields in domains of arbitrary geometry are described analytically [36,42,43,26]. This enables accurate assessment of the stress field near the crack tip without the use of nodal enrichment functions or local mesh refinement. The scaled boundary finite element method can also be formulated on polygons with arbitrary number of sides [44,45], which enables crack propagation problems to be handled more flexibly than the standard finite element method. The concept of the scaled boundary finite element method is presented in this section. The solution using generalized coordinates [36] is adopted to treat the cohesive tractions acting on the crack faces. Only the key equations are presented. For a detailed description of the methodology, interested readers are referred to the literature [36].

𝐮(𝜉, 𝜂) = 𝐍(𝜂)𝐮(𝜉)

(6)

where u(𝜉) are nodal displacement functions along the radial direction. The stress field 𝝈(𝜉,𝜂) is 𝛔(𝜉, 𝜂) = 𝐃𝐁1 (𝜂)𝐮(𝜉),𝜉 + 𝜉 −1 𝐃𝐁2 (𝜂)𝐮(𝜉)

(7)

where D is the elasticity constitutive matrix of the material. B1 (𝜂) and B2 (𝜂) are the standard SBFEM strain-displacement matrices [21]. 3.3. Equilibrium condition The equilibrium condition in a polygon can be derived from the principle of virtual work [46] or the method of weighted residuals [21]. Considering a polygon with non-zero side face tractions, the equilibrium condition is given by the scaled boundary finite element equation in displacement [36] [ ] 𝐄0 𝜉 2 𝐮(𝜉),𝜉𝜉 + 𝐄0 − 𝐄1 + 𝐄𝑇1 𝜉𝐮(𝜉),𝜉 − 𝐄2 𝐮(𝜉) + 𝐅(𝜉) = 0 (8)

3.1. Geometric representation on polygons Fig. 3 shows the scaled boundary finite representation of a cracked polygon for a crack. The geometry of the polygon is required to satisfy a scaling condition such that the entire boundary has to be visible from a point inside the domain. This point is called the scaling centre. The scaling condition can always be satisfied by subdividing a sub-domain to smaller ones. In an uncracked polygon, the scaling centre O is chosen at the geometric center. A radial coordinate 𝜉 with 𝜉 = 0 at the scaling centre and 𝜉 = 1 at the polygon boundary is defined. Standard one dimensional finite elements discretize the polygon boundary. On each element on the boundary, a local coordinate 𝜂 that varies from −1 to + 1 is defined.

where Ei , i = 0, 1, 2 are coefficient matrices that depend only on the geometry and material properties of an element sector in the polygon and F(𝜉) is the load vector due to the tractions on the side-faces. Eq. (8) is a system of linear nonhomogeneous second-order ordinary differential equations for the radial displacement functions u(𝜉). The solution of Eq. (8) in terms of generalized coordinates was reported by Song [36]. The variable X(𝜉) { } 𝐮(𝜉) 𝐗(𝜉) = (9) 𝐪(𝜉) 44

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

where q(𝜉) are the internal nodal forces along the radial lines 𝐪(𝜉) = 𝐄0 𝜉𝐮(𝜉),𝜉 +

𝐄𝑇1 𝐮(𝜉)

3.5. Solution for side-face tractions varying as a power function in 𝜉 (10)

When the traction distribution on the side face, tn and ts , can be expressed in the form of power functions in 𝜉, the nodal load vector F(𝜉) can be expressed as

is introduced to transform Eq. (8) into a first-order differential equation { 𝜉𝐗(𝜉),𝜉 = −𝐙𝐗(𝜉) −

} 0 𝐅(𝜉)

with the Hamiltonian matrix [ −1 𝑇 ] 𝐄 𝐄 𝐄−1 0 𝐙= 0 1 −𝐄2 + 𝐄1 𝐄−1 𝐄𝑇1 −𝐄1 𝐄−1 0 0

(11)

𝐅(𝜉) = 𝜉 𝑏 𝐅𝑡

(26)

Eq. (25) can be solved analytically. For the last N-1 diagonal blocks in S with positive real parts of eigenvalues, the condition for finiteness of W(𝜉) at 𝜉 = 0 is satisfied by setting the terms in brackets in Eq. (25) to zero. This determines the integration constants ci as

(12)

3.4. General solution procedure

0

( ( )) 𝜏 −𝐒𝑖 −𝐈 𝐅(𝑖𝑤) (𝜏)𝑑𝜏 𝑤ℎ𝑒𝑛 Re 𝜆 𝐒𝑖 > 0

A block diagonal Schur decomposition of Z in Eq. (12) is first performed. This results in for a polygon [36]

𝐜𝑖 =

𝐙𝚿 = 𝚿𝐒

Using Eq. (26), and substituting Eq. (24) into Eq. (25) and integrating analytically results in the particular solution [36] ( )−1 ( ( )) 𝐖𝑖 (𝜉) = 𝐖(𝑖𝑝) (𝜉) = − 𝐒𝑖 + 𝑏𝐈 𝐅̄ (𝑖𝑤) 𝜉 𝑏 𝑤ℎ𝑒𝑛 Re 𝜆 𝐒𝑖 > 0 (28)

(13)

where S is the Schur matrix and 𝚿 is a transformation matrix. S contains 2 N diagonal blocks of the form ( [ ] ) 𝐒 𝐈 𝐒 = 𝑑𝑖𝑎𝑔 𝐒1 , … , 𝐒𝑁−1 , 𝑁 , 𝐒𝑁+2 , … , 𝐒2𝑁 (14) 0 𝐒𝑁+1

(𝑤)

𝐅𝑖

The diagonal elements of S are equal to the real parts of the eigenvalues of Z with SN = SN +1 = 0. The transformation matrix 𝚿 is J2n -orthogonal [36] and also consists of 2 N blocks as [ ] 𝚿 = 𝚿1 , … , 𝚿𝑁−1 , 𝚿𝑁 , 𝚿𝑁+1 , … , 𝚿2𝑁 (16)

where [ 0 𝐉2𝑛 = 𝐈

] −𝐈 0

(17)

Eq. (31) can be used when b is sufficiently disjoint from 𝜆 (-Si ). However, when b is close to or equal to 𝜆 (-Si ), the matrix Si + bI becomes rank deficient. To overcome this problem, a stable approach to evaluate the particular solution is obtained by augmenting Eq. (23) as [36] { } [ ]{ } 𝐖 (𝜉) 𝐅̄ (𝑖𝑤) 𝐖𝑖 (𝜉) −𝐒𝑖 𝜉 𝑏𝑖 = (32) 𝜉 𝜉𝑏 0 𝑏 ,𝜉 The particular solution Wi (p) (𝜉) is obtained by evaluating the matrix function

The unknown functions X(𝜉) can be decomposed as [36]

⎡− 𝐒 𝑖 ⎢ ⎢ 0 𝜉⎣

(20)

where W(𝜉) are the generalized coordinate functions. 𝚿 can be partitioned into two row blocks with equal size, 𝚿(u) and 𝚿(q) , and further into four square matrices of size n × n as ] [ (𝑢) ] [ (𝑢) 𝚿𝑛 𝚿(𝑝𝑢) 𝚿 𝚿= = (21) 𝚿(𝑞) 𝚿(𝑛𝑞) 𝚿(𝑝𝑞)

] 𝐖(𝑖𝑝) (𝜉) 𝜉𝑏

𝜉 −𝐒𝑖 0

𝑝) 𝐖𝑁+1 (𝜉) = 𝐖(𝑁+1 (𝜉)

(33)

𝑤ℎ𝑒𝑛 𝑤ℎ𝑒𝑛

𝐒𝑁 = 𝐒𝑁+1 = 0 𝐒𝑁 = 𝐒𝑁+1 = 0

The particular solutions are obtained analytically as ( ) 𝑤) 𝑤ℎ𝑒𝑛 𝐒𝑁 = 𝐒𝑁+1 = 0 = −𝑏−1 𝐅̄ (𝑁𝑤) − 𝑏−1 𝐅̄ (𝑁+1

𝐖(𝑁𝑝) (𝜉)

𝑝) (𝑤) 𝐖(𝑁+1 (𝜉) = −𝑏−1 𝐅̄ 𝑁+1 𝜉 𝑏

(23)

where

The solution of Eq. (23) is ( ) 𝜉 𝐖𝑖 (𝜉) = 𝜉 −𝐒𝑖 𝐜𝑖 − 𝜏 −𝐒𝑖 −𝐈 𝐅(𝑖𝑤) (𝜉)𝑑𝜏 ∫1

[ =

𝐖𝑁 (𝜉) = 𝑐𝑁 + 𝐖(𝑁𝑝) (𝜉)

For the diagonal blocks with real parts of eigenvalues satisfying 𝜆(Si )≠0, substituting Eq. (20) into Eq. (21) pre-multiplied by -𝚿ī T J2n and using Eqs. (17) and (19 leads to [36]

𝑙

𝐅̄ (𝑖𝑤) ⎤⎥ 𝑏 ⎥⎦

For the diagonal blocks with zero eigenvalues in S, the solution for W(𝜉) is obtained following a similar procedure as for the diagonal blocks satisfying 𝜆 (Si )≠0, which leads to [36]

Using Eqs. (20) and (21), the radial displacement functions u(𝜉) and the nodal internal forces q(𝜉) can be expressed as { } [ (𝑢) ] 𝐮(𝜉) 𝚿 = 𝐖(𝜉) (22) 𝐪(𝜉) 𝚿(𝑞)

( )𝑇 𝑇 𝚿𝑢 𝐅(𝜉) 𝐅(𝑖𝑤) (𝜉) = −𝐇− 𝑖

(30)

Using Eq. (26) and substituting Eq. (24) into Eq. (25) and integrating analytically in 𝜉, results in the particular solution Wi (p) (𝜉) ( )−1 ( ( )) 𝐖(𝑖𝑝) (𝜉) = − 𝐒𝑖 + 𝑏𝐈 𝐅̄ (𝑖𝑤) 𝜉 𝑏 𝑤ℎ𝑒𝑛 Re 𝜆 𝐒𝑖 < 0 (31)

Pre-multiplying Eq. (13) with -𝚿ī T J2n and using Eq. (17) results in { 𝑇 𝐇𝑖 𝑤ℎ𝑒𝑛 𝑗=𝑖 −𝚿𝑇 𝐉2𝑛 𝚿𝑗 = (19) 𝑙 0 𝑤ℎ𝑒𝑛 𝑗≠𝑖

𝜉𝐖𝑖 (𝜉),𝜉 = −𝐒𝑖 𝐖𝑖 (𝜉) − 𝐅(𝑖𝑤) (𝜉)

(29)

𝑙

𝐖𝑖 (𝜉) = 𝜉 −𝐒𝑖 𝐜𝑖 + 𝐖(𝑖𝑝) (𝜉)

(18)

𝐗(𝜉) = 𝚿𝐖(𝜉)

( )𝑇 𝑇 = −𝐇− 𝚿(𝑢) 𝐅𝑡 𝑖

For the first N-1 diagonal blocks in S with negative real parts of eigenvalues, the condition of finiteness of W(𝜉) at 𝜉 = 0 is satisfied for any ci . The solution consists of two parts i.e. a general solution and a particular solution as

(15)

It satisfies the relation ( )𝑇 {𝐇 𝑤ℎ𝑒𝑛 𝑗 = 𝑖̄ 𝑖 𝚿𝑇𝑖 𝐉2𝑛 𝚿𝑗 = − 𝚿𝑇𝑗 𝐉2𝑛 𝚿𝑖 = 0 𝑤ℎ𝑒𝑛 𝑗 ≠ 𝑖̄

(27)

where

with N conjugate pairs of Si and Sī , and I is an identity matrix. The index of the block conjugate to the ith block is ī 𝑖̄ = 2𝑁 + 1 − 𝑖

∫1

𝑤ℎ𝑒𝑛

𝐒𝑁 = 𝐒𝑁+1 = 0

(34)

(35)

(36)

(37)

3.6. Stiffness matrix and equivalent load vector (24) Once the solutions for Wi (𝜉) have been determined, the displacement functions u(𝜉) and internal nodal forces q(𝜉) can be written as 𝐮(𝜉) = 𝚿(𝑛𝑢) 𝜉 −𝐒 𝐜 + 𝚿(𝑢) 𝐖(𝑝) (𝜉)

(25) 45

(38)

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

Fig. 4. Shadow domain procedure.

𝐪(𝜉) = 𝚿(𝑛𝑞) 𝜉 −𝐒 𝐜 + 𝚿(𝑞) 𝐖(𝑝) (𝜉)

enables the interface element to be directly coupled to the SBFEM just like the FEM. A nonlinear fracture analysis is then performed using the shadow domain. This determines the distribution of the cohesive traction at the crack faces tn and ts during the load step. Although the interface elements and interface elements can determine the distribution of the cohesive tractions on the crack faces, the asymptotic stress field in the vicinity of the crack, which is necessary for accurate evaluation of the stability of a cohesive crack, cannot be directly computed from the shadow domain. This information, however, can be determined from the cracked polygon. The distribution of the cohesive tractions is extracted from the interface elements to determine the nonhomogeneous term F(𝜉) in Eq. (8). The stress field in the vicinity of the crack is then determined considering the nonhomogeneous term F(𝜉) by Eq. (46).

(39)

The integration constants c can be evaluated from the solution at the boundary 𝜉 = 1 by the relations ( ) 𝐮(𝜉 = 1) = 𝚿(𝑛𝑢) 𝜉 −𝐒 𝐜 + 𝐖(𝑛𝑝) (𝜉 = 1) + 𝚿(𝑝𝑢) 𝐖(𝑝𝑝) (𝜉 = 1) (40) ( ) 𝐪(𝜉 = 1) = 𝚿(𝑛𝑞) 𝜉 −𝐒 𝐜 + 𝐖(𝑛𝑝) (𝜉 = 1) + 𝚿(𝑝𝑞) 𝐖(𝑝𝑝) (𝜉 = 1)

(41)

Substituting Eq. (40) into Eq. (41) results in 𝐪(𝜉 = 1) = 𝐊𝐮(𝜉 = 1) − 𝐑𝐹

(42)

where K is the stiffness matrix ( )−1 𝐊 = 𝚿(𝑛𝑞) 𝚿(𝑛𝑢)

(43)

and RF is the nodal load vector due to the side-face traction ( ) 𝐑𝐹 = − 𝚿(𝑛𝑞) − 𝐊𝚿(𝑝𝑢) 𝐖(𝑝𝑝) (𝜉 = 1)

4.2. Stress field at cohesive crack tip Once the distribution of the cohesive tractions on the crack face, tn (𝜉) and ts (𝜉), has been determined, they are assumed to vary as a power function in 𝜉 as

(44)

The integration constants ci can be evaluated using Eq. (38) once u(𝜉 = 1) has been determined. The complete solution for the displacement field is obtained by substituting Eq. (38) into Eq. (6) resulting in ( ) 𝐮(𝜉, 𝜂) = 𝐍(𝜂) 𝚿(𝑛𝑢) 𝜉 −𝐒 𝐜 + 𝚿(𝑢) 𝐖(𝑝) (𝜉) (45)

𝑡𝑛 (𝜉) = 𝑡(𝑛𝑢) 𝑡𝑠 (𝜉) = 𝑡(𝑠𝑢)

The stresses 𝝈(𝜉,𝜂) are obtained by substituting Eq. (38) into Eq. (7) resulting in 𝛔(𝜉, 𝜂) =

𝑁 ∑ 𝑖=1

𝚿𝜎𝑖 (𝜂)𝜉 −𝐒𝑖 −𝐈 𝐜𝑖 +

2𝑁 ∑ 𝑖=1

(

𝜉 −1 𝚿𝜎𝑖 (𝜂)𝐖(𝑖𝑝) (𝜉) −

where 𝚿𝜎𝑖 (𝜂) = [𝚿𝜎𝑥𝑥𝑖 (𝜂) 𝚿𝜎𝑦𝑦𝑖 (𝜂) 𝚿𝜎𝑥𝑦𝑖

𝑗=1 𝑚 ∑ 𝑗=1

𝑒𝑗 𝜉 𝑡𝑗

(47)

𝑔𝑗 𝜉 𝑡𝑗

(48)

with exponents tj = j-1, j = 1, 2, …, m. The coefficients ej and gj are determined by solving the system of m linear equations in Eqs. (47) and (48). The nodal side-face load vector resulting from normal and shear cohesive tractions are formulated by multiplying Eqs. (47) and (48) with the crack surface area A resulting in (𝑚 ) ∑ 𝑡 𝐅𝑡 (𝜉) = 𝑒𝑗 𝜉 𝑗 𝐅̄ 𝑡 (49)

)

𝐃𝐁1 (𝜂)𝚿(𝑖𝑢) 𝐅̄ (𝑖𝑤) (𝜉)

(46) (𝜂)]𝑇

𝑚 ∑

is the stress mode [36].

4. Cohesive crack propagation modelling with SBFEM

𝑛

(

When modeling cohesive fracture problems using the SBFEM, the cohesive tractions on the crack faces are the side-face loads resulting in the nonhomogeneous term F(𝜉). The distribution of F(𝜉) is nonlinear because it depends on the relative displacement of the crack faces during each load step. In order to use the equations derived in Section 3, a twostep approach known as the shadow domain procedure is adopted [29].

𝐅𝑡𝑠 (𝜉) =

𝑗=1

𝑚 ∑ 𝑗=1

𝑛

) 𝑔𝑗 𝜉 𝑡𝑗 𝐅̄ 𝑡𝑠

(50)

where 𝐅𝑡𝑛 and Fts are expressed in global Cartesian coordinates as [ ]𝑇 𝐅̄ 𝑡𝑛 = 𝐴𝑡(𝑛𝑢) − sin 𝛼 cos 𝛼 0 0 ⋯ 0 0 sin 𝛼 − cos 𝛼 (51)

4.1. Shadow domain procedure

]𝑇 [ 𝐅̄ 𝑡𝑠 = 𝐴𝑡(𝑛𝑢) − cos 𝛼 sin 𝛼 0 0 ⋯ 0 0 cos 𝛼 − sin 𝛼

In the first step, the distribution of the cohesive traction on the crack faces is determined using interface elements. To couple the interface elements to the SBFEM, the cracked polygon is partitioned into two normal polygons as shown in Fig. 4. The resulting mesh is called the shadow domain [29,45]. In shadow domain, the crack faces are discretized and

and 𝛼 is the crack angle shown in Fig. 4. Note that only the entries corresponding to the nodes on the crack faces are non-zero. The nonhomogeneous term F(𝜉) in Eq. (8) therefore, becomes (𝑚 ) (𝑚 ) ∑ 𝑡 +1 ∑ 𝑡𝑘 ̄ ̄ 𝑗 𝐅𝑡 + 𝐅𝑡 𝐅(𝜉) = 𝑒𝑗 𝜉 𝑔𝑘 𝜉 (53) 𝑗=1

46

𝑛

𝑘=1

𝑠

(52)

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

S(s) , it is possible to define the singular stress field as 𝛔(𝑠) (𝜉, 𝜂) =

𝑁 ∑ 𝑖=1

(𝑠)

𝚿𝜎 (𝑠) (𝜂)𝜉 −𝐒𝑖 𝑖

−𝐈 (𝑠) 𝐜𝑖

(55)

where 𝚿𝜎 (𝑠) (𝜂) and ci (s) are the stress modes and integration constants 𝑖

corresponding to S(s) . The stress intensity factors are defined by {√ } { } 𝐾𝐼 2𝜋𝑟𝜎𝑦𝑦 ||𝜃=0 = lim √ 𝐾𝐼𝐼 𝑟→0 2𝜋𝑟𝜏𝑥𝑦 ||𝜃=0

(56)

Substituting the stress components in Eq. (55) at angle 𝜃 = 0 into Eq. (56) and using the relation 𝜉 = r/LOA at 𝜃 = 0 in Fig. 4a results in { (𝑠) } { } √ 𝚿𝜎𝑦𝑦𝑖 (𝜂(𝜃 = 0))𝐜(𝑖𝑠) 𝐾𝐼 = 2𝜋𝐿𝑂𝐴 (57) ) (𝑠) 𝐾𝐼𝐼 𝚿(𝜏𝑠𝑥𝑦𝑖 (𝜂(𝜃 = 0))𝐜𝑖 For the hydraulic fracture analysis motivated by reservoir overflow of a typical gravity dam-foundation system (Fig. 5), the loads include self-weight G of the dam, water pressure FL1 on the upstream face and FL2 on the preset crack corresponding to a full reservoir, water pressure FO1 on the upstream face and FO2 on the preset crack corresponding to an elevated water level, and water pressure FC (wn ) in the newly formed crack due to the full reservoir and elevated water level. Among them, G, FL1 and FL2 are constant, FO1 and FO2 are linearly proportional to the water level elevation, FC (wn ) is dependent on the crack opening wn . In the first load step, the crack has not propagated, so no cohesive interface elements are present, and no overflow is considered. The following equation is solved:

Fig. 5. Schematic showing the external loads acting on a dam with a crack at the damfoundation interface.

where each j or k term has the same form as F(𝜉) in Eq. (26). Therefore, the stress field can be evaluated using Eq. (46).

4.3. Stability of cohesive crack

𝐾01 𝑢1 = 𝐺 + 𝐹𝐿1 + 𝐹𝐿2

There is limited documented literature on the stability of bi-material interface cracks subjected to cohesive tractions on the crack faces. A universal criterion to determine this condition is not available. In this paper, the stability of the bi-material interface crack is determined by the condition

here K01 is the stiffness matrix of the dam and foundation bulk for the current geometry, u1 is the displacement vector to be solved. In subsequent load steps, water pressure due to reservoir overflow is included, and cohesive interface elements are introduced along the newly-formed crack. The equilibrium equation for the (j + 1)th iteration of the ith load step becomes ( ( )) ( ) ( ) ( ) = Δ𝜆𝑗𝑖 𝐹𝑜1 + 𝐹𝑜2 + 𝐹𝑐 𝑢𝑗𝑖 − 𝐹𝑐 𝑢𝑖−1 𝐾0𝑖 + 𝐾𝐶𝑂𝐻 𝑢𝑗𝑖 Δ𝑢𝑗+1 (59) 𝑖

𝐾1 = 0, 𝐾2 = 0

(54)

This condition is similar to enforcing the zero-K condition [29,47] i.e. the stress at the crack has to be finite (i.e. there is no stress singularity). Unlike a crack in a homogeneous medium, both K1 and K2 have to be considered when determining the stability of a bi-material interface crack because the magnitudes of the tensile and shear tractions on the crack face are of the same order. The SBFEM enables accurate calculation of the stress intensity factors directly from their definitions. From Eq. (46), it can be identified that as 𝜉→0 the diagonal blocks in S with real parts of eigenvalues satisfying −1 < Re(𝜆(S)) < 0 becomes dominant. Denoting this diagonal block as

(58)

here K0 i is the stiffness matrix of the dam and foundation bulk for the current geometry, KCOH is the stiffness of the cohesive interface elements dependent on the crack profile, Δui j +1 is the displacement increment of the current iteration, Δ𝜆i j is the amplification factor to scale the water pressure due to reservoir overflow of the previous iteration, Fc is the water pressure inside the crack other than the preset crack, ui j is the total displacement vector of the previous iteration, ui-1 is the total displacement at the previous load step. A local arc-length procedure [29] is used for solution of this equation. Convergence of the nonlinear iteration is decided by the norm of the out-of-balance force vector.

Fig. 6. Embankment-foundation system.

47

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

Fig. 7. Overflow vs. crack mouth opening displacement.

Fig. 8. Overflow vs. horizontal displacement of embankment crest.

5. Numerical examples 5.1. An embankment on a rigid foundation An embankment resting on a rigid and impervious foundation shown in Fig. 6a is considered. This problem was previously investigated by Bolzon and Cocchetti [3] by using a friction-cohesive interface model and will be used to validate the developed formulation. An initial crack with the length of 1 m is introduced at the embankment-foundation interface (Fig. 6a). The embankment-foundation system is discretized by the polygon mesh shown in Fig. 6b. The mesh has 159 polygons and 796 degrees-of-freedom. The geometry of the embankment is L = 30 m and h = 10 m. The material properties used are the same as that reported in [3]. For the embankment, the Young’s modulus Ec = 24 GPa, Poisson’s ratio 𝜈 c = 0.15 and density 𝜌 = 2000 kg/m3 . For the foundation, the Young’s modulus Ef = 240 GPa and Poisson’s ratio 𝜈 f = 0. The bulk material is assumed to be linear elastic and the crack propagation takes place only along the interface. In the process zone, the relation between the normal cohesive traction tn and the crack opening displacement wn is 𝑡𝑛 = 𝑡(𝑛𝑢) 𝑒−𝛼𝑤𝑛

Fig. 9. Distribution of normal cohesive traction along the crack.

(60)

(u )

where tn = 0.3 MPa is the tensile strength of the interface and the coefficient 𝛼 is 3.333mm−1 . For this example, shear decohesion is not considered, which is implemented by assuming a very high shear strength for the cohesive interface elements. The external loads acting on the embankment-foundation system include the self-weight of the embankment, the water pressure corresponding to a full reservoir, the overloading of water pressure caused by an elevated reservoir and the hydrostatic pressure on the crack face. The distribution of the hydrostatic pressure on the crack face is assumed to follow from Eq. (1). A parametric study is carried out to determine the effect of 𝜌 on the crack propagation behavior of the interface crack. Five values of 𝜌 are considered, i.e., 𝜌 = 0, 𝛼/100, 𝛼/10, 𝛼, 3𝛼. For the case of 𝛼 = 0, no water pressure in the pre-set crack is considered whereas for other values of 𝛼, the water pressure corresponding to the hydrostatic pressure is considered at the preset crack. Fig. 7 shows the predicted overflow versus crack mouth opening displacement at the mouth of the initial crack for all the values of 𝜌 considered. For the case where there is no water pressure acting on the crack face (𝜌 = 0), the reservoir overloading is embodied only in the horizontal pressure on the upstream face. Therefore, the peak load of the system is the higher compared with the all the other cases where water pressure inside the crack is considered. From the equilibrium of forces at

Fig. 10. Distribution of water pressure along the crack.

the base point D, it is found that an overflow of 176.7 m is required to completely detach the embankment from the foundation. From Fig. 7, it is clear that both the present result and that of Bolzon and Cocchetti [3] converge to this value as the crack propagates. The main difference between the two models is that the friction at the interface is considered in [3], which contributes to the difference of the shape of the curves. With the increase of 𝜌, and hence the increase of water pressure, the 48

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53 Table 2 Material properties of dam and foundation.

Dam concrete Foundation rock

Elastic modulus (GPa)

Poisson’s ratio

Mass density (kg/m3 )

27.85 18.10

0.178 0.323

2354 2720

Fig. 11. Geometry and Scaled boundary polygon discretization of the dam-foundation system.

influence of friction at the interface becomes insignificant, which can be seen from the subsequent results. From Fig. 7, it is evident that as 𝜌 increases, the peak overflow decreases, and the consistency of the present result and that of Bolzon and Cocchetti [3] become better. For 𝜌 = 𝛼/10, 𝛼 and 3𝛼, excellent agreement is observed. For 𝜌 = 3𝛼, the peak overflow is 36.5 m. Compared with the case with no water pressure in the crack, the maximum overflow is reduced by 80% when the crack opening displacement reaches only 8.5% of the response when 𝜌 = 0. For the cases with 𝜌 = 𝛼/100, 𝛼/10, 𝛼, the maximum overflow were reduced by 44%, 58% and 74%, respectively. The significance of the water pressure on the load bearing capacity of a dam-foundation system is clearly demonstrated. The relation between the overflow and the horizontal displacement of the embankment crest is shown in Fig. 8. The distribution of normal cohesive traction along the crack corresponding to peak overflow is shown in Fig. 9. When 𝜌 = 0, peak overflow is attained when the crack is 21.3 m. As 𝜌 increases, the crack gets shorter dramatically, with the crack length equal to 16.9 m, 9.7 m, 5.33 m and 5.33 m corresponding to 𝜌 = 𝛼/100, 𝛼/10, 𝛼, 3𝛼, respectively. For 𝜌 = 0 and 𝛼/100, the cohesive

Fig. 12. Overflow vs. crack mouth opening displacement.

zone is fully developed. Fig. 10 shows the water pressure distribution inside the crack. Only for the case where 𝜌 = 3𝛼 is the water pressure at the tip of the preset crack is close to the hydrostatic pressure. For the other cases, the water pressure in the crack is much smaller than the hydrostatic pressure. 5.2. Dam-foundation system The concrete dam-rock foundation system shown in Fig. 11 proposed as a benchmark by the International Commission On Large Dams (ICOLD) [48] is considered. A preset crack with the length of 2 m is introduced at the dam-foundation interface. The dam-foundation system is discretized using 326 polygons with arbitrary number of sides (Fig. 11). The total number of degrees-of-freedom is 1540. The material properties obtained from experiments [40] on concreterock interfacial fracture are used for the dam-foundation system and

Fig. 13. Deformed dam-foundation system and distribution of cohesive traction.

49

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

Fig. 16. Overflow vs. horizontal displacement of dam crest. Fig. 14. Sketch map of water pressure distribution inside the crack.

5.2.2. Case with water pressure in the crack A parametric study is carried out considering the four types of water pressure distribution on the crack face shown in Fig. 14, where p0 is the full hydrostatic pressure corresponding to current water level. Four distribution schemes of water pressure inside the crack are considered (Fig. 14). For all the four schemes, the water pressure is a function of the crack opening displacement. In the water pressure Schemes 1, 2 and 3, no pressure is present from the crack tip to the knee point of the cohesive law (normal component as in Fig. 1a) to model the phenomenon of water lag. The water pressure starts from the knee point and follows an exponential relation (Eq. (1)) with the crack opening. Water pressure equals zero at the knee point and increases to full hydraulic pressure corresponding to crack openings of wnc , 0.5(wnc + wn1 ) and wn1 , respectively for the three schemes. The corresponding values for the parameter 𝜌 are 4797, 9594 and 479,700, respectively. In the fourth scheme, no water lag is present and it is assumed that the water penetrates to the crack tip. The pressure increases from zero at the crack tip to full hydrostatic pressure at the knee point. The predicted overflow versus crack mouth displacement relations are shown in Fig. 15. The relation between the overflow and horizontal displacement of the dam crest is shown in Fig. 16. From Fig. 15, it is observed that unlike the case where there is no water pressure, the overflow decreases after attaining a maximum value. The peak overflows for the first three types of water pressure distribution are very close (approximately 8.16 m). There is negligible difference in the predicted structural response. From this analysis, the distribution of water pressure from the knee point to the maximum opening of the cohesive

are shown in Table 2. The constitutive relations in the process zone are modified from the experimental measurements in [40] and are shown in Fig. 1. The parameters of the softening curves are listed in Table 1. 5.2.1. Case without water pressure in the crack For comparison, the case of no water pressure in the crack including the preset crack is considered first. Only the water pressure on the upstream face varies as the overflow varies. The predicted overflowcrack mouth opening displacement is depicted in Fig. 12. The overflow increases with the crack mouth opening displacement, which indicates that the load bearing capacity increases although the ligament on the interface becomes shorter due to cracking. This phenomenon is partly due to the stabilizing effect of the self-weight of the dam and the absence of the water pressure inside the crack. From the numerical simulation, the maximum predicted overflow is 42.79 m and the crack length is 54.9 m (Fig. 12). By considering the equilibrium of forces at the dam toe and assuming the dam is completely detached from the foundation, the maximum overflow can be predicted as 43.95 m. The slight difference is due to the remaining short ligament of interface when dam fails in numerical simulation (Fig. 13). The crack mouth opening displacement is 140.8 mm and is considerably large in the view of engineering practice. For a more realistic simulation, the hydraulic pressure acting on the crack faces has to be considered. Fig. 13 shows the polygon mesh of the dam along with the normal and shear traction distribution at maximum overflow. The deformation of the dam is amplified by a scaling factor of 50.

Fig. 15. Overflow vs. crack mouth displacement.

50

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

(a) Global view

(b) Partial view

Fig. 17. Distribution of water pressure along the crack.

Fig. 18. Distribution of cohesive traction along the crack.

Fig. 19. Crack deformation.

law does not affect the load bearing capacity of the dam. So for sake of convenience, rectangular distribution of water pressure can be assumed for this section. Compared with the water pressure distributions of Schemes 1−3, the predicted peak overflow for Scheme 4 is 7.22 m, which is approximately 11.5% lower. Although the additional water pressure is distributed over only a narrow area from the crack tip to the knee point, its influence is significant. However, it is also noticed that with the propagation of the crack, the crack opening near the crack tip increases, this influence gets less significant and the curves for all the four schemes coincide.

The distribution of water pressure, cohesive tractions along the crack and the relative displacements at maximum overflow are shown in Figs. 17–19, respectively. As can be seen from Fig. 17, the difference in the water pressure distributions of Schemes 1, 2 and 3 are distinct. The water pressure distribution of Schemes 3 and 4 are very similar, the only noticeable difference occurs at distances between 4.62 m and 4.73 m from the dam heel (Fig. 17b). For Scheme 3, the water pressure decreases from full hydrostatic pressure to zero. Whereas for Scheme 4, the water pressure decreases from full hydrostatic pressure to 0.22 MPa (equivalent to full hydrostatic pressure corresponding to a water level of 22 m). This difference in water pressure distribution, although exists

51

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

Fig. 20. Cracking process of ICOLD dam for Type 3 water pressure distribution.

over only a very short part of the crack, results in significantly different values of the peak overflow. From Figs. 18 and 19, the cohesive traction distribution and crack displacement are similar regardless of the type of water pressure distribution. Both the crack opening and sliding displacements decrease from the crack mouth to the crack tip. The magnitude of the crack sliding displacement is smaller than the crack opening displacement. Compared with the tangential component of the cohesive tractions, the normal component is much smaller for the whole newly-formed crack This phenomenon is different from the case where there is no water pressure inside the crack. In the latter case, the normal component is less prominent. This is due to the fact that water pressure inside the crack always forces the crack to open. In all the four schemes of water pressure distribution considered, the peak overflow is reached after the crack has propagated a distance of 2.4637 m. Fig. 20 shows the crack propagation process for the case with Scheme 3 water pressure distribution. The deformation is amplified by a scaling factor of 500. The crack propagation process when other types of water pressure distribution are considered is very similar.

sliding displacement since the water pressure inside the crack always has a tendency to force the crack to open. 2. Water pressure distributed from the fictitious crack tip to the knee point of the normal cohesive law significantly decreases the bearing capacity of the dam at the interface. But the detailed form of water pressure from the knee point to the maximum opening on the normal cohesive law has only negligible influence. For convenience, a rectangular distribution can be assumed for this regime. The water lag, although exists only in a very short part of the crack, has significant influence on the peak overflow that the dam can resist. For the case in the current research, a difference of 11.5% is witnessed. Acknowledgments This study was supported by the National Key Research and Development Program of China through Grant number 2017YFC0404903 and the National Science Foundation of China through Grant number 51579033. References

6. Conclusions

[1] Rescher OJ. Importance of cracking in concrete dams. Eng Fract Mech 1990;35:503–24. [2] Boone TJ, Ingraffea AR. A numerical procedure for simulation of hydraulically– driven fracture propagation in poroelastic media. Int J Numer Anal Methods Geomech 1990;14:27–47. [3] Bolzonn G, Cocchetti G. Direct assessment of structural resistance against pressurized fracture. Int J Numer Anal Methods Geomech 2003;27:353–78. [4] Alfano G, Marfia S, Sacco E. A cohesive damage-friction interface model accounting for water pressure on crack propagation. Comput Methods Appl Mech Eng 2006;196:192–209. [5] Federal Regulatory Commission, Engineering guidelines for the evaluation of hydropower projects, Federal Energy Regulatory Commission, Office of Hydropower Licensing (1991). [6] Dewey RR, Reich RW, Saouma VE. Uplift modeling for fracture mechanics analysis of concrete dams. ASCE J Struct Eng 1994;120:3025–44. [7] Plizzari GA. On the influence of uplift pressure in concrete gravity dams. Eng Fract Mech 1998;59:253–67. [8] Brühwiler E, Saouma VE. Water fracture interaction in concrete–Part II: hydrostatic pressure in cracks. ACI Mater J 1995;92:383–90. [9] Slowik V, Saouma VE. Water pressure in propagating concrete cracks. J Struct Eng 2000;126:235–42. [10] Dugdale D. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4. [11] Barenblatt G. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962;7:55–129. [12] Hillerborg A, Modeer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanisms and finite elements. Cement Concr Res 1976;6:773–82. [13] Barpi F, Valente S. Modeling water penetration at dam-foundation joint. Eng Fract Mech 2008;75:629–42. [14] Reich W, Bruhwiler E, Solwik V, Saouma VE. Experimental and computational aspects of a water/fracture interaction. Dam fracture and damage. Bourdarot E, Mazars J, Saouma VE, editors. The Netherlands: Balkema; 1994. [15] Bruhwiler E, Saouma VE. Water fracture interaction in concrete-Part 1: fracture properties. ACI Mater J 1995;92:296–303. [16] Barpi F, Valente S. The cohesive frictional crack model applied to the analysis of the dam-foundation joint. Eng Fract Mech 2010;77:2182–91. [17] Desroches J, Detournay E, Lenoach B, Papanastasiou P, Pearson JRA, Thiercelin M, Cheng A. The crack tip region in hydraulic fracturing. Proc R Soc 1994;447:39–48.

The scaled boundary finite element method (SBFEM) coupled with the cohesive crack model is extended to investigate the hydraulic fracture problem along the gravity dam-rock foundation interface. The crack is loaded with water pressure and cohesive traction, with the shear part as prominent as the normal part. The displacement and stress fields are semi-analytically obtained through the SBFEM. The stress intensity factors that are evaluated directly from this stress field enables the zeroK criterion that governs the stability of a cohesive crack to be conveniently and accurately determined. Therefore, a relatively coarse mesh compared to that required by the FEM can be used. The use of arbitrary n-side polygons provides great flexibility in mesh generation. The proposed procedure is validated through the hydraulic fracture of a rectangular embankment located on rigid foundation. As an application, the procedure is employed to study the hydraulic fracture of a benchmark dam-foundation system. The parameters for the cohesive crack model are obtained from previous lab tests on concrete-rock specimens, based on which the cracking behavior of the interface due to different mode mixity angles can be reproduced. Four schemes of water pressure distribution are investigated, and the case of no water pressure inside the crack is also studied for comparison. The following conclusions are arrived at: 1. Water pressure inside the crack greatly decreases the capacity of the dam at the interface. When no water pressure exists, the reservoir overflow keeps increasing as the crack grows due to the stabilizing effect of self-weight, and no softening section is observed. When water pressure inside the crack is considered, the peak overflow decreases to less than 20% of the case without water pressure. And in this case, the crack opening is much more prominent than the crack 52

H. Zhong et al.

Engineering Analysis with Boundary Elements 88 (2018) 41–53

[18] Alberto A, Valente S. Asymptotic fields at the tip of a cohesive frictional crack growing at the bi-material interface between a dam and the foundation rock. Eng Fract Mech 2013;108:152–61. [19] Xiao QZ, Karihaloo BL, Liu XY. Incremental-secant modulus interation scheme and stress recovery for simulating cracking process in quasi-brittle materials using XFEM. Int J Numer Methods Eng 2007;69:2606–35. [20] Karihaloo BL, Xiao QZ. Asymptotic fields ahead of mixed mode frictional cohesive cracks. J Appl Math Mech 2010;90:710–20. [21] Song C, Wolf JP. The scaled boundary finite-element method-alias consistent infinitesimal finite-element cell method-for elastodynamics. Comput Methods Appl Mech Eng 1997;147:329–55. [22] Doherty JP, Deeks AJ. Adaptive coupling of the finite-element and scaled boundary finite-element methods for non-linear analysis of unbounded media. Comput Geotech 2005;32:436–44. [23] Genes MC, Kocak S. Dynamic soil-structure interaction analysis of layered unbounded media via a coupled finite element/boundary element/ scaled boundary finite element model. Int J Numer Methods Eng 2005;62:798–823. [24] Man H, Song C, Gao W, Tin-Loi F. A unified 3D-based technique for plate bending analysis using scaled boundary finite element method. Int J Numer Methods Eng 2012;91:491–515. [25] Li QH, Chen SS, Luo XM. Steady heat conduction analyses using an interpolating element-free Galerkin scaled boundary method. Appl Math Comput 2017;300:103–15. [26] Song C, Tin-Loi F, Gao W. A definition and evaluation procedure of generalized stress intensity factors at cracks and multi-material wedges. Eng Fract Mech 2010;77:2316–36. [27] Goswami S, Becker W. Computation of 3-D stress singularities for multiple cracks and crack intersections by the scaled boundary finite element method. Int J Fract 2012;175:13–25. [28] Chen SS, Li QH, Liu YH, Xue ZQ. Mode III 2D fracture analysis by the scaled boundary finite element method. Acta Mech Solida Sin 2013;26:619–28. [29] Yang ZJ, Deeks AJ. Fully-automatic modeling of cohesive crack growth using a finite element-scaled boundary finite element coupled method. Eng Fract Mech 2007;74:2547–73. [30] Ooi ET, Song C, Tin-Loi F, Yang ZJ. Automatic modeling of cohesive crack propagation in concrete using polygon scaled boundary finite elements. Eng Fract Mech 2012;93:13–33. [31] Goodman RE, Taylor RL, Brekke TL. A model for the mechanics of jointed rock. J Soil Mech Found 1968;94(3):660–77. [32] Gerstle W, Xie M. FEM modeling of fictitious crack propagation in concrete. ASCE J Eng Mech 1992;118(2). [33] Shi M, Zhong H, Ooi ET, Zhang C, Song C. Modelling of crack propagation of gravity dams by scaled boundary polygons and cohesive crack model. Int J Fract 2013;183:29–48.

[34] Bouchard PO, Bay F, Chaastel Y, Tovena I. Crack propagation modelling using an advanced remeshing technique. Comput Methods Appl Mech Eng 2000;189:723–42. [35] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng 1999;46:131–50. [36] Song C. Analysis of singular stress fields at multi-material corners under thermal loading. Int J Numer Methods Eng 2006;65:620–52. [37] Alfano G. On the influence of the shape of the interface law on the application of cohesive-zone models. Compos Sci Technol 2006;66:723–30. [38] Mahalingam Sakethraman. Study of Interfacial Crack Propagation In Flip Chip Assemblies With Nano-Filled Underfill Materials Ph.D. thesis. Georgia Institute of Technology; August 2005. [39] Toftegaard H, Goutianos S, Sørensen BF. Interfacial crack growth in bimaterial composite specimens. In: Proceedings of the 17th International Conference on Composite Materials, Edinburgh, UK; 2009. [40] Zhong H, Ooi ET, Song C, Ding T, Lin G, Li H. Experimental and numerical study of the dependency of interface fracture in concrete–rock specimens on mode mixity. Eng Fract Mech 2014;124−125:287–309. [41] Bazant ZP, Yu Q, Size effect in concrete specimens and structures: new problems and progress, in Li VC, Leung KY, William KJ, Billington SL, ed., Fracture Mechanics of Concrete Structures (Vail, Colorado:, 2004), pp. 153–162. [42] Song C, Wolf JP. Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary finite-element method. Comput Struct 2002;80:183–97. [43] Song C. Evaluation of power-logarithmic singularities, t-stresses and higher order terms of in-plane singular stress fields at cracks and multi-material corners. Eng Fract Mech 2005;72:1498–530. [44] Ooi ET, Song C, Tin-Loi F, Yang Z. Automatic modelling of cohesive crack propagation in concrete using polygon scaled boundary finite elements. Eng Fract Mech 2012;93:13–33. [45] Ooi ET, Song C, Tin-Loi F, Yang Z. Polygon scaled boundary finite elements for crack propagation modeling. Int J Numer Methods Eng 2012;91(3):319–42. [46] Deeks AJ, Wolf JP. A virtual work derivation of the scaled boundary finite-element method for elastostatics. Comput Mech 2002;28(6):489–504. [47] Bazant ZP, Li NY. Stability of cohesive crack model: part 1 - energy principles. ASME J Appl Mech 1995;62:959–64. [48] CIGB/ICOLD Benchmark A2: Imminent Failure Flood for a Concrete Gravity Dam, 5th Benchmark Workshop on Numerical Analysis of Dams, June 2–5 Denver Colorado USA, 1999. [49] Dong W. Fracture mechanisms of rock-concrete interface: experimental and numerical. J Eng Mech ASCE 2016;142(7) 04016040.

53