under multiple loading with an outward revolution crack are confronted to a reference solution, and those for a double edge-crack specimen are compared to the theoretical solution.

## 2 PRESENTATION OF THE METHOD

An elastic body $\Omega$ with a crack is considered. We assume that the crack surfaces are traction-free. We consider the case of linear elasticity. The approach is limited to small deformations. This post-processing method uses the displacement field, solution of the variational formulation and two level sets in order to define a local basis to the crack front in 3D. We note that the solution displacement field can be obtained by different methods such as a standard FE method if the crack is meshed or the XFEM method [2-3,5] if the crack is not meshed.

### 2.1 Level sets representation and local basis

The crack is treated as a single surface $\Gamma_{cr}$. Moes et al. [5] described a crack geometry in 3D with two signed distance functions. The signed distance function $\ln s(x)$ defines the surface of the crack. It is given by the shortest distance of any point $x$ to the extended crack surface $\Gamma_{ext}$ ($\Gamma_{ext}$ is actually a smooth extension of $\Gamma_{cr}$ all across $\Omega$). The sign is given by the position of the point $x$: +1 if the point is “above” the surface, and -1 if the point is “below” the surface. Since the crack surface is well known, another signed function $\text{lst}(x)$ is needed to define the crack front. It is given in finding the nearest point $\underline{x}$ of the crack front $\Gamma_0$ from the current point $x$: $\text{lst}(x) = n \cdot (x - \underline{x})$ where $n$ is the outward normal to the crack front inside $\Omega$ at $\underline{x}$ (Figure 1). We notice that our choice for the level set functions may differ from other papers available in the literature [6,11] where the functions are given in terms of analytical expressions, and a special treatment must be carried out to ensure their orthogonality. In our case, the level sets are real distance functions, chosen to be orthogonal. We underline the fact that the crack surface $\Gamma_{cr}$ is given by $\ln s(x)=0$ and $\text{lst}(x)<0$, and that the crack front $\Gamma_0$ is given by $\ln s(x)=0$ and $\text{lst}(x)=0$ (Figure 1). Note that level sets can be introduced even if the crack is meshed. With this representation, it is easy to create a local basis on the crack front. Sukumar [6] uses the gradients of the level sets for such a basis. From the nodal values of the level sets ($\ln s^i$ and $\text{lst}^i$) the value of the gradients for each finite element $\Omega_e$ is computed:

$$ \nabla \ln s_j^{\text{elt}} = \sum_{i \in N_e} \ln s^i N_{j,i} \quad j = 1,2,3 \quad (1) $$

$$ \nabla \text{lst}_j^{\text{elt}} = \sum_{i \in N_e} \text{lst}^i N_{j,i} \quad j = 1,2,3 \quad (2) $$

where $N^i$ are the classical shape functions and $N_e$ the nodes of the finite element $\Omega_e$. Then we define the nodal value of the gradient at node i as the average of the gradients of the elements

<figure><img src="image_1.png" /><figcaption>Figure 1: The two iso-zero level sets defining the crack location</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2: Local basis on the front and cylindrical co-ordinates system</figcaption></figure>
connected to this node: $\nabla \mathbf{I} \mathbf{s}^1$ and $\nabla \mathbf{I} \mathbf{s}^1$. Note that the level sets and their gradients are approximated by the same shape functions as the displacement field.

Now, if we choose $\mathbf{e}_1=\nabla \mathbf{I} \mathbf{s}^1$, $\mathbf{e}_2=\nabla \mathbf{I} \mathbf{s}^1$ and $\mathbf{e}_3=\mathbf{e}_1 \wedge \mathbf{e}_2$, we have defined a local basis at any point of the crack front (even at any point of $\Omega$). In this basis (Figure 2), the expressions of the auxiliary fields are given in [5]. All the following quantities are expressed in this basis.

## 2.2 Domain integral and bilinear form

For linear thermo-elastostatics, in the absence of body forces and assuming traction-free crack surfaces, the volume form of the domain integral is given by [10]:

$$J(\mathbf{\theta}) = \int_{\Omega} \sigma_{ij} \mathbf{u}_{ip} \mathbf{\theta}_{pj} - \psi \mathbf{\theta}_{k,k} - \frac{\partial \psi}{\partial T} T_{,k} \mathbf{\theta}_k d\Omega \quad (3)$$

where $\sigma$ is the Cauchy stress tensor, $\mathbf{u}$ the displacement field, $T$ the temperature field, $\psi$ the energy density and $\mathbf{\theta}$ the virtual crack extension velocity field. $\mathbf{\theta}$ is tangent to the crack faces and normal to the front (i.e. $\mathbf{\theta}$ directed by $\mathbf{e}_1$).

It can be shown [9] that the energy release rate $G$ can be considered as a symmetric bilinear form of the displacement field $\mathbf{u}$. Let $g(\mathbf{u},\mathbf{v})$ be this form. In this case:

$$g(\mathbf{u}, \mathbf{v}) = \frac{1}{2} \int_{\Omega} \frac{\partial B}{\partial \nabla \mathbf{u}} \cdot (\nabla \mathbf{v} \cdot \nabla \mathbf{\theta}) + \frac{\partial B}{\partial \nabla \mathbf{v}} \cdot (\nabla \mathbf{u} \cdot \nabla \mathbf{\theta}) - B(\mathbf{u}, \mathbf{v}) \operatorname{div} \mathbf{\theta} - \frac{\partial B}{\partial T} T_{,k} \mathbf{\theta}_k d\Omega \quad (4)$$

where $B(\mathbf{u}, \mathbf{v}) = \varepsilon(\mathbf{u}) : \Lambda : \varepsilon(\mathbf{v})$ is a bilinear form of the energy density and $\Lambda$ the Hooke's tensor. And if the displacement field $\mathbf{u}$ is the solution of the elastic problem,

$$g(\mathbf{u}, \mathbf{u}) = J = \int_{\Omega} \sigma(\mathbf{u}) : (\nabla \mathbf{u} \nabla \mathbf{\theta}) - \psi(\varepsilon(\mathbf{u})) \operatorname{div} \mathbf{\theta} - \frac{\partial \psi}{\partial T} (\nabla T \cdot \mathbf{\theta}) d\Omega \quad (5)$$

in linear elasticity, it can be decomposed in a regular part and a singular part. The SIFs appear in the singular part:

$$\mathbf{u} = \mathbf{u}_R + K_I \mathbf{u}_S^1 + K_{II} \mathbf{u}_S^2 + K_{III} \mathbf{u}_S^3 \quad (6)$$

It can be proved than $\mathbf{u}_S^1$, $\mathbf{u}_S^2$ and $\mathbf{u}_S^3$ are orthogonal for the scalar product defined by $g(\mathbf{u},\mathbf{v})$ and that the contribution of the regular terms is equal to zero. Therefore, the SIFs are given by:

$$g(\mathbf{u}, \mathbf{u}_S^1) = K_I g(\mathbf{u}_S^1, \mathbf{u}_S^1) \quad (7)$$

$$g(\mathbf{u}, \mathbf{u}_S^2) = K_{II} g(\mathbf{u}_S^2, \mathbf{u}_S^2) \quad (8)$$

$$g(\mathbf{u}, \mathbf{u}_S^3) = K_{III} g(\mathbf{u}_S^3, \mathbf{u}_S^3) \quad (9)$$

where $g(\mathbf{u}_S^1, \mathbf{u}_S^1)$, $g(\mathbf{u}_S^2, \mathbf{u}_S^2)$ and $g(\mathbf{u}_S^3, \mathbf{u}_S^3)$ are determined using the Irwin formula.

To compute the SIFs, the three bilinear forms mentioned above are just needed, with $\mathbf{u}_S^1$, $\mathbf{u}_S^2$ and $\mathbf{u}_S^3$ the displacement fields for a plane crack in a infinite media. These expressions are to be found in [5]. We underline the fact that the latter expressions only use $\mathbf{u}_R$, $\mathbf{u}_S^1$, $\mathbf{u}_S^2$ and $\mathbf{u}_S^3$. That means that auxiliary strains are given in terms of displacement gradients and that stresses are derivated from the constitutive law. Consequently, these fields have a zero divergence, which means that locally the crack surface is a plane and the crack front is linear.

## 2.3 G-theta method in Code_Aster

The G-theta method is used to compute the energy release rate in Code_Aster [10]. The particularity of the method is to consider a crack front discretized with $nno$ nodes (or points), so
<figure><img src="image_3.png" /><figcaption>Figure 3: 3D tore around the crack front</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4: Section of the tore</figcaption></figure>

that a curvilinear coordinate system can be defined. The number of points of the crack front (nno) is chosen by the user. In linear elasticity the energy release rate verifies the variational equation:

$$ \int_{\Gamma_0} G(s)\boldsymbol{\theta}(s) \cdot \mathbf{n}(s)ds = J(\boldsymbol{\theta}), \forall \boldsymbol{\theta} \in \Theta \quad (10) $$

where $\mathbf{n}(s)$ denotes the normal at the crack front $\Gamma_0$ at the curvilinear coordinate $s$, and $J(\boldsymbol{\theta})$ is the value of the previous domain integral for all the structure. This is an important point of the method. Most of other methods (Moes et al. [5], Sukumar [6]) compute the SIFs using a local approach for the integration, i.e. compute $J$ on a box (virtual extension domain) around a point. Our method avoids the delicate choice of the box size by computing $J$ on the whole domain, for various virtual velocity fields $\boldsymbol{\theta}^i$ such as $\boldsymbol{\theta}^i(s) \cdot \mathbf{n}(s) = N_i(s)$, $i \in \{1, nno\}$ where $nno$ is the number of nodes (or points) of the crack front. The theta fields represent a virtual transformation of the domain in a crack propagation. They should satisfy $\boldsymbol{\theta}^i \cdot \mathbf{n} = 0$ on $\partial\Omega$. Moreover, the field must be regular. From a numerical point of view, it is interesting to use constant fields in the vicinity of the crack front $\Gamma_0$. We propose a convenient theta field considering a virtual tore around $\Gamma_0$ (Figure 3). $R_{\text{inf}}$ and $R_{\text{sup}}$ are the respective lower and upper radius of the tore. At each point of $\Gamma_0$, at the curvilinear co-ordinate $s$, a normal plane $P$ to $\Gamma_0$ can be defined (Figure 4), in which the theta field can be defined. Since $\boldsymbol{\theta}^i = \alpha \mathbf{e}_1$, the values of $\alpha$ are:

$$ \alpha = Ni(s) \quad \text{for } 0 \le r(s) \le R_{\text{inf}}(s) \quad (11) $$

$$ \alpha = 0 \quad \text{for } R_{\text{sup}}(s) \le r(s) \quad (12) $$

$$ \alpha \text{ is a linear function of } r(s) \text{ for } R_{\text{inf}}(s) \le r(s) \le R_{\text{sup}}(s) \quad (13) $$

The relations between the $J(\boldsymbol{\theta}^i)$ and the SIFs along the front are now detailed. One of the possibilities is to use the shape functions $N_j$ to discretize $\boldsymbol{\theta}^i$ and $G$, so that:

$$ \boldsymbol{\theta}^i(s) = \sum_{j \in nno} \boldsymbol{\theta}_j^i N_j(s) \quad (14) $$

$$ G(s) = \sum_{j \in nno} G_j N_j(s) \quad (15) $$

If we replace in eqn (10) the discretization of $G$ and $\boldsymbol{\theta}^i$, we have to solve the following linear system:

$$ \sum_{j=1}^{nno} G_j \int_{\Gamma_0} N_i(s) N_j(s) ds = J(\boldsymbol{\theta}^i), \quad i \in \{1, nno\} \quad (16) $$

In this system, the $J(\boldsymbol{\theta}^i)$ are known by eqn (5) and the unknowns are the $G_j, j \in \{1, nno\}$. Now, to compute the SIFs, the method is very similar. For modes I, II and III we discretize $K^{I,II,III}$ with the shape functions, so that:
$$K^{1,II,III}(s) = \sum_{j \in nno} K_j^{1,II,III} N_j(s) \qquad (17)$$

and the new system is:

$$\sum_{j=1}^{nno} K_j^{1,II,III} \int_{\Gamma_0} N_i(s) N_j(s) ds = g(\mathbf{u}, \mathbf{u}_s^{1,II,III})_{0i} , \quad i \in \{1, nno\} \qquad (18)$$

## 3 NUMERICAL EXAMPLES

Applications of the present formulation are shown in the following examples: a cylinder under multiple loading (traction and torsion) with an outward revolution crack and a double edge-crack specimen under an original loading (auxiliary fields imposed all over the specimen).

### 3.1 Cylinder under multiple loading

A cylinder with an outward revolution crack is considered. The crack is a circular crown in a plane orthogonal to the axis of the cylinder. The parameters *a* and *b* determine the radius of the inner sane cylinder and the radius of the outer cylinder (Figure 5). In order to simulate an infinite medium, the height of the cylinder is *h* = 10 *b*. The material is a standard steel with a Young Modulus E = 205000 MPa and a Poisson's ratio ν = 0.3. Mode I traction and mode III torsion are applied on the left edge of the cylinder. The right edge is fully constrained. The solutions of this benchmark are given in [12]. The mesh contains about 13000 linear elements (hexahedral and pentahedral elements). The SIFs are computed for each nodes of the front (20 nodes for mesh 1 and 40 nodes for mesh 2). We observe that the difference between the numerical result and the reference solution for mesh 2 are about 1.25 per cent for KI and 5.91 per cent for KIII (Table 1). The “displacement-jumps method” used in *Code_Aster* gives errors about 10 per cent for KI and KIII, which indicates that the new method is more accurate

### 3.2 Double edge-crack specimen

A double edge-crack specimen is analysed. The dimensions are shown at Figure 6, with w = 50 mm. The mesh is free and composed of about 10000 pentahedral and hexahedral finite

<figure><img src="image_5.png" /><figcaption>Figure 5: Outward revolution crack in a cylinder</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6: Double edge-crack specimen</figcaption></figure>

<table><caption>Table 1: Error between SIFs for two meshes for the revolution crack problem</caption><thead><tr><th rowspan="2"></th><th colspan="2">Mesh 1 (13000 elements)</th><th colspan="2">Mesh 2 (17000 elements)</th></tr><tr><th>KI</th><th>KIII</th><th>KI</th><th>KIII</th></tr></thead><tbody><tr><th>Error min (%)</th><td>1.39</td><td>6.11</td><td>0.84</td><td>5.84</td></tr><tr><th>Error max (%)</th><td>2.07</td><td>6.29</td><td>1.25</td><td>5.91</td></tr></tbody></table>
elements. We have chosen an original loading to test the method. Instead of imposing stresses on the top or bottom faces, the boundary conditions are nodal displacement imposed to be equal to the singular auxiliary fields [5]. In this case, if the mode I auxiliary displacement is imposed, we should obtain $KI = 1.0$ and $KII = KIII = 0$ (same process for modes II and III auxiliary fields). To illustrate this fact, the superposition of the three auxiliary modes is imposed. The SIFs are computed for each nodes on the crack front (20 nodes). The advantage of that kind of loading is that the solution displacement is also given as a data in eqn (18). As a matter of fact, we only test the ability of the method to recover correct SIFs. The results are compared with the reference solution which is $KI = KII = KIII =1.0$. The difference between the SIFs computed and the SIFs expected is about 1.0% for KI and KII and 0.4% for KIII.

## 4 CONCLUSIONS

The formulation and the implementation of the “G-theta” method for the computation of the stress intensity factors (SIFs) in linear elasticity is described. In this method, a bilinear form of the energy release rate is used to extract the SIFs separately. By considering several virtual crack extension velocity fields, the J-integrals are performed on the whole domain instead of what is usually done with an integration box and a virtual extension domain. The crack is represented by two level sets which permit to create a local co-ordinates system easily.
The performance of the “G-theta” method for 3D static planar cracks is studied. Benchmark mixed mode problems were solved for different crack geometries. This method appears to be more general and more accurate than the “displacement-jumps method” used before in *Code\_Aster*.
Improvements of the method are carried out. They concern the introduction of second-order finite elements and the adaptation of X-FEM to problems dealing with contact to simulate crack growth in 3D. Ultimately propagation criteria will be discussed in that framework.

## 5 REFERENCES

1. Henshell RD, Shaw KG. Crack tip finite elements are unnecessary. *International Journal For Numerical Methods In Engineering* 9:495-507, 1975

2. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. *International Journal For Numerical Methods In Engineering* 45:601-620, 1999

3. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal For Numerical Methods In Engineering* 46:131-150, 1999

4. Melenk JM, Babuška I. The partition of unity finite element method: Basic theory and applications. *Computer methods in applied mechanics and engineering* 139:389-314, 1996

5. 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* 53:2549-2568; 2002.

6. Sukumar N, Chopp DL, Moran B. Extended finite element method and fast marching method for three-dimensional fatigue crack propagation. *Engineering Fracture Mechanics* 70:29-48; 2003

7. Liu XY, Xiao QZ, Karihaloo BL. X-FEM for direct evaluation of mixed mode SIFs in homogeneous and bi-materials. *International Journal For Numerical Methods In Engineering* 59:1103-1118, 2004

8. Moës N. Utilisation des intégrales de domaines 3D pour l'extraction des facteurs d'intensité de contraintes de fissures en élastostatique représentées par « level sets ». *Internal report, École Centrale de Nantes* 2004

9. Visse E. Calcul des coefficients d'intensité de contraintes en thermo-élasticité linéaire plane. *EDF : Manuel de référence du Code\_Aster*, Document R7.02.05. 1995

10. Granet S, Boiteau O, Visse E. Taux de restitution de l'énergie en thermo-élasticité linéaire. *EDF : Manuel de référence du Code\_Aster*, Document R7.02.01. 2002

11. Gravouil A, Moës N, Belytschko T. 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* 53:2569-2588; 2002

12. Tada H, Harris PC, Irwin GR. The Stress Analysis Of Cracks Handbook. *Del Research Corporation, Hellertoxn, Pennsylvania* 1973