A finite element code for fracture propagation analysis within elasto-plastic continuum

~

Engineering Fracture Mechanics Vol. 53, No. 2, pp. t93-21 I, 1996

Pergamon

0013-7944(95)00105-0

A FINITE PROPAGATION

ELEMENT

CODE

ANALYSIS

FOR

WITHIN

Copyright :f; 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0013-7944/96 $15.00 + 0.00

FRACTURE ELASTO-PLASTIC

CONTINUUM I. L. L I M t CSIRO Division of Minerals, Pinjarra Hills, Brisbane 4067, Australia I. W. JOHNSTON Faculty of Engineering, Victoria University of Technology, Melbourne 3000, Australia and S. K. CHOI CSIRO Division of Petroleum Resources, Syndal, Melbourne 3150, Australia

AbstractA new code for the simulation of mixed-mode fracture propagation within a body deforming

elasto-plastically is presented in this paper. The analysis is fully automated in that it requires no user interaction. An automatic local remeshing technique is developed based on the window remeshing technique of Murti [Numerical fracture mechanics using finite element methods. Ph.D.Thesis, University of New South Wales (1986)]. This algorithm ensures that proper element geometries are maintained in the crack tip vicinity and that reasonable shaped elements are preserved further away from the crack tip. A new approach was developed to predict the fracture trajectory under load or displacement control. This new approach is particularly suited to applications where the applied loads or displacements are the controlling factors rather than the crack extension. The code incorporates recently developed profile reduction and mesh rezoning techniques for increased efficiency. The capability of this new code is demonstrated by application to a number of practical problems,

1. INTRODUCTION IN RECEYTyears there has been an increasing rate of development of codes for fracture propagation modelling. Two major approaches are possible: one being the use of discrete fractures and the other through the use of smeared cracking. Discrete fracture modelling is the preferred technique for simulating the propagation of a small number of individual fractures over the use of smeared cracking [1] since it does not suffer from stress locking, spurious modes and directional bias [2]. To date, there are over 15 codes for modelling mixed-mode fracture propagation. These codes employ boundary element [3], finite element [4], displacement discontinuity [5] or hybrid [6] methods. The development of a new code was motivated by the need to model fracture propagation in a soft rock such as Melbourne mudstone or its synthetic counterpart Johnstone [7], both of which possess uniaxial compressive strengths typically less than 10 MPa. Since the inelastic characteristic of this soft rock in the compressive regime plays an important role for many geomechanics problems [8], it was necessary that the code could model both fracture propagation and elasto-plastic behaviour simultaneously. At the time this work commenced in 1988, there were no codes which possessed these capabilities. Thus, it was decided to develop the new finite element code EPIC (Elasto-Plastic finite element code Incorporating Crack propagation) with these capabilities. Since then, several other codes with similar capabilities have been developed [9-11]. The purpose of this paper is to describe the capabilities and limitations of EPIC. 2. THE EPIC CODE

In essence, EPIC was designed to simulate the propagation of multiple discrete fractures in an elasto-plastic continuum in a completely automated mode without requiring any user tTo whom all correspondence should be addressed. 193

194

I.L. LIM et

al.

interaction. Complete automation was desirable as EPIC is intended to model elasto-plastic problems where extensive computational run time may be required. Moreover, there was a lack of appropriate hardware (such as a workstation) and software available to the authors at that time and consequently, interactive modelling was not an option. As extensive computation was expected, the code was developed with efficiency as an important criterion. The code itself has three basic components. The first component consists of the algorithm for the stress analysis of an e|asto-plastic continuum, while the second component is the automatic remeshing algorithm which modifies the topology of the mesh to propagate multiple fractures. The third component consists of algorithms for computing the stress intensity factors (SIF), the fracture trajectory, bandwidth reduction, rezoning of the relevant data to the new mesh and other data

l Atom°,c--0eneratorI for pile problems

Input data for new increment lin applied load or displacement I

I

l Stress analysis

I

No

Reapply present increment in applied tractions and employ mesh variables from previous iteration

l

I

Determine location and geometry of new cracks to be initiated

Do for all cracks Calculate KI,KII , Determine whether each crack propagates and its

trajectory

USe data from

~A

'storage' arrays for this part of analysis

o

Store present geometrical and state variables into 'storage'

arrays

1

utomaUcremeshing Profile reduction

I

I Rezone mesh I

,,

~

I"Rearrange data base I ~

Fig. 1. Flow chart of EPIC code.

Finite elementcode for fracture propagation analysis

195

manipulation. A flow chart is provided in Fig. 1 to illustrate the relationship between the primary components of EPIC. In many codes, the algorithm is crack length controlled (e.g. ref. [11]). In such approaches, the crack is extended by a predetermined length and the applied load or displacement is adjusted until the SIF at the crack tip falls just below some critical value assumed to be necessary for further fracture propagation. Such an approach is particularly suited to problems where crack growth is the primary controlling factor. In EPIC, the algorithm can be either load or displacement controlled. This approach was chosen to reflect the displacement controlled problems which are addressed elsewhere [12]. Consequently, a different approach to that outlined by Ingraffea [13] had to be employed. In this approach, it was desired that for a given load or displacement increment, the predicted crack extension length would be incrementally adjusted until it converged. It is assumed here that stable fracture propagation occurs. It should be noted that the crack is defined to be stable or unstable with respect to a given load or displacement increment. Stable fracture propagation is considered to be maintained if the computed crack extension length converges to a value less than that predefined by the user. Otherwise the crack is considered to be unstable. Only when the predicted crack extension length has converged will the next load/displacement increment be applied. Thus, with any load/displacement increment, the fracture path is sought by an iterative scheme. The primary steps for a displacement controlled analysis, and as schematically outlined in Fig. 1, may be given below: (1) For a given displacement increment i, stress analysis is performed on the fractured continuum and the SIF is computed for all existing crack tips. (2) With the computed SIFs, the algorithm decides whether any given crack is to propagate and its direction of propagation according to a selected mixed-mode fracture criterion. (3) Should any crack propagate, the crack is tentatively extended a small length predetermined by the user. (4) Remeshing is performed on the mesh corresponding to the previous displacement increment i - 1, to extend the crack along its predicted trajectory. (5) Nodal renumbering to minimise the profile of the global stiffness matrix is performed. Rezoning of the state variables from the old mesh, corresponding to the increment i - 1, to the new mesh is affected. (6) The displacement increment i is reapplied. (7) Steps 1 to 6 are repeated in an iterative manner and for each iteration, an improved estimate of the crack extension length is made using a predictor-corrector algorithm to be described in greater detail later. This continues until the predicted crack trajectory converges, i.e. it does not differ significantly from the previous prediction. The details of this procedure will be elaborated upon later in this paper. (8) Once convergence is obtained, the present geometrical and state variable data associated with the increment i are stored into "storage" arrays for use in subsequent displacement increments. (9) A new displacement increment i + 1 is applied and the entire process of steps 1 to 9 is repeated until the entire analysis is completed. Should unstable fracture propagation occur, steps 1 to 7 are performed iteratively with the exception that convergence of the predicted crack trajectory is not sought. Instead, the crack is incrementally extended by the user-predetermined length until the algorithm detects that stable crack propagation is once again attained. Once this occurs, the steps revert to that outlined previously to seek a converged solution. In the rest of this paper, the various components are described in greater detail. 3. FINITE ELEMENT FORMULATION The finite element formulation in EPIC was based upon ELASP [14, 15] which employs the FEABL profile solution scheme developed by Orringer and French [16]. Substantial modifications were undertaken, leading to efficiency improvements in excess of 25%. The code provides plane stress, plane strain and axisymmetric analysis capabilities. A combination of eight-noded

196

!. L. LIM et al.

quadrilateral and six-noded triangular quadratic serendipity elements can be employed. Full (3 × 3 point) and reduced (2 × 2 point) integration is available for the quadrilateral elements while three-point full integration is employed for the triangular elements. A Mohr-Coulomb elasto-plastic constitutive model with non-associated flow rule options [17] is available. The integration of the path dependent elasto-plastic constitutive relationship is performed with an iterative process based on the work of Nayak and Zienkiewicz [18]. Further details of the formulation employed in EPIC and the reasons for its choice may be found in ref. [1]. Verification of the elasto-plastic model in EPIC was performed on some benchmark problems, e.g. the thick-walled cylinder problem assuming a Tresca failure criterion [19] and Mohr-Coulomb material assuming both associated and non-associated flow rules [20]. Further verification was performed against theoretical collapse loads of rigid circular footings on purely cohesive and cohesive frictional materials [21].

4. AUTOMATIC WINDOW REMESHING 4.l. Basic approach One of the unique features of any discrete fracture propagation finite element code is its remeshing algorithm. To date, there are basically three remeshing strategies that have been reported in the literature. The first involves the distortion of the existing mesh in order to accommodate the mixed-mode fracture trajectory. This mesh distortion may be achieved either through slight modifications in the mesh topology or via mapping techniques [22]. The second approach is to remesh the entire mesh as advocated by Shephard et al. [23]. Similar approaches were employed by Sumi [24] and Taniguchi [25, 26]. The advantage of this approach is that well-shaped elements can usually be generated. However, its disadvantage is that a large number of state variables must be transferred from the old to the new mesh. The transference procedure tends to be error prone [27] and can be computationally expensive, particularly when an elasto-plastic continuum is involved. Remeshing over a large region will also naturally require more computational effort. The third strategy employs a local remeshing technique whereby only the region of elements within a certain distance of the crack tip is modified [13]. This approach has the advantage of requiring only the transference of state variables from a small region which is to be remeshed. It follows that any errors involved in the transference process are confined to the region and the computational effort is reduced. However, the generated elements may not be as well shaped as in the second approach since there tend to be greater constraints on the placement of the new elements due to the existing mesh surrounding the region to be remeshed. EPIC employs the third approach as it was desirable to reduce the computational effort and any errors during the transference process. Another important consideration in any such remeshing algorithm is the application of elements to adequately model the crack tip singularity [28]. For completeness, the relevant recommendations from Lira et al. [28] are recast here. (1) The smallest possible quarter point element (QPE) size should be utilised to minimise its limited non-singular stress field representation. (2) When fewer element layers are preferred, transition elements [29] can be employed as long as their radial size is within the range 1 < LT,,~:/LsRz,,~-tl <_ 10 for adequate singularity representation. Lr1~,,and LsRz,,~- v, are as defined by Lim et al. [28], where Lr,,~,, is the length of the transition element layer i away from the crack tip and LsRz,,i-~,, is the total length of elements with singular representation capability between the crack tip and element layer i. Transition elements beyond this size range may incur excessive errors due to large element distortion. (3) Normal elements should be employed ifa denser mesh around the crack tip is desired where LN'jl/LsRz',i- Ll ~ l. This arrangement should enable a good representation of both the singular and non-singular terms. (4) The optimal QPE shape and arrangement appear to be adequately realised by a rosette of equilateral triangular elements. (5) A smooth increasing progression of element sizes radially from the crack tip is preferred.

Finite element code for fracture propagation analysis

197

2 adjoining elements

-/

-

_

window Jsub-boundary

8 noded

element

/.

%

transition element

QPE

f/

window boundary window

\

~ - - - boundary node

f window sub-boundary crack

median

Fig. 2. Definitions of window region.

The remeshing algorithm employed in EPIC is based upon the work of Murti [30], and outlined by Valliappan and Murti [31]. The algorithm by Murti [30] works primarily for rectangular shaped meshes, and it was not efficient and not entirely robust. Therefore, substantial improvements were made to Murti's algorithm to achieve a more versatile, robust and efficient algorithm. Additional routines were added to increase the versatility of the code so that it could manage more diverse situations. Robustness of the code was improved by developing better algorithms wherever possible. Efficiency was obtained by employing optimised coding techniques, minimising the input-output and streamlining the code wherever possible. The algorithm employs a window which bounds a region of elements surrounding the crack tip (Fig. 2). This window is a region bounded by window nodes as illustrated in Fig. 2. The window boundary itself is divided into window sub-boundaries. In the original algorithm, the window must always remain a convex region, otherwise inadmissible elements may be generated. This requirement has been relaxed in this new algorithm to allow non-convex regions. However, convex regions are still preferred as they generally lead to the generation of better shaped elements. The definition of the window sub-boundaries are as depicted in Fig. 2, whereby they adjoin the edges of two adjacent elements. By its very definition, the window sub-boundary cannot include elements which are connected to the window region by only one window node, e.g. an element connected only to the corner of the window as depicted in Fig. 3. This figure also illustrates the location of window sub-boundaries around a window region. This definition was devised to promote the convexity of the window region. Basically, as the fracture extends, the window advances together with the crack tip to its new location (Fig. 4). When a new fracture extension is to be performed, all the elements and nodes within the window boundary are absorbed as shown in Fig. 4(a). A window sub-boundary that

®

~

,~

....

~

window sub-boundaries

f

crack

@ elements excluded from

.J,

i (~

being absorbed for the next windowadvancement.

Fig. 3. Definition of window sub-boundaries.

I.L. LIM et al.

198

windownodes ../ ~ \ _ ~___~ = - - ~'~. T ~l...... 4t, "~.,~ / ~

windowsub-boundary in directionof crack

all elementsand directionof nodeswithinwindow crackextension are absorbed backelements newwindow !(~

t

-~.J

"~'~-~ new window

(b)

~". . . . . . . _ (c)

!

-_-_-_~-_-_

" ~ ~ ~

.~

j generationof new elements withinwindow

Fig. 4. Steps in advancement of a window region.

lies in the direction of crack extension is identified as shown in Fig. 4(a). The window then moves by absorbing the two adjacent elements adjoining the window sub-boundary and its corresponding nodes as illustrated in Fig. 4(b). New elements are subsequently generated behind the redefined window to fill the gap left behind by the advancing window and are designated as back elements. The last phase involves the generation of new elements within the new window location as depicted in Fig. 4(c). Further details on the procedures of window remeshing are provided below for various situations. 4.2. Self-similar extension (mode I) When a crack extends in a self-similar mode, there are three possibilities catered for by the algorithm (Fig. 5). 4.2.1. The new crack tip is located within W - (the rear portion of the window). This occurs when the new crack tip has not advanced beyond the window median. The window median is as shown in Fig. 5(a). In this case, the window is not redefined. The elements and nodes within the window are absorbed and regenerated as pictured in Fig. 5(b) to advance the crack tip to its new position. This strategy enables good element shapes to be maintained within the window and avoids the unnecessary generation of back elements. 4.2.2. The new crack tip is located within W + (the front region of the window). The elements and nodes within the window are first absorbed and the window advanced by absorbing two adjacent elements and its corresponding nodes in the direction of the advancing crack as depicted in Fig. 5(c). Back elements are created behind the new window and new elements are generated within the redefined window.

Finite element code for fracture propagation analysis

199

4.2.3. The new crack tip extends beyond W ÷. The total extension is sequentially divided into smaller segments with the first segment approximately reaching the point where the fracture trajectory intercepts the original window boundary as shown in Fig, 5(d). The same procedure as in Section 4.2.2 is then invoked. The remaining extension is then progressively divided with the same procedure until no segments are left, as depicted in Fig. 5(d and e). 4.3 Mixed-mode extension There are three different procedures of window remeshing depending on the orientation of the trajectory with respect to the window region. Essentially, the basic approach of these procedures is similar to that in Sections 4.2.2 and 4.2.3. These are depicted in Fig. 6. Figure 6(a) depicts a non-self-similar crack extension towards the front of the window region. Figure 6(b and c) illustrates the crack extension to the left and right (viewed from the initial crack direction) regions of the window. Back elements are generated accordingly to fill the space left by the advancing window. The shape of these back elements are checked against a number of criteria to ensure they are not excessively distorted. Should that occur, the element may be divided into two triangular elements, if it leads to better shaped elements. Figure 7 depicts the simple criteria employed to achieve this. When the largest interior angle of a quadrilateral equals or exceeds 120° and its opposite corner has an interior angle equalling or exceeding 30 °, the quadrilateral is divided into two triangles (Fig. 7). The same applies should the largest interior angle and its opposite interior angle both equal or exceed 100°. Division into triangular elements is also performed if the largest interior angle equals or exceeds 150°. These simple criteria have been devised as they provide a

median

/.W +

W - .' (a)

median

~ I

'

~

\

(b) new ~ndow

(c) ~, ~

median

~

/ window median- . ~ - ~ /':/ ~ 'Y""'~/~ \ ti/~'~

/~// ~ ' ~

~

window ~ . sud-boundary indirection / of crack extension

7,, 1'1 - - 4 --adjacent "x,~',~ / elements

-i- ......... -i -! ....... 1st segment 2nd segment

(d) ,,

(e)

J1 last segment Fig. 5. Self-similar crack extension.

window

200

I. L. LIM et al.

/./new wind_ow __

(a)

(b)

(c)

Fig. 6. Mixed-mode crack extension. compromise between a simplistic criterion (such as employed by Valliappan and Murti [31]) and more complex criterion. Experience has shown these to be adequate. 4.4. Remeshing in the window region A totally different procedure from that of Murti [30] was employed in order to implement the recommendations of Section 4.1. As noted in Section 4.1, a critical parameter is the relative size of elements in the radial direction from the crack tip. The size of each element should grow progressively, as a function of its radial distance from the crack tip. The algorithm employed here seeks to achieve this objective. The algorithm consists of three stages as depicted in Fig. 8. In the first stage, a rosette of elements is generated. This rosette is similar to that employed by Wawrzynek [32]. The second stage encompasses the temporary generation of superelements that connect the rosette of elements to the window nodes. The last stage involves the division of the superelements into elements. The first layer in the rosette of elements consists of eight identical equilateral QPEs, surrounded by successive layers of transition elements as shown in Fig. 8(a). The singular element >_.120 °

~1 >--150°

>~ 30 °

>..100 °

~P1 >-120°

>-. 30 ° = largest interior angle

Fig. 7. Criteria for division of quadrilateral elements into triangular elements.

201

Finite element code for fracture propagation analysis

(a)

(b)

(c)

Fig. 8, Stages in generationof elements in windowregion. mid-nodes in the direction radiating from the crack tip are positioned according to the formula developed by Lynn and Ingraffea [33] to model the r ,..2 singularity in the radial direction. The maximum number of transition element layers may be defined by the user. The actual number of layers employed is dependent on the relative size of the rosette with respect to the window. It is desirable that the rosette of elements encompasses as large an area of the window region as possible, since in the rosette the elements are optimally (or near-optimally) shaped and sized for most practical applications. If the maximum number of transition element layers specified by the user is exceeded, then normal element layers are employed. Usually, transition elements are employed in preference to normal elements in the immediate vicinity of the crack tip, This is because for similar accuracy to be attained, larger size transition elements can be employed in place of normal elements due to their better singularity modelling capability [28]. This leads to the use of fewer elements and hence, better computational efficiency. However, this is at the expense of reducing the transition elements' ability to model linear strains in the radial direction. Thus, further from the crack tip where the singularity becomes less dominant, normal elements may be the more appropriate choice. A number of other variables may also be defined by the user, e.g. the size of the quarter-point elements (Lo) and the relative sizes of each surrounding layer of transition and normal elements (Ls). A simple technique was sought to generate the elements between the rosette and window boundaries such that a reasonable degree of control was possible on the element sizes in the radial direction from the crack tip. It was also preferred that iterative techniques, such as mesh smoothing, be avoided to reduce computational costs. This was achieved through the use of an intermediate step involving the generation of temporary superelements. These superelements are generated sequentially in a clockwise direction around the window region and consist of quadrilateral or triangular elements. Each superelement connects the window nodes to the rosette as shown in Fig. 8(b). The shapes of the superelements are important considerations as they would constrain the shapes of the final elements produced. A simple strategy was employed to obtain acceptably shaped superelements as illustrated schematically in Fig. 9 and outlined below: (1) Construct a potential superelement, such as a triangular or quadrilateral superelement, first. A second potential superelement adjacent to the first in the clockwise direction is then constructed. The second superelement may be either triangular or quadrilateral in shape. Figure 9 depicts some of the permutations considered by the algorithm. (2) A distortion shape factor D is calculated for each potential superelement. The factor D may range from 0 which represents an optimally shaped superelement to 1 which represents an unacceptably shaped superelement. The distortion shape factors for each combination of two adjoining potential superelements are added together. The combination with the lowest D, and hence the combination with the overall best shaped superelements, is chosen. (3) The first superelement in the combination is then generated. (4) The next combination of superelements consisting of the second and third potential superelements is then constructed and selected according to the procedures 1 to 3. This process continues until all superelements are generated within the window region. The distortion shape factor D referred to above is described by the interaction of various shape

1. L. LIM et al.

202

parameters similar to that developed by R o b i n s o n [34, 35]. F o r a quadrilateral, these shape parameters are the: aspect ratio (A) = max(edf~, ~/e2) skew (S) = {tan

L(e3/J~+fde2)}/(0.5rc)

taper in the x-direction (T~) = J~/.]~ (1)

taper in the y-direction (T,) = e4/e2, where e~ = 0.25(x, + x2 + x~ + x4),

~ = 0.25(yt + y2 q- y3 q- y4)

e2 = 0 . 2 5 ( - - x l

J~_ = 0.25(--y, + y: + y3 - - y4)

q- x2 q- x3 - - X4),

e3=O.25(-x,-x2+x3+x4),

[;=0.25(-y~-yz+y3+y4)

e4 = 0.25(x, - x2 + x3 - x4),

J4 = 0.25(y, -- y2 + y3 - y4)

rosette of elements

crack

(a)

(b)

(c)

(d)

(e)

(f )

(0)

(h)

inadmissible _ /element

(i)

(i)

Fig. 9. Strategy to obtain acceptably shaped superelements.

(2)

Finite element code for fracture propagation analysis

203 tentative division of superelements

clockwise direction of superelement division

~Y

tentative division of superelements

i

-7 3

superelement

actual

division of superelements

1

Y

/

~

transition

\ .,_~ element

"~QP'/ ~ "crack E global axes

division to form better shaped elements

(c)

.-X

Fig. 10. Axes employed in determining shape distortion factor D.

Fig. 11. Division of superelements into elements.

and x,, y, are local (transformed global) nodal cartesian coordinates of the quadrilaterals. A visual appreciation of the physical significance of the coefficients ej a n d f may be found in Lim et al. [27]. The basic difference of this definition from the work of Robinson [34, 35] is that he defined xj and y, as global nodal coordinates whereas here they are defined as local coordinates. The local axes are as defined in Fig. 10, with the y-axis passing through the crack tip. This modified definition allows skew and taper to be defined with respect to the radial direction emanating from the crack tip. This is of particular importance to controlling the element shape in the radial direction. The distortion shape factor D is now described by O = W + W,F~ + W,,.~,. + W,,FT, + W , ~

(3)

W + W~F, + W,,F,, + W,,F,, + W,F,, ' where the factors F,, F,,, F,, and F, are functions of S, T,., T,. and A, respectively. These factors fall within the range {0,1}. The weighting factors W,, IV,,, IV,, and W, correspond to the shape parameters S, T,., ~ and A, respectively. These factors allow varying emphasis (i.e. importance) to be associated with certain shape parameters. The weighting factor W serves as a toggle which may be employed to avoid quadrilaterals with internal angles exceeding a given angle. Presently, the critical angle is set at 150 °. In the last phase of the window remeshing, the superelements are tentatively divided into elements in the radial direction as depicted in Fig. 1 l(a). The final division of the superelements is performed sequentially in a clockwise direction as shown in Fig. 1 l(a). The actual manner of division is governed by the resulting element shapes as shown in Fig. l l(b). Whenever poorly shaped quadrilateral elements are present, they are divided into two better shaped elements as

204

I.L. LIM et al. free

/

surf~ (a) ;J" • principal/~~ stress / ] direction ~ ~ _ ~

(b)

,7

crack initiated

I

Gauss point with highest principal tensile stress free _ _7~-. surface ~___~~] I ..... (c) ,~~ i

!

new window region

crack initiated

windowregion Fig. 12. Fracture initiation at free surface. shown in Fig. 1 l(c), according to the criteria given in Section 4.3. This strategy generally leads to reasonably well-shaped elements and avoids an iterative process. 4.5. Fracture initiation and termination Cracks may be initiated or terminated at free surfaces. When the principal tensile stress in an element adjacent to a free surface exceeds the specified material tensile strength, EPIC is able to generate a crack extending from the surface (Fig. 12), At present, the criterion for calculating the position on the boundary where the crack first forms is based on the Gauss point with the highest principal tensile stress and in the principal tensile stress direction as depicted in Fig. 12(a). The point on the mesh boundary where the crack initiates may vary at any point along the element edge as illustrated in Fig. 12(b). The crack may also form from boundary corners as shown in Fig. 12(c). More refined criteria may be implemented but this basic model has proved to be adequate for the purposes of this present work, When a running fracture intersects a free surface, the fracture effectively terminates. Figure 13 illustrates an example of a fracture termination as modelled by EPIC. 5. FRACTURE TRAJECTORY Determination of the fracture trajectory is an essential component in discrete fracture propagation modelling. The first step requires the accurate computation of the SIF. Presently, five displacement-based SIF computation techniques have been implemented in EPIC. They are the: (1) quarter-point displacement technique (QPDT), (2) displacement correlation technique (DCT), (3) displacement extrapolation technique (DET),

I. . . . . . . . .

free ~-7 surface

i window 1,,~-~- region

i i

!1

'

crack terminated /

I

~

2V. °" (a)

(b)

Fig. 13. Fracture termination at free surface.

I

Finite elementcode for fracturepropagationanalysis

205

(4) reduced displacement extrapolation technique (RDET) and (5) limited displacement extrapolation technique (LDET). The details of these techniques and their relative performances can be found in refs [28, 36, 37] and will not be elaborated here. It should be sufficient to emphasise the points made in refs [28, 36, 37] that these techniques are easy to implement and they require little computational effort. In particular, the QPDT technique has been found to be capable of providing excellent accuracies especially when the guidelines given in Section 4.1 are adhered to. Once the SIF has been computed, some theoretical criterion is required to predict the direction of fracture propagation of the crack for a given material and loading condition. In EPIC, four mixed-mode fracture criteria are available to predict whether the crack will extend and its direction of propagation. These are the: (1) (2) (3) (4)

strain energy density criterion, Sm~,[38] maximum energy release rate criterion, Gm,x[39] maximum circumferential tensile stress criterion, a0.... [40] an empirical criterion [41] given by

+ \K,,cJ

1

for the fracture envelope combined with either the Smm, Gm,xor a0maxcriterion to predict the fracture extension direction. The parameters fl and 7 are material-related constants. In all the above criteria, the equations derived from these theories cannot be solved directly but require an iterative procedure to locate the appropriate maxima or minima. The Sm~°criterion, in particular, has been known [42] to possess some local minima which may be overlooked if the iterative step is not fine enough. The algorithm implemented in EPIC examines the variation of the pertinent quantity around the crack tip in small angular increments. The rate of variation (d/d0) and the rate of change (d2/d02) in the appropriate quantity are also determined as additional checks. Once the direction of crack propagation has been determined, an estimation of the length of crack extension is required. There is a number of approaches which has been employed to determine the crack extension length. A simple method is to calculate the extension length from the fracture criterion itself such as utilised by Murti [30] with the SED criterion. A similar approach was taken by Yehia and Shephard [43] with the NT criterion. However, such an approach inherently assumes that the crack extension does not alter the geometry of the problem nor the applied stress field significantly. This is only true if the incremental crack extension itself is very small with respect to the problem being analysed. The use of such small increments may be computationally expensive. This simplistic approach is unlikely to be applicable in many practical problems where the propagating crack significantly alters the geometry of the structure. A better approach is one advocated by Ingraffea [13] whereby the applied load or displacement is adjusted until the stress intensity factors at the crack tip fall below the critical fracture propagation threshold. This approach is able to account for significant changes in the geometry of the structure associated with the propagating crack as well as the changes in the stress field. Another approach employed by Saouma [44] was to compute the variation of energy release rate with crack extension length along a predicted trajectory. The final crack extension length occurs at the point where the energy release rate equals the crack surface energy. Both approaches require the crack length to be chosen first and may be defined as crack length controlled algorithms. EPIC employs a different approach, which is controlled by the applied load or displacement. Essentially, in this approach the applied load or displacement is incremented and the crack length is predicted by an iterative procedure as outlined in Section 2. The basic iterative algorithm for prediction of the crack extension length is now presented below in greater detail. A relatively simple algorithm was developed. To elucidate the functioning of the algorithm, the following example is provided. At a given displacement increment i, if the magnitude of the computed fracture quantity Q exceeds a critical value Q~, it is recorded in variable QPi • QP~ denotes that this quantity is more positive, i.e. larger than Q~. A predetermined fracture extension length dp is assumed initially as a starting point. This is followed by remeshing and associated data

206

1. L. LIM et

al.

Q QPi (a)

ec

QNi = Distance,d

dc

'~ dpn

QP i

(b)

Qe QNi I'

=" Distance,d

dc dcp

t"

dPn

Q QP i , ~ _ . . . . . QN i

(c)

de dc p

j~

~1 .

Qc

P

Distance,d

I 7dpn-I Fig. 14. Procedure in determining crack extension length.

Fig. 15. Deformed finite element mesh of simulated mixed-mode SCB test.

Finite element code for fracture propagation analysis

207

numerical rangeof experimental results

123

I 32

front

back

I

32

front

I

centre line ----,-

23

back

centre line

notch tip

notch tip

Fig. 16. Comparison between predicted and experimental fracture trajectories.

~XZ2

ring

crack

~ ~ j sheared .~'/~'/// ~ I ~ ' ~ zone // ~'

(//

/ ~ "<:-)

fan shaped wedge radial crack fan

Fig. 17. Development of failure mechanisms in embedded end-bearing piles in soft rock.

208

I.L. LIM et al.

manipulation. On reapplication of the displacement increment, if Qj remains larger than Qc, it is assumed that the fracture is unstable and the procedure outlined in Section 2 is taken. However, if Q, is now smaller than Q~, it is recorded in QN~, which denotes it is negative with respect to Q~. It is now assumed that the crack tip location where Q~ equals Q~ lies somewhere between QP~ and QNi, and may be approximated by a linear relationship as depicted in Fig. 14(a). In Fig. 14(a), the distance between QP, and QN~ is denoted by dr.... The predicted crack extension length is now given by d~ as shown in Fig. 14(a) where

(5) and d,7, is the distance from the crack tip to the location of QPj. In Fig. 14(a), d,r is effectively zero since QPj was computed at the crack tip. Remeshing is now performed with the new crack extension length of dc and the displacement increment i is reapplied, if the newly computed Q~ is now larger than Qc, it is recorded into QPi as illustrated in Fig. 14(b). The distance between QNi and QP~ is now considerably shorter than dr,,. Another estimate is now made on the predicted extension crack length dc. Remeshing is again performed on the new extension length and the displacement increment is reapplied. Now it is assumed that the computed Q, is less than Q~ and is recorded into QN~ as illustrated in Fig. 14(c). It is obvious from Fig. 14(a-c) that the distance between QP, and QN~ is being progressively shortened. Convergence occurs when the length dp,, falls within a specified tolerance. The value of dc at convergence would be the final predicted crack extension length corresponding to increment i of applied displacement. Experience indicates that this simple algorithm typically requires ca 3-6 iterations to achieve convergence. 6. PROFILE R E D U C T I O N

As pointed out in Section 3, EPIC employs a profile solution scheme. Thus, the size of the profile (and bandwidth) becomes a critical consideration for the computational effort associated with inversion of the global stiffness matrix and subsequent solution of the matrix equations. The size of the profile itself is governed by the sequence of the nodal numbering. On completion of the automatic rcmcshing, the nodal numbering will have been rearranged according to the nodes reemployed and/or generated. The resulting profile corresponding to the nodal sequence may have been adversely affected leading to greatly increased computational costs. Consequently, it is imperative that some form of profile reduction is applied, This led to a comparison of various existing algorithms for profile (and bandwidth) reduction. The technique by SIoan [45] was found to provide the best overall results in that it required very little computational effort and led to excellent profile reductions. Details of that study are contained in ref. [46]. Thus, after each remeshing, the Sloan algorithm is invoked to reduce the profile. It was found that this algorithm performed very well, leading to typical profile reductions in the range of 50%. 7. M E S H R E Z O N I N G

Once the new mesh has been generated, the state variables such as stresses and displacements from the old mesh must be transferred to the new mesh. This rezoning process can be error-prone. The techniques developed in refs [27, 47] were employed at this stage for determining the location of the new nodes with respect to the old mesh and for the mesh rezoning. However, only displacements may be rezoned with this technique. Other techniques [48-54] would be required to rezone the stresses. However, as far as the authors are aware, all these techniques are subject to accumulation of error since interpolation techniques are typically employed which require some empirical assumptions regarding the spatial variation of the stress field. In EPIC, it was assumed that the crack tip is embedded in a linear elastically deforming region. Thus the stresses may be easily computed from the nodal displacements [55]. This assumption is adequate for the purposes of the fracture propagation investigations described in refs [12, 56].

Finite elementcode for fracture propagation analysis

209

0.13 MP15 Fracture

ring crack 0.10

elastic region

E

.E

O~

-10.07 p 1 2 3 4 5 6

6

0.04

0

I 0.025 Radius (m)

(mm) 0.2 1.0 2.0 3.0 4.0 5.0 0.05

Fig. 18. Simulated developmentof failure mechanisms.

8. APPLICATION This code was applied to the simulation of fracture propagation in semi-circular specimens of soft rock subjected to three-point bending, known as the SCB fracture toughness test [7, 57]. Details of the experimental technique and results may be found in refs [7, 58]. Five SCB tests were modelled [1] assuming a linear elastic material, encompassing a wide range of mixed-mode loading conditions. Figure 15 shows two stages in the simulation of the fracture test I32 with an inclined notch at 54 °. The deformed meshes are shown with a 10-fold magnification. Figure 16 compares the simulated fracture trajectory against the range of experimental results. Clearly, good predictions were obtained. EPIC was subsequently applied to the modelling of several end-bearing pile problems in soft rock [1]. The typical final failure mechanisms observed experimentally are depicted in Fig. 17 for an embedment ratio of 2. The sheared region is considered to form due to inelastic deformations in the rock. Figure 18 depicts the result obtained from the finite element analysis which suggests that the simulated formation of ring crack and yielded regions corresponded very well with experimental observations. These results indicate the capability of EPIC in modelling of fracture propagation within an elasto-plastic continuum. 9. SUMMARY In this paper, an overall description of the development and capabilities of EPIC was discussed. EPIC is a finite element code for the modelling of multiple fractures propagating in a mixed-mode trajectory. It is able to model the fracture propagation in an elasto-plastic continuum while the region in the vicinity of the crack tip is assumed to remain linear elastic. The analyses are completely automated. The remeshing technique is a substantially improved version of the window remeshing technique originally developed by Murti [30]. This technique ensures that proper

210

I.L. LIM et al.

e l e m e n t g e o m e t r i e s a r e m a i n t a i n e d in the c r a c k tip v i c i n i t y f o r s i n g u l a r i t y m o d e l l i n g a n d t h a t r e a s o n a b l e e l e m e n t s h a p e s are p r e s e r v e d e v e r y w h e r e else. T h e f r a c t u r e t r a j e c t o r y is p r e d i c t e d w i t h a n e w a p p r o a c h w h i c h is l o a d / d i s p l a c e m e n t c o n t r o l l e d . T h i s a p p r o a c h is p a r t i c u l a r l y suited to p r o b l e m s w h e r e a p p l i e d l o a d / d i s p l a c e m e n t are the c o n t r o l l i n g f a c t o r s r a t h e r t h a n the c r a c k e x t e n s i o n . T h e a p p l i c a t i o n o f r e c e n t l y d e v e l o p e d profile r e d u c t i o n a n d m e s h r e z o n i n g t e c h n i q u e s was also d e s c r i b e d . T h i s p a p e r c o n c l u d e d w i t h a b r i e f d e s c r i p t i o n o f p r o b l e m s to w h i c h E P I C has b e e n successfully a p p l i e d , t h u s d e m o n s t r a t i n g its c a p a b i l i t y to p e r f o r m f r a c t u r e p r o p a g a t i o n a n a l y s i s w i t h i n a b o d y d e f o r m i n g e l a s t o - p l a s t i c a l l y . Acknowledgements--We wish to express our gratitude to Dr V. Murti for assistance in understanding his fracture

propagation code. I. L. Lim wishes to acknowledge the financial support provided by the Monash Graduate Scholarship.

REFERENCES [1] I. L. Lira, Fracture propagation in soft rock. Ph.D. Thesis, Monash University, Australia (1993). [2] J. G. Rots, Computational modelling of concrete fracture. Ph.D. Thesis, Delft University of Technology (1988). [3] M. Alcantud, M. Doblare, L. Gracia and J. Dominguez, Analysis of the contact problem between elastoplastic 2-D bodies by means of the BEM. Proc. 1st EUROBEM Conf., Niza (1989). [4] P. A. Wawrzynek and A. R. Ingraffea, Interactive finite element analysis of fracture processes: an integrated approach. Theor. appl. Fracture Mech. 8, 137-150 (1987). [5] G. X. Sun, B. N. Whittaker and D. J. Reddish, Computer simulation of crack propagation in proximity to a tunnel by displacement discontinuity method. Trans. Inst. Min. Metall. 10, Al14-122 (1992). [6] J. Nehme and W. Gerstle, Coupled finite and boundary element modelling of 2D cracks. Proc. 21st Midwestern Mech. Conf. (Edited by J. B. Ligon et al.), pp. 327-328 (1989). [7] I. L. Lim, I. W. Johnston, S. K. Choi and J. N. Boland, Fracture testing of a soft rock with semi-circular specimens under three-point bending. Part l--mode I. Int. J. Rock Mech. Min. Sci. and Geomech. Abstr. 31, 185-197 (1994). [8] C. M. Haberfield and I. W. Johnston, Numerical modelling for weak rock. Computer Methods and Advances in Geomechanics, 1159-1164 (1991). [9] T. N. Bittencourt, Computer simulation of linear and nonlinear crack propagation in cementitious materials. Ph.D. Thesis, Cornell University (1993). [10] D. V. Swenson, Private communication (1991). [11] A. Carpinteri and S. Valente, Analysis of mixed-mode fracture in concrete. ASCE J. Engng Mech. 115, 2344-2347 (1989). [12] I. L. Lim, S. K. Choi and I. W. Johnston, Modelling of failure mechanisms in an indentation problem. 1993 Conference on the Applications of Computers in the Mineral Industry (pp. 151 158), Wollongong, Australia (1993). [13] A. R. Ingraffea, Numerical modelling of fracture propagation. Rock Fracture Mechanics (Edited by H. P. Rossmanith), pp. 151-208. Springer, Wien (1983). [14] I. B. Donald, H. K. Chiu and S. W. Sloan, Theoretical analyses of rock socketed piles. Int. Conf. Foundations on Rocks (pp. 303-316), Sydney (1980). [15] H. K. Chiu and I. B. Donald, Finite element techniques to predict collapse loads of foundations. Proc. Int. Conj. Finite Element Meth. (pp. 341-347), Shanghai (1982). [16] O. Orringer and S. E. French, FEA BL (Finite Element Analysis Basic Library) User's Guide. Aeroelastic and Structures Research Lab., Dept of Aeronautics and Astronautics, MIT (1972). [17] E. H. Davis, Theories of plasticity and the failure of soil masses. Soil Mechanics (Edited by I. K. Lee), Butterworths, London (1968). [18] G. C. Nayak and O. C. Zienkiewicz, Elasto-plastic stress analysis. A generalization for various constitutive relations including strain softening. Int. J. numer. Meth. Engng 5, 113-135 (1972), [19] R. Hill, Mathematical Theory of Plasticity. Oxford University Press, London (1950). [20] C. M. Haberfield, The performance of the pressuremeter and socketed piles in weak rock. Ph.D. Thesis, Monash University, Australia (1987). [21] A D. Cox, G. Eason and H. G. Hopkins, Axially symmetric plastic deformations in soils. Royal Society of London, Phil Transactions Series ,4, 1-45 (1961). [22] T. Nishioka, R. Murakami and S. Matsuo, Finite element simulation of fast curving fracture tests. Proc. Int. Conf. Computational Engng Sci. (Edited by S. N. Atluri et al.), pp. 75(~755 (1991). [23] M. S. Shephard, N. A. B. Yehia, G. S. Burd and T. J. Weidner, Automatic crack propagation tracking. Compos. Structures 20, 211-223 (1985). [24] Y. Sumi, Computational crack path prediction. Theor. appl. Fracture Mech. 4, 149-156 (1985). [25] T. Taniguchi, Crack propagation analysis in civil engineering structures. Civil-Comp. 89, 371-380 (1989). [26] T. Taniguchi, Private communication (1991). [27] I. L. Lim, I. W. Johnston, S. K. Choi and V. Murti, An improved numerical inverse isoparametric mapping technique for 2D mesh rezoning. Engng Fracture Mech. 41, 417-435 (1992). [28] I. L. Lim, I. W. Johnston and S. K. Choi, Application of quadratic singular distorted isoparametric elements to linear fracture mechanics. Int. J. numer. Meth. Engng 36, 2473-2499 (1993). [29] I. L. Lim, 1. W. Johnston and S. K. Choi, The use of transition elements. Engng Fracture Mech. 40, 975-983 (1992). [30] V. Murti, Numerical fracture mechanics using finite element methods. Ph.D. Thesis, University of New South Wales (1986). [31] S. Valliappan and V. Murti, Automatic remeshing technique in quasistatic and dynamic crack propagation. Proc. of NUMETA'85 Conf. (Edited by J Middleton and G. N. Pande), pp. 107-116. Swansea (1985).

Finite element code for fracture propagation analysis

211

[32] P. A. Wawrzynek, Interactive finite element analysis of fracture processes: an integrated approach. Masters Dissertation, Cornell University (1987). [33] P. P. Lynn and A. R, Ingraffea, Transition elements to be used with quarter-point crack-tip elements. Int. J. numer. Meth. Engng 12, 1031-1036 (1978). [34] J. Robinson, CRE method of element testing and the Jacobian shape parameters. Engng Comp. 4, 113-118 (1987). [35] J. Robinson, Some new distortion measures for quadrilaterals. Finite Element Anal. Des. 3, 183-197 (1987). [36] I. L. Lim, I. W. Johnston and S. K. Choi, Comparison between various displacement-based stress intensity factor computation techniques. Int. J. Fracture 58, 193-210 (1992). [37] 1. L. Lira, l. W. Johnston and S. K. Choi, On stress intensity computation from the quarter-point element displacements. Comm. appl. Numer. Meth. 8, 291-300 0992). [38] G. C. Sih, Strain energy density factor applied to mixed-mode crack problems. Int. J. Fracture Mech. 10, 305-321 (1974). [39] M. A. Hussain, S. L. Pu and J. H. Underwood, Strain energy release rate for a crack under combined Mode I and Mode II. Fracture Analysis, A S T M S T P 560, 2-28 (1974). [40] F. Erdogan and G. C. Sih, On the crack extension path in plates under plane loading and transverse shear. A S M E J. has. Engng 85, 519-527 (1963). [41] N. M. Taha and S. E. Swartz, Crack propagation models for mixed-mode loading, in Fracture o/' Concrete and Rocks: Recent Developments (Edited by S. P. Shah, S. E. Swartz and B. Barr), pp. 5-17. Elsevier Applied Science (1989). [42] K. J. Chang, A further examination on the application of the strain energy density theory to the angled crack problem. A S M E J. appl. Mech. 49, 377-382 (1982). [43] N. A. B. Yehia and M. S. Shephard, The NT-criterion for predicting crack growth increments. Engng Fracture Mech. 26, 371-382 (1987). [44] V. E. Saouma. Finite element analysis of reinforced concrete: a fracture mechanics approach. Ph.D. Thesis, Cornell University (1981). [45] S. W. Sloan. A Fortran program for profile and wavefront reduction. Int. J. numer. Meth. Engng 28, 2651-2679 (1989). [46] I. L. Lim, I. W. Johnston and S. K. Choi, A comparison between various algorithms for profile reduction of sparse matrices. Comput. Structures 5'7, 297-302 (1995). [47] I. L. Lim, 1. W. Johnston and S. K. Choi, A mesh rezoning technique for discrete crack propagation analyses. 6th Int. Con[~ Aust. on Finite Element Meth., pp. 429~435 (1991). [48] E. Hinton and J. S. Campbell, Local and global smoothing of discontinuous finite element functions using a least square method. Int. J. numer. Meth. Engng 8, 461-480 (1974). [49] J. H. Cheng and N. Kikuchi, A mesh rezoning technique for finite element simulation of metal forming processes. Int. J. numer. Meth. Engng 23, 219-228 (1986). [50] R. H. Crawford, D. C. Anderson and W. N. Waggenspack, Mesh rezoning of 2D isoparametric elements by inversion. h~t. J. numer. Meth, Engng 28, 523-531 (1989). [51] C. J, M. Gelten and A. W. A. Konter, Application of mesh rezoning in the updated Langrangian method to metal forming analyses, in Numerical Methods in Industrial Forming Processes (Edited by J. F. T. Pittman et aLL pp. 511 521. Pineridge Press (1982). [52] H. T. Y. Yang, M. Heinstein and J. M. Shih, Adaptive 2D finite element simulation of metal forming process, lnt. J. numer. Meth. Engng 28, 1409 1428 (1989). [53] J. H. Cheng, Automatic adaptive remeshing for finite element simulation of forming processes. Int. J. numer. Meth. Engng 26, 1-18 (1988). [54] A. Habraken and S. Cescotto, An automatic technique for finite element simulation of forming processes. Int. J. numer. Meth. Engng 30, 1503 1525 (1990). [55] V. Murti, Y. Wang and S. Valliappan, Numerical inverse isoparametric mapping in 3D FEM. Compos. Structures 29, 611-622 (1988). [56] I. L. Lira, I. W. Johnston and S. K. Choi, A numerical model for simulation of fracture propagation in rock. Int. S.imp. on the Application ~/' Comp. Meth. in Rock Mech. and Engng (pp. 397-404), Xian, China (1993). [57] I. L. Lim, I. W. Johnston and S. K. Choi, Assessment of mixed-mode fracture toughness testing methods for rock. Int. J. Rock Mech. Min. Sei. Geomech. Abstr. 31, 265-272 (1994). [58] 1. L. Lira, I. W. Johnston, S. K. Choi and J. N. Boland, Fracture testing of a soft rock with semi-circular specimens under three-point bending. Part 2--mixed-mode. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 31, 199 212 (1994). (Received 26 August 1994)