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

Anthony Gravouil, Nicolas Moës, Ted Belytschko

► To cite this version:

Anthony Gravouil, Nicolas Moës, Ted Belytschko. Non-planar 3D crack growth by the extended finite element and level sets. Part II: Level set update. International Journal for Numerical Methods in Engineering, 2002, 53 (11), pp.2569-2586. 10.1002/nme.430. hal-01007111

HAL Id: hal-01007111

https://hal.science/hal-01007111

Submitted on 24 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Non-planar 3D crack growth by the extended finite element and level sets—Part II: Level set update

A. Gravouil†, N. Moës‡ and T. Belytschko*,§

Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, U.S.A.

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.

KEY WORDS: fracture; cracks; discontinuous approximation; finite elements; level sets

## 1. INTRODUCTION

Level set methods have been applied effectively 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 non-planar 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

<footnote>*Correspondence to: Ted Belytschko, Department of Mechanical Engineering, Northwestern University, 2145 N. Sheridan, Evanston, IL 60208, U.S.A.</footnote>

<footnote>†Laboratoire de Mécanique des Solides, Institut National des Sciences Appliquées de Lyon, 34 Av. des Arts, 69621 Villeurbanne, France</footnote>

<footnote>‡Laboratoire de Mécanique et Matériaux, Ecole Centrale de Nantes, 1 Rue de la Noe, 44321 Nantes, France</footnote>

<footnote>§E-mail: tedbelytschko@northwestern.edu</footnote>

<footnote>Contract/grant sponsor: Office of Naval Research</footnote>

<footnote>Contract/grant sponsor: Department of Energy</footnote>

<footnote>Contract/grant sponsor: Délégation générale pour l'Armement</footnote>
with the extended finite element method [6–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 modified to handle cracks. We describe the crack by two level sets: the first 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 differential 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 finite 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 differential equations which govern their evolution and hence the growth of the crack. The discretization of the partial differential 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.

## 2. LEVEL SETS: DESCRIPTION AND UPDATE

### 2.1. *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 $\phi$ level set, called the crack surface level set; its zero isosurface corresponds to the crack surface.

<figure><img src="image_2.png" /><figcaption>Figure 1. Representation of a non-planar crack in three-dimensions by two level sets $\phi$ and $\psi$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2. Level sets for a straight edge crack and an extension of the $\phi$ level set.</figcaption></figure>

2. The $\psi$ level set, called the front level set; the intersection of the crack surface zero level set ($\phi=0$) with the front zero level set ($\psi=0$) gives the crack front.

Throughout this paper we will assume that $\phi(\mathbf{x}, t)$ and $\psi(\mathbf{x}, t)$ are signed distance functions, i.e. the distance from $\mathbf{x}$ to $\phi(\mathbf{x}, t)=0$ is given by $\phi$, and similarly for $\psi$. The crack front and the crack surface are given by the following:

$$ \phi(\mathbf{x}, t) = 0, \quad (\mathbf{x}, t) < 0 \quad \text{defines the crack location} \quad (1) $$

$$ \phi(\mathbf{x}, t) = 0, \quad (\mathbf{x}, t) = 0 \quad \text{gives the crack front} \quad (2) $$

$$ (\mathbf{x}, t) > 0 \quad \text{does not intersect the crack} \quad (3) $$

Both functions are assumed to be Lipschitz continuous. They provide a convenient means of constructing the discontinuities in the displacement field in Reference [10].

The functions $\phi(\mathbf{x}, t)$ and $\psi(\mathbf{x}, t)$ are assumed to be orthogonal so that

$$ \nabla\phi \cdot \nabla\psi = 0 \quad \forall t \quad (4) $$

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 defined. The orthogonality of the two level sets also provides a simple way to define an orthogonal system of curvilinear co-ordinates intrinsic to the crack. The numerical approximation of $\phi$ and $\psi$ will not be exactly orthogonal, but will satisfy orthogonality closely enough to meet our needs.

The $\phi$ level set is not explicitly defined for $\psi>0$. However, for purposes of updating the level sets, the $\phi$ level set must be constructed in a small domain ahead of the crack front, i.e. for $\psi>0$. An extension of the $\phi$ 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].
<figure><img src="image_4.png" /><figcaption>Figure 3. Level set subdomain for a crack; level sets are updated only in this subdomain.</figcaption></figure>

The unit vector normal to the crack front is given by

$$ \mathbf{n}_{\psi} = \frac{\nabla \psi}{\|\nabla \psi\|} \qquad (5) $$

The unit vector normal to the crack surface is denoted by $\mathbf{n}_{\phi}$ and given by

$$ \mathbf{n}_{\phi} = \frac{\nabla \phi}{\|\nabla \phi\|} \qquad (6) $$

The tangent to the crack front $\mathbf{e}_t$ is then given by

$$ \mathbf{e}_t = \mathbf{n}_{\psi} \times \mathbf{n}_{\phi} \qquad (7) $$

Note that these base vectors (see Figure 4) are defined 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 defined by $\mathbf{n}_{\phi}$ and $\mathbf{n}_{\psi}$. This velocity depends on the crack growth law. An example of a crack growth law is given later. We denote this velocity by $\mathbf{V}$. This velocity is decomposed into components (see Figure 5) along the $\mathbf{n}_{\psi}$ and $\mathbf{n}_{\phi}$ directions by

$$ \mathbf{V} = V_{\psi}\mathbf{n}_{\psi} + V_{\phi}\mathbf{n}_{\phi} \qquad (8) $$

## 2.2. 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
<figure><img src="image_5.png" /><figcaption>Figure 4. Base vectors for the level set functions for a crack.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. The components of the crack front velocity $V_\phi$ and $V_\psi$.</figcaption></figure>

set functions are:

$$ \phi(\mathbf{x}, 0) = z \qquad (9) $$

$$ (\mathbf{x}, 0) = \frac{x^2}{a^2} + \frac{y^2}{b^2} - 1^2 \qquad (10) $$

The $\psi$ level set in this case is not a signed distance function. In order to make $\psi$ a signed distance function, the condition

$$ |\nabla \psi| = 1 \qquad (11) $$

needs to be imposed. We use the method proposed in Reference [5] of normalizing $\nabla\psi$ by solving the Hamilton–Jacobi equation:

$$ \frac{\partial \psi}{\partial \tau} + \operatorname{sign}(\psi)(|\nabla \psi| - 1) = 0 \qquad (12) $$

where $\operatorname{sign}(\psi)$ is the sign function and $\tau$ is a time like parameter; it does not correspond to the physical time of the problem but at any time $\tau$ is reset to zero and the equation is then updated in $\tau$ until steady state is reached. When Equation (12) is brought to steady state, Equation (11) is satisfied. Theoretical and numerical studies show that this method is efficient, 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 $\phi$ 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 first introduced by Chopp [15] to improve the stability of the level set method.

Before discussing the specific procedures for a crack, we briefly review some basic concepts in level set theory. The level sets are updated by solving Hamilton–Jacobi type equations. For any level set $\psi$, this equation can be obtained as follows. Let $\mathbf{x}(s, t)$ be the trajectory of a point on the crack front. By definition $\psi(\mathbf{x}(s, t), t) = 0$, so the derivative with respect to time
on any point of the crack front must vanish, which gives

$$ \frac{\partial \psi}{\partial t} + \mathbf{V} \cdot \nabla \psi = 0 \qquad (13) $$

Generally, the velocity is projected onto the normal to the zero isosurface of the level sets, so Equation (13) becomes

$$ \frac{\partial \psi}{\partial t} + V_{\psi} |\nabla \psi| = 0 \qquad (14) $$

### 2.3. Extension of velocity field

A key ingredient in applying level set methods to crack growth is the extension of the velocity field 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 defined by the intersection of $\phi=0$ and $\psi=0$. In order to update the level sets, the velocity must be known in the entire level set subdomain. Furthermore, the velocity field must be updated so that the $\phi$ level set is not changed for $\psi<0$, i.e. the crack surface must not be moved by the level set update. The extension process for the velocity field and update of the level sets that we have developed consists of three steps:

1. Preliminary extension of the 1D velocity field to the level set subdomain,

2. A modification of the $\phi$-component of the velocity, $V_{\phi}$,

3. The update of the $\phi$ and $\psi$ level sets by standard methods using the modified velocity.

To extend the 1D scalar velocity field $V_{\phi}$ from the crack front to the 3D subdomain, we first define this field at the nodes of the finite 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 $\phi$ and $\psi$ level sets, i.e.

$$ \nabla \phi \cdot \nabla V_{\phi} = 0 \quad \text{and} \quad \nabla \psi \cdot \nabla V_{\phi} = 0 \qquad (15) $$

We follow Peng et al. [5] by solving to steady state two Hamilton-Jacobi equations:

$$ \frac{\partial V_{\phi}}{\partial \tau} + \operatorname{sign}(\phi) \frac{\nabla \phi}{|\nabla \phi|} \cdot \nabla V_{\phi} = 0 \qquad (16a) $$

$$ \frac{\partial V_{\phi}}{\partial \tau} + \operatorname{sign}(\psi) \frac{\nabla \psi}{|\nabla \psi|} \cdot \nabla V_{\phi} = 0 \qquad (16b) $$

which allows us to satisfy the conditions (15). By solving to steady state these Hamilton-Jacobi equations, we extend the 1D velocity field $V_{\phi}$ to the entire 3D level set subdomain. A simple geometrical interpretation can be given to the above: each component of the field $V_{\phi}$ is extended along orthogonal lines from the interface. Indeed, the steady-state solutions of Equations (16a) and (16b) satisfy

$$ \mathbf{n}_{\phi} \cdot \nabla V_{\phi} = 0 \quad \text{and} \quad \mathbf{n}_{\psi} \cdot \nabla V_{\phi} = 0 \qquad (17) $$

where $\mathbf{n}_{\phi}$ and $\mathbf{n}_{\psi}$ are the normal vectors to the level sets $\phi$ and $\psi$, respectively. This process is used along the gradients of the two level sets $\phi$ and $\psi$ for each scalar velocity field $V_{\phi}$ and $V_{\psi}$.
Once the velocity has been extended into the level set subdomain, the component $V_\phi$ is
modified as follows:

$$
\bar{V}_{\phi} = H(\psi) \frac{V_{\phi} \psi}{V_{\psi} \Delta t} \quad \text{if } V_{\psi} > \varepsilon \|\mathbf{V}\| \quad (18)
$$

where $H$ is the Heaviside step function and $\varepsilon$ a small parameter. This modification serves the
following purposes:

1. It makes $\bar{V}_{\phi} = 0$ for $\psi < 0$, so on the crack surface, the velocity of the component normal to the crack vanishes,

2. It linearly increases $\bar{V}_{\phi}$ from the crack front to the $\psi$ surface where the crack front is expected to be at the end of the update.

The modified velocity is then used to update the level sets. If $V_\psi < \varepsilon \|\mathbf{V}\|$, the present method is not applicable since then the major non-zero component is $V_\phi$, 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.

2.4. Level set update

The entire procedure is summarized in Table I. Steps 1–4 have been discussed in the above.
The update of the $\phi$ level set is given in step 5, and is done by the customary level set

<caption>Table I. Scheme for level set update.</caption>

1—orthogonal extension of the crack level set $\phi$ in level set subdomain defined by $\psi \ge 0$ (Equation (20))

2—extend $V_\psi$ to the domain

$$
\frac{\partial V_{\psi}}{\partial \tau} + \mathrm{sign}(\psi) \frac{\nabla}{|\nabla \psi|} \cdot \nabla V_{\psi} = 0 \quad \frac{\partial V_{\psi}}{\partial \tau} + \mathrm{sign}(\phi) \frac{\nabla \phi}{|\nabla \phi|} \cdot \nabla V_{\psi} = 0
$$

3—extend $V_\phi$ to the domain

$$
\frac{\partial V_{\phi}}{\partial \tau} + \mathrm{sign}(\phi) \frac{\nabla \phi}{|\nabla \phi|} \cdot \nabla V_{\phi} = 0 \quad \frac{\partial V_{\phi}}{\partial \tau} + \mathrm{sign}(\psi) \frac{\nabla}{|\nabla \psi|} \cdot \nabla V_{\phi} = 0
$$

4—adjustment to prevent modification of previous crack surface

$$
\bar{V}_{\phi} = H(\psi) \frac{V_{\phi}}{\bar{V}_{\psi} \Delta t}
$$

5—update and reinitialize the $\phi$ level set

$$
\frac{\partial \phi}{\partial t} + \bar{V}_{\phi} |\nabla \phi| = 0 \quad \frac{\partial \phi}{\partial \tau} + \mathrm{sign}(\phi)(|\nabla \phi| - 1) = 0
$$

6—update the $\psi$ level set

$$
\frac{\partial \psi}{\partial t} + V_{\psi} |\nabla \psi| = 0
$$

7—orthogonalize and reinitialize the $\psi$ level set

$$
\frac{\partial \psi}{\partial \tau} + \mathrm{sign}(\phi) \frac{\nabla \phi}{|\nabla \phi|} \cdot \nabla = 0 \quad \frac{\partial \psi}{\partial \tau} + \mathrm{sign}(\psi)(|\nabla \psi| - 1) = 0
$$
update Equation (14). Immediately after its update, it is normalized by the counterpart of Equation (12), i.e. $\phi$ is reinitialized.

Steps 6 and 7 constitute the update of the $\psi$ level set. Prior to reinitializing $\psi$, it is orthogonalized to $\phi$. For this, we need to satisfy

$$ \nabla\phi \cdot \nabla\psi = 0 \qquad (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

$$ \frac{\partial\psi}{\partial\tau} + \text{sign}(\phi) \frac{\nabla\phi}{|\nabla\phi|} \cdot \nabla = 0 \qquad (20) $$

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 $\phi$ level set must continue to have a zero velocity field on the previous crack surface (i.e. for $\psi<0$) and has to exactly define the growing crack. This is the goal of the step 4 in Table I.

2. The extension of each scalar velocity field in step 3 of Table I is done in two steps: first we extend the value of the 1D scalar field 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 $\psi$. Finally, we extend these values in the normal direction to $\phi$. A very simple way to improve this method is to extend in the same steady-state process in the directions normal to $\psi$ and $\phi$.

### 3. 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

$$ \frac{\partial f}{\partial t} + H(\nabla f, \mathbf{x}, t) = 0 \qquad (21) \\ f(\mathbf{x}, 0) = f_0(\mathbf{x}) $$

where $H$ is the Hamiltonian. The level sets are approximated by finite elements so in the level set subdomain:

$$ \phi(\mathbf{x}, t) = \sum_{I} N_{I}(\mathbf{x})\phi_{I}(t) \qquad (22) \\ (\mathbf{x}, t) = \sum_{I} N_{I}(\mathbf{x})\psi_{I}(t) $$
where $N_I(\mathbf{x})$ are the shape functions (interpolants). The shape functions $N_I(\mathbf{x})$ are $C^0$, so the crack representation is piecewise continuously differentiable, i.e. there are kinks in the surface at element interfaces.

The velocities are approximated by the same shape function:

$$V_{\phi}(\mathbf{x}, t) = \sum_{I} N_{I}(\mathbf{x}) V_{\phi I}(t) \qquad (23)$$

$$V_{\psi}(\mathbf{x}, t) = \sum_{I} N_{I}(\mathbf{x}) V_{\psi I}(t)$$

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 first a finite element approach which verifies 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 $\text{sign}(f)$ is needed. The following expression is used [5]:

$$\text{sign}(f_i) = \frac{f_i}{\sqrt{f_i^2 + \Delta l^2 \nabla f_i^2}} \qquad (24)$$

where $\nabla f_i$ is the gradient of the level set on the considered vertex, and the coefficient $\Delta 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 $\Delta l$ is defined as a characteristic length of the smallest finite element of the level set subdomain. In the same way, for the Hamilton-Jacobi equations (16) and (20) the following definition of the sign function is used (see Reference [5]):

$$\text{sign}(f_i) = \frac{f_i}{\sqrt{f_i^2 + \Delta l^2}} \qquad (25)$$

The following second-order Runge-Kutta scheme is used for time integration:

$$\tilde{f}^{n+1} = f^n - \Delta t H(f^n) \qquad (26)$$

$$f^{n+1} = \frac{(f^n + \tilde{f}^{n+1})}{2} - \frac{\Delta t}{2} H(\tilde{f}^{n+1})$$

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

$$\Delta t = \min \frac{h(\mathbf{x})}{|V(\mathbf{x})|} \qquad (27)$$
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.

## 4. 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 figures 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 $\nu = 0.3$.

Given the two level sets $\phi$ and $\psi$, 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

$$ \frac{da}{dN} = CG^m \qquad (28) $$

where $C$ is a constant fit 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

$$ \mathbf{V} = CG^m(\cos \theta_c \mathbf{n}_\psi + \sin \theta_c \mathbf{n}_\phi) \qquad (29) $$

where $\theta_c$, the angle of the velocity to the plane tangent to the crack, is obtained by

$$ \theta_c = 2 \arctan \frac{1}{4} \left( \frac{K_1}{K_{II}} - \operatorname{sign}(K_{II}) \sqrt{\left(\frac{K_1}{K_{II}}\right)^2 + 8} \right) \qquad (30) $$

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$.

### 4.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 $\phi$-level set is horizontal and represents the crack, and the $\psi$-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 field is defined explicitly in each time step without the use of any mechanical law. In each time
<figure><img src="image_7.png" /><figcaption>Figure 6. The initial edge crack with the $\phi$ and $\psi$ level sets.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Evolution of the $\phi$ level set for circular growth; the initial crack is at the left hand bottom.</figcaption></figure>

<table><caption>Table II. Radial error compared to the exact circle.</caption><tr><td>Angle</td><td>22.5°</td><td>45°</td><td>67.5°</td><td>90°</td><td>112.5°</td><td>135°</td><td>157.5°</td><td>180°</td></tr><tr><td>Error (%)</td><td>0.027</td><td>0.037</td><td>0.051</td><td>0.124</td><td>0.193</td><td>0.207</td><td>0.312</td><td>0.361</td></tr></table>

step, the velocity components were prescribed by

$$ \left\{ \begin{aligned} V_{\phi} &= \frac{1}{\Delta t} 2R \sin \Delta\theta \sin \Delta\theta \\ V_{\psi} &= \frac{1}{\Delta t} 2R \sin \Delta\theta \cos \Delta\theta \end{aligned} \right. \qquad (31) $$

where $\Delta\theta$ 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

$$ \text{error} = \frac{1}{Nr_{\text{exact}}} \sum_{I}^{N} |r_{I}(\theta) - r_{\text{exact}}| \qquad (32) $$

Figure 7 shows the propagation of the $\phi$ level set and one instance of the $\psi$ level set. The front at various times is also shown. In Table II, we compare the numerical radius $r$ of the $\phi=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.
<figure><img src="image_9.png" /><figcaption>Figure 8. Side and top views of the beam with an initial crack.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Evolution of the crack front in the beam seen from the top.</figcaption></figure>

## 4.2. 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 $\alpha$ to the plane of the cross-section is used to initiate crack growth with $a=0.01$ m and $\alpha=45^\circ$ (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 field 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 specific 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 first grows to the lateral free surfaces, then normal to the surface. No changes are needed in the algorithm to account for these topological changes.

## 4.3. A lens-shaped crack

In this example, we consider a cube with a cusp crack subjected to hydrostatic tension as shown in Figure 12.
<figure><img src="image_11.png" /><figcaption>Figure 10. Crack and vector velocity field after 12 time steps.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Evolution of a semi-circular crack in a beam.</figcaption></figure>

The crack geometry is characterized by the radius $R$ and the azimuthal angle $\alpha$. 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 defined by $R = 0.005$ m and $\alpha = 45^\circ$ 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 field on the front.
<figure><img src="image_13.png" /><figcaption>Figure 12. Initial cusp crack subjected to hydrostatic tension.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Initial crack with the initial velocity field on the crack front.</figcaption></figure>

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.
<figure><img src="image_15.png" /><figcaption>Figure 14. Crack after 15 time steps.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. A cube with an inclined penny-shaped
crack subjected to a tensile loading.</figcaption></figure>

4.4. Inclined penny-shaped crack

In this last example, we consider the problem in Figure 15: a cube with an inclined penny-
shaped crack subjected to a tensile loading with h = 0.02 m and an initial crack defined 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.

5. CONCLUSIONS

A level set method for arbitrary non-planar cracks in three-dimensional bodies has been pre-
sented. A key feature that distinguishes this method from other level set methods is a technique
for the extension of the velocity field 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 surround-
ing the crack. This is an efficient 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
<figure><img src="image_17.png" /><figcaption>Figure 16. Evolution of the crack after 17 time steps for the inclined penny-shaped crack.</figcaption></figure>

adaptively refine the volume around the crack. This is particularly important for the discretization of the level set update, which requires a finer mesh to accurately capture the crack front. Methods with different refinements in the solid mechanics and level set discretization may also be desirable.

The level set technique couples naturally with the extended finite element method, wherein the discontinuous and near-tip asymptotic fields 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 simplifies the structure of the software and leads to great versatility in treating complex problems in crack growth. The extended and generalized finite element methods have recently seen rapid development [21–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 $\psi_I^* = w_I = 0$ for each node of the mesh.

*Step 2:* Loop on the finite elements.

We compute

$$\nabla \psi = \sum_{J=1}^{4} \nabla N_J \psi_J \quad \bar{V} = \int_{E} V \, dx$$
and the following coefficients for each node *I* of the tetrahedron:

$$
K_I = \bar{V} \frac{\nabla \cdot \nabla N_I}{|\nabla \psi|} \qquad \delta\psi = \sum_{J=1}^{4} K_J \psi_J
$$
$$
\alpha_I = \frac{H(\delta\psi_I / \delta\psi)}{\sum_{K=1}^{4} H(\delta\psi_K / \delta\psi)} \quad \text{with} \quad \delta\psi_I = H(K_I) \left(\sum_{L=1}^{4} H(K_L)\right)^{-1} \sum_{J=1}^{4} H(-K_J)(\psi_J - \psi_I)
$$

where *H* is the Heaviside function.

$$
\begin{align*}
\psi_I^* &= \psi_I^* + \alpha_I \delta\psi & (v \text{ is the volume of the finite element } E) \\
w_I &= w_I + \alpha_I v
\end{align*}
$$

Step 3: Loop on the nodes of the mesh for the time integration.

$$
\left\{
\begin{aligned}
\tilde{\psi}_I^{n+1} &= \psi_I^n - \Delta t \frac{(\psi_I^*)^n}{w_I^n} \\
\psi_I^{n+1} &= \frac{(\psi_I^n + \tilde{\psi}_I^{n+1})}{2} - \frac{\Delta t}{2} \frac{(\tilde{\psi}_I^*)^{n+1}}{\tilde{w}_I^{n+1}}
\end{aligned}
\right.
$$

ACKNOWLEDGEMENTS

The support of the Office of Naval Research and Department of Energy, to Northwestern University,
as well as the French support from Direction des Systèmes de Forces et de la Prospective, Délégation
Générale pour l'Armement, are gratefully acknowledged. We are also grateful to Jean-François Remacle
for his help in using the AOMD C++ library and Natarajan Sukumar for his comments.

REFERENCES

1. Sethian JA. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Material Science. Cambridge University Press: Cambridge, U.K., 1996.

2. Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. *Journal of Computational Physics* 1988; **79**(1):12–49.

3. Burchard P, Cheng L-T, Merriman B, Osher S. Motion of curves in three spatial dimensions using a level set approach. *Journal of Computational Physics* 2001; **170**(2):720–741.

4. Cheng L-T, Burchard P, Merriman B, Osher S. Motion of curves constrained on surfaces using a level set approach. UCLA CAM Report 00-32, University of California at Los Angeles (2000).

5. Peng D, Merriman B, Osher S, Zhao H, Kang M. A PDE-based fast local level set method. *Journal of Computational Physics* 1999; **155**:410–438.

6. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. *International Journal for Numerical Methods in Engineering* 1999; **45**(5):601–620.

7. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**:131–150.

8. Dolbow J, Moës N, Belytschko T. Discontinuous enrichment in finite elements with a partition of unity method. *Finite Elements in Analysis and Design* 2000; **36**:235–260.

9. Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. *International Journal for Numerical Methods in Engineering* 2000; **48**: 1741–1760.

10. Moës N, Gravouil A, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets—Part I: Mechanical model. *International Journal for Numerical Methods in Engineering* 2002; **53**:2549–2568.

11. Stolarska M, Chopp DL, Moës N, Belytschko T. Modeling crack growth by level sets in the extended finite element method. *International Journal for Numerical Methods in Engineering* 2001; **51**:943–960.
12. Sukumar N, Chopp DL, Moran B. Extended finite element method and fast marching method for three-dimensional fatigue crack propagation. *Engineering Fracture Mechanics* 2001, submitted for publication.

13. Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modeling. *International Journal for Numerical Methods in Engineering* 2000; **48**:1549–1570.

14. Sukumar N, Chopp DL, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite element method. *Computer Methods in Applied Mechanics and Engineering* 2001, accepted for publication.

15. Chopp DL. Computing minimal surfaces via level set curvature flow. *Journal of Computational Physics* 1993; **106**:77–91.

16. Chen S, Merriman B, Osher S, Smereka P. A simple level set method for solving Stefan problems. *Journal of Computational Physics* 1997; **135**:8–29.

17. Hou TY, Li Z, Osher S, Zhao H. A hybrid method for moving interface problems with application to the Hele–Shaw flow. *Journal of Computational Physics* 1997; **134**:236–252.

18. Barth TJ, Sethian JA. Numerical schemes for the Hamilton–Jacobi and level set equations on triangulated domains. *Journal of Computational Physics* 1998; **145**:1–40.

19. Remacle JF, Geuzaine C. Gmsh finite element grid generator. Available at www.geuz.org/gmsh, 1998.

20. Lazarus V, Leblond J-B, Mouchrif S-E. Crack front rotation and segmentation in mixed mode I+III or I+II+III—Part I: calculation of stress intensity factors. *Journal of the Mechanics and Physics of Solids* 2001; **49**(7):1399–1420.

21. Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. *International Journal for Numerical Methods in Engineering* 2001; **50**:993–1013.

22. Duarte C, Hamzeh O, Liska T. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**:2227–2262.

23. Melenk J, Babuska I. The partition of unity finite element method: basic theory and applications. *Computer Methods in Applied Mechanics and Engineering* 1996; **39**:289–314.

24. Strouboulis T, Copps K, Babuska I. The generalized finite element method: an example of its implementation and illustration of its performance. *International Journal for Numerical Methods in Engineering* 2000; **47**(8): 1401–1417.