# Energy release rate computation using material differentiation of elastic BIE in 3-D elastic fracture

Marc Bonnet, Haihong Xiao

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marc Bonnet, Haihong Xiao. Energy release rate computation using material differentiation of elastic BIE in 3-D elastic fracture. Computational Mechanics '95, 1995, Hawaii, United States. pp.2856-2861, 10.1007/978-3-642-79654-8_473. hal-00122018

HAL Id: hal-00122018

https://hal.science/hal-00122018v1

Submitted on 22 Oct 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Energy release rate computation using material differentiation of elastic BIE in 3-D elastic fracture.

Marc BONNET, Haihong XIAO

Laboratoire de Mécanique des Solides (URA CNRS 317)

Ecole Polytechnique, Palaiseau, France*

# 1 Introduction

The energy release rate $G(s)$, function of the arc length $s$ along the front $\partial\Gamma$ of a crack $\Gamma$:

$$ \int_{\partial\Gamma} G(s) \delta\ell(s) \, ds = -\delta W \quad (1) $$

where $\delta W$ is the perturbation of the elastic potential energy at equilibrium $W$ induced by a crack front normal extension $\delta \ell$, the load being kept fixed, is one of the basic quantities involved in elastic fracture mechanics. In linear fracture mechanics, $G$ is linked to the stress intensity factors $K_I(s), K_{II}(s), K_{III}(s)$ through Irwin formula:

$$ G(s) = \frac{1 + \nu}{2\mu} [K_I^2(s) + K_I^2(s)] + \frac{1}{2\mu} K_{III}(s) \quad (2) $$

($\mu$: shear modulus, $\nu$: Poisson ratio). In addition, $G$ has a clear thermodynamical meaning (Nguyen [6]) and plays a key role in Griffith-type crack extension criterions.

Thus the consideration of perturbations of $W$ under fictitious body changes associated to virtual crack extensions provides a computational tool for elastic crack analysis. In the first numerical applications derivatives of $W$ are evaluated using small finite crack perturbations and finite differences (Hellen [5]). In later works (e.g. Delorenzi [1], Destuynder et al. [2]) the concept of material differentiation is applied to $W$, leading to rigorous formulations for $G$. This approach has led to FEM implementations (Wadier & Malak [3]). The present paper deals with a BIE formulation of the virtual crack extension approach in 3D elasticity.

# 2 Material differentiation

Let us consider a three-dimensional elastic body $\Omega_p$ whose shape changes according to a geometrical transformation given in Lagrangian form:

$$ y \in \Omega_0 \rightarrow y^p = \Phi(Y; p) \in \Omega^p \quad \text{where} \quad (\forall y \in \Omega) \quad \Phi(y; 0) = y \quad (3) $$

The parameter $p$ acts as a fictitious time and the “initial” configuration $\Omega = \Omega_0$ is conventionally associated with $p = 0$. A given domain evolution considered as a whole admits many different representations (3). The initial transformation velocity $\theta(y)$ is defined by

$$ \theta(y) = \Phi_{,p}(y; 0) \quad (4) $$

The material derivative $\dot{f}(y)$ of a generic field $f(y, p)$ in the domain transformation, taken at $p = 0$, is defined as:

$$ \dot{f}(y) = \lim_{\delta p \to 0} [f(y^{\delta p}, \delta p) - f(y, 0)] \delta p^{-1} = f_{,p}(y) + \nabla f(y) \cdot \theta(y) \quad (5) $$

and the material derivative of a generic surface integral is then given by the formula [4]:

$$ \frac{d}{dp} \int_S f \, dS = \int_S \{ \dot{f} + f \operatorname{div}_S \theta \} \, dS \quad (6) $$

<footnote>*E-mail address: bonnet@athena.polytechnique.fr</footnote>
where $div_S()$ denotes the surface divergence of a vector field. The above definitions still hold when the geometrical transformation (3) depends on a finite number $p_1, \dots, p_n$ of parameters.

### 3 BEM formulation for the energy release rate

The external boundary $\partial\Omega$ of the elastic body $\Omega$ is split into $S_T$ (prescribed traction $t = t^D$) and $S_u$ (prescribed displacement $u: u = u^D$). No body forces are present. A crack $\Gamma$, with upper and lower traction-free faces $\Gamma^\pm$ and unit normal $n$ directed from $\Gamma^-$ to $\Gamma^+$, is embedded in $\Omega$. In order to use only usual displacement BIE formulations, the so-called “multiregion approach” [7] is considered: $\Omega$ is split into two subdomains $\Omega^+$, $\Omega^-$ separated by a surface $S$ containing the crack $\Gamma$ (figure 1): The cracked solid equilibrium is then formulated in terms of two elastic problems on $\Omega^+$, $\Omega^-$ coupled by perfect bonding conditions on $S - \Gamma$. The energy release rate $G$ associated with the cracked solid $\Omega$ and the loading $u^D, t^D$ is defined by (1), or, equivalently by:

$$ \int_{\partial\Gamma} G(\theta, \nu) \, ds = -\dot{W} \quad \forall \theta \in \Theta \qquad (7) $$

where $\nu$ is the unit normal to the crack front $\partial\Gamma$ exterior to $\Gamma$ and tangent to $\Gamma$. Also, $\Theta$ denotes the set of virtual crack extensions, that is, those transformation velocities $\theta$ associated with geometrical transformations $\Phi(\cdot; p)$ which describe a crack extension: one has

$$ \theta.n = 0 \quad \text{on } \Gamma \qquad \theta = 0 \quad \text{on } S_u, S_T \qquad (8) $$

Thus only regular virtual crack extensions (i.e. without kinking) are considered. Moreover, in eq. (7) the variation $\dot{W}$ of $W$ is taken for constant loading ($u^D, t^D$) so that one has

$$ \dot{u} = 0 \quad (\text{on } S_u^{\pm}) \qquad \dot{t} = 0 \quad (\text{on } S_T^{\pm}) \qquad \dot{t} = 0 \quad (\text{on } \Gamma^{\pm}) \qquad (9) $$

Now, for linear elastic problems, one has the well-known boundary-only expression of the potential energy at equilibrium $W$:

$$ W = \frac{1}{2} \int_{S_u} t \cdot u^D dS - \frac{1}{2} \int_{S_T} t^D \cdot u dS $$

The variation $\dot{W}$ of $W$ in a crack extension thus stems from application of formula (6) to the above equation, and is expressed in terms of the material derivatives ($\dot{u}, \dot{t}$) as follows:

$$ \dot{W} = \frac{1}{2} \int_{S_u} \dot{t} \cdot u^D dS - \frac{1}{2} \int_{S_T} t^D \cdot \dot{u} dS \qquad (10) $$

The primary elastic variables ($u, t$) and their material derivatives ($\dot{u}, \dot{t}$) are respectively governed by the following BIEs [8]:

$$ \int_{\partial\Omega} [u_i(y) - u_i(x)] \Sigma_{ij}^k(x, y) n_j(y) dS_y - \int_{\partial\Omega} t_i(y) U_i^k(x, y) dS_y = 0 \qquad (11) $$

$$ \int_{\partial\Omega} [\dot{u}_i(y) - \dot{u}_i(x)] \Sigma_{ij}^k(x, y) n_j(y) dS_y - \int_{\partial\Omega} \dot{t}_i(y) U_i^k(x, y) dS_y = \int_{\partial\Omega} [\theta_{\ell}(y) - \theta_{\ell}(x)] \{ U_{i,\ell}^k(x, y) - D_{\ell j} u_i(y) \Sigma_{ij}^k(x, y) \} + U_i^k(x, y) \mathrm{div}_S \theta(y) dS_y \qquad (12) $$

in terms of the Kelvin fundamental displacement $U_i^k$ and stress $\Sigma_{ij}^k$ and using the tangential differential operator $D_{\ell j}() = n_{\ell}(), j - n_j(), \ell$.

A BE discretization of (11) for each boundary $\partial\Omega^+$, $\partial\Omega^-$, using 8- or 9-noded surface elements, is set up. The three basic steps involved in the computation of $G$ using the present approach are:

1. SOLUTION OF THE PRIMARY BIE. The coupled elastostatic problems on $\Omega = \Omega^+ \cup \Omega^-$ are numerically solved: the usual BEM linear system

$$ [A] \{u\} + [B] \{t\} = \{0\} $$

is built and, after appropriate column switches, the governing system on the vector $\{v\}$ of elastostatic unknowns has the form:

$$ [K] \{v\} = \{g^0\} \qquad (13) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Left: multiregion modelling of cracked solids; right: crack surface (shaded area: $E(\partial\Gamma)$).</figcaption></figure>

## 2 - SOLUTION OF THE DERIVATIVE BIE.
This step involves the construction of a discrete set of admissible transformation velocity fields $\theta \in \Theta$. Denote by $E(\partial\Gamma)$ the set of boundary elements adjacent to the crack front $\partial\Gamma$ and let $A^1, \dots, A^{NC}$ be the NC mesh nodes located on $\partial\Gamma$ (figure 1). The local numbering of nodes on each element along $\partial\Gamma$ is arranged so that the curve $(\xi_2 = -1)$, associated with the nodes 1, 2, 3, is located on $\partial\Gamma$ ($\xi = (\xi_1, \xi_2)$ being the antecedent of $y$ in the parent element). In order to take into account the known fact that $\dot{W}$ ultimately depends only on the normal extension velocity ($\theta.\nu$) of the crack front, transformation velocities of the following form are introduced:

$$ \theta(y) = \sum_{k=1}^{NC} \theta_k B^k(\xi) \quad (14) $$

in terms of NC scalar nodal values $\theta_k = (\theta.\nu)(A^k)$ and vector interpolation functions $B^k$. The latter are built so that $\theta(y) = 0$ outside $E(\partial\Gamma)$ and, on any element $E_e \in E(\partial\Gamma)$:

$$ B^k(\xi_1, 1) = 0 \quad (\xi_1 \in [-1, 1]) \quad \text{and} \quad B^k(\eta_1^l, -1) = \delta_{kl}\nu(\eta_1^l, -1) \quad (l = 1, 2, 3) \quad (15) $$

where $\eta^l$ is the antecedent of the crack front node in the parent element $A^l \in E_e \subset E(\partial\Gamma)$ and the local numbering $k = 1, 2, 3$ is used. The definition of $B^k(\xi)$ (in local numbering) uses a continuation $\nu(\xi_1, \xi_2)$ of the unit normal $\nu(\xi_1, -1)$ to $\partial\Gamma$:

$$ B^k(\xi) = f(\xi_2)S_k(\xi_1)\nu(\xi) \qquad \nu(\xi) = \frac{1}{|a_1(\xi)|}(a_1(\xi) \wedge n(\xi)) \quad (16) $$

where $a_1(\xi) = y, \xi_1$ and $S_1, S_2, S_3$ are the classical one-dimensional quadratic shape functions; $f(\xi) = (3 - 2\xi - \xi^2)/4$ with quarter-node elements (allowing for a linear variation of $f$ in the physical space) or $f(\xi) = (1 - \xi)/2$ with ordinary elements. The interpolation of $(\theta.\nu)$ on $\Gamma$ takes the form of a standard one-dimensional interpolation:

$$ (\theta.\nu)(y) = \sum_{k=1}^{NC} \theta_k S_k(\xi_1) \quad (17) $$

The definition (16) is then substituted into the derivative BIE (12). Due to the linearity of the right hand side of (12) with respect to $\theta$, one has

$$ \dot{u} = \sum_{k=1}^{NC} \theta_k \dot{U}^k \qquad \dot{t} = \sum_{k=1}^{NC} \theta_k \dot{T}^k \quad (18) $$

where the pair $(\dot{U}^k, \dot{T}^k)$ satisfies the matrix equation:

$$ [A] \{\dot{U}^k\} + [B] \{\dot{T}^k\} = \{f^1(u, t; B^k)\} $$
where the right-hand side {$f^1$}(u, t; B^k) comes from the discretization of the right-hand side of eq. (12) with $\theta = B^k$. The above equation, together with homogeneous boundary conditions, leads to the governing linear systems of equations for the vector {$\dot{v}^k$} of unknown derivatives:

$$[K]\{\dot{v}^k\} = \{f^1(u, t; B^k)\} \qquad (19)$$

The same matrix [K] appears in (13) and (19), because the present construction of $\theta$ is such that the Dirichlet and Neumann parts $S_u^\pm, S_T^\pm$ remain fixed and thus are material surfaces.

### 3 - SOLUTION OF THE GOVERNING VARIATIONAL EQUATION FOR G.
The energy release rate is interpolated, along $\partial\Gamma$ using the quadratic shape functions $S_k$ and nodal values $G_k = G(A^k)$. Then eq. (18) is substituted into (10) so that a linear matrix equation for the NC unknowns $G_k$ is readily obtained from the following discretized form of the variational equation (7):

$$\text{find } G_k \quad \forall \theta_m \ (1 \le m \le \text{NC}) \\ \theta_m = \left\{ G_k \int_{\Gamma} S_k(s)S_m(s) \, ds + \frac{1}{2} \int_{S_u} \dot{T}^m \cdot \mathbf{u}^D \, dS - \frac{1}{2} \int_{S_T} t^D \cdot \dot{U}^m \, dS \right\} = 0 \qquad (20)$$

## 4 Numerical examples

**EXAMPLE 1 - ROUND BAR WITH A PENNY-SHAPED AXIAL CRACK.** An internal penny-shaped plane crack of radius $R_1$ is situated in the symmetry plane $y_3 = 0$ of a cylindrical bar (axis $Oy_3$, length $2H$, external radius $R > R_1$) subjected to a uniform tension $p$ along the axialdirection (mode I). Two different meshes (figure 2) were used for one-eighth of the structure. The numerical values obtained for $G$ at the crack front nodes (figure 2) show good agreement with the known semi-analytical solution for $K_I$ (Tada, Paris & Irwin [9]) reference solution (figure 2).

**EXAMPLE 2 - SEMI-ELLiptICAL SURFACE CRACK.** A semi-elliptical surface crack (see figure 3 for the geometrical notations) is situated in the symmetry plane $y_3 = 0$ of a rectangular parallelepiped subjected to a uniform tension $p$ along the axialdirection (mode I). Owing to geometrical symmetry, only one-quarter of the boundary is discretized (136 9-noded boundary elements); the crack front itself supports 13 nodal values of $G$ and $\theta$. Two variant meshes M1 and M2 (figure 3) were used, with respectively uniform ($\Delta\theta = \pi/24$, mesh M1) and non-uniform ($\Delta\theta = \pi/32$ (resp. $\pi/16$) for $\theta \in [0, \pi/4]$ (resp. $\theta \in [\pi/4, \pi/2]$), mesh M2) angular spacing between crack front nodes, the crack edge being located at $\theta = 0$.

Numerical values of the non-dimensional SIF $K_I^*Q/(p\sqrt{\pi a})$ were obtained from the values of $G$ computed with the present method using (2), with $Q = E(k)$ and $k^2 = 1 - (b/a)^2$ ($b \le a$) or $k^2 = 1 - (a/b)^2$ ($a \le b$); $E(k)$ denotes the complete elliptic integral of the second kind. They are compared (figure 4) to other numerical results from Newman and Raju [10] and Tanaka and Itoh [11]. The latter were obtained using a sophisticated special crack-front element which allows for the modelling of both the square-root crack front singularity and the crack edge singularity (whose exponent differs from $-1/2$ except for $\nu = 0$), and are thus expected to provide the better reference solution near the crack edge. Our results show an improvement in reproducing the small peak near the crack edge when mesh M2, which is finer than M1 near the crack edge, is used; generally speaking, they agree reasonably well with the reference ones. For comparison sake, figure 4 also shows the values of $K_I^*$ obtained by extrapolation of the crack opening displacement, which tend to be somewhat poorer than those obtained using the present approach, despite the use of quarter-node elements along the crack front.

## References

[1] DELORENZI H.G. - On the energy release rate and the J-integral for 3-D crack configurations. *Int. J. Fract.*, **19**, pp. 183-194, 1982.

[2] DESTUYNDER PH., DJAOUA M., LESCURE S. - Quelques remarques sur la mécanique de la rupture élastique. *J. Mécan. Théor. Appl.*, **2**, pp. 113-135, 1983.
[3] WADIER Y., MALAK O. — The theta method applied to the analysis of 3D-elastic-plastic cracked bodies. SMIRT, Los Angeles 1989.

[4] PETRYK H., MROZ Z. — Time derivatives of integrals and functionals defined on varying volume and surface domains. Arch. Mech., **38**, pp. 694-724, 1986.

[5] HELLEN T.K. — On the method of virtual crack extension. Int. J. Num. Meth. in Eng. 9, pp. 187-207, 1975.

[6] NGUYEN Q.S. — Bifurcation and stability in dissipative media (plasticity, friction, fracture). Appl. Mech. Rev., **37**, pp. 1-31, 1994.

[7] CRUSE T.A. — *Boundary Element Analysis in Computational Fracture Mechanics*. Kluwer Academic Publishers, 1988.

[8] BONNET M. — Regularized BIE formulations for first- and second-order shape sensitivity of elastic fields. *Special issue of Computers and Structures*, (S. Saigal, guest editor.), to appear, 1995.

[9] TADA, PARIS & IRWIN — *The stress analysis of crack handbook*. Del. Research Corporation, Hellertown, Pennsylvania, USA (1973).

[10] NEWMAN J.C., RAJU I.S. — An empirical stress-intensity factor equation for the surface crack. Engng. Frac. Mech., **15**, pp. 185-192, 1981.

[11] TANAKA M., ITOH H. — New crack elements for boundary element analysis of elastostatics considering arbitrary stress singularities. Appl. Math. Modelling, **11**, pp. 357-363, 1987.

<figure><img src="image_4.png" /><figcaption>Figure 2: Example 1: comparison between numerical results and reference solution (upper: coarse mesh, lower: fine mesh)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Example 2: Geometrical notation for the crack plane (upper left), boundary element subdivision of crack plane (lower left: mesh M2, $b = 0.6a$, right: mesh M1, $b = a$)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Example 2: comparison between results for $b = 0.4a$ (upper left), $b = 0.6a$ (upper right), $b = a$ (lower left), $b = 2a$ (lower right).</figcaption></figure>