# Development of a finite element enriched method adapted for the case of multiple cracked structure

Emilien Durif, Julien Réthoré, Alain Combescure

► To cite this version:

Emilien Durif, Julien Réthoré, Alain Combescure. Development of a finite element enriched method adapted for the case of multiple cracked structure. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2010, 19, pp.217-228. 10.3166/ejcm.19.217-228. hal-00938573

HAL Id: hal-00938573

https://hal.science/hal-00938573

Submitted on 3 Jun 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Development of a finite element enriched method adapted for the case of multiple cracked structure

Emilien Durif — Julien Réthoré — Alain Combescure

Université de Lyon, Laboratoire LaMCoS (UMR/CNRS)
INSA Lyon, Bât. J. d'Alembert
20 avenue Albert Einstein, F-69621 Villeurbanne cedex
emilien.durif@insa-lyon.fr

**ABSTRACT.** The present method is based on the eXtended Finite Element Method and the partition of unity properties which allows the addition of enriched functions in the displacement field discretization. Those correspond to the analytical displacement field near a straight crack. The additional degrees of freedom correspond directly to the Stress Intensity Factors. This method is validated in mode 1 and mixed mode. A parametric study has been performed to analyse the effect of different parameters on the accuracy of the method. The case of two shading cracks has also been investigated.

**RÉSUMÉ.** Cette méthode se base sur la méthode des éléments finis étendus et les propriétés de partition de l'unité qui autorisent l'ajout de fonctions d'enrichissement dans la base du champ de déplacement. Celles-ci correspondent aux champs de déplacements analytiques au voisinage d'une fissure droite. Les degrés de liberté supplémentaires correspondent directement aux facteurs d'intensité des contraintes. Cette méthode est validée en mode 1 et en mode mixte. Une étude paramétrique permet de déterminer les paramètres influant sur la précision des résultats. Enfin un exemple permet d'étudier l'interaction de deux fissures en écran.

**KEYWORDS:** eXtended finite element method, analytical enrichment, interacting cracks.

**MOTS-CLÉS** : méthode des éléments finis étendus, enrichissement analytique, interaction de fissures.
# 1. Introduction

The eXtended Finite Element Method (X-FEM) allows crack modeling independently of the mesh ((Moes *et al.*, 1999) and (Belytschko *et al.*, 1999)). Thanks to the partition of unity properties of finite element shape functions (Melenk *et al.*, 1996), enrichment functions can be added in the displacement formulation. Thus discontinuities and crack tip singularities can be captured more accurately. Heavyside enrichment function has been used to represent the discontinuity (HX-FEM (Moes *et al.*, 1999)). The crack tip singularity is given by the basis of the analytical displacement field of Westergaard near the crack tip. This formulation leads to additional degrees of freedom. The enrichment is applied at the nodes of the elements which are cut by the crack. To obtain a more accurate estimation of the stress intensity factors a wider enriched zone should be used around the crack tip (Laborde *et al.*, 2005). In the case of a structure with many cracks, the numbers of additional DOF can be high. In this kind of problem, the calculation of the stress intensity factors is very interesting to point out the critical cracks in the structure. The calculation of the stress intensity factors is usually performed using the interaction integral method (Yau *et al.*, 1980) or the least square identification (McNeill *et al.*, 1987). The use of these methods could be very difficult in the case of crack interacting closely or even intersecting. Specific discontinuous basis functions have been used to make calculation of intersecting cracks (Budyn *et al.*, 2004), (Zi *et al.*, 2004) and (Daux *et al.*, 2000). However if the number of crack is high and if the crack tips become too closer the calculation could become very complicated. Analytical studies had been made for interacting crack (Kachanov, 1987) and (Li *et al.*, 2003) or intersecting cracks (Cheung *et al.*, 1992) using superposition principle to estimate stress intensity factors or equivalent elastic properties of a multiple cracked structure (Helsing, 2000). These examples are made for semi-infinite crack in infinite medium. To combine the advantages of the X-FEM in making calculation of finite cracked structure and those of the analytical method in the accuracy of stress intensity factors evaluation, the enrichment functions can be taken as the analytical displacement fields (Liu *et al.*, 2004) such as the formalism of Williams' series (Williams, 1957). The method presented in this paper is based on a specific analytical enrichment for small straight cracks which provides directly the stress intensity factors with only two additional degrees of freedom per crack.

In the next section, the displacement formulation is detailed. In Section 3, the capabilities of the method are tested in mode 1 and mixed mode. The Section 4 discusses about the influence of different calculation parameters on the accuracy. Finally in Section 5, the interaction of two shading cracks is investigated.

# 2. Enrichment formulation

The enrichment developed herein is based on the analytical displacement field of Griffith-Inglis (Weertman, 1996). In this model cracks are accounted for such as a continuous sum of dislocations. The discontinuous displacement field ($Ud(x, y)$) in mode 1 and 2 linearly depends on the stress intensity factors at the crack tips ($K_i$)
Equation [1]). In this formulation cracks are symmetric. Thus the stress intensity factors are the same at both tips of the crack and then only one SIF is considered for each mode.

$$ \underline{Ud}(x, y) = \sum_{i=1}^{2} K_i \underline{F}_i(x, y) \qquad [1] $$

The total displacement field is decomposed into a discontinuous and a continuous part [2].

$$ \underline{U}(x, y) = \underline{Uc}(x, y) + \underline{Ud}(x, y) \qquad [2] $$

The continuous part is described using standard finite elements ($\underline{Uc}(x, y)$) and the discontinuous one ($\underline{Ud}(x, y)$) with the analytical enrichment functions ($\underline{F}_j(x, y)$) by using the partition of the unity properties (Melenk et al., 1996) (Equation [3]).

$$ \underline{U}^h(x, y) = \sum_{i=1}^{N} N_i(x, y) \underline{u}_i + \sum_{j=1}^{2} \sum_{i=1}^{N_H} K_j^i \cdot N_i(x, y) \underline{F}_j(x, y) \qquad [3] $$

$N_i(x, y)$ are the standard finite element shape functions. $N$ is the total number of nodes in the structure. A domain around the crack tip is chosen and $N_H$ defines the number of nodes inside this domain. Two additional degrees of freedom are considered ($K_1^i$ and $K_2^i$) at each enriched node $i$. The following condition $\forall i \in N_H$, $K_j^i = K_j$ is enforced to decrease the number of enriched degrees of freedom. Eventually there are only two additional degrees of freedom per crack. The corresponding eXtended Finite Element approximation is given by Equation [4]. $\Phi(x, y)$ is a cut-off function (Chahine et al., 2008). It is equal to 1 in all elements containing $N_H$ nodes and decreases linearly to 0 in adjacent elements. The displacement field is finally written as,

$$ \underline{U}^h(x, y) = \sum_{i=1}^{N} N_i(x, y) \underline{u}_i + \sum_{j=1}^{2} K_j \underline{F}_j(x, y) \Phi(x, y) \qquad [4] $$

The formulation leads to the calculation of the global elastic stiffness matrix $\underline{K}$ (Equation [5]). It is shared into a $2N \times 2N$ square matrix which is the classical finite element stiffness matrix $\underline{K}_c$, a $2 \times 2$ square matrix which components represent the enriched terms ($\underline{K}_c$) and a $2N \times 2$ which is the coupling matrix ($\underline{C}$) (Equation
[5]). The integration of enriched functions has ever been studied by several authors. The first developed idea was to introduce integration subcells into enriched elements. The construction of the subcells depends on the position of the crack into the element and it is often complicated to elaborate. Integration of discontinuous function using the contour integral transformation as proposed in (Ventura, 2006) and in (Ventura et al., 2009) has been used to circumvent this drawback. In our case the continuous components of the matrix $K_c$ are integrated with four gauss points. The enriched function are not given in a simple way and the integration cannot be evaluated exactly with usual quadrature rules. To integrate the coupled and the enriched components ($K_e$ and $C$) the elements are subdivided into a given number of square subcells in the reference frame. Inside each of these subcells, the integration is performed with four gauss points (Figure 1 and Equation [6]). This integration is an approximation but it gives accurate results with at least 10 × 10 subcells ($n_{cell}$).

$$ K = \begin{bmatrix} K_c & C \\ C^t & K_e \end{bmatrix} \begin{array}{l} \updownarrow N \\ \updownarrow 2 \end{array} \qquad [5] $$

$$ \int_{\Omega_e} f(\xi, \eta) d\xi d\eta = \sum_{i=1}^{n_{cell}} \sum_{j=1}^{4} \omega_j f(\xi_{ij}, \eta_{ij}) \qquad [6] $$

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Enriched and coupled components integration with square subcells in the reference frame of the element</figcaption></figure>
In this formulation the additional degrees of freedom correspond directly to stress intensity factors $K_j$ which are directly calculated during the resolution of the linear system without any postprocessing.

### 3. Study of the analytical extended finite element method in mode 1

The validity of the method is verified with a simple traction test on a sample containing a crack inclined by the angle $\beta$ (Figure 2). In mode 1, $\beta$ is equal to 0 and the relative crack length $\lambda = \frac{a}{W}$ increases from 0.02 to 0.25. a and W are respectively, the half of the crack length, and the width of the structure. The finite size correction on $K_1$ depends upon $\lambda$, it is given for example in (Lemaitre J., 1985) (Equation [7]). In mixed mode, stress intensity factors depends on $\beta$. They are given by Equation [8]. In this case, the crack length is fixed, $\lambda = 0.1$ and $\beta$ varies between 0 and $\pi$. The material is chosen with a Young's Modulus of $E = 210\,GPa$ and a Poisson's ratio of $\nu = 0.3$. A unit applied stress is used $\sigma_A = 1\,MPa$. A uniformed mesh of 45 × 45 quadrangles is used. Cut-off function support dimensions are defined by a radius equal to $r = 0.4W$. The evolution of the stress intensity factors with $\lambda$ is compared with those given by Equation [7] in Figure 3. A good agreement is observed between the two curves. However a bigger difference is observed when $\lambda$ is high. This is due to the fact that the crack length becomes similar to the cut-off function diameter. The accuracy depends on the relative size of the crack, the structure and the cut-off function radius. Figure 4 compares the stress intensity factors evaluated by the proposed method and by Equation [8]. No fluctuation is observed in mixed mode. $K_1$ is closed to 1 and $K_2$ closed to 0 when $\beta = 0$. With $\beta = \frac{\pi}{4}$, $K_2$ is maximal and equal to $K_1$. This study shows the validity of the method. However dimensional parameters seem to influence the accuracy of the stress intensity factors calculation.

$$K_1 \approx \sigma_A \sqrt{\pi a} [1 - 0,025\lambda^2 + 0,06\lambda^4] \left[\cos\left(\frac{\pi}{2}\lambda\right)\right]^{-1/2} \quad [7]$$

$$\left\{ \begin{array}{l} K_1 = \sigma_A \sqrt{\pi a} \cos^2(\beta) \\ K_2 = \sigma_A \sqrt{\pi a} \cos(\beta) \sin(\beta) \end{array} \right. \quad [8]$$

### 4. Parametric study of the enrichment

All the studied parameters are listed in the table 1. The convergence, as a function of the number of element along one direction h/W, where h is the mesh width and W is the width of the structure, is studied with a fixed crack length and a fixed cut-off function size ($\lambda = a/w = 0.1$ and $\gamma = r/w = 0.2$). The mesh is uniform and it
<figure><img src="image_3.png" /><figcaption>**Figure 2.** *Inclined centered crack in a traction loading*</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** *Stress intensity factor evolution in mode 1 with the relative size of the crack ($\lambda = a/W$): comparison with analytical models*</figcaption></figure>

consists in quadrangle elements. The convergence is studied with respect to the analytical displacement applied on the boundary nodes of the structure. This prescribed displacement is given by the same function as the one we use as an enrichment. The boundary conditions thus depends on prescribed stress intensity factors. A relative error is calculated between the prescribed stress intensity factors and the computed ones. The relative error ($\epsilon_r$) is plotted in Figure 5. The error is less than 10% even with a coarse mesh of 10 × 10 elements. The convergence rate is about 0.62 which corresponds to the order of magnitude for usual X-FEM approach with direct extraction of stress intensity factors (Nicaise *et al.*, 2009).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** *Calculation of the stress intensity factors in mixed mode with the inclination (β) ((a) mode 1 and (b) mode 2)*</figcaption></figure>

<table><caption>**Table 1.** *Parameters of the parametrical study*</caption><thead><tr><th>Parameters</th><th>λ</th><th>γ</th><th>h/W</th></tr></thead><tbody><tr><td>Relations</td><td>λ = <sup>a</sup>/<sub>W</sub></td><td>γ = <sup>r</sup>/<sub>W</sub></td><td>...</td></tr><tr><td>Intervals</td><td>λ ∈ [0.05; 0, 2]</td><td>γ ∈ [0.2; 0.5]</td><td>h/W ∈ [0.012; 0.067]</td></tr></tbody></table>

The study in traction mode emphasises boundary effects as the crack length and the cut-off function radius become similar to the structure width. The accuracy of the calculation is studied with the simultaneous variation of relative crack length (λ) and the relative cut-off function radius (γ). The values of interest are given in Table 1. When γ increases and λ decreases, The relative error (εr) on the evaluation of the stress intensity factors decreases (Figure 6). With a fixed value of λ, the precision increases with the size of the enriched zone. The cut-off function has to cover a sufficiently large area around the crack to ensure a correct accuracy. The ratio γ/λ has to be controlled in order to ensure a good precision. In Figure 6 a plane corresponding to εr = 0.05 indicates the values of λ and γ to provide a relative error lower than 5%.

## 5. Multiple cracks

A multiple cracked structure is considered. The material is supposed to be linear elastic and isotropic. The displacement field falls into the continuous part given by the standard finite element basis and the discontinuous part given by the analytical enrichment functions. Each of the discontinuous displacement is given by the corresponding
<figure><img src="image_6.png" /><figcaption>Figure 5. Evolution of the relative error $\epsilon_r$ with the relative mesh size (h/W)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Relative error in the calculation of $K_1$ with the relative variation of the cut-off function ($\gamma = r_B/W$) and the crack length ($\lambda = a/W$).</figcaption></figure>

enrichment function $F_j^k$. j corresponds to the mode and k is the number of the crack. Thus, the displacement field is searched based on the following decomposition:

$$ \underline{U}^h(x, y) = \sum_{i=1}^{N} N_i(x, y)\underline{u}_i + \sum_{k=1}^{N_{fiss}} \sum_{j=1}^{2} K_{kj} \underline{F}_j^k(x, y)\Phi_k(x, y) \quad [9] $$
The number of degrees of freedom is given by the standard finite elements part and two additional degrees of freedom per crack. Those correspond directly to the stress intensity factor in mode j (1 or 2) for the crack k ($K_{kj}$) as in the previous sections of the paper. Each crack has a cut-off function $\Phi_k(x, y)$. These functions can be separated or not. Two shading interacting cracks are considered (Figure 7). $\delta = \Delta/W$ is the parameter which controls the relative distance between the two cracks. Figure 7 also illustrates how two intersecting cut-off functions are built. Figure 8 gives the evolution of $K_1$ with $\delta$. The parameter $\lambda$ is chosen equal to 0.1 thus the crack length is five time smaller than the structure width. To evaluate the interaction between the cracks, two calculations have been performed. The first one considers one crack with a distance of $\Delta$ from the median line. A second one is carried out with two cracks such as described in Figure 7. In the first case $K_1$ is not constant with $\delta$ due to the boundary effects observed in section 3. When the crack is far from the boundaries $K_1$ remain almost constant. The calculation with two cracks shows that the second crack has a strong effect on the values of stress intensity factors. The closer the cracks, the lower the stress intensity factors are. Moreover the values tend to 0.5 which means that the two cracks tend a lonely crack with a $K_1$ equal to 1. However when $\delta$ is high the two curves are close because the interaction between the two cracks is small. This is confirmed by the amplified deformed shapes (pictures 9 (a) and (b)). When $\delta$ is small the two cracks are very close and it seems that the crack opening is lower. This observation is in agreement with the evolution of the stress intensity factor with $\delta$.

The cut-off functions defines the numerically influenced area of each crack. Despite the fact that the two cut-off functions intersect or not, the evolution of the stress intensity factor remains continuous. This observation shows that the cut-off function enables a good mechanical connection between the continuous and the discontinuous parts.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** *Study of two shadowed cracks*</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Evaluation of $K_1$ in the case of two shadowed cracks and the case of one crack</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Deformed mesh in the case of two shadowed cracks with an amplification of $50 \cdot 10^3$ with $\delta = 0.25$ (a) and $\delta = 0.025$ (b)</figcaption></figure>

**6. Conclusion and discussion**

This work shows that the X-FEM method is also a useful tool for multiple crack problems when analytical enrichment functions are adopted. The analytical enrichment used in the present paper corresponds to the displacement around a crack of finite size under symmetric loading (the values of the stress intensity factors are the same at both tips of the crack). The procedure has first been validated in a structure with one crack. Results show a good accuracy in the estimation of stress intensity factors. However calculation parameters such as crack length, cut-off function width, have to be controlled in order to ensure the quality of the numerical results. Cracks have to
stay relatively small compared to the structure length. However, the convergence rate
in term of stress intensity factors is the same as for other X-FEM based methods with
direct estimation of stress intensity factors. It is shown that crack interactions can be
efficiently investigated using the proposed method. As only two additional degrees
of freedom are required per crack, one can imagine to study structures with a large
number of cracks for a very limited computational cost.

7. References

Belytschko T., Black T., “Elastic crack growth in finite elements with minimal remeshing”, *International Journal for Numerical Methods in Engineering*, vol. 45, n° 5, p. 601-620, 1999.

Budyn E., Zi G., Moes N., Belytschko T., “A method for multiple crack growth in brittle materials without remeshing”, *Int. J. Numer. Meth. Engng*, vol. 61, p. 1741-1770, 2004.

Chahine E., Laborde P., Renard Y., “Crack tip enrichment in the XFEM method using a cut-off function”, *Int. J. Numer. Meth. Engng*, vol. 75, n° 6, p. 629-646, 2008.

Cheung Y., Woo C., Wang Y., “A general method for multiple crack problems in a finite plate”, *Computational Mechanics*, vol. 10, n° 5, p. 335-343, 1992.

Daux C., Moes N., Dolbow J., Sukumar N., Belytschko T., “Arbitrary branched and intersecting cracks with the extended ÿnite element method”, *Int. J. Numer. Meth. Engng*, vol. 48, p. 1741-1760, 2000.

Helsing J., “Fast and accurate numerical solution to an elastostatic problem involving ten thousand randomly oriented cracks”, *International Journal of Fracture*, vol. 100, n° 4, p. 321-327, 2000.

Kachanov M., “Elastic solids with many cracks- A simple method of analysis”, *International Journal of Solids and Structures*, vol. 23, n° 1, p. 23-43, 1987.

Laborde P., Pommier J., Renard Y., Salaun M., “High order extended finite element method for cracked domains”, *Int. J. Numer. Meth. Engng*, vol. 64, p. 354-381, 2005.

Lemaitre J. C. J., *Mecanique des materiaux solides*, Dunod, 1985.

Li Y., Tham L., Wang Y., Tsui Y., “A modified Kachanov method for analysis of solids with multiple cracks”, *Engineering Fracture Mechanics*, vol. 70, n° 9, p. 1115-1129, 2003.

Liu X., Xiao Q., Karihaloo B., “XFEM for direct evaluation of mixed mode SIFs in homogeneous and bi-materials”, *International Journal for Numerical Methods in Engineering*, vol. 59, n° 8, p. 1103-1118, 2004.

McNeill S., Peters W., Sutton M., “Estimation of stress intensity factor by digital image correlation”, *Engineering fracture mechanics*, vol. 28, n° 1, p. 101-112, 1987.

Melenk J., Babuška I., “The partition of unity finite element method: basic theory and applications”, *Computer methods in applied mechanics and engineering*, vol. 139, n° 1-4, p. 289-314, 1996.

Moes N., Dolbow J., Belytschko T., “A finite element method for crack growth without remeshing”, *International Journal for Numerical Methods in Engineering*, vol. 46, n° 1, p. 131-150, 1999.
Nicaise S., Renard Y., Chahine E., “Optimal convergence analysis for the eXtended Finite Element Method”, 2009.

Ventura G., “On the elimination of quadrature subcells for discontinuous functions in the eXtended Finite-Element Method”, *International Journal for Numerical Methods in Engineering*, vol. 66, n° 5, p. 761-795, 2006.

Ventura G., Gracie R., Belytschko T., “Fast integration and weight function blending in the extended finite element method”, *Int. J. Numer. Meth. Engng*, vol. 77, p. 1-29, 2009.

Weertman J., *Disclocation based on fracture mechanics*, World Scientific, 1996.

Williams M., “On the stress distribution at the base of a stationary crack”, *ASME Journal of Applied Mechanics*, vol. 24, p. 109-114, 1957.

Yau J., Wang S., Corten H., “A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity”, *Journal of Applied Mechanics*, vol. 47, p. 335, 1980.

Zi G., Song J., Budyn E., Lee S., Belytschko T., “A method for growing multiple cracks without remeshing”, *Modelling and Simulation in Materials Science and Engineering*, vol. 12, p. 901-915, 2004.