A three-dimensional crack growth simulator with displacement discontinuity method

Engineering Analysis with Boundary Elements 48 (2014) 73–86

Contents lists available at ScienceDirect

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

A three-dimensional crack growth simulator with displacement discontinuity method Jingyu Shi a,n, Baotang Shen a, Ove Stephansson b, Mikael Rinne c a

CSIRO Earth Science and Resource Engineering, QCAT, 1 Technology Court, Pullenvale, QLD 4069, Australia Helmholtz/Centre Potsdam, GFZ German Research Centre for Geosciences, Potsdam, Germany c Aalto University, School of Engineering, Department of Civil and Environmental Engineering, Geoengineering, Espoo, Finland b

art ic l e i nf o

a b s t r a c t

Article history: Received 31 January 2014 Received in revised form 7 July 2014 Accepted 7 July 2014

This paper first outlines the theory of a well established three dimensional boundary element method: displacement discontinuity method (DDM) and proposes to use a crack growth criterion based on maximum normal or shear stress for a three dimensional crack growth simulator, FRACOD3D. Triangular elements are used in the simulator code. A numerical scheme is used to overcome a difficulty associated with the evaluation of the basic solution for DDM in some special situations and another numerical scheme is used to calculate the stresses on the boundary elements where the stresses obtained from the normal DDM scheme have large errors. The crack growth is implemented incrementally in that new front elements are introduced at the crack front; thus no need to re-mesh the old part of the cracks. The effects of neighbouring front elements are taken into account in implementation of the crack growth to overcome severer twisting of the new front elements generated from the growth. The numerical results from FRACOD3D of two simple examples agree very well with analytical solutions, and propagation configuration of a circular disc crack in an infinite body under shear is close to that observed in an experiment in literature under similar loading condition. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Displacement discontinuity method Crack growth Maximum stress criteria

1. Introduction Numerical simulation becomes essential for engineering problems for which physical experiments are very difficult or even impossible, such as crack propagation around mining sites. There are several classes of numerical simulation methods, including Finite Difference Method (FDM), Finite Element Method (FEM), Boundary Element Method (BEM), Finite Volume Method (FVM) and Discrete Element Method (DEM). FEM may be the most popular method, but others have their advantages. One advantage of boundary element methods over finite element methods is that BEM uses formulations of one dimension less than others: onedimensional formulation for two-dimensional problems and twodimensional formulation for three-dimensional problems, because it only deals with the boundary of the problems. So generally the number of final equations in BEM is smaller than that of FEM for the same problem. Displacement discontinuity method (DDM), a BEM, for crack problems has at least one advantage over other BEMs in that the DDM treats the crack as one surface while other BEMs have to consider discretisation of both surfaces of the crack. Therefore the DDM further reduces the number of the final

n

Corresponding author. Tel.: þ 61 7 3327 4441. E-mail address: [email protected] (J. Shi).

http://dx.doi.org/10.1016/j.enganabound.2014.07.002 0955-7997/& 2014 Elsevier Ltd. All rights reserved.

equations. But one should note that the final system of equations from BEM has a dense coefficient matrix, or even a full matrix in most cases, different from the sparse matrix from FEM and thus the solution schemes for sparse system of equations cannot be used. DDM is based on the fundamental solution for displacements and stresses in infinite elastic bodies caused by displacement dislocation (displacement discontinuity – DD) on surfaces of a finite planar crack inside the body. The fundamental solution is valid for isotropic, homogeneous and linear elastic material. If the crack is curved and non-planar, then it is approximated as an assembly of a finite number of planar crack elements. The displacements and stresses at a position inside the body are then taken as the sum of the effects from all the planar elements. This superposition is possible since the material is assumed to be linear elastic. For problems of finite body or cavity in infinite body, the surface of the finite body or surface of the cavity can be treated as imaginary closed curved crack in an infinite body, with the region outside the finite body or the cavity being filled with the same material, allowing for the fundamental solution to be employed. Generally, in practice the displacement discontinuities at the elements are not known, so the fundamental solution cannot be used directly. If the displacements and/or stresses at a sufficient number of positions are specified as the boundary condition, then the displacement discontinuities at the planar elements that cause

74

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

the specified displacements and/or stresses can be found by solving a system of equations. The displacements and stresses at other positions in the body can then be calculated from the obtained displacement discontinuities on the elements. This is the basic principle of DDM. In two dimensional cases, Crouch and Starfield [1] outlined the theory and implemented the method in computer code, which is supplied in the book. Crouch and Starfield [1] also presented the governing equations for three dimensional situations. Shen [2] used the DDM to simulate crack growth with a newly proposed criterion for rock masses. The basic two dimensional DDM has been extended to couple hydraulic and thermal effects and to handle multiple material regions, see [3]. Various ways have been used to derive the governing equations of DDM. For example, Hong and Chen [4] showed that DDM is a special case of dual boundary element method. In three dimensional cases, Rongved [5] obtained the fundamental solution by using four special harmonic functions. There are other ways for the formulation [6]. Many papers on three dimensional DDM exist, see [7–12] for different aspects of the method. DDM has also been employed in three dimensional poroelastic and therm-elastic problem modelling [13,14]. The core of DDM is integral functions, integrals over the planar crack region of products of Green function and the displacement discontinuity components [5]. If the displacement discontinuity components are uniform over the crack region, then the integral functions depend on one single integral. In this case, to calculate displacements and stresses in the body, the derivatives of the integral function are needed. The analytical expression of the integral function has been obtained for a triangular element (see [15,7]) and its partial derivatives can be computed analytically, see [8] for those derivatives needed in DDM. However, if the perpendicular projection of the considered point on the integration plane is on a line along which a side of the triangular integration element lies, then the integral function and its derivatives become infinite. This makes the evaluation of these expressions difficult or impossible. This was reported by Kuriyama and Mizuta [8], but they did not provide an alternative scheme. Instead they tried to arrange elements in their examples so that this issue did not occur. It is noted that the whole analysis breaks down if this occurs on just any one element. This cannot be avoided for some of analysed domains, such as rectangular block and cylinder surfaces. In addition, when propagation of cracks is considered, it is very hard to avoid this occurring on elements on newly extended crack surface. In this paper, we propose to use a numerical integration scheme for such cases, or even for regular points. We also discuss various methods for computing the derivatives. In addition, we employ a numerical scheme to calculate stresses on the elements, where the normal DDM procedure yields large errors. Cracks exist in natural rock medium so a computational code used for simulating rock medium needs to be able to simulate crack growth in cases where crack front stresses exceed the rock strength. Three-dimensional crack growth has a very complicated mechanism. Various theories exist to describe crack growth, such as maximum principal stress criterion and minimum strain energy density factor criterion [16,17]. With these theories there are three basic modes in which a crack can grow: open mode (Mode I), shearing mode (Mode II) and tearing mode (Mode III). In practice, a crack could propagate in a way which combines these basic modes. A modified criterion based on maximum stress is proposed here to determine whether or not a crack will grow, and if it grows, in which direction it will go. DDM is based on the displacement discontinuity of cracks, so it should be suitable for crack growth simulation. Although there is a large number of literatures on crack growth simulation in twodimension using DDM, there are not many papers in three-

dimensional cases, as far as the authors are aware, possibly due to the difficulty of eliminating or reducing severe twisting of front elements generated from the growing crack surface after a few steps of propagation with general loading conditions [18]. This difficulty exists also in other numerical methods [19]. The difficulty is even worse when multiple cracks grow and their growths meet other cracks. This will be addressed separately in a future work. A paper published in Journal of Donghua University (2010) by Wang and Huang was brought into our attention recently. Wang and Huang gave some implementation details of 3D crack growth, together with basis of DDM and crack propagation criterion. In this paper we describe the theory and implementation of a simulation code FRACOD3D, based on three dimensional DDM, and the proposed criterion for crack growth. We attempt to make the code robust for general situations. The crack growth is implemented incrementally in that new front elements are introduced at the whole crack front step by step. In this way, mesh of the old part of the cracks is kept at each step of growth. The effects of neighbouring front elements are taken into account in implementation of the crack growth to overcome severer twisting of the new front elements generated from the growth. It is found in the present paper that this scheme works well. The DDM part of the code is verified with two examples for which analytical solutions exist. These include (1) circular planar crack under tension in the direction normal to the crack plane and (2) spherical cavity under uniform inner pressure on the cavity surface. Capability of the code for simulation of crack growth is shown with one example of circular planar crack under shear. The paper is arranged as follows. Section 2 outlines the fundamental solution due to displacement discontinuity on a crack and Section 3 gives the details of how to use the fundamental solution to solve practical problems. Various ways of evaluating the coefficients in DDM are outlined in Section 4. Section 5 formulates the criterion for crack growth. Implementations of DDM and the crack growth are described in Sections 6 and 7, respectively, before two verification examples and one crack growth example are shown in Section 8. Finally, some conclusions are given in Section 9.

2. Three-dimensional elasticity with displacement discontinuity The fundamental solution of DDM is for displacements and stresses at any point in an infinite body caused by known displacement discontinuity on a planar crack in the infinite body. To express the solution mathematically, a local coordinate system, oxyz, is set up such that the planar crack is located on the oxy plane, as shown in Fig. 1. Let the displacement components on the lower surface of the crack, which has outward normal in positive z direction, be uxþ ; uyþ ; uzþ , and let the displacement components

Fig. 1. Coordinate system and components of displacement discontinuity.

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

on the upper surface of the crack, which has outward normal in negative z direction, be ux ; uy ; uz . The displacement discontinuity (DD) components are defined as Dx ¼ uxþ  ux ; Dy ¼ uyþ  uy ; Dz ¼ uzþ  uz : It is noted that this definition yields negative Dz if the crack surfaces are pulled open and positive if the surfaces are compressed, or “over-closed”. The displacement components could be functions of the position on the surface of crack, so the DD components could depend on the position. Rongved [5] first expressed the displacement components in the elastic body that satisfy the general equilibrium equations in terms of four general harmonic functions and used the generalised Hooke's law to express the stress components. Then the harmonic functions were determined so that the displacement and stress components were continuous everywhere in the infinite body except that the displacements had discontinuity components Dx, Dy, Dz over the crack surfaces. Let the planar crack have domain area A. In general case of varying displacement discontinuity components over the domain area, the fundamental solutions for displacements and stresses depend on partial derivatives of integrals of the displacement discontinuity components. In this paper, we employ assumption that the displacement discontinuity components are uniform on the domain area A. Thus the displacements and stresses depend on one integral [1,7,8] 1 I ¼ ∬A dx0 dy0 r

ð1Þ

where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ ½xðQ Þ  x0 ðtÞ2 þ½yðQ Þ y0 ðtÞ2 þ½zðQ Þ2 0

being the distance from the integration point tðx ðtÞ; y ðtÞ; 0Þ to the collocating point Q(x,y,z) (see Fig. 1) inside the infinite elastic body, where the displacements and stresses are calculated. In this case, the displacement components and stress components inside the elastic body are given by ux ðQ Þ ¼

  1 ½2ð1  νÞI z  zI xx Dx zI xy Dy  ½ð1  2νÞI x þ zI xz Dz ; 8π ð1  νÞ

uy ðQ Þ ¼

  1  zI xy Dx þ ½2ð1  νÞI z  zI yy Dy  ½ð1  2νÞI y þ zI yz Dz ; 8π ð1  νÞ

 1 ½ð1  2νÞI x  zI xz Dx þ ½ð1  2νÞI y  zI yz Dy 8π ð1  νÞ þ ½2ð1  νÞI z  zI zz Dz g:

uz ðQ Þ ¼

ui ðQ Þ ¼ Bix Dx þ Biy Dy þBiz Dz ¼

∑ Bik Dk

i ¼ x; y; z;

ð4Þ

k ¼ x;y;z

σ ij ðQ Þ ¼ Aijx Dx þ Aijy Dy þAijz Dz ¼ ∑ Aijk Dk

i; j ¼ x; y; z

ð5Þ

k ¼ x;y;z

If the crack has a curvature such that the conditions for the above expressions are invalid, then an approximation can be made. The approximation is made by subdividing the crack into small elements with simple planar shapes, such as triangles or rectangles, approximating the displacement discontinuity components on each element as uniform so that the above formula can be used, and then adding up the effects from all elements to calculate the approximate value for the whole crack. It should be noted that the above solutions are expressed in the local coordinate system on each element. Before the superposition, the solutions from each element must be transformed so that they refer to a universal coordinate system. Let B0ij e and A0jik e (i,j,k¼ x,y,z) be the coefficients of the DD components on element e in the displacement and stress component expressions referred to the universal coordinate system o0 x0 y0 z0 that is fixed relative to the DD element. The displacement and stress components at point Q, referred to the system o0 x0 y0 z0 , produced by the DD components on all the elements, are given by u0i ðQ Þ ¼ ∑

3

∑

e ¼ 1 k ¼ x;y;z Ne

σ 0ij ðQ Þ ¼ ∑

B0ik eDke

3

∑

e ¼ 1 k ¼ x;y;z

A0ijk eDke

i ¼ x; y; z

i; j ¼ x; y; z

ð6Þ

ð7Þ

where Ne is the total number of the elements in the discretisation. The above expressions show how to compute the displacement and stresses inside an infinite body that are caused by the displacement discontinuity on complicated cracks.

3. Three-dimensional DD boundary element method (3D DDM) ð2Þ

and  G ½2I xz  zI xxx Dx þ ½2νI yz  zI xxy Dy þ ½2νI zz ð1  2νÞI xx 4π ð1  νÞ   zI xxz Dz ;  G ½2νI xz  zI xyy Dx þ ½2I yz  zI yyy Dy þ½2νI zz  ð1  2νÞI yy σ yy ¼ 4π ð1  νÞ   zI yyz Dz ;   G  zI xzz Dx  zI yzz Dy þ½I zz zI zzz Dz ; σ zz ¼ 4π ð1  νÞ  G ½ð1  νÞI yz  zI xxy Dx þ ½ð1  νÞI xz zI xyy Dy  ½ð1  2νÞI xy σ xy ¼ 4π ð1  νÞ  þ zI xyz Dz ;   G ½ð1  νÞI zz  νI xx zI xxz Dx  ½νI xy þ zI xyz Dy  zI xzz Dz ; σ xz ¼ 4π ð1  νÞ

σ xx ¼

σ yz ¼

where G and v are the shear modulus and the Poisson ratio of the elastic material, respectively. The subscripts of I denote the partial derivatives. Using notations Bij (i,j¼ x,y,z) for the coefficients of the displacement discontinuity components Dx, Dy, Dz in the displacement expressions and Aijk (i,j,k¼ x,y,z) for the coefficients in the stress components, the above solutions can be expressed as

Ne

0

75

  G  ½νI xy þ zI xyz Dx þ ½ð1  νÞI zz  νI yy  zI yyz Dy  zI yzz Dz : 4π ð1  νÞ ð3Þ

The above expressions are for displacement and stress components when the displacement discontinuity components on a crack in an infinite elastic body are known. However, in most practical problems related to cracks, the known variables in advance are displacements or tractions on the crack, not the DD components. In some other practical problems, the elastic body is finite, not infinite, in size. For finite bodies, the displacement and/ or traction components on an exterior boundary are generally given. The problems associated with cracks can be solved using various numerical methods and the displacement discontinuity method (DDM) is considered to be most suitable. With DDM, such problems are approximately solved with the following steps: 1. Discretise the crack into small planar elements. In the case of a finite body and an infinite body with a cavity, an infinite body with an imaginary crack is assumed. The shape of the imaginary crack is identical to the boundary of the finite body or the surface of the cavity in the infinite body. The imaginary crack is discretised into small planar elements. Assume that the number of the elements is Ne.

76

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

2. Approximate the known continuously distributed tractions and/or displacements on the whole crack with a piece-wise continuous distribution, such as the uniform distribution, over each element. 3. Assume that the known tractions and/or displacements on the crack are caused by some DD components over each planar element, and then try to determine the correct values of the DD components that produce the known values of tractions and/or displacements. This is done by (a) expressing the specified tractions and/or displacements at each element in terms of the assumed DD components on all elements, (b) setting up a system of governing equations for the DD components and (c) solving the system of equations for the DD components. In the process of setting up the equations, variables must be transformed into a universal coordinate system in which the tractions and/or displacements are given before the contributions from all the elements can be added together. 4. Calculate the displacements and the stresses inside the body with the found values of the DD components according to the formulae given above. When the displacement components on the surface of an exterior boundary of a finite body, or of the cavity in an infinite body, are given, the governing equations for DD components are obtained by applying the expression (6) to the collocating points Q at the centroid of each element where displacements are specified. Note that although the specified displacement components can refer to a global universal coordinate system, the displacement discontinuity components DDs refer to the local coordinate systems on the elements. For traction boundary value problems, we need to express tractions in terms of the displacement discontinuity components DD. With the stress components given in Eq. (7), we can obtain the traction components T x ; T y ; T z on a plane element at the point Q as Ne

T i ðQ Þ ¼ ∑

3

∑

e ¼ 1 k ¼ ðx;y;zÞ

H eik Dke

ði ¼ x; y; zÞ

ð8Þ

where the coefficients H eij ði; j ¼ x; y; zÞ are related to the coefficients A0jik e. Applying these expressions to collocating points at the centroid of each element where tractions are given, we can obtain the governing equations for DD components. If displacement is given on part of the boundary surface and traction is specified on the remaining part, then the problem has a type of mixed boundary conditions. In this case, the discretisation should make the boundary between the traction part and displacement part lie on the sides of the small elements. Employing the expression (6) for specified displacement on the displacement part and Eq. (8) for specified traction on the traction part will produce a system of equations for the DD components. There is another type of mixed boundary conditions: on an element, conditions are not specified for all three displacement components or not specified for all three traction components. It is that one or two displacement component/s and two or one traction components are given on some or all elements. Sliding of a body on another hard body is one example of this. In this case, expressions for the corresponding component in Eqs. (6) and (8) should be used to establish the governing equations.

4. Integrations for 3D DDM The coefficients in the expressions for the displacement and stress (or traction) components contain partial derivatives of the integral, I, over the discrete elements, with the integrand being inversely proportional to the distance between the integration point and the collocating point. In this section, we consider the

calculation of these coefficients, i.e. calculation of the integral itself and its derivatives. When displacements and stresses at points inside the body are calculated, the integrand is finite and the integration is regular. In the formation of the governing equations for the DD components, besides regular integrations in which the collocating points are on other elements rather than the integration element, there are integrations in which the collocating point is on the integration element. In the latter case, the integrand becomes infinite at the collocating point so the integral is singular. We will discuss methods of calculating the integrations for each of these cases in the following sections.

4.1. Regular integrations When the collocating point is not on the integration element, the integrand is finite. In this case, the derivatives of the integral can be calculated in two ways. (1) Differentiate the integration result of the integral directly with a numerical scheme, or if possible, analytically. (2) Exchange the order of differentiation and integration, and then integrate the derivatives of the integrand, numerically, or if possible analytically, if certain conditions are satisfied. The integration of 1/r, denoted as I in Eq. (4), has been obtained analytically for triangular elements [7,8,15]. Therefore the coefficients can be evaluated according to the first method by differentiating the analytical integration expression analytically. Certainly, numerical differentiation can be employed too. With the second method, i.e. exchanging the order of differentiation and integration and then integrating the derivatives of 1/r, after the exchange, the integrands can be as high as 1/r7 (for third derivatives) and have various forms, such as x(t)y(t)/r3 or x(t)2/r5. Medina and Liggett [15] reported analytical integrals up to x(t)2y(t)/r5 on triangular elements. There are no analytical integrations of order 1/r7 found and a numerical method may be needed for these higher orders. For simplicity, it may be better to do all the integrations numerically, although the analytical ones are exact. Although the analytical expression of the integral, I, has been obtained for a triangular element [7,15] and its partial derivatives can be done analytically [8], if the perpendicular projection of the collocating point on the integration plane is on a line along which a side of the triangular integration element lies, then the expressions of the integral and its derivatives become infinite. This was reported by Kuriyama and Mizuta [8], but they did not provide an option for dealing with such cases, except for arranging elements in their examples so that no such problem arises. For problems with boundary surfaces mutually perpendicular to one another, such as rectangular blocks or cylinders, this situation cannot be avoided. Furthermore, in problems with crack propagation, it is very hard to avoid adding a new element with its centroid projection falling on a side of other elements. In this case, we propose to use the second method and numerically integrate the derivatives of the 1/r. We assume uniform DDs on the elements in our code. But here we outline some schemes that can be used to handle cases where the DD components are not uniform on the integration element. In this case, the integrations of Dx/r, Dy/r and Dz/r depend on the form of the variations of the DD components and no papers on the analytic expressions of the integrations have been found. Most likely the integrations cannot be carried out analytically even for simple variations of the DD components and therefore the derivatives need to be obtained with numerical schemes. With the first method, numerical differentiation of the numerical

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

integrations can be used. For example, for the integral Ix Dx I x ¼ ∬A qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffidx0 ðtÞdy0 ðtÞ 2 0 ½xðQ Þ  x ðtÞ þ½yðQ Þ  y0 ðtÞ2 þ ½zðQ Þ2

ð9Þ

the derivative with respect to x (i.e. x(Q)) of Ix can be done by firstly numerically integrating two integrations, one with x(Q) and the other with x(Q)þ h (h being a very small value), and then dividing the difference of the two integrations by the small value h ∂I x =∂x  fI nx ½xðQ Þ þ h  I nx ½xðQ Þg=h; where I nx ½n being the numerical integration result. With the second method of exchanging the order of differentiation and integration, the derivatives of Dx/r, Dy/r and Dz/r can be integrated numerically. Note that Dx, Dy and Dz are functions of integration variables only, and do not depend on the differentiation variables. So the integrands will be the products of Dx, Dy and Dz and derivatives of 1/r. These derivatives are the integrands for the case of uniform DD components with the second method of order exchange. 4.2. Singular integrations When the collocating point is on the integration element, the integral becomes singular since at the collocating point the integrand becomes infinite. This occurs when traction and/or displacement components on the crack are calculated and the contribution from the element that contains the collocating point is needed. This situation cannot be avoided in formation of the governing equations for DD components. We follow the procedure of Cayol and Cornet [7] for this situation. The traction and/or displacement components should take the values when the collocating point approaches the crack surface from the interior of the body, that is, the limit is taken as zðQ Þ-0: The limiting procedure is after the integration and differentiation procedures, so before the limiting procedure is taken, the integrands are continuous and the integration and differentiation procedures are the same as in the regular case, for example, the order of integration and differentiation can be exchanged. If analytical results for the integration and differentiation procedures can be obtained, then the limiting procedure can be carried out. Unfortunately the limiting procedure cannot use numerical results of the integration and differentiation. It is noted that the second derivatives of the integral in displacement expressions, Eq. (1), and the third derivatives of the integral in stress expressions, Eq. (2), disappear after the limitation is taken, since the factor z (i.e. z(Q)) in the front of these terms becomes zero. In the case of uniform DD components on the element, the analytic expressions for the integration and differentiation are available and the limits can be carried out. The limits of the needed first and second derivatives of the integral of 1/r have been given in the paper of Cayol and Cornet[7], with the collocating point being at the centroid of triangular element. It is noted that when the DD components are not uniform, integrations of Dx/r, Dy/r and Dz/r depend on the variation of the DD components and generally it is hard to evaluate the integrations and their derivatives analytically. Therefore it is hard to find the limits.

and tearing Mode III, depending on the stress state around the crack front [16]. A crack growing with open mode is caused by tensile stress in the direction normal to the crack plane. The sliding and tearing modes are caused by shear stresses in different directions on the crack plane. Expressions for dominant parts of three dimensional stresses and displacements near the crack front have been obtained in literatures and they are expressed in terms of three stress intensity factors. The stress intensity factors are quantities that describe the strength of the stress singularity on the crack front. There are three basic stress intensity factors, KI, KII, KIII, corresponding to the three basic crack modes. The basic stress intensity factors are used in various criteria for the crack growth. They are calculated by either limits of stresses near the front or the differences of the displacements on the two surfaces of the crack near the front. Since the differences of the displacements on the two surfaces of the crack are the DD components of DDM, we reproduce the expressions for the three stress intensity factors in terms of the displacement differences [16] rffiffiffiffiffi E π pffiffiffiffiffiffiffiffiffiffiffiffiffiffi KI ¼ ½ub ðr; θ ¼ π ; ϕÞ ub ðr; θ ¼  π ; ϕÞ; 4ð1  ν2 Þ cos ϕ 2r rffiffiffiffiffi E π pffiffiffiffiffiffiffiffiffiffiffiffiffiffi K II ¼ ½un ðr; θ ¼ π ; ϕÞ  un ðr; θ ¼  π ; ϕÞ; 4ð1  ν2 Þ cos ϕ 2r rffiffiffiffiffi E π pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð10Þ K III ¼ ½ut ðr; θ ¼ π ; ϕÞ  ut ðr; θ ¼  π ; ϕÞ 4ð1 þ νÞ cos ϕ 2r The notations in these expressions are shown in Fig. 2 and E is Young's modulus of the material. There are three aspects in crack growth that need to be considered: (1) when a crack starts to grow at a point on the crack front; (2) in which direction the crack will go and (3) how far the crack will go. For each of these aspects, there are various criteria. For example, for the aspect of (1) when a crack starts to grow, there is an energy release rate criterion [16] and minimum strain energy factor criterion [17]. For aspect (2) in which direction the crack grows, there is a modified maximum principal stress criterion, minimum strain energy density factor criterion and maximum energy release rate criterion. Paris law and minimum strain energy density factor criterion are two examples for aspect (3) how far the crack will grow. It can be seen that some criteria cover all the three aspects and some criteria determine two aspects simultaneously. With three dimensional stress states, crack growth will most likely not follow the old crack orientation. In this paper, we adopt a quasi-plane criterion to determine when a front element will grow and in which direction it grows. We assume that if a front element will grow, then it will grow in a direction which is on the plane normal to the crack front line, that is the plane given by ϕ ¼0 in Fig. 2. We propose that on this plane, crack will propagate in the direction along which normal stress or shear stress takes a

5. Three-dimensional crack growth with FRACOD3D Cracks exist in natural rock medium so we need to consider crack growth for engineering problems related to rock. A major advantage of the DDM against other BEM lies in its simple presentation of cracks. It is well known that there are three basic modes at a crack front in which the crack can grow: open Mode I, sliding Mode II

77

Fig. 2. Failure element on a sphere centred on the crack front.

78

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

maximum value. In the following we will consider exactly in which direction, represented by angle θ, the crack will go. On this plane, the three basic stress intensity factors can be calculated in terms of the DDs as rffiffiffiffiffiffi E 2π ð Db Þ KI ¼ 8ð1  υ2 Þ d rffiffiffiffiffiffi E 2π Dn K II ¼ 8ð1  υ2 Þ d rffiffiffiffiffiffi E 2π Dt K III ¼ ð11Þ 8ð1 þ υÞ d

with φ being the internal friction angle and c being the cohesion of the material. Note that tensile normal stress is positive and compressive normal stress is negative. p Since ffiffiffiffiffiffiffiffi the term c does not have a singularity, its product with 2π r is neglected and we have pffiffiffiffiffiffiffiffi K IIce ¼ K IIc þ 2π r ð  σ θ jθ ¼ θIIe Þ tan φ    θIIe 3θIIe þ cos ¼ K IIc  K I 3 cos 2 2   θIIe 3θIIe  K II 3 sin þ 3 sin tan φ ð16Þ 2 2

where d is the distance from centroid of the crack front element to the crack front and the displacement discontinuity components take their values at the centroid of the crack front element (see [16]. As defined earlier, Db is negative for displacements related to open crack surface.) The stress state near the crack front on this plane due to three dimensional deformations has the following dominate terms [20]:     KI θ 3θ K II θ 3θ þ pffiffiffiffiffiffiffiffi  5 sin þ 3 sin þ OðrÞ; 5 cos  cos σ r ¼ pffiffiffiffiffiffiffiffi 2 2 2 2 4 2π r 4 2π r     KI θ 3θ K II θ 3θ þ pffiffiffiffiffiffiffiffi  3 sin  3 sin þ OðrÞ; σ θ ¼ pffiffiffiffiffiffiffiffi 3 cos þ cos 2 2 2 2 4 2π r 4 2π r     KI θ 3θ K II θ 3θ r  1=2 sin þ sin þ pffiffiffiffiffiffiffiffi cos þ 3 cos þ OðrÞ; τrθ ¼ pffiffiffiffiffiffiffiffi 2 2 2 2 4 2π r 4 2π r

Similarly, we can define the combined Mode III stress intensity factor θ θIIIe ð17Þ ¼ K III K IIIe ¼ max K III cos ¼ K III cos 2 2 θ A ½  π ;π 

K III θ sin þ OðrÞ; τrϕ ¼  pffiffiffiffiffiffiffiffi 2 2π r K III θ τθϕ ¼ pffiffiffiffiffiffiffiffi cos þOðrÞ: 2 2π r

ð12Þ

First we look at the possibility for the opening mode. This mode of crack growth is due to the tensile stress σ θ in the θ-direction. From the definition of the basic Mode I stress intensity factor KI, it can be seen that it is proportional to the maximum normal pffiffiffiffiffiffiffiffistress from the pure Mode I deformation, with coefficient of 2π r . The material toughness for Mode I, KIc, is a property that measures the material ability to resist against the Mode I crack growth. It is accepted that if KI 4KIc, the crack will grow with Mode I. If we define the quantity that is proportional to the maximum normal stress from the threepffiffiffiffiffiffiffiffi dimensional combined opening deformation, with coefficient 2π r , as combined Mode I stress intensity factor KIe, we have      θ 3θ θ 3θ  K II 3 sin þ3 sin K I 3 cos þ cos K Ie ¼ max 2 2 2 2 θ A ½  π ;π      θIe 3θIe θIe 3θIe þ cos  K II 3 sin þ 3 sin : ð13Þ ¼ K I 3 cos 2 2 2 2 The angle θIe defines the direction in which the maximum value KIe occurs. Cracks may propagate in shear mode. Similar to the opening mode of crack surfaces, we define the quantity proportional to the maximum shear stress as combined Mode II stress intensity factor     θ 3θ θ 3θ K IIe ¼ max K I sin þ sin þ K II cos þ3 cos 2 2 2 2 θ A ½  π ;π      θ 3 θ θ 3 θ IIe IIe IIe IIe : ð14Þ ¼ K I sin þ sin þ K II cos þ 3 cos 2 2 2 2 The toughness KIIc measures the material limit of crack growth under pure shearing deformation Mode II. The internal friction and cohesion add some more strength to resist the growth of the crack when there is compressive normal stress near the crack front. So we can set up a combined toughness for shear deformation mode pffiffiffiffiffiffiffiffi K IIce ¼ K IIc þ 2π r fð  σ θ jθ ¼ θIIe Þ tan φ þ cg; ð15Þ

and combined toughness for tearing deformation mode pffiffiffiffiffiffiffiffi K IIIce ¼ K IIIc þ 2π r ð  σ θ jθ ¼ θIIIe Þ tan φ    θIIIe 3θIIIe þ cos ¼ K IIIc  K I 3 cos 2 2   θIIIe 3θIIIe þ 3 sin tan φ  K II 3 sin 2 2

ð18Þ

Clearly, KIIIe occurs at θ ¼ θIIIe ¼ 0. In this direction, K IIIce ¼ K IIIc  4K I tan φ: To determine the direction and mode for the crack propagation, we define three ratios RI ¼ K Ie =K Ic ;

RII ¼ K IIe =K IIce ;

RIII ¼ K IIIe =K IIIce :

ð19Þ

Let RC be the maximum of RI, RII and RIII. If Rc o 1, then the crack does not propagate. If Rc Z 1 and Rc ¼ RI , then it is regarded as Mode I propagation and the direction is given by the angle θIe . If Rc Z 1 and Rc ¼ RII , then it is regarded as Mode II propagation and the direction is given by the angle θIIe . If Rc Z 1 and Rc ¼ RIII , then it is regarded as Mode III propagation and the direction is given by the angle θIIIe . The matter of how far the crack would propagate is handled as follows. For the steady-state problem, it is considered that the propagation distance is “one element length”. If some element propagates, one or two or more new elements with similar size are added to its front, depending on the length of the front edge of the element and propagating status of its neighbours. The distance of the new elements in the growth direction is the average of side lengths of the element. If there are many elements propagating, new elements are added to each of the propagating elements. With the newly created area from the new elements, the crack propagation criteria are checked again and new elements are added if propagation occurs. The procedure continues until there are no more elements propagating. Although, with the newly added elements, the growth test shows that no more elements propagate, the crack might grow if smaller elements are used. So strictly speaking for more accurate results, the newly added elements at last step should be taken off and smaller elements should be used, instead, for further analysis and growth test until some criterion of accuracy on the element size is satisfied.

6. Implementation of the basic of 3D DD method The displacement and stress components due to a displacement discontinuity on a crack element are related to derivatives of a function, the so-called integral function, which is an integral over a crack element. These derivatives of the integral function are involved in the influence coefficients for the displacements and stresses and, as stated earlier, can be carried out in one of two methods. These are: (1) first find the analytical expression of the

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

integral function and then differentiate the resulting function, analytically or numerically; and (2) exchange the order of integration and differentiation by first differentiating the integrand analytically and then integrating the derivatives of the integrand analytically or numerically. Cayol and Cornet [7] give the analytical expression for an integral function; Kuriyama and Mizuta [8] give analytical expressions for the integral function and its derivatives. So the influence coefficients could be evaluated analytically with the first method. But at a point whose perpendicular projection on the plane of an integration element is close to prolonging the line of a side of the integration element, the coefficients become very large. The coefficients become infinite if the projection is on the line. In such cases, the second method is employed in this paper and integrations of the derivatives of the integrand are done numerically. The numerical procedure can also be employed for other points. In the numerical integration scheme, the integration can be approximated by multiplying the area of the whole element by the value of the integrand evaluated at the centroid of the element. The accuracy of the numerical integration can be improved by subdividing the boundary element into smaller regions and the integration on the boundary element being taken as the sum of the integrations on each of the small regions. On a small region, the integration is approximated by the product of the area of the small region and value of the integrand evaluated at a central point inside the small region. In FRACOD3D, the small regions are triangles. Another difficulty of DDM lies in the fact that integrations become singular with the integrand being infinite at the centres of crack elements, over which the integrations are evaluated, in the calculation of influence coefficients for displacements and/or stresses at the centres. The displacement discontinuity boundary element method must overcome this difficulty to obtain the basic variable, i.e. the displacement discontinuity on the crack surfaces, since displacement or stresses at the centres of crack elements are needed in the set up of the basic equations. This is also true if displacements and/or stresses at the centre of elements are to be calculated for other purposes. This difficulty can be overcome by employing the limits of Cayol and Cornet [7]. In the two dimensional case, Crouch and Starfield [1] reported that if stresses at points on the crack surface are calculated with normal procedure from displacement discontinuities after the DDs are obtained from the governing equations, then they have large errors. They suggested the use of numerical procedures to calculate the derivatives of displacement, i.e. the strain, numerically and then use the generalised Hook's law and the given boundary traction to calculate the stress field on the crack surface. The strains are calculated numerically with displacements at the neighbouring elements. Large errors also occur in three dimensional DDM. In the three dimensional case, the strains contain partial derivatives of the displacement components in the local x and y directions. In general, neighbouring elements do not locate at the positions with the same coordinates and the partial derivatives, with respect to the coordinates, cannot be calculated directly. However, the directional derivatives of the displacement components, with respect to direction related to a neighbouring element, can be calculated numerically. Then the relations between the directional derivatives and the partial derivatives with respect to the coordinates are employed to obtain the partial derivatives. This procedure is outlined by Shi and Shen [21] and implemented in the code. It gives relatively accurate results. Here we outline a basic scheme for general case of triangular element with three neighbours sharing common edges. For more details and some special situations, readers are referred to Shi and Shen [21]. Assume that triangular element P has three neighbours Q1, Q2 and Q3 on its three sides and we want to calculate the stresses at

79

element P. The displacements at neighbouring elements have been transformed into the local coordinate system at element P. We can numerically calculate the directional derivatives of the displace

!

!

! ment components along the three directions PQ 1 , PQ 2 and PQ 3 .

! For instance, the directional derivative of ux in the direction PQ 1 can be calculated numerically with ∂ux =∂n1  ½ux ðQ 1 Þ  ux ðPÞ= ΔlPQ 1 , where ΔlPQ 1 is the distance from the centre of element P to the centre of element Q1. With the relation between directional derivative and partial derivatives with respect to coordinates, we can establish ∂ux ∂ux ∂ux ux ðQ 1 Þ  ux ðPÞ cos ðQ 1 xÞ þ cos ðQ 1 yÞ þ cos ðQ 1 zÞ  ; ∂x ∂y ∂z ΔlPQ 1 where (cos(Q1x), cos(Q1y), cos(Q1z) ) are the directional cosines of

! the direction PQ 1 in the local coordinate system at element P. Similarly, we can obtain two more equations for the three partial

!

! derivatives with directional derivatives along PQ 2 and PQ 3 ∂ux ∂ux ∂ux ux ðQ 2 Þ  ux ðPÞ cos ðQ 2 xÞ þ cos ðQ 2 yÞ þ cos ðQ 2 zÞ  ; ∂x ∂y ∂z ΔlPQ 2 ∂ux ∂ux ∂ux ux ðQ 3 Þ  ux ðPÞ cos ðQ 3 xÞ þ cos ðQ 3 yÞ þ cos ðQ 3 zÞ  : ∂x ∂y ∂z ΔlPQ 3 Generally, the three partial derivatives ∂ux =∂x; ∂ux =∂y; ∂ux =∂z can be obtained from the three equations. In the same way, partial derivatives of uy and uz can be calculated. With these partial derivatives, we can compute the strains and then the stresses from the generalised Hook's law.

7. Implementation of the crack propagation Under certain conditions, a crack inside a solid body can propagate to enlarge its size. The criterion for a crack to propagate has been outlined in Section 5 and here we outline the implementation procedure in FRACROD3D for the crack propagation. After the displacement discontinuities on the cracks (as well as the displacement and stress at interior points) are obtained, the crack front elements are iterated for a propagation criterion check, including the determination of propagation orientations of those front elements that would grow. If the crack does not propagate at a front element, then the front element remains as a front element. If the crack propagates at a front element, new elements are added in front of the front element and the old front element becomes an interior element. The new elements contain one or more new front elements. Orientation of the new elements is determined according to the procedure from the state of stress near the front edge of the front element as outlined in Section 5. The size of the new elements depends on the size of the old front element. From simulations of several examples, it has been found that use of the orientation determined from only one front element could lead to the propagation being very unsmooth: propagations of neighbouring elements could go in very different orientations and the new elements could possibly overlap, or severely twisted. This can be avoided if the propagation orientation of a front element is modified by taking into account the propagation orientations of its neighbouring front elements. In FRACOD3D, the modified propagation orientation of a front element combines the propagation orientations of three front elements: the considered front element and its immediate neighbours at both sides. In the combination, the importance of each element is controlled by weighting factors. The weighting factors are chosen as 0.6 for the considered front element and 0.2 for its two neighbours. Similarly, the size of the new elements can contain effects from the neighbouring front elements.

80

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

The number and layout of new elements depend on the size of the front element and propagation status of the neighbouring elements. For example, if the crack propagates at an element, but does not propagate at its neighbours on both sides, then four new elements, with three new front elements, are added. If the previous neighbour (the iteration direction is anticlockwise) propagates, two or four new elements, with one or two new front elements, are added depending on the size of the old front element and a potential front element is formed, which will be a new front element if the next neighbour does not propagate. Two new elements are added for small front elements, and four for large ones. The critical value for distinguishing small and large front elements is the maximum side length of elements in the initial discretisation. If the previous neighbour does not propagate, i.e. the next neighbour does propagate since the case of both neighbours not propagating has been considered, then one or two front elements are added. The edge in the propagation direction of a newly added element is set perpendicular to the front edge of the propagating front element. After the new elements and new nodes are added and proper boundary conditions are formed for the new system of the whole problem, a new round calculation for the displacement discontinuities is performed. The propagation criterion is checked again with the new results. The calculation is terminated if no further propagation is found. Strictly speaking, for more accurate results, if after the new elements are added, propagation test shows no further propagation, then the newly added elements at last step should be taken off and smaller elements should be employed, instead, for another round calculation and test. If the crack still does not grow, then elements with even smaller size should be used within certain limitation. It is noted that calculation of the combined stress intensity factors and combined toughnesses in Section 5 and using them in the crack growth check has taken into account of mixed mode loading. With the incremental step by step growth scheme, the crack at a front element grows in one mode at this step and it may grow in another mode at next step. Thus overall, the crack may propagate in a mixed mode at the front element.

8. Examples In this section, we present three examples simulated with FRACOD3D. The first two are used to verify the code with analytical solutions, while the third shows that the code can simulate crack growth. 8.1. Circular disc crack under tension As the first example for verification of the code, we investigate the problem with a planar circular crack, or penny-shaped crack, in an infinite elastic body, but we do not consider crack propagation. The circular crack lies in the xy-plane and the infinite body is under uniform tension in the z direction at far field. Young's modulus and Poisson ratio of the material is 5 GPa and 0.2, respectively. The tensile stress at far field is 10 MPa perpendicular to the crack plane. Radius of the crack is 3 m. The crack is divided into 144 triangular elements with 81 nodes, see Fig. 3. The induced normal displacement on the crack surface has the analytical expression [22] qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1  νÞT 0 a uz ðr; z ¼ 0Þ ¼ 1  ðr=aÞ2 ; πG where a is the radius of the crack and T0 is the value of tensile stress at far field. This expression is also true for the displacement due to the uniform internal pressure on the crack surfaces. Table 1

Fig. 3. Mesh for simulation of circular crack in tensile stress field.

Table 1 Normal displacement at five elements near the x-axis. r/a uz =a(10  3 ) (F) uz =a(10  3 ) (A) Difference (%)

0.1307 2.4157

0.3272 2.3198

0.5237 2.1136

0.7204 1.7527

2.4236

2.3101

2.0826

1.6958

 0.33

0.42

1.49

3.36

0.9169 1.1016 0.9758 12.89

shows the normal displacements on the crack surface at centroids of 5 elements from FRACOD3D (F) and those from the analytical solution (A) at the same distances from the centre of the crack to the centroids of the elements. The FRACOD3D underestimates the displacement near the centre of the crack, with very limited error, and overestimates the displacement away from the centre. The difference between results of the FRACOD3D and those of the analytical solution increases with the distance from the centre of the crack. Overestimation of displacement on the whole crack by DDM is observed by Kuriyama and Mizuta [8] and Ishijima et al. [23]. Ishijima et al. [23] observed that the displacement difference increases with the distance from the centre of the crack, while graphs in [8] show the opposite. Figs. 4 and 5 show the profiles of the displacement and the stress field around the crack on the xz plane. We have not been able to find a simple analytical expression for stresses around the crack. A careful check of the analytical values shown in [8,23] shows that the analytical values of σ zz near the crack on the x-axis are very close to the values from FRACOD3D.

8.2. Spherical cavity in an infinite body under inner pressure The second example for verification of the code considers deformation and stress state due to inner uniform pressure on the spherical cavity surface in an infinite body. In this case, the analytic radial displacement given in [24] is uR ¼

P i R a 3 4G R

and the normal stresses are  a 3 P  a 3 1 σ RR ¼  P i ; σ θθ ¼ σ φφ ¼ i ¼  σ RR : R 2 2 R

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

Fig. 4. Distribution of displacement on the xz plane around the crack.

Fig. 5. Distribution of principal stresses on the xz plane around the crack.

81

82

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

where a is the radius of the cavity, R is the radial distance from the centre of the cavity and Pi is the internal pressure on the surface of the cavity. In the simulation with FRACOD3D, the material parameters used are E ¼ 1 GPa and ν ¼ 0:25. The radius of cavity a¼3 m and internal pressure is P i ¼ 1 MPa. The mesh for the simulation is shown in Fig. 6, with 224 triangular elements and 114 nodes. With these values, the analytical radial displacement on the surface of the cavity is 1.875E-3 m (outward). Element #1 is close to the pole of the sphere and its centre is at (0.7362, 0.1464,2.8477). The radial direction from the centre of the cavity to this element is close to

the negative z direction, therefore we can compare the radial components of displacement and stress from the analytical solution to the z-components from FRACOD3D. The displacement and stress components at element #1 from FRACOD3D are ux ¼ 4.5543E-4, uy ¼ 9.1093E-5, uz ¼  1.9177E-3 and σxx ¼705589.06, σyy ¼576087.81, σzz ¼  929628.28, σxy ¼26736.66, σxz ¼ 340320.57, σyz ¼67677.98. It can be seen that there are differences between the analytical solution and the results from FRACOD3D. This is partly because the displacement and stresses from FRACOD3D are calculated at the centres of the elements, which are not exactly on the spherical surface due to the mesh discretisation. In addition, the stresses on the elements are computed with a numerical scheme [21]. The distance from the centre of the cavity to the centroid of the element #1 is 2.9450 m and with this value for R, the value of the analytical radial displacement will be 1.9820E-3 m. The magnitude of the displacement from FRACOD3D is 1.9731E-3 m. These two values are very close. Figs. 7 and 8 show the displacement field and principal stresses around the spherical cavity from code FRACOD3D, respectively. These figures show that the symmetry feature of the problem holds. To verify the code, we plotted the variations of normal stresses along the x-axis in the body from FRACOD3D and the analytical solution in Figs. 9 and 10. The x-axial direction is along the radial direction, so we have ux ¼ uR , σ xx ¼ σ RR ; σ yy ¼ σ θθ ;σ zz ¼ σ φφ . It can be seen that the stresses are very close to the analytical values. The radial displacement also agrees well with the analytical value (Fig. 11).

8.3. Propagation of circular disc crack under shear

Fig. 6. Mesh on spherical cavity surface.

FRACOD3D is based on the fundamental solution from a small planar crack and uses meshes on a crack. When a crack grows, new elements can be added on the old meshes and no re-meshing of the old part of crack is required. Furthermore, with some mathematical manipulations beforehand, the propagation criterion test

Fig. 7. Displacement field around the xz plane and spherical cavity.

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

83

Fig. 8. Principal stresses around the xz plane and spherical cavity.

Fig. 9. Variation of analytical and DDM radial normal stress along the x-axis. (The normal stress at x/a¼1 from FRACOD3D is calculated at centre of an element, which is close to, but not on, the x-axis).

can be done with DDs on the front edge elements. Thus it is suitable and efficient for the simulation of crack propagation. In the next example, we have a circular crack in an infinite body under shear stress at the far field which is parallel to the crack surface. Young's modulus and Poisson ratio are the same as those for the first example of the circular crack under tension. The toughnesses are given in Table 2. The radius of the crack is also 3 m. The crack is discretised into 216 triangular elements with 121 nodes. The mesh has the similar layout as that shown in Fig. 3. The far field shear stress is 1.0E8 Pa along the fracture plane and Fig. 12

Fig. 10. Variation of analytical and DDM circumferential normal stresses along the x-axis. (The normal stresses at x/a ¼1 from FRACOD3D are calculated at centre of an element, which is close to, but not on, the x-axis).

shows the crack configuration after a few steps of propagation with a displacement field around the crack, while the principal stresses around the crack are shown in Fig. 13. Generally, the variation of stresses on the crack surface is smooth. The propagation angle at the symmetrical fracture front in shear direction is 711, close to the value of the 2D case. The configuration of growing crack is similar to those observed in experiments [25] and other numerical simulations [26,27] with inclined penny-shaped crack under uni-axial tension or compression.

84

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

9. Conclusion We have outlined the fundamental solution for an infinite isotropic homogeneous linear elastic three-dimensional body subject to displacement dislocation or discontinuity (DD) over a

finite planar crack. We also discussed how the fundamental solution can be used to solve practical three dimensional (3D) engineering problems using the boundary element method. This scheme is called three-dimensional displacement discontinuity method (3D DDM). The basic formulations and equations for 3D DDM have been presented and they form the theoretical foundation for the computational code FRACOD3D. Three-dimensional crack growth has a very complicated mechanism. We proposed that a crack will propagate at its front where maximum tensile stress or maximum shear stress reaches some critical value. The propagation orientation of the crack follows the direction in which the stress is maximum. These crack growth features have been implemented in the three dimensional simulation code FRACOD3D. The basic aspect of the code is verified with two examples for which analytic solutions exist. Growth of a circular planar crack under shear is also simulated and the propagation angle at the symmetric line is close to that in the two dimensional case and growth configuration is similar to those observed in experiments and other numerical simulations under similar loading condition.

Fig. 11. Variation of analytical and DDM radial displacement along the x-axis.

Table 2 Material toughness and elasticity for crack propagation under shear. E

υ

KIc

KIIc

KIIIc

5.0e9 Pa

0.2

1.0e7 Pa m0.5

1.0e13 Pa m0.5

1.0e13 Pa m0.5

Acknowledgements We are very grateful to Hazama Corp., Japan and TEKES, Finland for the financial support towards the early development of FRACOD3D and the International Collaboration Project on Coupled Fracture Mechanics Modelling, (project team consisting of CSIRO,

Fig. 12. Configuration, with displacement field, of the crack after a few steps of growth.

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

85

Fig. 13. Configuration, with principal stresses, of the crack after a few steps of growth.

Posiva, KIGAM, SKEC, Geodynamics Ltd., Geomecon, FRACOM, SNU and ETH) for the financial support towards further development of the code. We would also like to thank Dr Binzhong Zhou for developing the graphic functions and interfaces for FRACOD3D. J. Shi also likes to thank Yang Shi for checking and editing the text.

References [1] Crouch SL, Starfield AM. Boundary element methods in solid mechanics. London: George Allen & Unwin; 1983. [2] Shen B. Mechanics of fracture and intervening bridges in hard rocks. (Ph.D. Thesis). Stockholm, Sweden: The Royal Institute of Technology; 1993. [3] Shen B, Stephansson O, Rinne M. Modelling rock fracturing processes. A fracture mechanics approach using FRACOD. Dordrecht: Springer; 2014. [4] Hong HK, Chen JT. Generality and special cases of dual integral equations of elasticity. J Chin Soc Mech Eng 1988;9:1–9. [5] Rongved L. Dislocation over a bounded plane area in an infinite solid. J Appl Mech 1957;24:252–4. [6] Wiles T.D., Curran J.H.. A general 3-D displacement discontinuity method. In: Eisenstein E., editor. Proceedings of the 4th international conference on numerical methods in geomechanics. Edmoton, Canada; June 1982. P. 103–111. [7] Cayol V, Cornet FH. 3D mixed boundary elements for elastostatic deformation field analysis. Int J Rock Mech Min Sci 1997;34:275–87. [8] 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:111–23. [9] 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:317–22.

[10] Li H, Liu CL, Mizuta Y, Kayupov MA. Crack edge element of three-dimensional displacement discontinuity method with boundary division into triangular leaf elements. Commun Numer Meth Eng 2001;17:365–78. [11] Vijayakumar S, Yacoub TE, Curran JH. A node-centric indirect boundary element method: three-dimensional displacement discontinuities. Comput Struc 2000;74:687–7003. [12] Wang F, Huang XC, Deng JH. Numerical analysis of crack under compression by 3D displacement discontinuity method. J Shanghai Jiaotong Univ 2010;15: 297–300. [13] Cheng AH-D, Detournay E. On singular integral equations and fundamental solutions of poroelasticity. Int J Solids Struct 1998;35:4521–55. [14] Ghassemi A, Tarasovs S, Cheng AH-D. A 3-D study of the effects of thermomechanical loads on fracture slip in enhanced geothermal reservoirs. Int J Rock Mech Min Sci 2007;44:1132–48. [15] Medina DE, Liggett JA. Exact integrals for three-dimensional boundary element potential problems. Commun Appl Numer Methods 1989;5: 555–61. [16] Wilde A, Dual A. Boundary element formulation for three-dimensional fracture analysis. Southampton: WIT Press; 2000. [17] Kassir MK, Sih GC. Three-dimensional crack problems. Leyden: Noordhoff International Publishing; 1974. [18] Vandamme L, Curran JH. A three-dimensional hydraulic fracturing simulator. Int J Numer Meth Eng 1989;28:909–27. [19] Rabczuk T, Bordas a, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Comput Mech 2007;40:473–95. [20] Bower AF. Applied mechanics of solids. Florida: CRC Press; 2010. [21] Shi J, Shen B. Approximation schemes of stresses on elements for the threedimensional displacement discontinuity method. Eng Anal Boundary Elem 2014 (Accepted for publication). [22] Sneddon IN, Lowengrub M. Crack problems in the classical theory of elasticity. New York: John Wiley & Sons; 1969; 139. [23] Ishijima Y, Sato K, Kinoshita S. Application of the displacement discontinuity method to some crack problems. Theor Appl Mech 1980;28:167–86.

86

J. Shi et al. / Engineering Analysis with Boundary Elements 48 (2014) 73–86

[24] Saada AS. Elasticity: theory and application. Malabar, Florida: Krieger Publishing Company; 1993; 331–4. [25] Dyskin AV, Sahouryeh E, Jewell RJ, Joer H, Ustinov KB. Influence of shape and locations of initial 3-D cracks on their growth in uniaxial compression. Eng Frac Mech 2003;70:2115–36.

[26] Mi Y, Aliabadi MH. Three-dimensional crack growth simulation using BEM. Comput Struct 1994;52:871–8. [27] Mi Y, Aliabadi MH. Dual boundary element method for three-dimensional fracture mechanics analysis. Eng Anal Boundary Elem 1992;10:161–71.