# Boundary element analysis of an integral equation for three-dimensional crack problems

Anh Le Van, Bernard Peseux

► To cite this version:

Anh Le Van, Bernard Peseux. Boundary element analysis of an integral equation for three-dimensional
crack problems. International Journal for Numerical Methods in Engineering, 1988, 26 (11), pp.2383
- 2402. 10.1002/nme.1620261104. hal-01005293

HAL Id: hal-01005293

https://hal.science/hal-01005293

Submitted on 28 Oct 2016

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>
BOUNDARY ELEMENT ANALYSIS OF AN INTEGRAL
EQUATION FOR THREE-DMENSIONAL CRACK
PROBLEMS

A. LE VAN* AND B. PESEUX

Laboratoire de Mécanique des Structures, Ecole Nationale Superieure de Mecanique, 1, rue de la Noe, Nantes 44072 Cedex,
France

In another paper, the authors proposed an integral equation for arbitrary shaped three-dimensional cracks. In the present paper, a discretization of this equation using a tensor formalism is formulated. This approach has the advantage of providing the displacement discontinuity vector in the local basis which varies as a function of the point of the crack surface. This also facilitates the computation of the stress intensity factors along the crack edge. Numerical examples reported for a circular crack and a semi-elliptical surface crack in a cylindrical bar show that one can obtain good results, using few Gaussian points and no singular elements.

INTRODUCTION

The integral equation method has been thoroughly experienced in structure analysis. Integral
equations derived from the Somigliana representation were successfully applied to three-
dimensional problems, such as in References 1 and 2. In the latter reference, a particular
symmetrical problem of cracked bodies was also investigated. Integral equations with kernels
containing singular solutions, usually referred to as Kupradze elastic potentials,³ are particularly
well suited to crack analysis. The use of these potentials in linear fracture mechanics is
advantageous as it permits one to obtain directly stress intensity factors from the computed vector
density. Embedded cracks in an infinite medium were studied in References 4–7. In particular,
Bui⁴ showed that for plane cracks the mode I is entirely uncoupled from modes II and III. Later
on, integral equations for three-dimensional cracks were proposed;⁸,⁹ in the last reference use was
made of a Kupradze double layer potential, and the unknown was the vector density directly
related to the displacement discontinuity through the crack surface. In this paper, we carry out the
discretization of the integral equation proposed in Reference 9, with a view to studying imbedded
or surface crack problems.

For numerical purposes, the main difficulty is the reckoning of singular integrals defined in the
sense of the principal value. Cruse¹ succeeded in giving the closed form of the principal value
integral, using planar boundary elements. A more general way to evaluate two-dimensional
singular integrals was investigated by Kazantzakis and Theocaris,¹⁰ who also reduced these to
one-dimensional finite-part integrals.¹¹ Among the most recent works, one can find References 12
and 13 where exact expressions for some specific two-dimensional finite-part integrals are derived.
As for Kupradze potentials, Bui⁴ took a small polygon instead of a circle to perform principal value integrals, it seems, however, that this procedure cannot be used with ease. Lastly, one can find in Reference 14 an improvement of the evaluation of integrals of order 1/r, 1/r², 1/r³, by means of a minimum number of Gaussian points. As the considered distance r, small though it may be, remains always finite, these integrals are actually not singular ones. In any case, difficulties must certainly be expected when dealing with principal value integrals defined on curved surfaces.

In our analysis, we limit consideration to classical Gaussian quadrature, with no actual principal value computations. Nevertheless, two numerical examples given in the last section show that one can obtain fair results, using few Gaussian points, no singular integral computations and in particular, no singular elements around the crack edge. First, we deal with the embedded penny-shaped crack. The numerical results obtained proved to fit well to the analytical solution of Sneddon.¹⁵ Then we treat the problem of a semi-elliptical crack in a cylindrical bar under combined loads. In the opening mode, the numerical results can be compared to the existing results¹⁶⁻¹⁸ which involved analogous crack geometries.

## THEORETICAL BACKGROUND

Let us consider the problem of an elastic solid $\mathcal{D}$ containing an arbitrary shaped, imbedded or surface crack $S$, with the mixed boundary conditions:

$$y_0 \in S \cup D: \quad \lim_{\mathcal{D} \setminus S \ni x \to y_0^\pm, \mathbf{n}_x = \mathbf{n}_{y_0}} \mathbf{t}(x, \mathbf{n}_x) = \mathbf{t}^g(y_0, \mathbf{n}_{y_0}) \quad (1a)$$

$$\text{or: } y_0 \in S \cup D: \quad \lim_{\mathcal{D} \setminus S \ni x \to y_0^\pm} \mathbf{u}(x) = \mathbf{u}^g(y_0) \quad (1b)$$

where the superscript g stands for given, and the double sign $\pm$ is related to the surface orientation, locally defined by the normal $\mathbf{n}_{y_0}$ at each point $y_0$.

### 1. Displacement field expression

Let the displacement field be expressed by means of the Kupradze double-layer potential of the first kind³

$$\forall x \in E_3 \setminus S: \quad \mathbf{u}(x) = \int_S T^T(y-x, \mathbf{n}_y) \varphi(y) d_y S \quad (2)$$

where the density $\varphi$ is a vector function defined in the crack $S$, and $T$ the Kupradze tensor

$$T(y-x, \mathbf{n}_y) = \frac{2\mu}{8\pi(\lambda+2\mu)\cdot r^2} [\mathbf{n}_y \otimes \mathbf{e}_r - \mathbf{e}_r \otimes \mathbf{n}_y - (\mathbf{e}_r \cdot \mathbf{n}_y)I] - \frac{6(\lambda+\mu)\mathbf{e}_r \cdot \mathbf{n}_y}{8\pi(\lambda+2\mu)\cdot r^2} \mathbf{e}_r \otimes \mathbf{e}_r$$

$I$ is the unit tensor, $r = \|y-x\|$, $\mathbf{e}_r = (y-x)/\|y-x\|$.

The displacement discontinuity on the crack is obtained from

$$\forall y_0 \in S, \quad [\mathbf{u}(y_0)] \equiv \mathbf{u}(y_0^+) - \mathbf{u}(y_0^-) = \varphi(y_0) \quad (3)$$

Using (3), one can calculate the stress intensity factors which are now directly related to the density $\varphi$.
## 2. Integral equations

It was shown in Reference 9 that, under the assumptions $S \in C^{1,\alpha}$, $0 < \alpha \le 1$, $\phi \in C^{1,\beta}(S)$, $0 < \beta \le 1$, the integral equation for a general crack is expressed by

$$ 
\begin{aligned}
\forall y_0 \in S, \quad \mathbf{t}(\mathbf{y}_0, \mathbf{n}_{y0}) = & \frac{E}{16\pi(1-v^2)} \mathrm{pv} \int_{\Delta} \left(2(\phi_{,u}, \mathbf{e}_r, \mathbf{F}_{,v})\mathbf{n}_{y0} - (1-2v)((\phi_{,u}, \mathbf{e}_r, \mathbf{n}_{y0})\mathbf{F}_{,v} \right. \\
& \left. \qquad + (\mathbf{F}_{,v} \cdot \mathbf{n}_{y0})\phi_{,u} \wedge \mathbf{e}_r \right) + 3(\mathbf{e}_r \cdot \phi_{,u})((\mathbf{F}_{,v}, \mathbf{n}_{y0}, \mathbf{e}_r)\mathbf{e}_r + (\mathbf{n}_{y0} \cdot \mathbf{e}_r)\mathbf{e}_r \wedge \mathbf{F}_{,v}) \\
& - 2(\phi_{,v}, \mathbf{e}_r, \mathbf{F}_{,u})\mathbf{n}_{y0} + (1-2v)((\phi_{,v}, \mathbf{e}_r, \mathbf{n}_{y0})\mathbf{F}_{,u} + (\mathbf{F}_{,u} \cdot \mathbf{n}_{y0})\phi_{,v} \wedge \mathbf{e}_r) \\
& - 3(\mathbf{e}_r \cdot \phi_{,v})((\mathbf{F}_{,u}, \mathbf{n}_{y0}, \mathbf{e}_r)\mathbf{e}_r + (\mathbf{n}_{y0} \cdot \mathbf{e}_r)\mathbf{e}_r \wedge \mathbf{F}_{,u}) \mathrm{d}u \mathrm{d}v \\
& \qquad \qquad \qquad \qquad \qquad \qquad - \frac{E}{16\pi(1-v^2)} \int_{\partial S} \mathcal{C}(\mathbf{n}_{y0}, y-y_0, \tau)\phi(y)\mathrm{d}y 
\end{aligned} 
\quad (4) $$

with now (Figure 1): $r = \|y - y_0\|$, $\mathbf{e}_r = (y - y_0)/\|y - y_0\|$

$F$: a parametrization of $S$, defined on a domain $\Delta$ in $\mathbb{R}^2: \Delta \ni (u, v) \to y = F(u, v) \in S$

$\phi(u, v) \equiv \phi(F(u, v)) = \phi(y)$

$\tau$: the unit vector tangent to the crack edge $\partial S$, oriented according to the orientation of $S$.

One should notice that the symbol $\mathrm{pv}$ before the two-dimensional integral relates to the principal value on $S$ (not on $\Delta$), and that the expression still holds in the general case when $\mathbf{F}_{,u}$ is not

<figure><img src="image_2.png" /><figcaption>Figure 1. General crack configuration with respectively a particular and an ordinary point, $y_0$ and $y$</figcaption></figure>
perpendicular to $F, v$. The kernel of the line integral has been expressed for brevity as the product of a linear operator $\mathcal{C}$ and the density $\phi$. Such line integrals appear when surface crack problems are involved.

In the next section, we shall concern ourselves with the discretization of the two-dimensional integral. As previously mentioned, this discretization will not include the study of the principal value.

# DISCRETIZATION OF THE INTEGRAL EQUATION

## 1. Shape functions

The co-ordinates transformation is given by

$$ \tau: \xi \to (u, v) = (\langle N(\xi) \rangle \{u\}^e, \langle N(\xi) \rangle \{v\}^e) \quad (5) $$

where $\xi$ is a point of the reference element and the superscript e relates to an ordinary element. In this paper, we are interested to the most commonly encountered co-ordinate systems: the Cartesian, the polar and the cylindrical ones. In Cartesian co-ordinates $(u, v)=(y_1, y_2)$, in polar co-ordinates $(u, v)=(\rho, \theta)$ (Figure 2) and in cylindrical co-ordinates, with $\rho$ constant, $(u, v)=(\theta, y_3)$ (Figure 3).

## 2. Interpolation of $\phi$

In practice, it is advantageous to obtain the components of the density $\phi$ in the local basis which varies with the point $y$ of $S$. For this reason, we attach much importance to the choice of the co-ordinates system and we will consider different ways to interpolate the components of $\phi$ in the local basis. Use will be made of isoparametric elements.

In Cartesian co-ordinates, for instance, the local basis $b$ is invariant; it coincides with the global, fixed basis $e$ of the space $E_3$:

$$ b \equiv (b_1; b_2; b_3) = (e_1; e_2; e_3) $$

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

a)

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

b)

<caption>Figure 2. Interpolation with polar co-ordinates: (a) local co-ordinates $\xi=(\xi_1, \xi_2)$; (b) global coordinates $(\rho, \theta)$</caption>
<figure><img src="image_5.png" /><figcaption>Figure 3. Interpolation with cylindrical co-ordinates: (a) local co-ordinates $\xi=(\xi_1, \xi_2)$; (b) global co-ordinates $(\theta, y_3)$</figcaption></figure>

The density $\phi$ is considered in the fixed basis:

$$ \Phi = \phi_1 e_1 + \phi_2 e_2 + \phi_3 e_3 $$

Each of the components of $\Phi$ is then interpolated classically by

$$ \phi_i = \langle N(\xi) \rangle \{\phi_i\}^e, \quad i \in \{1, 2, 3\} \qquad (6a) $$

In cylindrical co-ordinates, the local basis varies with the point y:

$$ b = (b_1; b_2; b_3) = (e_\rho; e_\theta; e_3) $$

The density will be considered in the variable basis:

$$ \Phi = \phi_{\rho} e_{\rho} + \phi_{\theta} e_{\theta} + \phi_{3} e_{3} $$

and its components are now interpolated by

$$ \phi_i = \langle N(\xi) \rangle \{\phi_i\}^e, \quad i \in \{\rho, \theta, 3\} \qquad (6b) $$

Relations (6) in turn yield

$$ b\{\phi(\xi)\} = b[\mathcal{N}^{c_n}(\xi)] \cdot c_n\{\phi\}^e \qquad (7) $$

where $c_n$ denotes the canonical basis of $\mathbb{R}^n$. The letters $b, c_n$ added in small characters recall the bases under consideration. In the case of four-node elements, the matrix $[\mathcal{N}]$ can be explicitly as
follows:

$$b[\mathcal{N}^{c_{12}}(\xi)] = \begin{bmatrix} N_1 & 0 & 0 & | & N_2 & 0 & 0 & | & N_3 & 0 & 0 & | & N_4 & 0 & 0 \\ 0 & N_1 & 0 & | & 0 & N_2 & 0 & | & 0 & N_3 & 0 & | & 0 & N_4 & 0 \\ 0 & 0 & N_1 & | & 0 & 0 & N_2 & | & 0 & 0 & N_3 & | & 0 & 0 & N_4 \end{bmatrix}$$

and the density $\phi$ in the case of cylindrical co-ordinates can be expressed by

$$\langle \phi \rangle^e = \langle (\phi_\rho, \phi_\theta, \phi_3)^1, (\phi_\rho, \phi_\theta, \phi_3)^2, (\phi_\rho, \phi_\theta, \phi_3)^3, (\phi_\rho, \phi_\theta, \phi_3)^4 \rangle^e$$

One can then express the derivatives of $\phi$ with respect to the parameters $u$ and $v$:

$$\phi_{,u} = \mathcal{M}_u \cdot \phi^e \quad (8a)$$

$$\phi_{,v} = \mathcal{M}_v \cdot \phi^e \quad (8b)$$

where

$$\mathcal{M}_u = \begin{cases} \mathcal{N}_{,u} & \text{if } u \neq \theta \\ \mathcal{N}_{,u} + \mathcal{N}_1 & \text{if } u = \theta \end{cases}; \quad \mathcal{M}_v = \begin{cases} \mathcal{N}_{,v} & \text{if } v \neq \theta \\ \mathcal{N}_{,v} + \mathcal{N}_1 & \text{if } v = \theta \end{cases}$$

The matrical representation of $\mathcal{N}_1$ in the bases $b$ and $c_{12}$ (case of four-node elements) is

$$b[\mathcal{N}_1^{c_{12}}(\xi)] = \begin{bmatrix} 0 & -N_1 & 0 & | & 0 & -N_2 & 0 & | & 0 & -N_3 & 0 & | & 0 & -N_4 & 0 \\ N_1 & 0 & 0 & | & N_2 & 0 & 0 & | & N_3 & 0 & 0 & | & N_4 & 0 & 0 \\ 0 & 0 & 0 & | & 0 & 0 & 0 & | & 0 & 0 & 0 & | & 0 & 0 & 0 \end{bmatrix} \quad (9)$$

### 3. Discretization of the integral equation

We study next how to transform each of the terms of the right-hand side of (4). First, let us consider the terms containing $\phi_{,u}$. We have

$$\begin{align}
(\phi_{,u}, \mathbf{e}_r, \mathbf{F}_{,v}) \mathbf{n}_{y0} &= [\phi_{,u} \cdot (\mathbf{e}_r \wedge \mathbf{F}_{,v})] \mathbf{n}_{y0} \nonumber \\
&= [(\mathcal{M}_u \phi^e) \cdot (\mathbf{e}_r \wedge \mathbf{F}_{,v})] \mathbf{n}_{y0} \quad (\text{using (8a)}) \nonumber \\
&= [\mathbf{n}_{y0} \otimes \mathcal{M}_u^T (\mathbf{e}_r \wedge \mathbf{F}_{,v})] \phi^e \tag{10a}
\end{align}$$

Similarly

$$(\phi_{,u}, \mathbf{e}_r, \mathbf{n}_{y0}) \mathbf{F}_{,v} = [\mathbf{F}_{,v} \otimes \mathcal{M}_u^T (\mathbf{e}_r \wedge \mathbf{n}_{y0})] \phi^e \quad (10b)$$

As for the term $\phi_{,u} \wedge \mathbf{e}_r$, we have

$$\phi_{,u} \wedge \mathbf{e}_r = (\phi_{,u}, \mathbf{e}_r, \mathbf{b}_{0i}) \mathbf{b}_{0i} \quad (11)$$

where an implicit sum is implied on the repeated subscript $i \in \{1, 2, 3\}$ or $\{\rho, \theta, 3\}$, $b_0 = (\mathbf{b}_{0i})$ is the basis related to $y_0$. It should be noted that the components of the given stress $\mathbf{t}(y_0, \mathbf{n}_{y0})$ in the left-hand side of (4) are obtained in the basis $b_0$. Furthermore, in interacting or surface crack problems, $b_0$ may differ or not from $b$, according to the relative position of the particular point $y_0$ and the ordinary one $y$. Now in virtue of (11):

$$(\mathbf{F}_{,v} \cdot \mathbf{n}_{y0}) \phi_{,u} \wedge \mathbf{e}_r = (\mathbf{F}_{,v} \cdot \mathbf{n}_{y0}) [\mathbf{b}_{0i} \otimes \mathcal{M}_u^T (\mathbf{e}_r \wedge \mathbf{b}_{0i})] \phi^e \quad (10c)$$
In the same manner, the fourth terms of (4) can be expressed by

$$ (e_r \cdot \phi_{,u})[(F_{,v}, n_{y0}, e_r)e_r + (n_{y0} \cdot e_r)e_r \wedge F_{,v}] \\ = [((F_{,v}, n_{y0}, e_r)e_r + (n_{y0} \cdot e_r)e_r \wedge F_{,v}) \otimes \mathcal{M}_u^T e_r] \phi^e \quad (10d) $$

Eventually, relations (10) are expressed in matrical form by

$$ (\phi_{,u}, e_r, F_{,v}) n_{y0} = [\{n_{y0}\} (\langle e_r \wedge F_{,v} \rangle [\mathcal{M}_u])] \{\phi\}^e \quad (12a) $$

$$ (\phi_{,u}, e_r, n_{y0}) F_{,v} = [\{F_{,v}\} (\langle e_r \wedge n_{y0} \rangle [\mathcal{M}_u])] \{\phi\}^e \quad (12b) $$

$$ (F_{,v} \cdot n_{y0}) \phi_{,u} \wedge e_r = (F_{,v} \cdot n_{y0}) [\{b_{0i}\} (\langle e_r \wedge b_{0i} \rangle [\mathcal{M}_u])] \{\phi\}^e \quad (12c) $$

$$ (e_r \cdot \phi_{,u})((F_{,v}, n_{y0}, e_r)e_r + (n_{y0} \cdot e_r)e_r \wedge F_{,v}) \\ = [\{(F_{,v}, n_{y0}, e_r)e_r + (n_{y0} \cdot e_r)e_r \wedge F_{,v}\} (\langle e_r \rangle [\mathcal{M}_u])] \{\phi\}^e \quad (12d) $$

It should be noted that, in (12d) for instance, as the given stress in the left-hand side is expressed in the basis $b_0$, the vector $e_r$ between braces $\{ \}$ (column-matrix, see Appendix) must be expressed in the basis $b_0$, whereas the same $e_r$ between angle brackets $\langle \rangle$ (row-matrix) must be expressed in the basis $b$.

The terms containing $\phi_{,v}$ are discretized in the same manner; they yield relations analogous to (12) where $\phi_{,v}$ takes the place of $\phi_{,u}$, and $F_{,u}$ of $F_{,v}$.

Thus, for each particular point $y_0$ of the crack, we have discretized the integral equation using a tensor formalism, equations (10). This equation being discretized as a *vectorial* equation, and not as three separate scalar equations, this formalism allows us a compact programmation for numerical purposes. One has to solve an algebraic system in the form

$$ [A] \{\phi\} = \{B\} $$

where $[A]$ is a full, non-symmetrical matrix, $\{\phi\}$ the unknown vector and $\{B\}$ contains the components of the given stresses. In interacting or surface crack problems, $\{\phi\}$ contains different vector densities corresponding respectively to different cracks or boundaries of the solid $\mathcal{D}$.

N.B.: From the dimensional point of view, the matrices appearing in equations (12) (terms between brackets $[ ]$) are expressed in the inverse of a length or are dimensionless according as Cartesian, polar or cylindrical co-ordinates are involved.

## NUMERICAL PARTS

### 1. Embedded circular crack

First, the program has been tested on the simple case of a circular crack centred at O, with radius $a$, under a uniform pressure $p$, the analytical solution of which is well known:¹⁵

$$ \Phi(\rho) = \frac{8(1 - v^2) \cdot pa}{\pi E} \sqrt{1 - (\rho/a)^2} \cdot e_3, \quad \rho = \|Oy\| $$

As for any plane crack imbedded in the infinite medium, the three modes are uncoupled. The stress intensity factor in the opening mode is given by

$$ K_1 = 2p \sqrt{a/\pi} $$

The natural choice of parameters for the penny-shaped crack is $(u, v)=(\rho, \theta)$. One should notice that this choice yields a co-ordinates transformation which is well suited to the crack geometry:
the crack edge is especially better shaped than with the usual grid pattern, as shown in Figure 4.

Moreover, this choice of parameters facilitates the computation of the stress intensity factors which require to be taken as a limit when moving along a direction *perpendicular* to the crack edge. The given stress vector is expressed at each point $y_0$ in a local basis $b_0 = (e_\rho 0; e_\theta 0; e_3)$, whereas the vector density $\phi$ at a point $y$ in $b = (e_\rho; e_\theta; e_3)$; this situation is shown in Figure 5.

Figure 6 represents a quarter of a mesh with 181 interior nodal points (7 nodes along the radius and 36 along the polar angle $\theta$, each sector being equal to $10^\circ$). The numerical integration was performed with 2 × 2 Gaussian points.

As expected, the numerical values of the radial and tangential displacement discontinuities $\phi_\rho$ and $\phi_\theta$ are zero; moreover, all the nodes of the same radius have identical $\phi_3$ values. It is interesting to note that, in mixed modes problems, the formulation herein proposed permits us to obtain directly the radial and tangential components of $\phi$ ($\phi_\rho$ and $\phi_\theta$) without additional calculations, because $\phi$ is expressed in the local polar basis. Table I and Figure 7 give respectively the values of $\phi_3$ normalized by the theoretical maximum value $(\phi_{3max})_{theor} = 8(1 - v^2)pa/\pi E$, and its curve compared to Sneddon's solution. It should be noticed that one obtains in fact

<figure><img src="image_6.png" /><figcaption>Figure 4. Circular crack: (a) use of polar co-ordinates $(u, v) = (\rho, \theta)$; (b) usual grid pattern, use of Cartesian co-ordinates</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Variation of the local basis between the points $y$ and $y_0$ on the crack</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. Circular crack: first mesh</figcaption></figure>

<table><caption>Table I. Circular crack. First mesh. $\phi_3$ values normalized by the theoretical maximum value</caption><thead><tr><th>ρ/a</th><th>(φ<sub>3</sub>)<sub>theor</sub>/(φ<sub>3 max</sub>)<sub>theor</sub> = √<span style="text-decoration:overline;">1 - (ρ/a)<sup>2</sup></span></th><th>φ<sub>3</sub>/(φ<sub>3 max</sub>)<sub>theor</sub></th></tr></thead><tbody><tr><td>0</td><td>1</td><td>0.91244</td></tr><tr><td>0.3</td><td>0.95393</td><td>0.87363</td></tr><tr><td>0.6</td><td>0.8</td><td>0.72500</td></tr><tr><td>0.8</td><td>0.6</td><td>0.56512</td></tr><tr><td>13/15</td><td>0.49889</td><td>0.47032</td></tr><tr><td>14/15</td><td>0.35901</td><td>0.35041</td></tr><tr><td>1</td><td>0</td><td>0</td></tr></tbody></table>

$E/(16\pi(1 - v^2)) \cdot \phi_3$ directly from the program, so that the Young's modulus $E$ datum is not necessary. The relative error on $\phi_3$ is maximum at the crack centre and equals 8.75 per cent; that of the computed stress intensity factor $K_1$ is 3.97 per cent.

Figure 8 represents a finer mesh with 529 interior nodes (13 along the radius and 48 along the polar angle $\theta$, each sector being equal to 7.5°). Use is always made of 2 × 2 Gaussian points. Table II and Figure 7 give the values of $\phi_3$ normalized by the theoretical maximum $\phi_3$-value. The relative error on $\phi_3$ at the centre is 2.97 per cent; that of the computed stress intensity factor $K_1$ is 0.25 per cent.

## 2. Semi-elliptical crack in a cylindrical bar under combined loads

Let us consider a cylindrical bar with radius $R$, containing a semi-elliptical crack with semi-axes $a$ and $b$. The crack occurs perpendicularly in the lateral surface of the cylinder, along a line which will be referred to as the surface-line of the crack (Figure 9). Three types of loads will be applied uniquely on the crack, the external surface of the bar being left free: a uniform pressure, a linear pressure along the $y_1$-axis, and a shear distribution corresponding to the torsion of the bar (Figure 10).

According to the discussion in Reference 9, the surface crack problem can be investigated by considering a two crack problem: the first crack is the studied one, the second is a cylinder-shaped crack high enough to simulate a free bar. We introduce the notations:
<figure><img src="image_9.png" /><figcaption>Figure 7. Circular crack: $\phi_3$ normalized by the theoretical maximum value: analytical results; — analytical results; — x — computed values with the first mesh; — ○ — computed values with the second mesh</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Circular crack: second mesh</figcaption></figure>

$S_{cr}$: the semi-elliptical crack, admitting the Cartesian parametrization:

$$F_{cr}: \Delta_{cr} \ni (y_1, y_2) \to y \in S_{cr}$$

$S_{cyl}$: the cylinder-shaped crack, admitting the parametrization:

$$F_{cyl}: \Delta_{cyl} = [0, 2\pi[x[-h, h] \ni (\theta, z)] \to y \in S_{cyl}$$

$L$: the surface-line, which is the intersection of $S_{cyl}$ and the closure of $S_{cr}$

$$S_{cr} \ni y \to \phi_{cr}(y) \in \mathbb{R}^3, \quad S_{cyl} \ni y \to \phi_{cyl}(y) \in \mathbb{R}^3,$$
<table><caption>Table II. Circular crack. Second mesh. $\phi_3$ values normalized by the theoretical maximum value</caption><thead><tr><th>ρ/a</th><th>(φ<sub>3</sub>)<sub>theor</sub>/(φ<sub>3max</sub>)<sub>theor</sub> = √<span style="text-decoration:overline;">1 - (ρ/a)<sup>2</sup></span></th><th>φ<sub>3</sub>/(φ<sub>3max</sub>)<sub>theor</sub></th></tr></thead><tbody><tr><td>0</td><td>1</td><td>0.97028</td></tr><tr><td>0.13</td><td>0.99151</td><td>0.96311</td></tr><tr><td>0.26</td><td>0.96561</td><td>0.93761</td></tr><tr><td>0.38</td><td>0.92499</td><td>0.89805</td></tr><tr><td>0.49</td><td>0.87172</td><td>0.84761</td></tr><tr><td>0.59</td><td>0.80740</td><td>0.78810</td></tr><tr><td>0.67</td><td>0.74236</td><td>0.72839</td></tr><tr><td>0.74</td><td>0.67261</td><td>0.66313</td></tr><tr><td>0.8</td><td>0.6</td><td>0.59344</td></tr><tr><td>0.85</td><td>0.52678</td><td>0.52047</td></tr><tr><td>0.9</td><td>0.43589</td><td>0.42999</td></tr><tr><td>0.95</td><td>0.31225</td><td>0.31752</td></tr><tr><td>1</td><td>0</td><td>0</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Figure 9. Semi-elliptical surface crack in a cylindrical bar. L the surface-line of the crack</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10. Three loading types on the surface crack: (a) uniform pressure; (b) linear pressure; (c) shear distribution resulting from the torsion of the bar ($\sigma_{max}$, $\tau_{max}$: maximum stress values)</figcaption></figure>

respectively the unknown density functions over $S_{cr}$ and $S_{cyl}$

$$ \Phi_{cr} = \Phi_{cr} \circ F_{cr}, \quad \Phi_{cyl} = \Phi_{cyl} \circ F_{cyl} $$

The edge of $S_{cyl}$ is the union of the upper and the lower crack edges of the cylinder, plus a vanishing, counterclockwise contour surrounding $L$. The set of equations of the problem can be written as follows:

$$ \forall y_0 \in S_{cr}, \qquad \mathbf{t}(y_0, \mathbf{n}_{y_0}) = vp \int_{\Delta_{cr}} (\text{kernel containing } \Phi_{cr}) du dv \\ \qquad \qquad + \int_{\Delta_{cyl}} (\text{kernel containing } \Phi_{cyl}) du dv \quad (13a) $$

$$ \forall y_0 \in S_{cyl}, \qquad \mathbf{t}(y_0, \mathbf{n}_{y_0}) = \mathbf{0} = \int_{\Delta_{cr}} (\text{kernel containing } \Phi_{cr}) du dv \\ \qquad \qquad + vp \int_{\Delta_{cyl}} (\text{kernel containing } \Phi_{cyl}) du dv \quad (13b) $$

$$ \forall \tilde{y} \in \partial S_{cr} \setminus L, \qquad \lim_{S_{cr} \ni y \to \tilde{y}} \Phi_{cr}(y) = \mathbf{0} \quad (13c) $$

$$ \forall \tilde{y} \in \partial S_{cyl} \setminus L, \qquad \lim_{S_{cyl} \ni y \to \tilde{y}} \Phi_{cyl}(y) = \mathbf{0} \quad (13d) $$

$$ \forall \tilde{y} \in L, \qquad \lim_{S_{cr} \ni y \to \tilde{y}} \Phi_{cr}(y) + \lim_{S_{cyl}^{sup} \ni y \to \tilde{y}} \Phi_{cyl}(y) - \lim_{S_{cyl}^{inf} \ni y \to \tilde{y}} \Phi_{cyl}(y) = \mathbf{0} \quad (13e) $$

where the kernels in (13a) and (13b), given by (4), are not made explicit because of their length; $S_{cyl}^{sup}$ and $S_{cyl}^{inf}$ denote respectively the upper and lower part of the cylindrical crack, separated by the section containing $S_{cr}$. The relation (13e) holds if we choose the normal to $S_{cr}$ upward with respect to the crack and the normal to $S_{cyl}$ outward with respect to the cylinder. With another choice of surfaces orientations, however, this relation still holds but the signs before each term may change. Furthermore, one can easily verify that, in various surface crack problems, one or more relations of type (13e) may be obtained. Relations of type (13e) allow one to verify that, though every contour integral is different from zero, there always appears a set of these which cancel one another.
For numerical purposes, a crack depth-radius ratio $b/R=0.4$ and a semi-axes ratio $b/a=5/7 \approx 0.714$ are considered. The cylinder lateral surface is divided into 748 elements with 731 interior nodes, while the semi-elliptical crack contains but 58 elements and 60 interior nodes (Figure 11). So large a node number on the free surface is due to the use of four-node elements with their sides along the parallel of the cylinder or along its axis. It should be wiser to use curved cylindrical elements (i.e. elements with distinct cylindrical nodal values) so as to realize a relatively refined mesh around the surface-line together with a coarse mesh far away, and to reduce notably the node number on the free surface. We have, however, not chosen such elements because of the difficulties caused to automatic mesh generations. In fact, the sole aim of this paper is to show the feasibility and the validity of the proposed approach. Anyway, the interesting feature is the little node number on the crack itself.

Like the previous numerical example, use is made of four-node elements and the numerical quadrature is performed with $2 \times 2$ Gaussian points. Tables III, IV and V give the $\phi_i$-values, $i \in \{1, 2, 3\}$, resulting from three loading types. From Tables III and IV, one can note that, under uniform or linear pressure loadings, the crack is seen to be in the opening mode. Table V shows on the contrary that, under torsion loadings, the crack deforms in a mixed mode II+III. With the same maximum stress values (i.e. $\sigma=\sigma_{max}$, see notation, Figure 10), $\phi_3$-values due to uniform pressures are found to be greater than those due to linear pressure.

It is important to note that one obtains Cartesian components of $\phi$ over the crack, and cylindrical components on the cylinder. Some interesting results can be deduced from the $\phi_i$-values, $i \in \{\rho, \theta, 3\}$, over the free surface. In the opening mode for instance, the point belongs to the upper crack face and located at the centre of the ellipse has, besides a positive $\phi_3$, a negative $\phi_\rho$. This means that it moves upward with respect to the crack and outward with respect to the cylinder, as to some extent predicted by physical considerations.

Figures 12, 13 and 14 represent the crack in non-deformed and deformed states, where for a better representation deformations are emphasized by a multiplier equal to 5.

It is rather attractive to represent the stress intensity factor variations along the crack edge, but in surface crack problems, the stress intensity factor definition is almost difficult to establish.

<figure><img src="image_13.png" /><figcaption>Figure 11. Surface crack mesh: 58 elements, 60 interior nodes. The elements are shaped into the elliptic co-ordinates defined on the ellipse</figcaption></figure>
<table><caption>Table III. Components of $E/(16\pi(1-v^2)) \cdot \phi_{cr}$ resulting from a uniform pressure: $\phi_1=\phi_2=0$, $\phi_3 \neq 0$; ($v=0.3$, $b/R=0.4$, $b/a=5/7$, $b=5$, $\sigma=1$), the unit of a, b, R is [L], that of $\sigma$ is [F \cdot L$^{-2}$], that of the above components is [F \cdot L$^{-1}$])</caption><thead><tr><th>Node</th><th>E / (16π(1 - v²))·φ₃</th><th>Node</th><th>E / (16π(1 - v²))·φ₃</th><th>Node</th><th>E / (16π(1 - v²))·φ₃</th><th>Node</th><th>E / (16π(1 - v²))·φ₃</th></tr></thead><tbody><tr><td>1</td><td>0.4658</td><td>2</td><td>0.4583</td><td>3</td><td>0.4583</td><td>4</td><td>0.4351</td></tr><tr><td>5</td><td>0.4351</td><td>6</td><td>0.3955</td><td>7</td><td>0.3955</td><td>8</td><td>0.3445</td></tr><tr><td>9</td><td>0.3540</td><td>10</td><td>0.3745</td><td>11</td><td>0.3857</td><td>12</td><td>0.3894</td></tr><tr><td>13</td><td>0.3857</td><td>14</td><td>0.3745</td><td>15</td><td>0.3540</td><td>16</td><td>0.3445</td></tr><tr><td>17</td><td>0.2777</td><td>18</td><td>0.2838</td><td>19</td><td>0.2981</td><td>20</td><td>0.3110</td></tr><tr><td>21</td><td>0.3184</td><td>22</td><td>0.3211</td><td>23</td><td>0.3184</td><td>24</td><td>0.3110</td></tr><tr><td>25</td><td>0.2981</td><td>26</td><td>0.2838</td><td>27</td><td>0.2777</td><td>28</td><td>0.2179</td></tr><tr><td>29</td><td>0.2358</td><td>30</td><td>0.2503</td><td>31</td><td>0.2620</td><td>32</td><td>0.2681</td></tr><tr><td>33</td><td>0.2702</td><td>34</td><td>0.2681</td><td>35</td><td>0.2620</td><td>36</td><td>0.2503</td></tr><tr><td>37</td><td>0.2358</td><td>38</td><td>0.2179</td><td>39</td><td>0.1554</td><td>40</td><td>0.1866</td></tr><tr><td>41</td><td>0.2014</td><td>42</td><td>0.2110</td><td>43</td><td>0.2157</td><td>44</td><td>0.2173</td></tr><tr><td>45</td><td>0.2157</td><td>46</td><td>0.2110</td><td>47</td><td>0.2014</td><td>48</td><td>0.1866</td></tr><tr><td>49</td><td>0.1554</td><td>50</td><td>0.0876</td><td>51</td><td>0.1358</td><td>52</td><td>0.1444</td></tr><tr><td>53</td><td>0.1489</td><td>54</td><td>0.1500</td><td>55</td><td>0.1511</td><td>56</td><td>0.1500</td></tr><tr><td>57</td><td>0.1489</td><td>58</td><td>0.1444</td><td>59</td><td>0.1358</td><td>60</td><td>0.0876</td></tr></tbody></table>

<table><caption>Table IV. Components of $E/(16\pi(1-v^2)) \cdot \phi_{cr}$ resulting from a linear pressure with $\sigma_{max}=1$: $\phi_1=\phi_2=0$, $\phi_3 \neq 0$</caption><thead><tr><th>Node</th><th>E / (16π(1 - v²))·φ₃</th><th>Node</th><th>E / (16π(1 - v²))·φ₃</th><th>Node</th><th>E / (16π(1 - v²))·φ₃</th><th>Node</th><th>E / (16π(1 - v²))·φ₃</th></tr></thead><tbody><tr><td>1</td><td>0.3994</td><td>2</td><td>0.3927</td><td>3</td><td>0.3927</td><td>4</td><td>0.3722</td></tr><tr><td>5</td><td>0.3722</td><td>6</td><td>0.3371</td><td>7</td><td>0.3371</td><td>8</td><td>0.2921</td></tr><tr><td>9</td><td>0.3018</td><td>10</td><td>0.3193</td><td>11</td><td>0.3287</td><td>12</td><td>0.3317</td></tr><tr><td>13</td><td>0.3287</td><td>14</td><td>0.3193</td><td>15</td><td>0.3018</td><td>16</td><td>0.2921</td></tr><tr><td>17</td><td>0.2341</td><td>18</td><td>0.2394</td><td>19</td><td>0.2500</td><td>20</td><td>0.2592</td></tr><tr><td>21</td><td>0.2643</td><td>22</td><td>0.2661</td><td>23</td><td>0.2643</td><td>24</td><td>0.2592</td></tr><tr><td>25</td><td>0.2500</td><td>26</td><td>0.2394</td><td>27</td><td>0.2341</td><td>28</td><td>0.1829</td></tr><tr><td>29</td><td>0.1967</td><td>30</td><td>0.2061</td><td>31</td><td>0.2129</td><td>32</td><td>0.2163</td></tr><tr><td>33</td><td>0.2175</td><td>34</td><td>0.2163</td><td>35</td><td>0.2129</td><td>36</td><td>0.2061</td></tr><tr><td>37</td><td>0.1967</td><td>38</td><td>0.1829</td><td>39</td><td>0.1300</td><td>40</td><td>0.1537</td></tr><tr><td>41</td><td>0.1623</td><td>42</td><td>0.1669</td><td>43</td><td>0.1688</td><td>44</td><td>0.1695</td></tr><tr><td>45</td><td>0.1688</td><td>46</td><td>0.1669</td><td>47</td><td>0.1623</td><td>48</td><td>0.1537</td></tr><tr><td>49</td><td>0.1300</td><td>50</td><td>0.0732</td><td>51</td><td>0.1099</td><td>52</td><td>0.1134</td></tr><tr><td>53</td><td>0.1139</td><td>54</td><td>0.1133</td><td>55</td><td>0.1135</td><td>56</td><td>0.1133</td></tr><tr><td>57</td><td>0.1139</td><td>58</td><td>0.1134</td><td>59</td><td>0.1099</td><td>60</td><td>0.0732</td></tr></tbody></table>
<table><caption>Table V. Components of $E/(16\pi(1-v^2))\cdot\phi_{cr}$ resulting from a torsion shear distribution with $\tau_{max}=1: \phi_3=0, \phi_1, \phi_2 \neq 0$</caption><thead><tr><th>Node</th><th>E/(16π(1-v²))·φ₁</th><th>E/(16π(1-v²))·φ₂</th><th>Node</th><th>E/(16π(1-v²))·φ₁</th><th>E/(16π(1-v²))·φ₂</th></tr></thead><tbody><tr><td>1</td><td>0.0000</td><td>-0.2913</td><td>2</td><td>0.0074</td><td>-0.2858</td></tr><tr><td>3</td><td>-0.0074</td><td>-0.2858</td><td>4</td><td>0.0125</td><td>-0.2698</td></tr><tr><td>5</td><td>-0.0125</td><td>-0.2698</td><td>6</td><td>0.0135</td><td>-0.2431</td></tr><tr><td>7</td><td>-0.0135</td><td>-0.2431</td><td>8</td><td>0.0099</td><td>-0.2107</td></tr><tr><td>9</td><td>0.0265</td><td>-0.2554</td><td>10</td><td>0.0234</td><td>-0.2758</td></tr><tr><td>11</td><td>0.0132</td><td>-0.2871</td><td>12</td><td>0.0000</td><td>-0.2907</td></tr><tr><td>13</td><td>-0.0132</td><td>-0.2871</td><td>14</td><td>-0.0234</td><td>-0.2758</td></tr><tr><td>15</td><td>-0.0265</td><td>-0.2554</td><td>16</td><td>-0.0099</td><td>-0.2107</td></tr><tr><td>17</td><td>0.0044</td><td>-0.1761</td><td>18</td><td>0.0258</td><td>-0.2104</td></tr><tr><td>19</td><td>0.0352</td><td>-0.2334</td><td>20</td><td>0.0302</td><td>-0.2456</td></tr><tr><td>21</td><td>0.0165</td><td>-0.2510</td><td>22</td><td>0.0000</td><td>-0.2528</td></tr><tr><td>23</td><td>-0.0165</td><td>-0.2510</td><td>24</td><td>-0.0302</td><td>-0.2456</td></tr><tr><td>25</td><td>-0.0352</td><td>-0.2334</td><td>26</td><td>-0.0258</td><td>-0.2104</td></tr><tr><td>27</td><td>-0.0044</td><td>-0.1761</td><td>28</td><td>-0.0002</td><td>-0.1488</td></tr><tr><td>29</td><td>0.0363</td><td>-0.1927</td><td>30</td><td>0.0413</td><td>-0.2039</td></tr><tr><td>31</td><td>0.0345</td><td>-0.2104</td><td>32</td><td>0.0185</td><td>-0.2120</td></tr><tr><td>33</td><td>0.0000</td><td>-0.2126</td><td>34</td><td>-0.0185</td><td>-0.2120</td></tr><tr><td>35</td><td>-0.0345</td><td>-0.2104</td><td>36</td><td>-0.0413</td><td>-0.2039</td></tr><tr><td>37</td><td>-0.0363</td><td>-0.1927</td><td>38</td><td>0.0002</td><td>-0.1488</td></tr><tr><td>39</td><td>0.0050</td><td>-0.1378</td><td>40</td><td>0.0320</td><td>-0.1548</td></tr><tr><td>41</td><td>0.0394</td><td>-0.1640</td><td>42</td><td>0.0336</td><td>-0.1661</td></tr><tr><td>43</td><td>0.0184</td><td>-0.1652</td><td>44</td><td>0.0000</td><td>-0.1648</td></tr><tr><td>45</td><td>-0.0184</td><td>-0.1652</td><td>46</td><td>-0.0336</td><td>-0.1661</td></tr><tr><td>47</td><td>-0.0394</td><td>-0.1640</td><td>48</td><td>-0.0320</td><td>-0.1548</td></tr><tr><td>49</td><td>-0.0050</td><td>-0.1378</td><td>50</td><td>0.0016</td><td>-0.0972</td></tr><tr><td>51</td><td>0.0611</td><td>-0.1335</td><td>52</td><td>0.0580</td><td>-0.1193</td></tr><tr><td>53</td><td>0.0455</td><td>-0.1103</td><td>54</td><td>0.0222</td><td>-0.1024</td></tr><tr><td>55</td><td>0.0000</td><td>-0.1006</td><td>56</td><td>-0.0222</td><td>-0.1024</td></tr><tr><td>57</td><td>-0.0455</td><td>-0.1103</td><td>58</td><td>-0.0580</td><td>-0.1193</td></tr><tr><td>59</td><td>-0.0611</td><td>-0.1335</td><td>60</td><td>-0.0016</td><td>-0.0972</td></tr></tbody></table>

Generally speaking, the elastic state varies from plane strain in the interior to plane stress at the free surface. In the following, we shall limit our consideration to the stress intensity factors at the deepest point of the crack. There exists at this point a plane strain state; in the case of $v=0.3$, $b/R=0.4$, $b/a=5/7$, the stress intensity factors are given by:

$$ \begin{aligned} \text{for a uniform pressure, } \quad K_1/(\sigma\sqrt{\pi b}) &= 1/(\sigma\sqrt{\pi b}) \times 2\pi\sqrt{2\pi} \lim_{r\to 0} \{E/(16\pi(1-v^2))\cdot\phi_3/\sqrt{r}\} \\ &= 0.820 \end{aligned} $$

$$ \begin{aligned} \text{for linear pressure, } \quad K_1/(\sigma_{\max}\sqrt{\pi b}) &= 1/(\sigma_{\max}\sqrt{\pi b}) \times 2\pi\sqrt{2\pi} \lim_{r\to 0} \{E/(16\pi(1-v^2))\cdot\phi_3/\sqrt{r}\} \\ &= 0.595 \end{aligned} $$
<figure><img src="image_14.png" /><figcaption>Figure 12. Surface crack under a uniform pressure: $\phi_1 = \phi_2 = 0$, $\phi_3 \neq 0$. Crack shapes before and after deformation: (a) front view; (b) spatial view</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13. Surface crack under a linear pressure: $\phi_1 = \phi_2 = 0$, $\phi_3 \neq 0$. The deformation shape is analogous to that due to a uniform pressure: (a) front view; (b) spatial view</figcaption></figure>

for a torsion shear distribution,

$$ 
\begin{aligned}
K_{III}/(\tau_{max}\sqrt{\pi b}) &= 1/(\tau_{max}\sqrt{\pi b}) \times 2\pi\sqrt{2\pi(1-v)} \times \lim_{r \to 0} \{E/(16\pi(1-v^2)) \cdot \phi_2/\sqrt{r}\} \\
&= -0.298
\end{aligned}
$$

(As previously mentioned, $E/(16\pi(1-v^2)) \cdot \phi_i$, $i=2$ or 3, are obtained directly from the program). The stress intensity factors are normalized by $\sigma_{max}\sqrt{\pi b}$ or $\tau_{max}\sqrt{\pi b}$ related to a slit crack with length $2b$ in the infinite medium, under a remote tensile or shear loading.
<figure><img src="image_16.png" /><figcaption>Figure 14. Surface crack under a torsion shear distribution: upper view</figcaption></figure>

In Reference 16, stress intensity factors are computed by means of the boundary integral equation method for semi-elliptical cracks in cylindrical bars of diameter equal to 12 mm, in tension or in bending. The crack geometry corresponding to that considered herein is $b/a = 2.4$ mm/3.36 mm = 5/7, $b/2R = 2.4$ mm/12 mm = 0.2, for which the stress intensity factor can be obtained by interpolation:

$$ \text{in tensile loading with } \sigma = 1000 \text{ MPa, } \quad K_1 = 72.62 \text{ MPa}\sqrt{m} \rightarrow K_1/(\sigma\sqrt{\pi b}) = 0.836 $$

$$ \text{in pure bending with } \sigma_{\max} = 953.4 \text{ MPa, } \quad K_1 = 47.10 \text{ MPa}\sqrt{m} \rightarrow K_1/(\sigma_{\max}\sqrt{\pi b}) = 0.569 $$

One can find in Reference 17 a list of stress intensity factors for various crack and solid geometries. The so-called surface crack in solid cylinder is studied under tensile loading and pure bending (Figure 15). These actually are almond-shaped cracks which occur perpendicularly to the lateral surface of the cylinder. For lack of matter for comparisons, we consider, however, the crack of this type with $b/R = 0.4$; at the deepest point of the crack, the computed values are expected to be close to the above ones. The stress intensity factor is then expressed in the additive form¹⁷

$$ K_1 = [\sigma \cdot F_0 + \sigma_{\max} \cdot F_1] \sqrt{\pi b} $$

$$ F_0 = G[0.752 + 1.286\beta + 0.37Y^3] $$

$$ F_1 = G[0.923 + 0.199Y^4] $$

$$ G = 0.92(2/\pi) \cdot 1/\cos \beta \cdot [\tan \beta/\beta]^{1/2} $$

$$ Y = 1 - \sin \beta $$

$$ \beta = \pi/2 \cdot b/2R $$

$$ \text{which yields in tensile loading: } \quad K_1/(\sigma\sqrt{\pi b}) = 0.800 $$

$$ \text{and in pure bending: } \quad (K_1/(\sigma_{\max}\sqrt{\pi b})) = 0.606. $$

Recently, in Reference 18, use is made of the finite element method to study semi-elliptical cracks in cylindrical bars under tensile loadings. When computing the stress intensity factor $K_1$, the assumption that a plane strain state holds in the neighbourhood of the crack is made in order to apply the virtual crack extension method. Charts for $K_1$ at the centre of the crack are given for $b/a = 0.01, 0.2, 0.5, 1.0, 2.0$ and for various $b/R$ ratios. An extrapolation from these $K_1$-charts then gives us approximately the interval containing the $K_1/(\sigma\sqrt{\pi b})$ ratio: for $b/R = 0.4, b/a = 5/7$, with a uniform pressure $\sigma$, $K_1/(\sigma\sqrt{\pi b}) \in [0.93, 1.00]$.

The discrepancy between the above results is mainly due to the crack geometry differences; nevertheless, the values obtained are of the same order. For precise comparisons to be made,
<figure><img src="image_17.png" /></figure>

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

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

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

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

<figure><img src="image_22.png" /><figcaption>Figure 15. Surface crack in a solid cylinder (figure from Reference 17)</figcaption></figure>
investigations are being pursued to obtain more complete abacuses of stress intensity factors, reported for various crack configurations.

## CONCLUDING REMARKS

Equation (4) has been discretized and written in a tensor form, equation (10). Whereas lengthy expressions should be expected if (4) had been dealt with as three separate scalar equations, this approach yields a concise form of the discretized equation and thus allows a compact way of programming. Moreover, a crack parametrization being chosen; this enables one to obtain directly the components of the displacement discontinuity vector at each point, as well as the stress intensity factors along the crack edge, in a *local* basis. Unlike the usual finite element method, no singular element is required near the crack edge. The method herein proposed can also easily be applied to non-symmetrical crack configuration problems for which the boundary integral equation method becomes inefficient. Moreover, since no *a priori* assumptions of a plain strain or plain stress are to be made, the integral equation method can be used to investigate the stress singularity in the neighbourhood of the crack.

The example on the circular crack shows that one can obtain good results, using only four-node elements and 2 × 2 Gaussian points. The problem of the semi-elliptical crack in a cylindrical bar was also investigated with combined loadings. The numerical results obtained in opening mode were found to be in the range of different existing results for analogous crack geometries. The proposed method also allows one to obtain results for more complex loadings, in particular for torsion loading to which few works in the literature are devoted.

# APPENDIX

## Notation

$b=(\mathbf{b}_1; \mathbf{b}_2; \mathbf{b}_3)$ a local basis at a point $y$ in $E_3$

$c_n$ the canonical basis of $\mathbb{R}^n$

$\mathcal{D}$ the elastic body characterized by constants $(E, \nu)$

$\partial \mathcal{D}$ the boundary of $\mathcal{D}$

$e=(\mathbf{e}_1; \mathbf{e}_2; \mathbf{e}_3)$ the fixed basis of $E_3$

$E_3$ the three-dimensional Euclidean space

$F$ a parametrization of $S$, defined on a domain $\Delta$ in $\mathbb{R}^2$: $\Delta\ni(u, v)\rightarrow y=F(u, v)\in S$. In general, $F_{,u}$ may not be perpendicular to $F_{,v}$

$I$ unit tensor

$\mathbf{n}_y$ normal to $S$ at $y\in S$

$S$ the crack surface

$\partial S$ the crack edge

$\mathbf{t}(y, \mathbf{n}_y)$ stress vector at $y\in S$, with respect to normal $\mathbf{n}_y$

$T$ Kupradze tensor

$y, y_0$ respectively an ordinary, a particular point of $S$

$\mathrm{pv}\int_S$ surface integral on $S$, defined in the principal value sense

$\{a\}$ column-matrix of components of a vector **a** in a basis. One may encounter $b\{a\}$ which emphasizes the basis **b** in question

$\langle a \rangle$ row-matrix, $=\{a\}^T$, transpose of $\{a\}$
[$\mathcal{N}$] matrical representation of a linear transformation $\mathcal{N}$ from $E_3$ into $\mathbb{R}^n$, also denoted more precisely by $e[\mathcal{N}]$

$$ (a, b, c) = (a \wedge b) \cdot c $$

$a \otimes b$ tensor product defined by $\forall c \in \mathbb{R}^3$, $(a \otimes b)c = a(b \cdot c)$

REFERENCES

1. T. A. Cruse, *Int. J. Solids Struct.*, **5**, 1259-1274 (1969).

2. J. M. Boissenot, J. C. Lachat and J. Watson, *Revue Physique Appliquée*, **9**, 611 (1974).

3. V. D. Kupradze, T. G. Gegelia, M. O. Bashelishvili and T. V. Burchuladze, *Three-Dimensional Problems of the Mathematical Theory of Elasticity and Thermoelasticity*, North-Holland, Amsterdam, 1979.

4. H. D. Bui, *Compt. Rend. Acad. Sci. (Paris)*, Serie A, **280**, 1157 (1975).

5. T. A. Cruse, *AFOSR-TR-0813*, 1975, pp. 13-20.

6. P. S. Theocaris and J. G. Kazantzakis, *Int. J. Fract.*, **R117**-R119 (1978).

7. J. Weaver, *Int. J. Solids Struct.*, 321-330 (1977).

8. V. Sladek and J. Sladek, *Int. J. Solids Struct.*, **19**, 425-436 (1983).

9. A. Levan and J. Royer, *Int. J. Fract.*, **31**, 125-142 (1986).

10. J. G. Kazantzakis and P. S. Theocaris, *Int. J. Solids Struct.*, **15**, 203-207 (1979).

11. P. S. Theocaris and J. G. Kazantzakis, *Int. J. Fract.*, R165-R167 (1979).

12. P. S. Theocaris, N. I. Ioakimidis and J. G. Kazantzakis, *Int. J. Numer. Methods Eng.*, **629** (1980).

13. N. I. Ioakimidis, *Appl. Numer. Methods*, **1**, 183-189 (1985).

14. Liu Jun, G. Beer and J. L. Meek, *Eng. Analysis*, **2**, (3) (1985).

15. I. N. Sneddon, *Proc. Roy. Soc. London*, **187**, 229 (1946).

16. A. Athanassiadis, J. M. Boissenot, P. Brevet, D. Francois and A. Raharinaivo, *Int. J. Fract.*, 553-566 (1981).

17. R. C. Forman, V. Shivakumar, J. C. Newman, L. Williams and S. Piotrowski, *Fatigue Crack Growth Computer Program*. NASA/FLAGRO, Version June 1985.

18. M. A. Astis, M. Elices and A. Valiente, ECF6, June 1986.