- Email: [email protected]
A novel linear triangular element of a three-dimensional displacement discontinuity method
Engineering Analysis with Boundary Elements 59 (2015) 89–96
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A novel linear triangular element of a three-dimensional displacement discontinuity method Wan Cheng a,b,c, Yan Jin a,b,n, Hong Li d, Mian Chen a,b a
College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China State Key Laboratory of Petroleum Resources and Prospecting, Beijing 102249, China c School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA d Dalian University of Technology, Liaoning, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 2 March 2015 Received in revised form 27 April 2015 Accepted 30 April 2015
Since only the boundary of the domain requires discretization, the boundary element method (BEM) is very efficient for the semi-infinite or infinite rock-related engineering problems, e.g., hydraulic fracturing in reservoir stimulation and rock cutting during excavation. A real fracture in the solid is usually of an arbitrary geometry in three dimensions, which usually requires a three-dimensional displacement discontinuity method (3D DDM) to determine the deformation and stress field in order to achieve reliable results. However, the use of 3D DDM with triangular elements is limited by the singularities of the integral either within or nearby the domain. In this paper, a novel linear triangular element with three nodes on its vertices is proposed. The analytical integral expressions of this linear triangular element are also theoretically derived. A solution procedure is also described which can be applied to determine the displacement and stress field around a three-dimensional fracture inside the infinite solid. The accuracy of these results are compared with the analytical solutions of the displacements and stresses induced by a pressurized penny-shaped. This procedure takes a shorter time and requires less elements than the usual constant DDM when achieving the same accuracy. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Linear triangular element Boundary element method Displacement discontinuity method Three-dimensional fracture Hydraulic fracturing
1. Introduction Many important practical problems in science and engineering can be reduced into mathematical models that belong to a class of problems known as boundary value problems [5,6]. These are all characterized by a region of interest R enclosed within a boundary C. If R is a two dimensional region, C will be a bounding contour. If R is a three dimensional region then C is a bounding surface. For a fracture or cavity inside an infinite region, C will be the fracture or cavity surface. Only the boundary C is divided into elements in BEM, which is totally different from the network of elements for the entire region R in a finite element method (Fig. 1). The displacement discontinuity in BEM can be considered as a relative displacement between the two surfaces of a plane crack [24,2,23]. The 2D DDM of constant displacement element in linear elastic mechanics was systematically developed by Crouch [5] who proposed a general numerical algorithm to solve various boundary value problems. The 3D DDM is based on the analytical solution to the problem of a discontinuous displacement over a finite area in
n Corresponding author at: College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China. E-mail address: [email protected] (W. Cheng).
http://dx.doi.org/10.1016/j.enganabound.2015.04.020 0955-7997/& 2015 Elsevier Ltd. All rights reserved.
an infinite region. For a crack of arbitrary non-planar shape in an infinite three-dimensional space, the crack is usually required to be divided into triangular elements. Even though the strong singularities in boundary integral equation [7,29,1,21,22], the 3D DDM on the triangular element is improved from constant element to higher order element. Kuriyama et al. [9,10] proposed the 3D DDM with the boundary surface division into constant triangular elements and briefly presented the solution procedure. Davey and Hinduja [8] established the integral solutions for triangular elements applied to potential problems. Based on Davey's idea, Milroy et al. [14] obtained the general solution by a linear combination of the three integral solutions on the nodes of a triangular element. Napier and Malan [16] and Shou et al. [25] described a higher order displacement discontinuity method in conjunction with shape functions. Nikolski et al. [17,18] presented a boundary element by decoupling the complex variables for the shear tractions and real variable equation for the normal traction for solving 3D elastostatic problems of fractured rock under gravity. Then the complex variables BEM is extended to incorporate higher order approximations of the displacement discontinuities Nikolski [19,20,15]. The recent development of multi-stage hydraulic fracturing technology in a horizontal well in petroleum reservoir requires hydraulic fracturing simulators which can consider the three-
90
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
To numerically implement the displacement discontinuity element in a boundary element method, we need the analytical solution for the displacement discontinuity over an area which is usually defined as a boundary element. The general form solution for a displacement discontinuity element can be expressed as follows [6,23,24,26]:
Fig. 1. The element in BEM and FEM.
8 u ¼ 2ð1 νÞΦx;z zΦx;xx zΦy;xy ð1 2νÞΦz;x þ zΦz;xz > < x uy ¼ zΦx;xy þ 2ð1 νÞΦy;z zΦy;yy ð1 2νÞΦz;y þ zΦz;yz > : u ¼ ð1 2νÞΦ zΦ z x;x x;xz þ ð1 2νÞΦy;y zΦy;yz þ 2ð1 νÞΦz;z zΦz;zz
ð2Þ 8 σ xx ¼ 2G 2Φx;xz zΦx;xxx þ 2νΦy;yz zΦy;xxy þ Φz;zz þ ð1 2νÞΦz;yy zΦz;xxz > > > > > σ yy ¼ 2G 2vΦx;xz zΦx;xyy þ 2Φy;yz zΦy;yyy þ Φz;zz þ ð1 2νÞΦz;xx zΦz;yyz > > > > < σ zz ¼ 2G zΦx;xzz zΦy;yzz þ Φz;zz zΦz;zzz σ xy ¼ 2G ð1 νÞΦx;yz zΦx;xxy þ ð1 νÞΦy;xz zΦy;xyy ð1 2νÞΦz;xy þ zΦz;xyz > > > > > σ yz ¼ 2G νΦx;xy zΦx;xyz þ Φy;zz þ νΦy;xx zΦy;yyz zΦz;yzz > > > > : σ xz ¼ 2GΦx;zz þ νΦx;yy zΦx;xxz νΦy;xy þ zΦy;xyz zΦz;xzz
ð3Þ
Fig. 2. Triangular element in the 3D DDM.
dimensional formation and stress field around the fracture [3,4]. Physically, any fracture in an infinite three-dimensional space should be in three dimensions. The 3D DDM can be applied to simulate this kind of engineering problem numerically by assuming that the formation is a linear elastic medium. For the strong singularities in boundary integral equation, researchers hope to improve both the computational efficiency and accuracy utilizing the analytical integration instead of the numerical integration. Therefore, a novel linear displacement triangular element with three nodes on its vertices is proposed in this paper. If the displacement discontinuities on these three nodes keep the same, this linear displacement triangular element will be simplified as the constant displacement triangular. The idea of this linear triangular element can also be applied to determine linear triangular element in a fictitious stress method [26]. The analytical integral solutions of this linear triangular element are also theoretically derived. A solution procedure is also briefly described which can be applied to determine the stresses field around a three-dimensional fracture inside the solid. The accuracy of these results is compared with the analytical solutions of displacements and stresses induced by a pressurized penny-shaped crack [27].
2. Three-dimensional displacement discontinuity element As shown in Fig. 2, the two surfaces of a crack are separated with relative displacements. The two surfaces can be distinguished by the positive side (z¼ 0 þ ) and negative side (z¼ 0 ). Crossing from one side to the other, the displacement discontinuities are given by the relative displacement in three orthogonal directions respectively [23,24]. 8 > < Dx ðx; yÞ ¼ ux ðx; y; 0 Þ ux ðx; y; 0 þ Þ Dy ðx; yÞ ¼ uy ðx; y; 0 Þ uy ðx; y; 0 þ Þ ð1Þ > : D ðx; yÞ ¼ u ðx; y; 0 Þ u ðx; y; 0 Þ z z z þ where Di (i¼ x,y,z) are the displacement discontinuities in x,y,z direction respectively. ui (i¼ x,y,z) are the displacements in x,y,z direction respectively.
where σ ij ; ui ði; j ¼ x; y; zÞ represent the stress and displacement of Kelvin's problem in three dimensional space respectively. G, v are the shear modulus and Possion's ratio of the isotropic medium, respectively. Φi;j ; Φi;jk ; Φi;jkl ði; j; k; l ¼ x; y; zÞ are the partial derivatives with respect to each corresponding suffix of the following function,Φi . The function, Φi , from a triangular element against a point can be represented by [12,15] Φi ðx; y; zÞ ¼
1 Di ðξ; ηÞ ffidξ dη ∬Δ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8πð1 νÞ ðx ξÞ2 þ ðy ηÞ2 þ z2
i ¼ x; y; z
ð4Þ
where (x,y,z) are the coordinates of an arbitrary point in an infinite 3D space. (ξ, η, 0) are the local coordinates of the loading point. Di ðξ; ηÞ represents the discontinuous displacement on point (ξ, η, 0). Because of the strong singularity of Eq. (4) when (x,y,z) goes to (ξ, η, 0), it is extremely difficult to get the analytical integration [7]. The constant displacement element (Fig. 3a), which means the same displacement over the entire element, will give lower accuracy, compared with the linear displacement element (Fig. 3b). For an opened fracture, the fracture width varies continuously from the crack tip to the crack center. The normal displacement Dz ðξ; ηÞ of each triangular element represents the fracture width physically. For constant displacement element, the triangular element only has one node, which is usually the gravity point (Fig. 4a). Hence, Di ðξ; ηÞ is not related to the coordinates (ξ, η, 0) in the element, i.e., constant displacement element, then the Eq. (4) can be expressed by [9,10,26]. Φi ðx; y; zÞ ¼
Di 1 ∬Δ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffidξ dη 8πð1 νÞ ðx ξÞ2 þ ðy ηÞ2 þ z2
i ¼ x; y; z
ð5Þ
By contrast with constant displacement element, a linear triangular element (Fig. 4b) has three nodes at its vertexes, P1 (x1,y1,0), P2 (x2,y2,0), P3 (x3,y3,0). The displacements in each element form a plane in the three dimensional space, and the equation of this plane reads ξx η y1 Di ðξ; ηÞ Di ðP 1 Þ 1 x21 y21 Di ðP 2 Þ Di ðP 1 Þ ¼ 0 i ¼ x; y; z ð6Þ x31 y31 Di ðP 3 Þ Di ðP 1 Þ Di ðξ; ηÞ in Eq. (6) is solved as h i þ ðη y1 Þðx21 x31 Þ Di ðP 1 Þ Di ðξ; ηÞ ¼ 1 ðξ x1 Þðy31 x21yy21 Þ x31 y 31
þ
21
ðξ x1 Þy31 ðη y1 Þx31 ðP2 Þ ðη y1 Þx21 ðξ x1 Þy21 ðP 3 Þ Di þ Di x21 y31 x31 y21 x21 y31 x31 y21
ð7Þ
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
91
Fig. 3. The side views of a 3D fracture in constant and variable displacement element. (a) Constant displacement element. (b) Linear displacement element.
Fig. 4. Three-dimensional fracture surface divided into triangular element (red dot is the node). (a) Constant displacement triangular element. (b) Linear displacement triangular element. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
where xmn ¼ xm xn ; ymn ¼ ym yn m; n ¼ 1; 2; 3. Di ðP 1 Þ ; Di ðP2 Þ ; Di ðP 3 Þ represent the displacement discontinuities of the three nodes in a triangular element respectively. Also, the Eq. (7) can be expressed by the linear combinations of [1, ξ, η] h Di ðξ; ηÞ ¼ Di ðP 1 Þ
Di ðP 2 Þ
i Di ðP3 Þ U ½C 33 U 1
ξ
η
T
i ¼ x; y; z ð8Þ
The coefficient matrix 2
½C 33 ¼
2x1 y3 þ x2 y3 x3 y2 1 6 x y y x 4 31 1 31 1 ðx21 y31 x31 y21 Þ y21 x1 x21 y1
y21 y31
x31 x21
y31 y21
x31 x21
3 7 5
ð9Þ When a modeled fracture boundary is divided into triangular elements, the positions, shapes, areas and directions of those elements are generally different from each other. Therefore, the displacement discontinuity components of each element must be represented with respect to its own coordinate system. The local coordinate system corresponding to each triangular element is defined as shown in Fig. 5. Note that P1, P2 and P3 are in a
! counterclockwise direction with P1 on origin and the vector P 1 P 2 along the ξ-axis. The triangular element locates inside the plane ðξ; o; ηÞ. Hence, x1 ¼y1 ¼y2 ¼ 0 in this local coordinate system and the coefficient matrix [C] is simplified as 2 x2 y3 1 6 ½C 33 ¼ 40 x2 y 3 0
y3 y3 0
x32
3
x3 7 5 x2
ð10Þ
Fig. 5. The local coordinates of a triangular element [9,10].
Substituting Eq. (8) into Eq. (4), the double integral depends on the basic integral as follows: 8 f ðx; y; zÞ ¼ ∬Δ 1ρdξ dη > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < gðx; y; zÞ ¼ ∬Δ ρξdξ dη ρ ¼ ðx ξÞ2 þ ðy ηÞ2 þ z2 ð11Þ > > : hðx; y; zÞ ¼ ∬Δ ηdξ dη ρ
The analytical expression of f(x,y,z) in Eq. (11) was given by Kuriyama et al. [9,10], which is not listed in this paper. When x1 ¼y1 ¼ y2 ¼0, the analytical expression of g(x,y,z) and h(x,y,z) in Eq. (11) are given in Appendix A.
3. Solution procedure 3.1. Basic equations in 3D-DDM The displacement and stress components, acting on an arbitrary point in a three dimensional space, induced by a linear displacement discontinuity element can be determined by using Eqs. (2) and (3). Linear superposition of the displacement and
92
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
stress components is must when the number of displacement discontinuity elements is more than one. The procedure to determine these components is stated in detail as follows. Let the modeled fracture surface be divided into N triangular elements, M nodes in total and L nodes on the crack-edge, as shown in Fig. 4b. The displacement and stress components induced at the gravity of the ith element by the displacement discontinuity of the jth element can be given by Eqs. (2) and (3). But they are represented by their local coordinates of the jth element. Stress and displacement superposition must be conducted after coordinate transformations from the jth local coordinate to the ith local coordinate. Therefore, the z-stress component applied on the plane of the ith element is given by σ zz ðiÞ ¼
N X h T ij Sji Djx ðP 1 Þ
Djx ðP 2 Þ
Djx ðP 3 Þ
j¼1
Djy ðP 1 Þ Djy ðP 2 Þ Djy ðP3 Þ Djz ðP 1 Þ Djz ðP2 Þ Djz ðP 3 Þ
iT
i; j ¼ 1; 2; :::; N
ð12Þ
where [Tij] is the matrix for coordinate transformations from the jth local coordinate to the ith local coordinate and [Sji] are the coefficients in Eq. (3) with respect to the displacement discontinuity components on the jth element. In general, Djk ðP1 Þ ; Djk ðP 2 Þ ; Djk ðP 3 Þ k ¼ x; y; z are not independent variables because of the adjacent triangular elements will share one or two nodes (Fig. 4b), but the total number of nodes is M. Therefore, there are 3M variables and N equations in Eq. (12). By contrast, the total number of linear equations and variables both equal to the number of triangular in constant triangular element [9,10]. 3.2. Boundary conditions To simulate the fracture by using the boundary element method, the boundary conditions usually include both external (far-field stress) and internal (acting on the fracture surface directly) conditions [6]. For instance, the internal boundary condition of a hydraulic fracture in the subsurface is the fluid pressure, while the external boundary condition is the in-situ stress field, which is represented by the global coordinate. Coordinate transformation from the global coordinate to each local coordinate on the boundary should be conducted when using the boundary conditions [9,10]. 3.3. Crack tip The real displacements at the crack tip are zero. By contrast, the displacements of elements at the crack-edge in constant triangular element are usually not zeros because its nodes are on the gravity point of crack-tip triangular element, not on the crack-tip [9,10], as shown in Fig. 6a. This fact makes the constant displacement
discontinuity element hard to get reliable numerical results on the crack tip element. Hence, Li et al. [11] extended the constant triangular element to the square-root crack-tip triangular element in accordance with the square-root displacement behavior near the crack-tip and obtained a high accuracy in estimating the stress intensity factors (SIFs) at the crack-tip. In this article, we do not employ special tip element as what Li et al. [11] did. However, the linear displacement discontinuity element is made of three vertices of each triangular. The discontinous displacements (Dx, Dy, Dz) of the nodes on the crack-edge are zero, as shown in Fig. 6b. In addition, the total number of the variables in Eq. (12) decreases into 3(M–L). In Fig. 6a, (Dx, Dy, Dz) on the gravity point of triangular represent the discontinous displacements in constant element. (Dx, Dy, Dz) can be transformed into the displacements in crack-tip coordinate system by 2 tip 3 2 3 32 Dx Dx cos α sin α 0 6 tip 7 6 6 7 6 Dy 7 ¼ 4 sin α cos α 0 7 ð13Þ 54 D y 5 4 5 Dtip z
0
0
1
Dz
tip tip where ðDtip x ; Dy ; Dz Þ represent the discontinous displacements in crack-tip coordinate syetem; α represents the angle which starts
!
! anticlockwise from P 1 P 2 to P G P V . In linear elastic mechanics, SIFs of crack-tip are given by [28] 8 qffiffiffiffi 0 G 2π tip > > > K I ¼ 4ð1 νÞ γ Dz > > q ffiffiffiffi < tip K 0II ¼ 4ð1G νÞ 2π ð14Þ γ Dx > > qffiffiffiffi > > 0 G 2π tip > : K III ¼ 4 γ Dy
where K 0i ði ¼ I; II; III Þ represent the stress intensity factors in mode
! I, mode-II, mode-III respectively; γ represents the length of P G P V . Due to the special shape of crack-tip, a modified coefficient was proposed by [13] rffiffiffi 2 0 Ki ¼ ð15Þ K ði ¼ I; II; III Þ π i If the gap between the line P2–P3 and the real crack-tip is ignored, γ represents the radial distance from the node to crack-tip approximately in constant triangular element. Therefore, stress intensity factor in constant element can be calculated
triangular ! easily by substituting γ ¼ P G P V and Eq. (14) to Eq. (15). Note that this hypothesis is acceptable when the size of triangular is infinite small. By contrast to the constant element, the displacement discontinuities in linear triangular element are represented by the linear combination of the displacements on the three vertices of each triangular. In geometry, the displacement discontinuities on the
Fig. 6. Triangular element at the crack-tip. (a) Constant displacement element. (b) Linear displacement element. (PG represents the gravity point of crack-tip triangular element; the line PG–PV is perpendicular to the line P2–P3; ðxtip ; ytip ; ztip Þ represents the local coordinate system at the crack-tip.)
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
gravity of the linear triangular are the average values of those on the three vertices. Hence, the displacement discontinuities on the gravity point in crack-tip coordinate system of linear element are given by 2
3 2 Dtip cos α x 6 tip 7 16 6 Dy 7 ¼ 4 sin α 4 5 3 0 Dtip z
sin α cos α 0
3 32 ðP1 Þ Dx þ Dx ðP2 Þ þDx ðP 3 Þ 6 7 D ðP 1 Þ þ D ðP 2 Þ þ D ðP3 Þ 7 7 0 56 y y 4 y 5 ðP 1 Þ ðP 2 Þ ðP 3 Þ 1 Dz þ Dz þ Dz 0
ð16Þ
! Substituting the Eqs. (14), (16) and γ ¼ P G P V into Eq. (15), SIFs at crack-tip in linear displacement element can be easily obtained. SIFs at crack-tip in constant and linear displacement elements will be compared with the analytical solution of a pressurized pennyshape crack in Section 4.2.
93
The analytical integral of Eq. (19) depends on the basic integral as follows: 8 2 > lim∂ f ðx;y;zÞ ¼ ∬Δ ℜ13 dξ dη > > z-0 ∂z2 > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < 2 lim∂ gðx;y;zÞ ¼ ∬Δ ℜξ3 dξ dη ℜ ¼ ðx ξÞ2 þ ðy ηÞ2 2 ð20Þ ∂z z-0 > > > > ∂2 hðx;y;zÞ η > : lim ∂z2 ¼ ∬Δ ℜ3 dξ dη z-0
The analytical expression of the first equation in Eq. (20) was given by Kuriyama et al. [9,10], which is not listed in this paper. The analytical expression of the rest equations in Eq. (20) are given in Appendix B. Substituting Dx ¼ Dy ¼0 and the boundary conditions into Eq. (12) yields P net ðiÞ ¼
N X
h ðP Þ T ij Sji Djz 1
Djz ðP 2 Þ
Djz ðP 3 Þ
iT
i; j ¼ 1; 2; …; N
j¼1
ð21Þ
4. Validation for a pressurized penny-shape crack 4.1. Width validation The radius of a penny-shape crack is R and the internal fluid pressure is Pnet, as shown in Fig. 7. This crack is buried inside an infinite isotropically linear elastic medium. Young's modulus and Possion's ratio of this medium are 40 GPa and 0.25, respectively. The remote stress field is zero, i.e., the traction of external boundary is free. The analytical solution of the crack width is given by Sneddon [27]. 8 1 v2 P net pffiffiffiffiffiffiffiffiffiffiffiffiffiffi WðrÞ ¼ R2 r 2 ; ð0 or oRÞ ð17Þ πE where Pnet is the fluid pressure, MPa; R is the radius of the circle crack, m; r is the distance from (x,y) to circle center, m. E and v are Young's modulus and Possion's ratio, respectively. The penny-shaped crack model is divided into triangular mesh. Two examples of triangular mesh on a penny-shaped crack are given in Fig. 8. In the penny-shaped crack model, the shear displacements (Dx and Dy) and the shear stresses acting on the crack surface are zero. The normal stress in Eq. (3) can be simplified as σ zz ¼ 2G zΦx;xzz zΦy;yzz þ Φz;zz zΦz;zzz ð18Þ Since this crack is a planar crack, the coordinate z ¼0. Substituting z¼ 0 and Eq. (4) into Eq. (18) gives σ zz ¼
G ∂2 Dz ðξ; ηÞ lim ∬Δ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffidξdη 4πð1 νÞz-0∂z2 ðx ξÞ2 þ ðy ηÞ2 þ z2
i ¼ x; y; z
ð19Þ
Solving Eq. (21), we can obtain the normal displacement Dz (fracture width in physical meaning) at each node. The displacements by 3D DDM in this paper were compared with the analytical solutions, as shown in Fig. 9. 4.2. SIF validation According to the linear elastic fracture mechanics, the analytical expression of KI in a pressurized penny-shaped crack is given by rffiffiffi R K I ¼ 2P net ð22Þ π Eqs. (13)–(16) give the method in detail to determined SIFs in constant and linear displacement element of the boundary element method. The analytical and numerical KI are compared in Fig. 10. In Fig. 10a, the average relative error is 1.9167% for constant element, while 1.4226% for linear element. In Fig. 10b, the average value of the relative errors is 3.2042% for constant element, while 1.5158% for linear element. It is easy to find that the numerical solution on KI in linear triangular element is more precise than that in constant triangular element. 4.3. Stress validation The analytical solution of stress field in the neighborhood of a penny-shaped crack in an infinite elastic solid was given by Sneddon [27].
8 1=2 2 > σ > v þ 12 arc sin Rr > Pnetr ¼ 2π Rr 2 1 < r 4 R; z ¼ 0; ð23Þ 1=2 > σ zz 2 r2 R > > ¼ 1 arc sin : Pnet π R2 r where σr is the radial stress, MPa; σθ is the hoop stress, MPa; σzz is the z-stress in cylindrical coordinate system, MPa. An arbitrary point P0 (x0, y0, z0) in the region adjacent to the penny-shaped crack, stress induced by a single triangular element on point P0 can be expressed by the following equation, according to Eq. (3). ( σ xx ¼ 2G Φz;zz þ ð1 2νÞΦz;yy zΦz;xxz ð24Þ σ zz ¼ 2G Φz;zz zΦz;zzz After the linear superposition of stress, the total stress acting on point P0 by the entire penny-shaped crack can be given as σ zz ðP0 Þ ¼
N h ih ih X ðP Þ ðP 0 Þ Djz ðP 1 Þ Sj Tj 0
Djz ðP2 Þ
Djz ðP3 Þ
iT
i ¼ 1; 2; …; N
j¼1
Fig. 7. Penny-shaped crack model.
ð25Þ
94
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
Fig. 8. Example of the triangular mesh on a penny-shaped crack (N, M and L represent the number of triangular, vertex and edge-vertex, respectively). (a) N ¼258, M ¼146, L ¼ 32 and (b) N ¼488, M ¼269, L ¼ 48.
Fig. 10. Relative error of KI on a pressurized penny-shaped crack. (a) N ¼258, M¼ 146, G ¼ 32 and (b) N ¼488, M¼ 269, G ¼ 48.
Fig. 9. Comparisons of analytical and numerical width. (a) N ¼ 258, M ¼146, L ¼32 and (b) N ¼488, M¼ 269, L ¼ 48.
where σ zz ðP0 Þ is the normal stress on point P0 in global coordinate h i ðP Þ system; T j 0 is the matrix for coordinate transformation from the h i ðP Þ jth local coordinate to global coordinate; Sj 0 are the coefficients
of P0 in Eq. (3) with respect to the displacement discontinuity components on the jth element. Note that the other stress components acting on point P0 have the similar equations as Eq. (25). With the axial symmetry of this crack and its boundary conditions, σxx and σzz on x-axis in Eq. (24) should be equal to σr and σzz in Eq. (23), respectively, even though they are in different coordinate system. Their comparisons prove a high accuracy between the numerical solutions with analytical solutions, as
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
∬Δ ξ ρ xdξ dη ¼ 1
AI ¼
2
R y3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R y pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AI η2 þ BI η þ C I dη 0 3 AII η2 þ BII η þ C II dη 0
þ 1;
BI ¼
k2 1 AII ¼ 2 þ1; k1
hðx; y; zÞ ¼ ∬Ω
95
2ðx2 xÞ 2y; k2
BII ¼
2x 2y; k1
C I ¼ ðx2 xÞ2 þ y2 þ z2 C II ¼ x2 þ y2 þ z2
η y dξ dη þ yf ðx; y; zÞ ρ
ρ¼
ðA2Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx ξÞ2 þ ðy ηÞ2 þ z2 ðA3Þ
∬Δ η ρ ydξ dη ¼
ffi R x pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R x3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AIII ξ2 þ BIII ξ þ C III dξ þ x32 AIV ξ2 þ BIV ξ þ C IV dξ x1 Z
þ 2
AIII ¼ k1 þ 1;
x1 x2
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AV ξ2 þBV ξ þC V dξ
BIII ¼ 2x 2k1 y;
C III ¼ x2 þ y2 þ z2
2
AIV ¼ k2 þ 1; BIV ¼ 2x 2k2 ðk2 x2 þ yÞ; 2
2
Fig. 11. Comparisons of the analytical solutions with numerical solutions (x¼ ρ, y¼ 0, z ¼0).
shown in Fig. 11. Note that the error between analytical and numerical solution becomes larger when the point gets closer to pffiffiffi crack tip. This is because 1= r stress singularity exists at the adjacent domain of crack-tip. The stress in crack-tip domain is dominated by the SIF, not dominated by the linear superposition of the stresses induced by each element in the boundary. When r/R1 þ , σxx- þ1 and σzz-þ 1 according to Eq. (23). However, the distance between the point P0 and the gravity point PG of triangular element is finite, i.e., the integrals in Eqs. (4), (5) and (11) are non-singular. Hence, the numerical solutions of σr and σzz close to crack-tip are both smaller than their analytical solutions.
C IV ¼ x þ ðk2 x2 þ yÞ þ z2 AV ¼ 1; BV ¼ 2x; C V ¼ x2 þ y2 þ z2
ðA4Þ
k1 ¼ y3 =x3 ;
ðA5Þ
Eqs. A2 and A4 have the following form: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ax2 þ bx þ cdx ¼ 2ax4aþ b ax2 þ bx þ c þ pffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4ac b pffiffiffiffiffi ln2axþ b þ 2 a ax2 þ bx þ c þ Const 3 8 a 2
(1) A novel linear displacement triangular element with three nodes on its vertices is proposed. If the displacement discontinuities on these three nodes keep the same, this linear displacement triangular element will be simplified as the constant displacement triangular. (2) The solution procedure in this paper can be applied to determine the displacement and stress field around a three-dimensional fracture inside the infinite linear elastic solid. This procedure takes a shorter time and requires less elements than the constant DDM when achieving the same accuracy. If element subdivision is required for a reasonable accuracy then the analytical integration is computationally more efficient than the numerical integration.
∂2 gðx; y; zÞ xξ ∂2 f ðx; y; zÞ ¼ ∬Ω 3 dξ dη þ xlim 2 z-0 z-0 ∂z ∂z2 ℜ
ℜ¼
lim
AV I ¼
R y3 0
AVI η2 þBVI η þC VI
1 þ 1; k2 2
AV II ¼
1 2 þ 1; k1
ðA6Þ
1=2
BVI ¼ 2 x2k2 x y ; BVII ¼
2x 2y; k1
dη
R y3 0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx ξÞ2 þ ðy ηÞ2
ðB1Þ AVII η2 þ BV II η þC VII
C V I ¼ ðx2 xÞ2 þ y2 C V II ¼ x2 þ y2
∂2 hðx; y; zÞ yη ∂2 f ðx; y; zÞ ¼ ∬Ω 3 dξ dηþ ylim z-0 z-0 ∂z2 ∂z2 ℜ lim
ℜ¼
1=2
dη
ðB2Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx ξÞ2 þ ðy ηÞ2
ðB3Þ
1=2 Rx dξ ∬Ω yℜ3ηdη dξ ¼ x13 AIII ξ2 þ BIII ξ þ C III Z x2 Z x1 þ AIX ξ2 þ BIX ξ þ C IX dξ þ AX ξ2 þ BX ξ þ C X dξ x3
x2
2
AV III ¼ k1 þ 1; 2
AIX ¼ k2 þ 1; AX ¼ 1;
Acknowledgments This research was supported by the National Natural Science Foundation of China (51234006; 51325402; 51221003) and the China Scholarship Council (201406440009). Thanks to China Scholarship Council for supporting the first author of this paper to research in Georgia Institute of Technology, USA.
ða 4 0Þ
Appendix B. 2nd order partial derivatives
∬Ω xℜ3ξdξdη ¼
5. Conclusions
k2 ¼ y3 =ðx3 x2 Þ
BVIII ¼ 2x 2k1 y;
C V III ¼ x2 þ y2
BIX ¼ 2x 2k2 ðk2 x2 þ yÞ; BX ¼ 2x;
C IX ¼ x2 þ ðk2 x2 þ yÞ2 C X ¼ x2 þ y2 ðB4Þ
Eqs. B2 and B4 have the following form: Z
pffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1=2 1 dx ¼ pffiffiffi ln2ax þ b þ 2 a ax2 þ bx þ c þ Const ax þ bxþ c a
ða 4 0Þ
ðB5Þ References
Appendix A. Double integral of linear triangular element qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ξx gðx; y; zÞ ¼ ∬Ω dξ dη þ xf ðx; y; zÞ ρ ¼ ðx ξÞ2 þ ðy ηÞ2 þ z2 ρ ðA1Þ
[1] Ariza MP, Saez A, Dominguez J. A singular element for three-dimensional fracture mechanics analysis. Eng Anal Bound Elem 1997;20:275–85. [2] Berry DS. An elastic treatment of ground movement due to mining-I. Isotropic ground. J Mech Phys Solids 1960;8:280–92.
96
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
[3] Cheng W, Jin Y, Chen M, et al. A criterion for identifying hydraulic fractures crossing natural in 3D space. Pet Explor Dev 2014;41(3):336–40. [4] Cheng W, Jin Y, Lin Q, et al. Experimental investigation about influence of preexisting fracture on hydraulic fracture propagation under tri-axial stresses. Geotech Geol Eng 2015;37(3):413–8. [5] Crouch SL. Solution of plane elasticity problems by the displacement discontinuity method. Int J Numer Methods Eng 1976;10:301–43. [6] Crouch SL, Starfield AM. Boundary element methods in solid mechanics. London: Goerge Allen and Unwin Publishers; 1983. [7] Dyka CT, Remondi Millwater HR. 3D boundary elements and the integration of strong singularities. Comput Struct 1991;39(5):513–23. [8] Davey K, Hinduja S. Analytical integration of linear three-dimensional triangular elements in BEM. Appl Math Model 1989;13:450–61. [9] Kuriyama K, Mizuta Y. Three-dimensional elastic analysis by the displacement discontinuity method with boundary division into triangular leaf elements. Int J Rock Mech Min Sci Geomech Abstr 1993;30(2):111–23. [10] Kuriyama K, Mizuta Y, Mozumi H, Watanabe T. Three-dimensional elastic analysis by the boundary element method with analytical integrations over triangular leaf elements. Int J Rock Mech Min Sci Geomech Abstr 1995;32 (1):77–83. [11] Li H, Liu CL, Mizuta Y, et al. Crack edge element of three-dimensional displacement discontinuity method with boundary division into triangular leaf elements. Commun Numer Methods Eng 2001;17:365–78. [12] Lin’kov AM, Zubkov VV, Kheib MA. A method of solving three-dimensional problems of seam workings and geological faults. J Min Sci 1997;33 (4):295–315. [13] Liu G, Tu J, Zhang J. Solution for crack propagation in multi-crack body by displacement–discontinuous boundary element method. J Tsinghua Univ (Sci Technol) 1996;36:59–64. [14] Milroy J, Hinduja S, Davey K. The elastostatic three dimensional boundary element method: analytical integration for linear isoparametric triangular elements. Appl Math Model 1997;21:763–82. [15] Mogilevskaya SG, Nikolskiy DV. The use of complex integral representations for analytical evaluation of three-dimensional BEM integrals-potential and elasticity problems. Q J Mech Appl Math 2014;67(3):505–23.
[16] Napier JAL, Malan DF. The computational analysis of shallow depth tabular mining problems. J S Afr Inst Min Metall 2007;107:725–42. [17] Nikolski DV. Three-dimensional boundary element modeling of fractures under gravity load. In: Proceedings of the 46th American Rock Mechanics Association. 12-450. Chicago, USA; 24–27 June 2012. [18] Nikolski DV, Mogilevskaya SG, Labuz JF. Complex variables boundary element analysis of the three-dimensional crack problems. Eng Anal Bound Elem 2013;37:1532–44. [19] Nikolski DV and Mogilevskaya SG. Complex variables BEM modeling of three -dimensional crack problems. In: Proceedings of the 48th American Rock Mechanics Association symposium on US rock mechanics/geomechanics. Minneapolis, MN, USA; 1–4 June 2014. [20] Nikolskiy DV, Mogilevskaya SG, Labuz JF. Boundary element analysis of nonplanar three-dimensional cracks using complex variables. Int J Rock Mech Min Sci 2015;76:44–54. [21] Fata SN. Explicit expressions for the three-dimensional boundary integrals in linear elasticity. J Comput Appl Math 2011;235:4480–95. [22] Qin TY, Tang RJ. Three dimensional crack problem analysis using boundary element method with finite-part integrals. Int J Fract 1997;84:191–202. [23] Salamon MG. Elastic analysis of displacements and stresses induced by mining of seam or roof deposits. Part IV. J S Afr Inst Min Metall 1964;65:319–38. [24] Sokolnikoff. Mathematical theory of elasticity. New York: Mcgraw-hill book company; 1956. [25] Shou KJ, Siebrits E, Crouch SL. A higher order displacement discontinuity method for three-dimensional elastostatic problems. Int J Rock Mech Min Sci 1997;34(2):317–22. [26] Shou KJ. A three-dimensional hydrid boundary element method for non-linear analysis of a weak plane near an underground excavition. Tunn Undergr Space Technol 2000;15(2):215–26. [27] Sneddon IN. The distribution of stress in the neighborhood of a crack in an elastic solid. Proc R Soc Lond Ser: A Math Phys Sci 1946;46(1):187–229. [28] Xiao HT, Yue ZQ. A three-dimentional displacement discontinuity method for crack problems in layered rocks. Int J Rock Mech Min Sci 2011;48:412–20. [29] Wang C, Khoo BC. An indirect boundary element method for threedimensional explosion bubbles. J Comput Phys 2004;194:451–80.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A novel linear triangular element of a three-dimensional displacement discontinuity method Wan Cheng a,b,c, Yan Jin a,b,n, Hong Li d, Mian Chen a,b a
College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China State Key Laboratory of Petroleum Resources and Prospecting, Beijing 102249, China c School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA d Dalian University of Technology, Liaoning, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 2 March 2015 Received in revised form 27 April 2015 Accepted 30 April 2015
Since only the boundary of the domain requires discretization, the boundary element method (BEM) is very efficient for the semi-infinite or infinite rock-related engineering problems, e.g., hydraulic fracturing in reservoir stimulation and rock cutting during excavation. A real fracture in the solid is usually of an arbitrary geometry in three dimensions, which usually requires a three-dimensional displacement discontinuity method (3D DDM) to determine the deformation and stress field in order to achieve reliable results. However, the use of 3D DDM with triangular elements is limited by the singularities of the integral either within or nearby the domain. In this paper, a novel linear triangular element with three nodes on its vertices is proposed. The analytical integral expressions of this linear triangular element are also theoretically derived. A solution procedure is also described which can be applied to determine the displacement and stress field around a three-dimensional fracture inside the infinite solid. The accuracy of these results are compared with the analytical solutions of the displacements and stresses induced by a pressurized penny-shaped. This procedure takes a shorter time and requires less elements than the usual constant DDM when achieving the same accuracy. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Linear triangular element Boundary element method Displacement discontinuity method Three-dimensional fracture Hydraulic fracturing
1. Introduction Many important practical problems in science and engineering can be reduced into mathematical models that belong to a class of problems known as boundary value problems [5,6]. These are all characterized by a region of interest R enclosed within a boundary C. If R is a two dimensional region, C will be a bounding contour. If R is a three dimensional region then C is a bounding surface. For a fracture or cavity inside an infinite region, C will be the fracture or cavity surface. Only the boundary C is divided into elements in BEM, which is totally different from the network of elements for the entire region R in a finite element method (Fig. 1). The displacement discontinuity in BEM can be considered as a relative displacement between the two surfaces of a plane crack [24,2,23]. The 2D DDM of constant displacement element in linear elastic mechanics was systematically developed by Crouch [5] who proposed a general numerical algorithm to solve various boundary value problems. The 3D DDM is based on the analytical solution to the problem of a discontinuous displacement over a finite area in
n Corresponding author at: College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China. E-mail address: [email protected] (W. Cheng).
http://dx.doi.org/10.1016/j.enganabound.2015.04.020 0955-7997/& 2015 Elsevier Ltd. All rights reserved.
an infinite region. For a crack of arbitrary non-planar shape in an infinite three-dimensional space, the crack is usually required to be divided into triangular elements. Even though the strong singularities in boundary integral equation [7,29,1,21,22], the 3D DDM on the triangular element is improved from constant element to higher order element. Kuriyama et al. [9,10] proposed the 3D DDM with the boundary surface division into constant triangular elements and briefly presented the solution procedure. Davey and Hinduja [8] established the integral solutions for triangular elements applied to potential problems. Based on Davey's idea, Milroy et al. [14] obtained the general solution by a linear combination of the three integral solutions on the nodes of a triangular element. Napier and Malan [16] and Shou et al. [25] described a higher order displacement discontinuity method in conjunction with shape functions. Nikolski et al. [17,18] presented a boundary element by decoupling the complex variables for the shear tractions and real variable equation for the normal traction for solving 3D elastostatic problems of fractured rock under gravity. Then the complex variables BEM is extended to incorporate higher order approximations of the displacement discontinuities Nikolski [19,20,15]. The recent development of multi-stage hydraulic fracturing technology in a horizontal well in petroleum reservoir requires hydraulic fracturing simulators which can consider the three-
90
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
To numerically implement the displacement discontinuity element in a boundary element method, we need the analytical solution for the displacement discontinuity over an area which is usually defined as a boundary element. The general form solution for a displacement discontinuity element can be expressed as follows [6,23,24,26]:
Fig. 1. The element in BEM and FEM.
8 u ¼ 2ð1 νÞΦx;z zΦx;xx zΦy;xy ð1 2νÞΦz;x þ zΦz;xz > < x uy ¼ zΦx;xy þ 2ð1 νÞΦy;z zΦy;yy ð1 2νÞΦz;y þ zΦz;yz > : u ¼ ð1 2νÞΦ zΦ z x;x x;xz þ ð1 2νÞΦy;y zΦy;yz þ 2ð1 νÞΦz;z zΦz;zz
ð2Þ 8 σ xx ¼ 2G 2Φx;xz zΦx;xxx þ 2νΦy;yz zΦy;xxy þ Φz;zz þ ð1 2νÞΦz;yy zΦz;xxz > > > > > σ yy ¼ 2G 2vΦx;xz zΦx;xyy þ 2Φy;yz zΦy;yyy þ Φz;zz þ ð1 2νÞΦz;xx zΦz;yyz > > > > < σ zz ¼ 2G zΦx;xzz zΦy;yzz þ Φz;zz zΦz;zzz σ xy ¼ 2G ð1 νÞΦx;yz zΦx;xxy þ ð1 νÞΦy;xz zΦy;xyy ð1 2νÞΦz;xy þ zΦz;xyz > > > > > σ yz ¼ 2G νΦx;xy zΦx;xyz þ Φy;zz þ νΦy;xx zΦy;yyz zΦz;yzz > > > > : σ xz ¼ 2GΦx;zz þ νΦx;yy zΦx;xxz νΦy;xy þ zΦy;xyz zΦz;xzz
ð3Þ
Fig. 2. Triangular element in the 3D DDM.
dimensional formation and stress field around the fracture [3,4]. Physically, any fracture in an infinite three-dimensional space should be in three dimensions. The 3D DDM can be applied to simulate this kind of engineering problem numerically by assuming that the formation is a linear elastic medium. For the strong singularities in boundary integral equation, researchers hope to improve both the computational efficiency and accuracy utilizing the analytical integration instead of the numerical integration. Therefore, a novel linear displacement triangular element with three nodes on its vertices is proposed in this paper. If the displacement discontinuities on these three nodes keep the same, this linear displacement triangular element will be simplified as the constant displacement triangular. The idea of this linear triangular element can also be applied to determine linear triangular element in a fictitious stress method [26]. The analytical integral solutions of this linear triangular element are also theoretically derived. A solution procedure is also briefly described which can be applied to determine the stresses field around a three-dimensional fracture inside the solid. The accuracy of these results is compared with the analytical solutions of displacements and stresses induced by a pressurized penny-shaped crack [27].
2. Three-dimensional displacement discontinuity element As shown in Fig. 2, the two surfaces of a crack are separated with relative displacements. The two surfaces can be distinguished by the positive side (z¼ 0 þ ) and negative side (z¼ 0 ). Crossing from one side to the other, the displacement discontinuities are given by the relative displacement in three orthogonal directions respectively [23,24]. 8 > < Dx ðx; yÞ ¼ ux ðx; y; 0 Þ ux ðx; y; 0 þ Þ Dy ðx; yÞ ¼ uy ðx; y; 0 Þ uy ðx; y; 0 þ Þ ð1Þ > : D ðx; yÞ ¼ u ðx; y; 0 Þ u ðx; y; 0 Þ z z z þ where Di (i¼ x,y,z) are the displacement discontinuities in x,y,z direction respectively. ui (i¼ x,y,z) are the displacements in x,y,z direction respectively.
where σ ij ; ui ði; j ¼ x; y; zÞ represent the stress and displacement of Kelvin's problem in three dimensional space respectively. G, v are the shear modulus and Possion's ratio of the isotropic medium, respectively. Φi;j ; Φi;jk ; Φi;jkl ði; j; k; l ¼ x; y; zÞ are the partial derivatives with respect to each corresponding suffix of the following function,Φi . The function, Φi , from a triangular element against a point can be represented by [12,15] Φi ðx; y; zÞ ¼
1 Di ðξ; ηÞ ffidξ dη ∬Δ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8πð1 νÞ ðx ξÞ2 þ ðy ηÞ2 þ z2
i ¼ x; y; z
ð4Þ
where (x,y,z) are the coordinates of an arbitrary point in an infinite 3D space. (ξ, η, 0) are the local coordinates of the loading point. Di ðξ; ηÞ represents the discontinuous displacement on point (ξ, η, 0). Because of the strong singularity of Eq. (4) when (x,y,z) goes to (ξ, η, 0), it is extremely difficult to get the analytical integration [7]. The constant displacement element (Fig. 3a), which means the same displacement over the entire element, will give lower accuracy, compared with the linear displacement element (Fig. 3b). For an opened fracture, the fracture width varies continuously from the crack tip to the crack center. The normal displacement Dz ðξ; ηÞ of each triangular element represents the fracture width physically. For constant displacement element, the triangular element only has one node, which is usually the gravity point (Fig. 4a). Hence, Di ðξ; ηÞ is not related to the coordinates (ξ, η, 0) in the element, i.e., constant displacement element, then the Eq. (4) can be expressed by [9,10,26]. Φi ðx; y; zÞ ¼
Di 1 ∬Δ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffidξ dη 8πð1 νÞ ðx ξÞ2 þ ðy ηÞ2 þ z2
i ¼ x; y; z
ð5Þ
By contrast with constant displacement element, a linear triangular element (Fig. 4b) has three nodes at its vertexes, P1 (x1,y1,0), P2 (x2,y2,0), P3 (x3,y3,0). The displacements in each element form a plane in the three dimensional space, and the equation of this plane reads ξx η y1 Di ðξ; ηÞ Di ðP 1 Þ 1 x21 y21 Di ðP 2 Þ Di ðP 1 Þ ¼ 0 i ¼ x; y; z ð6Þ x31 y31 Di ðP 3 Þ Di ðP 1 Þ Di ðξ; ηÞ in Eq. (6) is solved as h i þ ðη y1 Þðx21 x31 Þ Di ðP 1 Þ Di ðξ; ηÞ ¼ 1 ðξ x1 Þðy31 x21yy21 Þ x31 y 31
þ
21
ðξ x1 Þy31 ðη y1 Þx31 ðP2 Þ ðη y1 Þx21 ðξ x1 Þy21 ðP 3 Þ Di þ Di x21 y31 x31 y21 x21 y31 x31 y21
ð7Þ
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
91
Fig. 3. The side views of a 3D fracture in constant and variable displacement element. (a) Constant displacement element. (b) Linear displacement element.
Fig. 4. Three-dimensional fracture surface divided into triangular element (red dot is the node). (a) Constant displacement triangular element. (b) Linear displacement triangular element. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
where xmn ¼ xm xn ; ymn ¼ ym yn m; n ¼ 1; 2; 3. Di ðP 1 Þ ; Di ðP2 Þ ; Di ðP 3 Þ represent the displacement discontinuities of the three nodes in a triangular element respectively. Also, the Eq. (7) can be expressed by the linear combinations of [1, ξ, η] h Di ðξ; ηÞ ¼ Di ðP 1 Þ
Di ðP 2 Þ
i Di ðP3 Þ U ½C 33 U 1
ξ
η
T
i ¼ x; y; z ð8Þ
The coefficient matrix 2
½C 33 ¼
2x1 y3 þ x2 y3 x3 y2 1 6 x y y x 4 31 1 31 1 ðx21 y31 x31 y21 Þ y21 x1 x21 y1
y21 y31
x31 x21
y31 y21
x31 x21
3 7 5
ð9Þ When a modeled fracture boundary is divided into triangular elements, the positions, shapes, areas and directions of those elements are generally different from each other. Therefore, the displacement discontinuity components of each element must be represented with respect to its own coordinate system. The local coordinate system corresponding to each triangular element is defined as shown in Fig. 5. Note that P1, P2 and P3 are in a
! counterclockwise direction with P1 on origin and the vector P 1 P 2 along the ξ-axis. The triangular element locates inside the plane ðξ; o; ηÞ. Hence, x1 ¼y1 ¼y2 ¼ 0 in this local coordinate system and the coefficient matrix [C] is simplified as 2 x2 y3 1 6 ½C 33 ¼ 40 x2 y 3 0
y3 y3 0
x32
3
x3 7 5 x2
ð10Þ
Fig. 5. The local coordinates of a triangular element [9,10].
Substituting Eq. (8) into Eq. (4), the double integral depends on the basic integral as follows: 8 f ðx; y; zÞ ¼ ∬Δ 1ρdξ dη > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < gðx; y; zÞ ¼ ∬Δ ρξdξ dη ρ ¼ ðx ξÞ2 þ ðy ηÞ2 þ z2 ð11Þ > > : hðx; y; zÞ ¼ ∬Δ ηdξ dη ρ
The analytical expression of f(x,y,z) in Eq. (11) was given by Kuriyama et al. [9,10], which is not listed in this paper. When x1 ¼y1 ¼ y2 ¼0, the analytical expression of g(x,y,z) and h(x,y,z) in Eq. (11) are given in Appendix A.
3. Solution procedure 3.1. Basic equations in 3D-DDM The displacement and stress components, acting on an arbitrary point in a three dimensional space, induced by a linear displacement discontinuity element can be determined by using Eqs. (2) and (3). Linear superposition of the displacement and
92
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
stress components is must when the number of displacement discontinuity elements is more than one. The procedure to determine these components is stated in detail as follows. Let the modeled fracture surface be divided into N triangular elements, M nodes in total and L nodes on the crack-edge, as shown in Fig. 4b. The displacement and stress components induced at the gravity of the ith element by the displacement discontinuity of the jth element can be given by Eqs. (2) and (3). But they are represented by their local coordinates of the jth element. Stress and displacement superposition must be conducted after coordinate transformations from the jth local coordinate to the ith local coordinate. Therefore, the z-stress component applied on the plane of the ith element is given by σ zz ðiÞ ¼
N X h T ij Sji Djx ðP 1 Þ
Djx ðP 2 Þ
Djx ðP 3 Þ
j¼1
Djy ðP 1 Þ Djy ðP 2 Þ Djy ðP3 Þ Djz ðP 1 Þ Djz ðP2 Þ Djz ðP 3 Þ
iT
i; j ¼ 1; 2; :::; N
ð12Þ
where [Tij] is the matrix for coordinate transformations from the jth local coordinate to the ith local coordinate and [Sji] are the coefficients in Eq. (3) with respect to the displacement discontinuity components on the jth element. In general, Djk ðP1 Þ ; Djk ðP 2 Þ ; Djk ðP 3 Þ k ¼ x; y; z are not independent variables because of the adjacent triangular elements will share one or two nodes (Fig. 4b), but the total number of nodes is M. Therefore, there are 3M variables and N equations in Eq. (12). By contrast, the total number of linear equations and variables both equal to the number of triangular in constant triangular element [9,10]. 3.2. Boundary conditions To simulate the fracture by using the boundary element method, the boundary conditions usually include both external (far-field stress) and internal (acting on the fracture surface directly) conditions [6]. For instance, the internal boundary condition of a hydraulic fracture in the subsurface is the fluid pressure, while the external boundary condition is the in-situ stress field, which is represented by the global coordinate. Coordinate transformation from the global coordinate to each local coordinate on the boundary should be conducted when using the boundary conditions [9,10]. 3.3. Crack tip The real displacements at the crack tip are zero. By contrast, the displacements of elements at the crack-edge in constant triangular element are usually not zeros because its nodes are on the gravity point of crack-tip triangular element, not on the crack-tip [9,10], as shown in Fig. 6a. This fact makes the constant displacement
discontinuity element hard to get reliable numerical results on the crack tip element. Hence, Li et al. [11] extended the constant triangular element to the square-root crack-tip triangular element in accordance with the square-root displacement behavior near the crack-tip and obtained a high accuracy in estimating the stress intensity factors (SIFs) at the crack-tip. In this article, we do not employ special tip element as what Li et al. [11] did. However, the linear displacement discontinuity element is made of three vertices of each triangular. The discontinous displacements (Dx, Dy, Dz) of the nodes on the crack-edge are zero, as shown in Fig. 6b. In addition, the total number of the variables in Eq. (12) decreases into 3(M–L). In Fig. 6a, (Dx, Dy, Dz) on the gravity point of triangular represent the discontinous displacements in constant element. (Dx, Dy, Dz) can be transformed into the displacements in crack-tip coordinate system by 2 tip 3 2 3 32 Dx Dx cos α sin α 0 6 tip 7 6 6 7 6 Dy 7 ¼ 4 sin α cos α 0 7 ð13Þ 54 D y 5 4 5 Dtip z
0
0
1
Dz
tip tip where ðDtip x ; Dy ; Dz Þ represent the discontinous displacements in crack-tip coordinate syetem; α represents the angle which starts
!
! anticlockwise from P 1 P 2 to P G P V . In linear elastic mechanics, SIFs of crack-tip are given by [28] 8 qffiffiffiffi 0 G 2π tip > > > K I ¼ 4ð1 νÞ γ Dz > > q ffiffiffiffi < tip K 0II ¼ 4ð1G νÞ 2π ð14Þ γ Dx > > qffiffiffiffi > > 0 G 2π tip > : K III ¼ 4 γ Dy
where K 0i ði ¼ I; II; III Þ represent the stress intensity factors in mode
! I, mode-II, mode-III respectively; γ represents the length of P G P V . Due to the special shape of crack-tip, a modified coefficient was proposed by [13] rffiffiffi 2 0 Ki ¼ ð15Þ K ði ¼ I; II; III Þ π i If the gap between the line P2–P3 and the real crack-tip is ignored, γ represents the radial distance from the node to crack-tip approximately in constant triangular element. Therefore, stress intensity factor in constant element can be calculated
triangular ! easily by substituting γ ¼ P G P V and Eq. (14) to Eq. (15). Note that this hypothesis is acceptable when the size of triangular is infinite small. By contrast to the constant element, the displacement discontinuities in linear triangular element are represented by the linear combination of the displacements on the three vertices of each triangular. In geometry, the displacement discontinuities on the
Fig. 6. Triangular element at the crack-tip. (a) Constant displacement element. (b) Linear displacement element. (PG represents the gravity point of crack-tip triangular element; the line PG–PV is perpendicular to the line P2–P3; ðxtip ; ytip ; ztip Þ represents the local coordinate system at the crack-tip.)
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
gravity of the linear triangular are the average values of those on the three vertices. Hence, the displacement discontinuities on the gravity point in crack-tip coordinate system of linear element are given by 2
3 2 Dtip cos α x 6 tip 7 16 6 Dy 7 ¼ 4 sin α 4 5 3 0 Dtip z
sin α cos α 0
3 32 ðP1 Þ Dx þ Dx ðP2 Þ þDx ðP 3 Þ 6 7 D ðP 1 Þ þ D ðP 2 Þ þ D ðP3 Þ 7 7 0 56 y y 4 y 5 ðP 1 Þ ðP 2 Þ ðP 3 Þ 1 Dz þ Dz þ Dz 0
ð16Þ
! Substituting the Eqs. (14), (16) and γ ¼ P G P V into Eq. (15), SIFs at crack-tip in linear displacement element can be easily obtained. SIFs at crack-tip in constant and linear displacement elements will be compared with the analytical solution of a pressurized pennyshape crack in Section 4.2.
93
The analytical integral of Eq. (19) depends on the basic integral as follows: 8 2 > lim∂ f ðx;y;zÞ ¼ ∬Δ ℜ13 dξ dη > > z-0 ∂z2 > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < 2 lim∂ gðx;y;zÞ ¼ ∬Δ ℜξ3 dξ dη ℜ ¼ ðx ξÞ2 þ ðy ηÞ2 2 ð20Þ ∂z z-0 > > > > ∂2 hðx;y;zÞ η > : lim ∂z2 ¼ ∬Δ ℜ3 dξ dη z-0
The analytical expression of the first equation in Eq. (20) was given by Kuriyama et al. [9,10], which is not listed in this paper. The analytical expression of the rest equations in Eq. (20) are given in Appendix B. Substituting Dx ¼ Dy ¼0 and the boundary conditions into Eq. (12) yields P net ðiÞ ¼
N X
h ðP Þ T ij Sji Djz 1
Djz ðP 2 Þ
Djz ðP 3 Þ
iT
i; j ¼ 1; 2; …; N
j¼1
ð21Þ
4. Validation for a pressurized penny-shape crack 4.1. Width validation The radius of a penny-shape crack is R and the internal fluid pressure is Pnet, as shown in Fig. 7. This crack is buried inside an infinite isotropically linear elastic medium. Young's modulus and Possion's ratio of this medium are 40 GPa and 0.25, respectively. The remote stress field is zero, i.e., the traction of external boundary is free. The analytical solution of the crack width is given by Sneddon [27]. 8 1 v2 P net pffiffiffiffiffiffiffiffiffiffiffiffiffiffi WðrÞ ¼ R2 r 2 ; ð0 or oRÞ ð17Þ πE where Pnet is the fluid pressure, MPa; R is the radius of the circle crack, m; r is the distance from (x,y) to circle center, m. E and v are Young's modulus and Possion's ratio, respectively. The penny-shaped crack model is divided into triangular mesh. Two examples of triangular mesh on a penny-shaped crack are given in Fig. 8. In the penny-shaped crack model, the shear displacements (Dx and Dy) and the shear stresses acting on the crack surface are zero. The normal stress in Eq. (3) can be simplified as σ zz ¼ 2G zΦx;xzz zΦy;yzz þ Φz;zz zΦz;zzz ð18Þ Since this crack is a planar crack, the coordinate z ¼0. Substituting z¼ 0 and Eq. (4) into Eq. (18) gives σ zz ¼
G ∂2 Dz ðξ; ηÞ lim ∬Δ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffidξdη 4πð1 νÞz-0∂z2 ðx ξÞ2 þ ðy ηÞ2 þ z2
i ¼ x; y; z
ð19Þ
Solving Eq. (21), we can obtain the normal displacement Dz (fracture width in physical meaning) at each node. The displacements by 3D DDM in this paper were compared with the analytical solutions, as shown in Fig. 9. 4.2. SIF validation According to the linear elastic fracture mechanics, the analytical expression of KI in a pressurized penny-shaped crack is given by rffiffiffi R K I ¼ 2P net ð22Þ π Eqs. (13)–(16) give the method in detail to determined SIFs in constant and linear displacement element of the boundary element method. The analytical and numerical KI are compared in Fig. 10. In Fig. 10a, the average relative error is 1.9167% for constant element, while 1.4226% for linear element. In Fig. 10b, the average value of the relative errors is 3.2042% for constant element, while 1.5158% for linear element. It is easy to find that the numerical solution on KI in linear triangular element is more precise than that in constant triangular element. 4.3. Stress validation The analytical solution of stress field in the neighborhood of a penny-shaped crack in an infinite elastic solid was given by Sneddon [27].
8 1=2 2 > σ > v þ 12 arc sin Rr > Pnetr ¼ 2π Rr 2 1 < r 4 R; z ¼ 0; ð23Þ 1=2 > σ zz 2 r2 R > > ¼ 1 arc sin : Pnet π R2 r where σr is the radial stress, MPa; σθ is the hoop stress, MPa; σzz is the z-stress in cylindrical coordinate system, MPa. An arbitrary point P0 (x0, y0, z0) in the region adjacent to the penny-shaped crack, stress induced by a single triangular element on point P0 can be expressed by the following equation, according to Eq. (3). ( σ xx ¼ 2G Φz;zz þ ð1 2νÞΦz;yy zΦz;xxz ð24Þ σ zz ¼ 2G Φz;zz zΦz;zzz After the linear superposition of stress, the total stress acting on point P0 by the entire penny-shaped crack can be given as σ zz ðP0 Þ ¼
N h ih ih X ðP Þ ðP 0 Þ Djz ðP 1 Þ Sj Tj 0
Djz ðP2 Þ
Djz ðP3 Þ
iT
i ¼ 1; 2; …; N
j¼1
Fig. 7. Penny-shaped crack model.
ð25Þ
94
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
Fig. 8. Example of the triangular mesh on a penny-shaped crack (N, M and L represent the number of triangular, vertex and edge-vertex, respectively). (a) N ¼258, M ¼146, L ¼ 32 and (b) N ¼488, M ¼269, L ¼ 48.
Fig. 10. Relative error of KI on a pressurized penny-shaped crack. (a) N ¼258, M¼ 146, G ¼ 32 and (b) N ¼488, M¼ 269, G ¼ 48.
Fig. 9. Comparisons of analytical and numerical width. (a) N ¼ 258, M ¼146, L ¼32 and (b) N ¼488, M¼ 269, L ¼ 48.
where σ zz ðP0 Þ is the normal stress on point P0 in global coordinate h i ðP Þ system; T j 0 is the matrix for coordinate transformation from the h i ðP Þ jth local coordinate to global coordinate; Sj 0 are the coefficients
of P0 in Eq. (3) with respect to the displacement discontinuity components on the jth element. Note that the other stress components acting on point P0 have the similar equations as Eq. (25). With the axial symmetry of this crack and its boundary conditions, σxx and σzz on x-axis in Eq. (24) should be equal to σr and σzz in Eq. (23), respectively, even though they are in different coordinate system. Their comparisons prove a high accuracy between the numerical solutions with analytical solutions, as
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
∬Δ ξ ρ xdξ dη ¼ 1
AI ¼
2
R y3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R y pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AI η2 þ BI η þ C I dη 0 3 AII η2 þ BII η þ C II dη 0
þ 1;
BI ¼
k2 1 AII ¼ 2 þ1; k1
hðx; y; zÞ ¼ ∬Ω
95
2ðx2 xÞ 2y; k2
BII ¼
2x 2y; k1
C I ¼ ðx2 xÞ2 þ y2 þ z2 C II ¼ x2 þ y2 þ z2
η y dξ dη þ yf ðx; y; zÞ ρ
ρ¼
ðA2Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx ξÞ2 þ ðy ηÞ2 þ z2 ðA3Þ
∬Δ η ρ ydξ dη ¼
ffi R x pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R x3 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AIII ξ2 þ BIII ξ þ C III dξ þ x32 AIV ξ2 þ BIV ξ þ C IV dξ x1 Z
þ 2
AIII ¼ k1 þ 1;
x1 x2
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AV ξ2 þBV ξ þC V dξ
BIII ¼ 2x 2k1 y;
C III ¼ x2 þ y2 þ z2
2
AIV ¼ k2 þ 1; BIV ¼ 2x 2k2 ðk2 x2 þ yÞ; 2
2
Fig. 11. Comparisons of the analytical solutions with numerical solutions (x¼ ρ, y¼ 0, z ¼0).
shown in Fig. 11. Note that the error between analytical and numerical solution becomes larger when the point gets closer to pffiffiffi crack tip. This is because 1= r stress singularity exists at the adjacent domain of crack-tip. The stress in crack-tip domain is dominated by the SIF, not dominated by the linear superposition of the stresses induced by each element in the boundary. When r/R1 þ , σxx- þ1 and σzz-þ 1 according to Eq. (23). However, the distance between the point P0 and the gravity point PG of triangular element is finite, i.e., the integrals in Eqs. (4), (5) and (11) are non-singular. Hence, the numerical solutions of σr and σzz close to crack-tip are both smaller than their analytical solutions.
C IV ¼ x þ ðk2 x2 þ yÞ þ z2 AV ¼ 1; BV ¼ 2x; C V ¼ x2 þ y2 þ z2
ðA4Þ
k1 ¼ y3 =x3 ;
ðA5Þ
Eqs. A2 and A4 have the following form: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ax2 þ bx þ cdx ¼ 2ax4aþ b ax2 þ bx þ c þ pffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4ac b pffiffiffiffiffi ln2axþ b þ 2 a ax2 þ bx þ c þ Const 3 8 a 2
(1) A novel linear displacement triangular element with three nodes on its vertices is proposed. If the displacement discontinuities on these three nodes keep the same, this linear displacement triangular element will be simplified as the constant displacement triangular. (2) The solution procedure in this paper can be applied to determine the displacement and stress field around a three-dimensional fracture inside the infinite linear elastic solid. This procedure takes a shorter time and requires less elements than the constant DDM when achieving the same accuracy. If element subdivision is required for a reasonable accuracy then the analytical integration is computationally more efficient than the numerical integration.
∂2 gðx; y; zÞ xξ ∂2 f ðx; y; zÞ ¼ ∬Ω 3 dξ dη þ xlim 2 z-0 z-0 ∂z ∂z2 ℜ
ℜ¼
lim
AV I ¼
R y3 0
AVI η2 þBVI η þC VI
1 þ 1; k2 2
AV II ¼
1 2 þ 1; k1
ðA6Þ
1=2
BVI ¼ 2 x2k2 x y ; BVII ¼
2x 2y; k1
dη
R y3 0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx ξÞ2 þ ðy ηÞ2
ðB1Þ AVII η2 þ BV II η þC VII
C V I ¼ ðx2 xÞ2 þ y2 C V II ¼ x2 þ y2
∂2 hðx; y; zÞ yη ∂2 f ðx; y; zÞ ¼ ∬Ω 3 dξ dηþ ylim z-0 z-0 ∂z2 ∂z2 ℜ lim
ℜ¼
1=2
dη
ðB2Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx ξÞ2 þ ðy ηÞ2
ðB3Þ
1=2 Rx dξ ∬Ω yℜ3ηdη dξ ¼ x13 AIII ξ2 þ BIII ξ þ C III Z x2 Z x1 þ AIX ξ2 þ BIX ξ þ C IX dξ þ AX ξ2 þ BX ξ þ C X dξ x3
x2
2
AV III ¼ k1 þ 1; 2
AIX ¼ k2 þ 1; AX ¼ 1;
Acknowledgments This research was supported by the National Natural Science Foundation of China (51234006; 51325402; 51221003) and the China Scholarship Council (201406440009). Thanks to China Scholarship Council for supporting the first author of this paper to research in Georgia Institute of Technology, USA.
ða 4 0Þ
Appendix B. 2nd order partial derivatives
∬Ω xℜ3ξdξdη ¼
5. Conclusions
k2 ¼ y3 =ðx3 x2 Þ
BVIII ¼ 2x 2k1 y;
C V III ¼ x2 þ y2
BIX ¼ 2x 2k2 ðk2 x2 þ yÞ; BX ¼ 2x;
C IX ¼ x2 þ ðk2 x2 þ yÞ2 C X ¼ x2 þ y2 ðB4Þ
Eqs. B2 and B4 have the following form: Z
pffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1=2 1 dx ¼ pffiffiffi ln2ax þ b þ 2 a ax2 þ bx þ c þ Const ax þ bxþ c a
ða 4 0Þ
ðB5Þ References
Appendix A. Double integral of linear triangular element qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ξx gðx; y; zÞ ¼ ∬Ω dξ dη þ xf ðx; y; zÞ ρ ¼ ðx ξÞ2 þ ðy ηÞ2 þ z2 ρ ðA1Þ
[1] Ariza MP, Saez A, Dominguez J. A singular element for three-dimensional fracture mechanics analysis. Eng Anal Bound Elem 1997;20:275–85. [2] Berry DS. An elastic treatment of ground movement due to mining-I. Isotropic ground. J Mech Phys Solids 1960;8:280–92.
96
W. Cheng et al. / Engineering Analysis with Boundary Elements 59 (2015) 89–96
[3] Cheng W, Jin Y, Chen M, et al. A criterion for identifying hydraulic fractures crossing natural in 3D space. Pet Explor Dev 2014;41(3):336–40. [4] Cheng W, Jin Y, Lin Q, et al. Experimental investigation about influence of preexisting fracture on hydraulic fracture propagation under tri-axial stresses. Geotech Geol Eng 2015;37(3):413–8. [5] Crouch SL. Solution of plane elasticity problems by the displacement discontinuity method. Int J Numer Methods Eng 1976;10:301–43. [6] Crouch SL, Starfield AM. Boundary element methods in solid mechanics. London: Goerge Allen and Unwin Publishers; 1983. [7] Dyka CT, Remondi Millwater HR. 3D boundary elements and the integration of strong singularities. Comput Struct 1991;39(5):513–23. [8] Davey K, Hinduja S. Analytical integration of linear three-dimensional triangular elements in BEM. Appl Math Model 1989;13:450–61. [9] Kuriyama K, Mizuta Y. Three-dimensional elastic analysis by the displacement discontinuity method with boundary division into triangular leaf elements. Int J Rock Mech Min Sci Geomech Abstr 1993;30(2):111–23. [10] Kuriyama K, Mizuta Y, Mozumi H, Watanabe T. Three-dimensional elastic analysis by the boundary element method with analytical integrations over triangular leaf elements. Int J Rock Mech Min Sci Geomech Abstr 1995;32 (1):77–83. [11] Li H, Liu CL, Mizuta Y, et al. Crack edge element of three-dimensional displacement discontinuity method with boundary division into triangular leaf elements. Commun Numer Methods Eng 2001;17:365–78. [12] Lin’kov AM, Zubkov VV, Kheib MA. A method of solving three-dimensional problems of seam workings and geological faults. J Min Sci 1997;33 (4):295–315. [13] Liu G, Tu J, Zhang J. Solution for crack propagation in multi-crack body by displacement–discontinuous boundary element method. J Tsinghua Univ (Sci Technol) 1996;36:59–64. [14] Milroy J, Hinduja S, Davey K. The elastostatic three dimensional boundary element method: analytical integration for linear isoparametric triangular elements. Appl Math Model 1997;21:763–82. [15] Mogilevskaya SG, Nikolskiy DV. The use of complex integral representations for analytical evaluation of three-dimensional BEM integrals-potential and elasticity problems. Q J Mech Appl Math 2014;67(3):505–23.
[16] Napier JAL, Malan DF. The computational analysis of shallow depth tabular mining problems. J S Afr Inst Min Metall 2007;107:725–42. [17] Nikolski DV. Three-dimensional boundary element modeling of fractures under gravity load. In: Proceedings of the 46th American Rock Mechanics Association. 12-450. Chicago, USA; 24–27 June 2012. [18] Nikolski DV, Mogilevskaya SG, Labuz JF. Complex variables boundary element analysis of the three-dimensional crack problems. Eng Anal Bound Elem 2013;37:1532–44. [19] Nikolski DV and Mogilevskaya SG. Complex variables BEM modeling of three -dimensional crack problems. In: Proceedings of the 48th American Rock Mechanics Association symposium on US rock mechanics/geomechanics. Minneapolis, MN, USA; 1–4 June 2014. [20] Nikolskiy DV, Mogilevskaya SG, Labuz JF. Boundary element analysis of nonplanar three-dimensional cracks using complex variables. Int J Rock Mech Min Sci 2015;76:44–54. [21] Fata SN. Explicit expressions for the three-dimensional boundary integrals in linear elasticity. J Comput Appl Math 2011;235:4480–95. [22] Qin TY, Tang RJ. Three dimensional crack problem analysis using boundary element method with finite-part integrals. Int J Fract 1997;84:191–202. [23] Salamon MG. Elastic analysis of displacements and stresses induced by mining of seam or roof deposits. Part IV. J S Afr Inst Min Metall 1964;65:319–38. [24] Sokolnikoff. Mathematical theory of elasticity. New York: Mcgraw-hill book company; 1956. [25] Shou KJ, Siebrits E, Crouch SL. A higher order displacement discontinuity method for three-dimensional elastostatic problems. Int J Rock Mech Min Sci 1997;34(2):317–22. [26] Shou KJ. A three-dimensional hydrid boundary element method for non-linear analysis of a weak plane near an underground excavition. Tunn Undergr Space Technol 2000;15(2):215–26. [27] Sneddon IN. The distribution of stress in the neighborhood of a crack in an elastic solid. Proc R Soc Lond Ser: A Math Phys Sci 1946;46(1):187–229. [28] Xiao HT, Yue ZQ. A three-dimentional displacement discontinuity method for crack problems in layered rocks. Int J Rock Mech Min Sci 2011;48:412–20. [29] Wang C, Khoo BC. An indirect boundary element method for threedimensional explosion bubbles. J Comput Phys 2004;194:451–80.