Non‐planar 3D crack growth by the extended finite element and level sets—Part II: Level set update

We present a level set method for treating the growth of non‐planar three‐dimensional cracks.The crack is defined by two almost‐orthogonal level sets (signed distance functions). One of them describes the crack as a two‐dimensional surface in a three‐dimensional space, and the second is used to describe the one‐dimensional crack front, which is the intersection of the two level sets. A Hamilton–Jacobi equation is used to update the level sets. A velocity extension is developed that preserves the old crack surface and can accurately generate the growing surface. The technique is coupled with the extended finite element method which approximates the displacement field with a discontinuous partition of unity. This displacement field is constructed directly in terms of the level sets, so the discretization by finite elements requires no explicit representation of the crack surface. Numerical experiments show the robustness of the method, both in accuracy and in treating cracks with significant changes in topology. Copyright © 2002 John Wiley & Sons, Ltd.


INTRODUCTION
Level set methods have been applied eectively to many problems of moving fronts (see References [1; 2]). Methods for curves in three dimensions are given by Burchard et al. [3] and Cheng et al. [4]. However, these methods cannot be applied directly to the growth of nonplanar three-dimensional cracks. In this paper, an extension of the method of Peng et al. [5] applicable to three-dimensional crack growth is described. The method couples very naturally 1 with the extended nite element method [6][7][8][9], and makes possible the modeling of arbitrary crack growth in three-dimensional bodies without remeshing (see the companion paper [10]).
Level set methods for crack growth in two-dimensional models have been described by Stolarska et al. [11]. Sukumar et al. [12] developed a fast marching method for the growth of planar cracks in three dimensions. Because the crack was planar, they were able to represent the crack front by a single level set, which was updated by the fast marching method. To obtain accurate solutions, they used a separate mesh for the level set update [13; 14].
Our level set update procedures are based on Peng et al. [5] and Burchard et al. [3] but are modied to handle cracks. We describe the crack by two level sets: the rst describes the crack surface, whereas the second is constructed so that the intersection of two level sets gives the crack front. Signed distance functions are used for the level sets. The level sets are updated by solving hyperbolic partial dierential equations. These need only be solved in a small domain around the crack front.
The method is applicable to arbitrary three-dimensional solids. The crack surface need not bear any relationship to the mesh, and in the examples, tetrahedral meshes constructed by automatic mesh generators are used. Therefore, the method can easily be implemented in standard nite elements programs.
An outline of the paper is as follows. In Section 2, we describe our level set representation of the crack and the partial dierential equations which govern their evolution and hence the growth of the crack. The discretization of the partial dierential equations is given in Section 3. Section 4 describes examples we have devised to check the evolution of the level sets and applications of the method to three-dimensional fatigue crack growth problems are given.

Description of crack by level sets
Consider a crack within a three-dimensional body as shown in Figure 1. The crack may either be completely interior to the body or intersect an edge. Two level sets are used to represent the crack: 1. The level set, called the crack surface level set; its zero isosurface corresponds to the crack surface.   2. The level set, called the front level set; the intersection of the crack surface zero level set ( = 0) with the front zero level set ( = 0) gives the crack front.
Throughout this paper we will assume that (x;t) and (x;t) are signed distance functions, i.e. the distance from x to (x;t) = 0 is given by , and similarly for . The crack front and the crack surface are given by the following: does not intersect the crack Both functions are assumed to be Lipschitz continuous. They provide a convenient means of constructing the discontinuities in the displacement eld in Reference [10].
The functions (x;t) and (x;t) are assumed to be orthogonal so that Since the two level sets are orthogonal, they will never have the same tangent plane on the crack front. Indeed, if the two level sets have the same tangent plane, then the crack front is not dened. The orthogonality of the two level sets also provides a simple way to dene an orthogonal system of curvilinear co-ordinates intrinsic to the crack. The numerical approximation of and will not be exactly orthogonal, but will satisfy orthogonality closely enough to meet our needs. The level set is not explicitly dened for ¿0. However, for purposes of updating the level sets, the level set must be constructed in a small domain ahead of the crack front, i.e. for ¿0. An extension of the level set is shown in Figure 2. It is not unique but a procedure is described later that generates a useful extension. We will show more complicated examples of extensions and how they are generated later. Extensions of other variables, such as the crack front velocity, will also be needed.
The level sets are only updated in a small subdomain around the crack, which we call the level set subdomain. An example of a level set subdomain is shown in Figure 3. While only a tube around the crack front is needed for the level set update, we need the level sets in the entire subdomain that encloses the crack for the stress analysis in Reference [10].  The unit vector normal to the crack front is given by The unit vector normal to the crack surface is denoted by n and given by The tangent to the crack front e t is then given by Note that these base vectors (see Figure 4) are dened everywhere in the level set domain. For the purpose of updating the level sets, we will consider the propagation of the crack front by a velocity in the plane dened by n and n . This velocity depends on the crack growth law. An example of a crack growth law is given later. We denote this velocity by V. This velocity is decomposed into components (see Figure 5) along the n and n directions by

Initialization and reinitialization
The initial level sets are given as data. For many cracks, level set functions which are not signed distance functions can easily be constructed. For example, for an elliptical planar crack normal to the z-axis with a and b the semi-major and semi-minor axes, respectively, the level The level set in this case is not a signed distance function. In order to make a signed distance function, the condition |∇ | = 1 (11) needs to be imposed. We use the method proposed in Reference [5] of normalizing ∇ by solving the Hamilton-Jacobi equation: where sign( ) is the sign function and is a time like parameter; it does not correspond to the physical time of the problem but at any time is reset to zero and the equation is then updated in until steady state is reached. When Equation (12) is brought to steady state, Equation (11) is satised. Theoretical and numerical studies show that this method is ecient, does not change the zero isosurface, only needs three or four iterations in a narrow band and has a complexity of O(N ) where N is the number of unknowns [5]. The level set is similarly initialized when it is not a signed distance function. This process is called initialization. When data are available only on the zero level sets, this process is called extension. In updates, the level sets often deviate from signed distance functions, so the same process is applied at every level set update; this is then called reinitialization which was rst introduced by Chopp [15] to improve the stability of the level set method.
Before discussing the specic procedures for a crack, we briey review some basic concepts in level set theory. The level sets are updated by solving Hamilton-Jacobi type equations. For any level set , this equation can be obtained as follows. Let x(s; t) be the trajectory of a point on the crack front. By denition (x(s; t);t) = 0, so the derivative with respect to time on any point of the crack front must vanish, which gives @ @t Generally, the velocity is projected onto the normal to the zero isosurface of the level sets, so Equation (13) becomes @ @t

Extension of velocity eld
A key ingredient in applying level set methods to crack growth is the extension of the velocity eld from the crack front to the entire level set subdomain. Recall that the velocity of the crack front is only given along the crack front, i.e. the curve dened by the intersection of = 0 and = 0. In order to update the level sets, the velocity must be known in the entire level set subdomain. Furthermore, the velocity eld must be updated so that the level set is not changed for ¡0, i.e. the crack surface must not be moved by the level set update. The extension process for the velocity eld and update of the level sets that we have developed consists of three steps: 1. Preliminary extension of the 1D velocity eld to the level set subdomain, 2. A modication of the -component of the velocity, V , 3. The update of the and level sets by standard methods using the modied velocity.
To extend the 1D scalar velocity eld V from the crack front to the 3D subdomain, we rst dene this eld at the nodes of the nite elements cut by the crack front. Subsequently, an interpolation process [16; 17] is used. Then we use the extension process, where the velocity is extended from the nodes into the level set subdomain, so that it is constant along the normal to the and level sets, i.e.
We follow Peng et al. [5] by solving to steady state two Hamilton-Jacobi equations: which allows us to satisfy the conditions (15). By solving to steady state these Hamilton-Jacobi equations, we extend the 1D velocity eld V to the entire 3D level set subdomain. A simple geometrical interpretation can be given to the above: each component of the eld V is extended along orthogonal lines from the interface. Indeed, the steady-state solutions of Equations (16a) and (16b) satisfy where n and n are the normal vectors to the level sets and , respectively. This process is used along the gradients of the two level sets and for each scalar velocity eld V and V .
Once the velocity has been extended into the level set subdomain, the component V is modied as follows: where H is the Heaviside step function and a small parameter. This modication serves the following purposes: 1. It makes V = 0 for ¡0, so on the crack surface, the velocity of the component normal to the crack vanishes, 2. It linearly increases V from the crack front to the surface where the crack front is expected to be at the end of the update.
The modied velocity is then used to update the level sets. If V ¡ V , the present method is not applicable since then the major non-zero component is V , which indicates the development of a kink, i.e. a discontinuity in the crack surface. In this paper, we restrict our attention to smooth cracks, so this situation should not arise.

Level set update
The entire procedure is summarized in Table I. Steps 1-4 have been discussed in the above. The update of the level set is given in step 5, and is done by the customary level set Table I. Scheme for level set update.
1-orthogonal extension of the crack level set in level set subdomain dened by ¿0 (Equation (20)) 2-extend V to the domain 4-adjustment to prevent modication of previous crack surface V = H ( ) V V t 5-update and reinitialize the level set 7-orthogonalize and reinitialize the level set update Equation (14). Immediately after its update, it is normalized by the counterpart of Equation (12), i.e. is reinitialized.
Steps 6 and 7 constitute the update of the level set. Prior to reinitializing ,i ti s orthogonalized to . For this, we need to satisfy ∇ ·∇ = 0 (19) which is imposed by procedures similar to those used for the extension. As we have explained before, Equation (16a) extends a variable from an interface along orthogonal lines. In the same way, this process is able to reorthogonalize one level set to another one. Orthogonalization (often called reorthogonalization because the level sets were orthogonal in the previous time step) is performed in step 7 by solving to steady state. This equation is equivalent to Equation (16a). This process has been used for the motion of curves in three dimensions in Reference [3]. We make the following remarks.
1. The crack must not be changed by the update, so the new level set must continue to have a zero velocity eld on the previous crack surface (i.e. for ¡0) and has to exactly dene the growing crack. This is the goal of the step 4 in Table I. 2. The extension of each scalar velocity eld in step 3 of Table I is done in two steps: rst we extend the value of the 1D scalar eld from the crack front to the nodes of the elements cut by the crack front [17]; then we extend these values in the normal direction to . Finally, we extend these values in the normal direction to . A very simple way to improve this method is to extend in the same steady-state process in the directions normal to and .

DISCRETIZATION IN TIME AND SPACE FOR UNSTRUCTURED MESHES
In this part, we will describe the discretization in time and space. All the equations discussed so far (Equations (12)- (20) and Table I) are Hamilton-Jacobi equations of the form @f @t where H is the Hamiltonian. The level sets are approximated by nite elements so in the level set subdomain: where N I (x) are the shape functions (interpolants). The shape functions N I (x) are C 0 ,s ot h e crack representation is piecewise continuously dierentiable, i.e. there are kinks in the surface at element interfaces. The velocities are approximated by the same shape function: When the level set method is used in conjunction with X-FEM, we have used the same mesh and shape functions for the stress analysis. The level set method has previously been applied to unstructured meshes by Barth and Sethian [18], and we used the same procedures. They proposed at rst a nite element approach which veries a monotonicity condition. However, this condition is too strong because it does not preserve the Lipshitz continuity of the numerical Hamiltonian. One of the main consequences of this is to reduce the accuracy of the numerical scheme (in practice, bad accuracy arises from obtuse tetrahedrons). To gain Lipshitz continuity of the numerical Hamiltonian, they proposed that the monotonicity condition be relaxed to a positivity condition. This provides a robust scheme which is stable, accurate and convergent (see Appendix A for an example of the algorithm in the case of Equation (14)).
In the Hamilton-Jacobi equations such as (12), a discrete sign function sign(f) is needed. The following expression is used [5]: where ∇f i is the gradient of the level set on the considered vertex, and the coecient l is a small parameter which is useful close to the iso-zero level set. The sign function is smooth and permits the zero level set to be preserved during the reinitialization process. The value of l is dened as a characteristic length of the smallest nite element of the level set subdomain. In the same way, for the Hamilton-Jacobi equations (16) and (20) the following denition of the sign function is used (see Reference [5]): The following second-order Runge-Kutta scheme is used for time integration: The level set update is explicit and conditionally stable, so the time step in this scheme has to meet the CFL condition. The critical time step can be estimated by where h is the smallest element size. Since the time step for the solid mechanics model (which is static and not restricted by any stability requirements) is generally several times the above critical time step, several updates of the level sets are required for each update of the solid mechanics model.

NUMERICAL EXPERIMENTS
In this section, numerical experiments are performed for an uncoupled level set update and three crack growth problems. Most of the meshes and gures were made with gmsh [19].
In the examples, the material properties are elastic and isotropic with Young's modulus E = 280 GPa and Poisson's ratio =0:3. Given the two level sets and , the crack front location is extracted as a set of 1D segments. The stress intensity factors at the mid-point of each of these segments is evaluated by the domain integral described in Reference [10]. We consider fatigue crack growth governed by the Paris law, which gives the rate crack of growth in mode I in terms of load cycles N by where C is a constant t to experimental results and G is the maximum energy release rate. We consider the cycles N a time-like variable so the expression for the crack front velocity is where c , the angle of the velocity to the plane tangent to the crack, is obtained by Thus, the crack growth direction depends on modes I and II stress intensity factors, whereas the crack speed depends on all three through the energy release rate G. In the numerical studies, for simplicity we chose m = 1 and C =1.

Level set experiments
We begin with a simple example to demonstrate the robustness of the discretization scheme for the level sets. Consider the geometry shown in Figure 6: a cubic box with an edge crack of length a. The initial level sets are two orthogonal planes. The -level set is horizontal and represents the crack, and the -level set is vertical and represents the crack front ( Figure 6).
The goal of this study is to update the crack front following a circle. The velocity eld is dened explicitly in each time step without the use of any mechanical law. In each time where is the angle between the planes tangent to the current and next crack surfaces. The above was chosen so that an exact solution gives a cylindrical surface of radius R. The radius was taken to be unity. The error was measured by error = 1 Nr exact N I |r I () − r exact | (32) Figure 7 shows the propagation of the level set and one instance of the level set. The front at various times is also shown. In Table II, we compare the numerical radius r of the = 0 isoset with the exact value for various angles. In all cases, the error is less than 1 per cent. This example shows that the algorithm is able to describe the evolution of non-planar 3D surfaces with good accuracy.

A beam under bending
In this example, we consider the beam shown in Figure 8 subjected to a bending load (see Reference [20]) for a theoretical study of this problem). The beam dimensions are h =0:02 m, l =0:1 m and d =0:01 m: An initial edge crack along a plane at an angle to the plane of the cross-section is used to initiate crack growth with a =0:01 m and =45 • (see Figure 8). The initial geometry is discretized with a unstructured mesh of 2553 nodes and 13 393 tetrahedrons. It is emphasized that the mesh does not conform to the crack geometry, and that the same mesh is used throughout the simulation. The crack is driven by a Paris fatigue law with the maximum circumferential stress hypothesis for the direction of crack propagation. From 15 to 25 points were used on the crack front to compute the velocity and the stress intensity factors. In this study, the beam is completely cut into two after 17 time steps. Figure 9 shows the evolution of the crack front from a top view. We can observe that the crack grows asymptotically to a plane orthogonal to the axis of the beam. Figure 10 shows the crack after 12 time steps, and the vector velocity eld on the crack front. An important remark concerns the topological properties of the method: with the level set approach, the crack does not need any specic treatment on the boundary, and the crack can completely cut the beam. Furthermore, as we will show in the next examples, this representation is very versatile in describing the topological evolution of the crack front: for instance a front cut into two or more subfronts. Figure 11 shows the evolution of a semi-circular crack. Again, we observe that the crack grows asymptotically to a plane orthogonal to the axis of the beam. Here, the crack front rst grows to the lateral free surfaces, then normal to the surface. No changes are needed in the algorithm to account for these topological changes.

A lens-shaped crack
In this example, we consider a cube with a cusp crack subjected to hydrostatic tension as shown in Figure 12.  The crack geometry is characterized by the radius R and the azimuthal angle . This problem was studied in the companion paper [10]. A simulation of the evolution of the crack in a cube with h =0:01 m and an initial crack dened by R =0:005 m and =45 • was considered. The initial geometry is discretized with a unstructured mesh of 1767 nodes and 8895 tetrahedrons. As in the last example, the mesh does not conform to the crack geometry, and the same mesh is used throughout the simulation. Figure 13 represents the initial crack with the initial vector velocity eld on the front. Figure 12. Initial cusp crack subjected to hydrostatic tension. From 55 to 67 points were used on the front to compute the crack front velocity and the stress intensity factors. Figure 14 shows the crack after 15 time steps. Note that the initial front has involved into four subfronts in each corner. Furthermore, in this example, one can notice that the convexity of the front has changed: this may be due to faster growth of the crack near the boundary of the cube.

Inclined penny-shaped crack
In this last example, we consider the problem in Figure 15: a cube with an inclined pennyshaped crack subjected to a tensile loading with h =0:02 m and an initial crack dened by a = b =0:005 m and =45 • . The mesh consists of 1747 nodes and 8847 tetrahedrons. Again, the mesh does not conform to the crack geometry, and the same mesh is used throughout the simulation. Figure 16 shows the evolution of the crack after 17 time steps. From 46 to 58 points were used on the front to compute the stress intensity factors. At the end of the computation, the box is completely cut by the crack. Furthermore, one can notice that the crack front has a complex path. Like the previous problem, at the beginning the front is one entity, then four entities, until the structure is completely cut by the crack.

CONCLUSIONS
A level set method for arbitrary non-planar cracks in three-dimensional bodies has been presented. A key feature that distinguishes this method from other level set methods is a technique for the extension of the velocity eld and a reinitialization process that preserves the shape of the crack but allows arbitrary growth of the crack front. The update of the level sets is performed through the solution of Hamilton-Jacobi equations in a small subdomain surrounding the crack. This is an ecient technique which generally requires only about 10 per cent of the total computation time.
These are our initial studies into this coupled fracture mechanics-level set method, and several weaknesses remain in the method. For very small cracks, it may be necessary to adaptively rene the volume around the crack. This is particularly important for the discretization of the level set update, which requires a ner mesh to accurately capture the crack front. Methods with dierent renements in the solid mechanics and level set discretization may also be desirable.
The level set technique couples naturally with the extended nite element method, wherein the discontinuous and near-tip asymptotic elds are constructed through a partition of unity. The resulting combined method requires no explicit representation of the crack except in its visualization. Instead, the crack and its growth are described entirely in forms of nodal data. This simplies the structure of the software and leads to great versatility in treating complex problems in crack growth. The extended and generalized nite element methods have recently seen rapid development [21][22][23][24], so the combination with level sets is very promising.
Although the formulation was applied here to elastostatic problems, and in particular, fatigue crack growth, the method is applicable to many other types of fracture problems. For example, the level set update is applicable to non-linear fracture and dynamic fracture. Only the constitutive equations and the solver in the solid mechanics software would need to be altered for these applications: the level set update and the displacement approximation would be identical.
APPENDIX A: THE EXPLICIT POSITIVE COEFFICIENT SCHEME FOR EQUATION (14) Step 1: Initialize * I = w I = 0 for each node of the mesh. where H is the Heaviside function. * I = * I + I w I = w I + I v (v is the volume of the nite element E) Step 3: Loop on the nodes of the mesh for the time integration.