# Regular Boundary Integral Equations for Three-Dimensional Finite or Infinite Bodies With and Without Curved Cracks in Elastodynamics

Marc Bonnet

Laboratoire de Mécanique des Solides, (UA CNRS 317),
Ecole Polytechnique, 91128 Palaiseau cedex, France

## Summary.

In this paper, a regular BIE is proposed for the time-domain elastodynamic analysis of 3D cracked solids. Since the displacement BIE applied to cracks leads to an ill-posed problem, we need a separate treatment using a traction BIE which is originally hypersingular. The CPV exclusion technique is first applied to the displacement BIE. Then the hypersingular traction BIE is reformulated in terms of singular CPV integrals and finally regularized.

The regularization of the steady-state elastodynamic BIE formulation (or, equivalently, of the Fourier- or Laplace-transformed domain BIE) is obtained as a special case of the previous results.

Finally some features of the boundary element discretization process related to the regularization technique are outlined.

## Introduction.

The boundary element method (BEM) is now well-established and applied in numerous engineering fields (structural analysis, geomechanics, acoustics, among others) [2][3][10][13][21]. However the integral equations of elastodynamics involve strongly singular kernels, integrable only in the Cauchy principal value sense. This seems to be the main obstacle to a satisfactory numerical implementation of BEM, mainly because these singular integrals correspond to near-diagonal terms in the matrix of the associated discrete problem. As pointed out in [4] [5], the usual techniques of integral calculus (differentiation under the integral sign, change of variable) applied to CPV integrals may give false results. On the other side, the implementation of BEM which handles correctly the CPV integrals becomes too difficult a task for general interpolations of the geometry and unknowns.

To overcome this difficulty, reformulations of BIE which use systematically integrable kernels, which we call regularization techniques, appear to be a good approach. Rizzo & Shippy [18] were the first to state a regularized BIE for bounded uncracked bodies in elastostatics. No extension to dynamics and/or infinite media was made until the simultaneous appearance of two similar works [7] [19] giving regularized BIE for finite or infinite bodies in steady-state elastodynamics. One year later Bonnet [1] gave the seemingly first regular BIE formulation for a general 3D curved crack in an infinite medium, based on the
CPV formulation of Sladek & Sladek [20]. A review of regularization methods can be found in Bui [6]. An up-to-date review of elastodynamic BIE can also be found in Beskos & Manolis [15].

The present work proposes a unified regularization approach for the BIE formulation of time-domain linear elastodynamic problems (including elastostatics and steady-state elastodynamics as special cases), applicable to three-dimensional finite or infinite bodies, with or without curved cracks. This approach includes all the situations treated in the above-mentioned works and extends them to the time-domain elastodynamic case.

# 1 - Statement of the problem and integral representation of the displacement and stress fields

We consider the dynamic response of a 3D cracked elastic body $\Omega$, with Lamé constants $(\lambda, \mu)$, subjected to surface loadings. Without loss of generality regarding the regularization problem we assume, for the simplicity of presentation, the body forces and the initial conditions to be zero. Furthermore we assume symmetrical loads for the crack, ie the tractions applied on the two crack faces to be opposite. $\Omega$ can be either infinite, semi-infinite or finite (see figure 1). The boundary of $\Omega$ is split in two parts $\Gamma$ and $\Gamma_c$:

(a) The surface $\Gamma$ represents the boundary of either the finite body and/or the internal cavities or inclusions, according to the geometry in consideration. $\Gamma$ is a sectionally Lyapunov [9][11] surface, with possibly edges or corners at discrete locations. Whatever the geometry considered, the unit normal vector $n$, where it exists, is directed towards the exterior of the material body (see figure 1). The surface $\Gamma$ can be the union of $M$ disjoint surfaces $\Gamma_1, ..., \Gamma_M$; this is the case if $\Omega$ is a multiply-connected domain.

(b) The surface $\Gamma_c$ represents a perfect crack embedded in $\Omega$. $\Gamma_c$ is a sectionally twice continuously differentiable surface, and coincides with the lower face $\Gamma_c$ of the crack. The unit normal $n$ to the crack is directed from $\Gamma_c$ to $\Gamma_c^+$.

The motion (or equilibrium) of $\Omega$ is then governed by the homogeneous Navier equation for the displacement field $u(y,t)$ ($y \in \Omega$ and $t \in T^+ = [0, +\infty)$:

$$ \left\{ \frac{\partial u}{\partial \rho} \Delta + \frac{\lambda + \mu}{\rho} \operatorname{grad} \operatorname{div} u - \frac{\partial^2 u}{\partial t^2} = 0 \right. \quad (1a) $$

with zero initial ($t=0$) conditions:

$$ u(x,0) = 0, \quad \dot{u}(x,0) = 0 \quad (1b) $$

Furthermore if $\Omega$ is unbounded, the radiation conditions, expressing the absence of any source or reflector at infinity, must be satisfied [11].

# 1.1 - Stokes tensors $U$ and $\Sigma$ [8]

Let $x \in \mathbb{R}^3$ and $f$ be a twice continuously differentiable function of time that vanishes for $t<0$. Let $U_{ik}[x,t,y|f]$ be the $i^{th}$ component of the (singular) displacement field at $y$ due to the time-dependent point force of magnitude $f(t)$ acting on $x$ along $e_k$-direction in an infinite elastic medium. The displacement tensor $U$ possesses the symmetry properties:

$$ U[x,t,y|f] = U[y,t,x|f] = {}^t U[x,t,y|f] \quad (2a) $$
$$ \mathrm{grad}_y U = - \mathrm{grad}_x U \qquad (2b) $$

The application of Hooke's law yields the corresponding Stokes stress tensor $\Sigma$:

$$ \Sigma_{is;k}[x,t,y|f] = C_{ispq} \frac{\partial}{\partial y_q} U_{p;k}(x,t,y|f) \qquad (3a) $$

$$ C_{ispq} = \lambda \delta_{is} \delta_{pq} + \mu (\delta_{ip} \delta_{sq} + \delta_{iq} \delta_{sp}) \qquad (3b) $$

Let $r = \|x-y\|$ be the euclidean distance between $x$ and $y$. A key feature of the Stokes tensors is their singular behaviour for $r$ arbitrary small:

$$ U[x,t,y|f] = O(1/r) \qquad (4a) $$

$$ \Sigma[x,t,y|f] = O(1/r^2) \qquad (4b) $$

Another crucial property of the Stokes tensors is the following one:

$$ U[x,t,y|f] - f(t)U^0(x,y) = \dot{f}(t)O(1) \qquad (5a) $$

$$ \Sigma[x,t,y|f] - f(t)\Sigma^0(x,y) = \ddot{f}(t)O(1) \qquad (5b) $$

in which we see that the singularity (4a-4b) of the dynamic tensors is entirely included into their static counterparts $U^0$ and $\Sigma^0$. (5a-5b) is an immediate consequence of the Taylor expansions (A.10) or (A.11) of $U$ and $\Sigma$ for small $r$. This property is not surprising because the dynamic tensors written for a constant force function $f(t) = 1$ must obviously be the static tensors.

## 1.2 - Integral representation of the displacement and stress fields.

Let $z$ be a fixed interior point of $\Omega$. The application of Maxwell-Betti reciprocal identity to the unknown displacement field $u(y,t)$ and the Stokes impulsive displacement tensor $U[z,t,y|\delta(t)]$, where $\delta$ denotes the Dirac delta distribution, leads to the well-known representation formula of $u(z,t)$:

$$ \begin{aligned} u_k(z,t) + \int_{\Gamma} n_s(y) \Sigma_{is;k}[z,t,y|u_i(y,t)] \, dS_y - \int_{\Gamma} U_{i;k}[z,t,y|t_i(y,t)] \, dS_y \\ = \int_{\Gamma_c} n_s(y) \Sigma_{is;k}[z,t,y|\phi_i(y,t)] \, dS_y \end{aligned} \qquad (6) $$

where $\phi_i(y,t) = u(y^+,t) - u(y^-,t)$ is the crack opening displacement (COD).

To obtain the contribution of $\Gamma_c$ in (6), the crack is viewed as a cavity bounded by two identical surfaces $\Gamma_c^+$ and $\Gamma_c^-$ (ie $y^+ \in \Gamma_c^+$, $y^- \in \Gamma_c^-$) of opposite normals $n(y^+) = -n(y^-)$ and the assumed symmetry of the crack loading (ie $t^+ = -t^-$) is taken into account.

The given crack tractions do not appear explicitly in the representation formula (6). The consequence is that all symmetric loadings of the crack lead to the same integral equation for the unknown COD. In order to circumvent this ill-posedness, we shall need the following representation of the stress tensor $\sigma$ associated with $u$ by Hooke's law, in order to get separately the tractions on the two crack faces. Taking into account (2a) (2b) (3a) one has:
$$
\sigma_{lp}(z,t) = C_{lpkr} \int_{\Gamma} n_s(y) \frac{\partial}{\partial y_r} \Sigma_{is;k}[z,t,y|u_i(y,t)] dS_y - \int_{\Gamma} \Sigma_{lp;\dot{i}}[z,t,y|t_i(y,t)] dS_y \\
- C_{lpkr} \int_{\Gamma_c} n_s(y) \frac{\partial}{\partial y_r} \Sigma_{is;k}[z,t,y|\phi_i(y,t)] dS_y \qquad (7)
$$

ie

$$
\sigma_{lp}(z,t) = \sigma_{lp}^{1}(z,t) + \sigma_{lp}^{2}(z,t) + \sigma_{lp}^{3}(z,t)
$$

2 - Regularized displacement BIE .

Letting $z \to x \in \Gamma$ in (6) and using the limit behaviour of single and double
layer elastic potentials [11], one obtains the well-known singular displacement
BIE:

$$
\begin{align}
c_{ki}(x)u_i(x,t) &+ \int_{\Gamma} n_s(y) \Sigma_{is;k}[x,t,y|u_i(y,t)] dS_y - \int_{\Gamma} U_{i;k}[x,t,y|t_i(y,t)] dS_y \nonumber \\
&= \int_{\Gamma_c} n_s(y) \Sigma_{is;k}[x,t,y|\phi_i(y,t)] dS_y \tag{8}
\end{align}
$$

The tensor c(x) is the so-called free term, which depends only on the local
geometry of Ω near the fixed point x. The value of c for an interior point of Ω
and for a regular boundary point is I and (1/2)I respectively (I being the
identity tensor). The general expression of c, valid for edges and corner points
x, is:

$$
c_{ki}(x) = \lim_{\epsilon \to 0} \int_{\partial B_{\epsilon}(x) \cap \Omega} n_s(y) \Sigma_{u_i k}^{0}(x,y) dS_y \quad (9)
$$

where $B_{\epsilon}(x)$ denotes the exclusion ball of radius $\epsilon$ centered at $x$ (see figure 2).
If $c^c(x)$ denotes the free term associated with the complementary set $\Omega^c = \mathbb{R}^3 - \Omega$,
(9) yields the following relation, valid for an arbitrary point $x$ in $\mathbb{R}^3$:

$$
c_{ki}(x) + c_{ki}^{c}(x) = \delta_{ki} \quad (10)
$$

The first two integrals of (8) are singular : U is integrable on Γ in an ordi-
nary sense, but the integral of Σ on Γ exists only in the sense of a Cauchy prin-
cipal value (CPV), denoted by the symbol (*).

Without loss of generality we assume in this part Ω to be uncracked, as we
focus our attention to the CPV singularity on Γ.

2.1 - Regular BIE for a bounded body under static loading

Let us restrict ourselves to the simplest case of Ω bounded under static
loading, for which the displacement BIE (8) reduces to :

$$
c_{ki}(x)u_i(x) + \int_{\Gamma} n_s(y) \Sigma_{is;k}^{0}(x,y) u_i(y) dS_y - \int_{\Gamma} U_{i;k}^{0}(x,y) t_i(y) dS_y = 0 \quad (11)
$$

Let us consider, as an auxiliary solution, a rigid displacement of Ω defined
by :
$$ u(y) = u^0 (\text{constant}) \qquad t = 0 \text{ on } \Gamma \qquad (12) $$

This rigid-body motion satisfies the elastostatic equations. Its associated BIE reads:

$$ c_{kl}(\hat{x}) u_i^0 + u_i^0 \int_{\Gamma}^{*} n_s(y) \Sigma_{is;kl}^0(x,y) dS_y = 0 \qquad (13) $$

The rigid-body identity (13) remains in particular true for the choice $u^0 = u(x)$, ie the unknown displacement at the fixed point $x$ considered as a constant field on $\Omega$. Subtracting (13) from (11) with $u^0 = u(x)$ gives:

$$ \int_{\Gamma} n_s(y) \Sigma_{is;kl}^{0}(x,y) (u_i(y) - u_i(x)) dS_y - \int_{\Gamma} U_{i;kl}^{0}(x,y) t_i(y) dS_y = 0 \qquad (14) $$

The kernel $\Sigma^0$ is singular like $\|x-y\|^{-2}$. Therefore, assuming $u(y)$ to be a Hölder function, ie the existence of $\alpha > 0$ and $C > 0$ such that:

$$ \|u(x)-u(y)\| \le C \|x-y\|^\alpha \qquad (15) $$

the singularity of the first integrand of (14) reduces to $\|x-y\|^{-2+\alpha}$, which is integrable on $\Gamma$, while the second integrand, unchanged, is weakly singular. In fact, the use of the rigid-body identity (13) allowed us to *transfer* the Cauchy principal value into the known term $c(x)$, all the remaining integrals being ordinary improper integrals.

Equation (14), first proposed by Rizzo & Shippy [18], gives the regular BIE for a bounded body $\Omega$ in elastostatics. The numerical implementation of (14), completely free of Cauchy principal values, can be done with all standard boundary elements, allowing high-degree interpolation of $\Gamma$, $u$ and $t$. This is a significant improvement of the classical boundary element formulations, in which the CPV integrals are satisfactorily evaluated only for plane boundary elements and low-degree interpolations. An additional advantage of (14) compared to (11) is the disappearance of the free term $c(x)$; this exempts one from computing its value for a corner point $x$.

However, the identity (13) is no more valid if we want to apply the previous method to an unbounded body and/or a dynamic loading because:

(a) If $\Omega$ is not bounded, a rigid-body motion has not a finite energy and thus violates the radiation condition.

(b) A static rigid-body displacement does not satisfy the elastodynamic equation (1a).

## 2.2 - Generalisation of the rigid-body identity to unbounded domains.

First we examine the point (a) by giving an adequate generalized form of the rigid-body identity (13) for certain unbounded domains.

(i) $\Omega$ is infinite: $\Omega = \mathbb{R}^3 - \Omega^c$, where $\Omega^c$ is a bounded open subset of $\mathbb{R}^3$ of boundary $\Gamma$. Consider a rigid displacement $u^0$ of $\Omega^c$. The identity (12) holds for the free-term $c^0(x)$. Recalling that in our convention the unit normals $n$ and $n^c$ of $\Omega$ and $\Omega^c$ are opposite, we use the relation (10) between $c(x)$ and $c^0(x)$, obtaining:
$$ (c_{ki}(\hat{x}) - \delta_{ki}) u_i^0 + u_i^0 \int_{\Gamma} n_s(y) \Sigma_{is;k}^0(x,y) dS_y = 0 \quad (16) $$

(ii) $\Omega$ is a half-space bounded by the plane $P$. In this case we simply have, for any constant vector $u^0$:

$$ u_i^0 \int_P n_s(y) \Sigma_{is;k}^0(x,y) dS_y = 0 \quad (17) $$

(17) can be established either by a direct evaluation of the integral in polar coordinates centered at $x \in \Pi$ or by using the antisymmetry property of $\Sigma^0$ with respect to the vector $x-y$ (see equation (A.9)).

(iii) $\Omega$ is an indented half-space, ie has a plane surface $\Gamma_2$, except for a finite region (see figure 2) bounded by a surface $\Gamma_1$. This geometrical configuration can be used to describe the embedded foundations of a structure in an infinite elastic soil. The following identity holds [1] for any point $x \in \Gamma_1 \cup \Gamma_2$, in particular for the edge points located on $\Gamma_1 \cap \Gamma_2$, and for any constant vector $u^0$:

$$ (c_{ki}(x) - \frac{1}{2}\delta_{ki}) u_i^0 + u_i^0 \int_{\Gamma} n_s(y) \Sigma_{is;k}^0(x,y) dS_y = 0 \quad (18) $$

To summarise (i)-(i)-(iii), it is easy to see that the identities (13), (16), (17) and (18) can be gathered into the following form:

$$ (c_{ki}(x) - \gamma \delta_{ki}) u_i^0 + u_i^0 \int_{\Gamma} n_s(y) \Sigma_{is;k}^0(x,y) dS_y = 0 \quad (19) $$

$$ \text{with } \begin{cases} \gamma = 1: \text{unbounded domain} \\ \gamma = 1/2: \text{half-space (possibly indented)} \\ \gamma = 0: \text{bounded domain} \end{cases} $$

### Regularization of the BIE (8) in the general case

Let us recall the displacement BIE (8) and rewrite it in the following form:

$$ c_{ki}(x)u_i(x,t) + \int_{\Gamma} n_s(y) \left\{ \Sigma_{is;k}[x,t,y|u_i(y,t)] - \Sigma_{is;k}^0(x,y) u_i(y,t) \right\} dS_y + \\ + \int_{\Gamma} n_s(y) \Sigma_{is;k}^0(x,y) u_i(y,t) dS_y - \int_{\Gamma} U_{i;k}[x,t,y|t_i(y,t)] dS_y = 0 \quad (20) $$

Thanks to the key property (5b) of the second Stokes tensor, the kernel in the first integral is regular. The CPV integral is then shifted into the static kernel $\Sigma^0$. Now use the generalized rigid-body identity (19), giving $u^0$ the value $u(x,t)$ for fixed $x$ and $t$, and subtract it from (20). This yields:
$$
\begin{equation}
\begin{aligned}
& \gamma \delta_{ki} u_i(x,t) + \int_{\Gamma} n_s(y) \left[ \Sigma_{ijk}[x,t,y|u_i(y,t)] - \Sigma_{ijk}^{0}(x,y) u_i(y,t) \right] dS_y + \\
& \qquad + \int_{\Gamma} n_s(y) \Sigma_{ijk}^{0}(x,y)(u_i(y,t) - u_i(x,t)) dS_y - \int_{\Gamma} U_{i,k}[x,t,y|t_i(y,t)] dS_y = 0
\end{aligned}
\tag{21}
\end{equation}
$$

Equation (21) is the general regularized BEM for uncracked domains we look for : the three integrals are weakly singular thanks to (5b), (15) and (4a) respectively . This modified BIE formulation, free of CPV integrals, is fully suitable for boundary element discretization . Although the tensor c(x) does not appear explicitly any more, (21) remains valid for x located at an edge or corner point of Γ . The boundedness or unboundedness of Ω in (21) is entirely taken into account, given the coefficient γ and the orientation choice for n .

Finally, if one considers a cracked body, the complete displacement BIE on
Γ is obtained by adding the two (regular) integrals over Γc of (8) to the right-
hand side of equation (21) .

3 - Regularized traction BIE for curved cracks

It has been pointed out in part 1 that the displacement BIE leads to an ill-posed problem for a prescribed crack surface traction. Thus the traction BIE is needed. The latter should be obtained by letting $z \to x_c \in \Gamma_c$ in the integral representation (7) of the stress tensor. However the third term $\sigma_3(z)$ is written in terms of the space derivatives of $\Sigma$, thus involving a $||x_c-y||^{-3}$ singularity on $\Gamma_c$ which is not integrable on $\Gamma_c$ even in a CPV sense. Since classical results of elastic potential theory do not apply to $\sigma_3(z)$, it is necessary, before any regularization process, to reduce the singularity of the kernel of $\sigma^2$.

3.1 - Tangential operators .

In order to do so, we are led to perform an integration by parts over curved surfaces, thus to introduce several mathematical identities using tangential operators. Consider first a field $f(y)$ function of $y \in \Gamma_c$. The function $f$ is not defined outside $\Gamma_c$, as illustrated by the examples $f = \phi_i(y)$ (crack opening displacement) and $f = n_i(y)$. Therefore the cartesian derivatives $f_j$ are meaningless unless $e_j$ are the basis vectors of the plane of a (planar) crack. For curved cracks, we have to extend in some way the domain of definition of $f$, for example by introducing the function $\hat{f}$

$$
(y \in V) \quad \hat{f}(y) = f(P(y)) \tag{22}
$$

where V is a neighbourhood of Γ_c, P(y) is the orthogonal projection of y onto
Γ_c. Clearly the restriction of f to Γ_c is equal to f. Moreover the normal deri-
vative of f is equal to zero, i.e. the vector grad f is tangent to Γ: therefore it
may be used to define the tangential gradient grad Γ_c of the function f. For the
consistency of notations we put :

$$
\mathrm{grad}_{\Gamma_c} = \mathrm{grad} - n \frac{\partial}{\partial n} = \epsilon_r D_r^y \quad (23a)
$$

$$
\text{Thus} \qquad \mathrm{grad}_{\Gamma_c} f = \mathrm{grad}_{\Gamma_c} \hat{f} = \mathrm{grad} \hat{f} \qquad (23b)
$$
In the following, unless stated otherwise, we shall omit the symbol (^), keeping in mind when necessary the extension (22). Let f and v be respectively a scalar and a vector field on Γ<sub>c</sub>. Application of Stokes identity to the vector field f(nx v) yields the following identity :

$$
\int_{\Gamma} ( \mathrm{f.n.rot}(\mathrm{nx v}) + \mathrm{v.grad}_{\Gamma} \mathrm{f} ) \, \mathrm{dS} = \int_{\partial \Gamma} \mathrm{f.} v. \, \mathrm{d}l \quad (24)
$$

Furthermore the tangential operator

$$
\mathbf{n} \cdot \mathbf{grad} = \mathbf{e}_j \epsilon_{jrs} (n_r D_s^y - n_s D_r^y) = \mathbf{e}_j \epsilon_{jrs} D_{rs}^y \quad (25)
$$

satisfies the identity :

$$
\int_{\Gamma_c} \mathbf{n} \cdot \mathbf{grad} \, \mathbf{f} \, \mathrm{dS} = \int_{\partial \Gamma_c} \mathbf{f} \, \mathrm{d}l \quad (26)
$$

where $\tau$ is the unit tangent vector to $\partial\Gamma_c$ and $\nu = \tau \cdot \mathbf{n}$

3.2 - Reformulation of σ³ with CPV integrals.

This task has been performed by Sladek & Sladek [20] for elastostatics and steady-state elastodynamics (see also Le Van [14]). Our aim here is not to repeat in full detail the derivation of the CPV traction BIE but rather to outline some key steps of the method of derivation used in [14][20], which we generalize for transient elastodynamics and a complex geometry of cracks in elastodynamics.

First, setting $\mathbf{f} = \sum_{i=1}^{N} [z_i, t, y_i] \phi_i(y, t)$, we integrate by parts the expression (7) of $\sigma^3$ using identity (26). This leads to:

$$
\sigma_{lp}^3(z,t) = -C_{lpkr} \left\{ \int_{\Gamma_c} \sum_{i=1}^{N} [z_i, t, y_i] D_{rs}^y \phi_i(y,t) dS_y + \rho \int_{\Gamma_c} U_{ijk} [z, t, y] \ddot{\phi}_i(y,t) n_r(y) dS_y \right\} (27)
$$

The above equation (27) holds a priori for a smooth crack only. If one considers a kinked crack with N branches : $\Gamma_c = \Gamma_c^1 \cup \dots \cup \Gamma_c^N$, a line integral term appears in the right hand side of (27) and is given by :

$$
\text{[line integral]} = \frac{1}{2} C_{lpkr} \sum_{\alpha=1}^{N} \int_{\partial \Gamma_c^{\alpha}} \sum_{i=1}^{N} [z_i, t, y] \phi_i(y,t) \epsilon_{jrs} \tau_{j[\alpha]}(y) \, dS_y \quad (28)
$$

In fact, the line integral (28) vanishes because the CODs $\phi^\alpha$ of each branch
$\alpha$ satisfy :

$$
\sum_{\alpha=1}^{N} \phi^{\alpha}(y) = 0 \quad (y \in \mathbb{R}) \qquad (29)
$$

with the respective orientations of the branches chosen as shown on figure 3.
Thus expression (27) holds for a smooth crack as well as a kinked or crossed
crack. This is pointed out by Le Van [14].

Once the integration by parts is done, the kernel of the first integral in σ³
leads to a CPV integral (27) as z → x ∈ Γ_c (the second one being an ordinary
improper integral). In order to relate the (given) surface tractions on both sides
of the crack to the (unknown) COD, we investigate the limit of
$t^3(z, n^±(x)) = σ^3(z) n^±(x)$ as $z \to x^± \in \Gamma_c^±$, $x$ being an interior point of the surface $\Gamma_c$. Using results of elastic potential theory [12], the free term arising from this limit processing is shown [14],[20] to be vanishing. Thus the surface traction associated to $σ^3$ is given by:

$$ t_1^{2±}(x,t) = - C_{ipkr} n_p^±(x) \left\{ \int_{\Gamma_c} \sum_{i,j,k} [x,t,y] D_{r\mu}^Y \phi_i(y,t) \right\} dS_y \\ + \rho \int_{\Gamma_c} n_r(y) U_{i,j,k} [x,t,y] \bar{\phi}_i(y,t) dS_y \quad (30) $$

### 3.3 - Exclusion of the CPV.

As the interior of the (mathematical) crack surface is empty, we cannot derive a rigid-body identity to perform the singularity exclusion as in part 2. Nevertheless, let us rewrite (30) in order to reject the CPV singularity in an integral of the static kernel $\Sigma^0$ alone:

$$ t_1^{2±}(x,t) = - C_{ipkr} n_p^±(x) \left\{ \int_{\Gamma_c} \left[ \sum_{i,j,k} [x,t,y] D_{r\mu}^Y \phi_i(y,t) - D_{r\mu}^Y \phi_i(y,t) \sum_{i,j,k} \Sigma^0_{i,j,k}(x,y) \right] dS_y + \right. \\ \left. + \int_{\Gamma_c} (D_{r\mu}^Y \phi_i(y,t) - D_{r\mu}^X \phi_i(x,t)) \Sigma^0_{i,j,k}(x,y) dS_y + \rho \int_{\Gamma_c} n_r(y) U_{i,j,k} [x,t,y] \bar{\phi}_i(y,t) dS_y + \right. \\ \left. + D_{r\mu}^X \phi_i(x,t) \int_{\Gamma_c} \Sigma^0_{i,j,k}(x,y) dS_y \right\} \quad (31) $$

By virtue of the property (5b) the first integrand is regular, while the second one, if we assume a Hölder condition analogous to (15) on the tangential derivatives of the COD, has a weak integrable singularity. The third integral, unchanged, is convergent in the ordinary sense. Then the regularization is complete if we are able to rewrite the remaining CPV integral

$$ A_{i,j,k}(x) = \int_{\Gamma_c} \Sigma_{i,j,k}^0(x,y) dS_y \\ = - \frac{1}{4\pi} \left\{ \beta^2 (\delta_{ik} I_s(x) + \delta_{jk} I_l(x) - \delta_{il} I_k(x)) + 3(1-\beta^2) J_{isk}(x) \right\} \quad (32) $$

$$ \text{with } I_p(x) = \int_{\Gamma_c} \frac{\partial r}{\partial y_p} \frac{dS_y}{r} ; \quad J_{isk}(x) = \int_{\Gamma_c} \frac{\partial r}{\partial y_i} \frac{\partial r}{\partial y_s} \frac{\partial r}{\partial y_k} \frac{dS_y}{r^2} \quad \text{and} \quad \beta^2 = \frac{\mu}{\lambda+2\mu} \quad (33) $$

in terms of regular integrals.

The integrals (33) can be integrated by parts using identity (24) [1]. Taking care of the singular nature of the integrals when we perform this transformation, we arrive at an expression using regular integrals $I_p$ and $J_{isk}$:
$$I_p(x) = \int_{\Gamma_c} \left[ \frac{1}{r} \frac{\partial r}{\partial n} - K(y) \right] n_i(y) \frac{dS_y}{r} - \int_{\partial \Gamma_c} \nu_i(y) \frac{d\Gamma_y}{r} \quad (34)$$

$$3J_{ik}(x) = \delta_{ik}I_s(x) + \delta_{sk}I_i(x) + \int_{\Gamma_c} \left[ \frac{2\partial r}{\partial y_i} \frac{\partial r}{\partial y_s} - 2n_i(y)n_s(y) \right] n_k(y) \frac{1}{r} \frac{\partial r}{\partial n} \frac{dS_y}{r} + \\ + \int_{\Gamma_c} \left[ 2n_i(y)n_s(y) - \frac{\partial r}{\partial y_i} \frac{\partial r}{\partial y_s} \right] n_k(y)K(y) \frac{dS_y}{r} - \int_{\Gamma_c} \left[ D_s^y(n_i n_k) + D_i^y(n_s n_k) \right] \frac{dS_y}{r} + \\ + \int_{\partial \Gamma_c} \left[ \nu_i(y)n_k(y)n_s(y) + \nu_s(y)n_k(y)n_i(y) - \frac{\partial r}{\partial y_i} \frac{\partial r}{\partial y_s} \nu_k(y) \right] \frac{d\Gamma_y}{r} \quad (35)$$

In the above, $K(y)$ stands for the sum of the principal curvatures $K = 1/R_1 + 1/R_3$. All integrands of (10), including the term $1/r^2 \partial r / \partial n$, are integrable (of order $O(1/r)$), assuming $\Gamma_c$ to be twice continuously differentiable (this hypothesis is necessary since (34)-(35) involve curvatures and derivatives of $n$). Finally, it should be noticed that the line integrals are regular for interior points $x \in \Gamma_c$.

The regularization (31)-(35) of the hypersingular term of (7) is complete. It remains valid without modification for a system of kinked cracks, the line integral (28) being zero. In this case, $\Gamma_c$ should be replaced by the branch $\Gamma_c^a$ on which $x$ is located, the integrals over other branches being regular.

## 4 - Regular BIE for steady-state elastodynamics

Consider the special case in which the boundary conditions (and body forces, if any) are harmonic in time with an angular frequency $\omega$, e.g.

$$t^{\pm}(x,t) = T^{\pm}(x;\omega) e^{-i\omega t} \quad (36)$$

where $(T, \text{etc...})$ are complex-valued fields. We keep in mind that whenever a formula with complex quantities appears, the physically relevant quantity is (say) its real part. We assume that the displacement may be represented as:

$$u(y,t) = u^t(y,t) + u^s(y;\omega)e^{-i\omega t} \quad (37)$$

in which $u^t$ is called the *transient part* of the solution, $u^s$ is determined so as to satisfy the boundary conditions whereas the transient part is obtained as an homogeneous solution of the elastodynamic equations under homogeneous boundary conditions with the initial conditions:

$$u^t(x,0) = u^o(x) - \text{Re}(u^s(x;\omega)) \qquad \dot{u}^t(x,0) = v^o(x) + \text{Re}(i\omega u^s(x;\omega)) \quad (38)$$

We assume that the transient part vanishes after a "sufficient" time, thus we are left with the steady-state solution. We will rewrite the regular BIEs obtained in the previous parts 2 and 3 with the time-dependence (36) (which is infinitely differentiable). The Stokes tensors become the Helmholtz tensors $U^\omega(x,y)$ (A.5) and $\Sigma^\omega(x,y)$ (A.6), with singularity and symmetry properties analogous to (2)(3)(4)(5) (see (A.11)). The regular displacement (for $x \in \Gamma$) BIE for a cracked body becomes:
$$
\begin{equation}
\begin{aligned}
& \gamma_{\delta k} u_i(x) + \int_{\Gamma} u_i(y) n_s(y) \left[ \Sigma_{i; k}^{\omega}(x, y) - \Sigma_{i; k}^{0}(x, y) \right] dS_y + \\
& \qquad + \int_{\Gamma} (u_i(y) - u_i(x)) n_s(y) \Sigma_{i; k}^{0}(x, y) dS_y \\
& \qquad - \int_{\Gamma} t_i(y) U_{i; k}^{w}(x, y) dS_y - \int_{\Gamma_c} \phi_i(y) n_s(y) \Sigma_{i; k}^{\omega}(x, y) dS_y = 0
\end{aligned}
\tag{39}
\end{equation}
$$

The regular steady-state traction BIE can be written in an equally straight-forward manner as a special case of the transient regular traction BIE (31)-(35).

## 5 - Numerical implementation of the regular boundary element method.

Numerical results concerning (a) time-domain elastodynamics and (b) curved crack analysis for our regularized BIE are not yet available. However the solution of three-dimensional elastostatic and steady-state elastodynamic problems using the regularized BIE (39) has been implemented. Several numerical tests are given in detail in [1]; they appear to be very satisfactory, giving a good accuracy with few degrees of freedom when compared to analytical solutions. Our corresponding computer code ASTRID can also deal with interior or exterior problems governed by the scalar Helmholtz equation.

We will not discuss here the whole discretization procedure of the time-domain regular BEM, which are beyond the scope of this paper, but only focus ourselves on two salient aspects of the spatial interpolation procedure closely related to the regularization techniques:

(a) The numerical counterpart of the CPV exclusion process: how to perform properly the numerical evaluation of the remaining singular integrals.

(b) the interpolation of the tangential derivatives of the COD on a general curved crack.

### 5.1 - Discretization of u and numerical evaluation of singular integrals.

The time interval of interest $[0, T]$ and the surface $\Gamma$ are divided into NT time intervals $[t_{i-1}, t_i]$ ($1 \le i \le NT$) and NEL boundary elements ($E_1, \dots, E_{NEL}$) respectively. The $E_m$ are mapped, in a standard fashion, on a reference element $E_0$, which is generally the square $(\xi_1, \xi_2) \in [-1, 1]^2$ or the triangle $0 \le \xi_1+\xi_2 \le 1$. Each component of the resulting approximate surface is continuous. The location of a point $y$ on the approximate surface is written in terms of NU shape functions $N_k$ and geometrical nodes $A_k$, NU being the total number of geometrical nodes:

$$ Oy = N_k(\xi_1, \xi_2) OA_k \quad (\xi_1, \xi_2) \in E_0 \quad (40) $$

and we choose for $u(y,t)$ an interpolation of the following type, which is isoparametric (with respect to $\Gamma$):

$$ u(y,t) = N_k(\xi_1, \xi_2) u_k(t) \quad (\xi_1, \xi_2) \in E_0; \ t \in [t_{i-1}, t_i] \quad (41) $$

In (41) $u_k(t)$ is the nodal value $u(A_k, t)$.
Evaluation of singular surface integrals During the assembly process, the sur-
face integrals are evaluated element by element for each fixed collocation point
x. Singular integrals occur for the boundary elements containing x. The CPV
exclusion procedure leads to element integrals of the form :

$$
I^2 = \int_E n_s(y) \Sigma_{i,j;k}^0 (x,y) (N_m(\hat{y}) - N_m(\hat{x})) J(\xi_1, \xi_2) d\xi_1 d\xi_2 \quad (42)
$$

where $\hat{x} = (\eta_1, \eta_2)$ is the antecedent of $x$ on the reference element $E_0$. Let $(\rho, \alpha)$ be a polar coordinate system centered at $\hat{x}$ on $E_0$ (figure ). Then for every standard polynomial shape function $N_m$ a reduced shape function $N_m^\alpha$ exists such that:

$$
N_m(\hat{y}) - N_m(\hat{x}) = \rho \hat{N}_m^\alpha(\rho, \alpha) \quad (43)
$$

The derivation of $\hat{N}_m^\alpha$ is straightforward, using for example a Taylor for-
mula (which is exact for polynomials). Furthermore $r = \|y-x\|$ can be rewrit-
ten, using (40) and (43), as :

$$
r = \rho r_0 = \rho \| \hat{N}_m^\alpha(\rho, \alpha) OA_m \| \quad (\text{summation over the nodes of E}) \quad (44)
$$

The factor 1/r₀ is always nonsingular at ρ = 0. The kernel Σº is transformed
into:

$$
\Sigma_{i,j;k}^{0}(x,y) = \frac{1}{\rho^2} \hat{\Sigma}_{i,j;k}^{0}(\rho, \alpha) \quad (\hat{\Sigma}^{0} \text{ regular at } \rho = 0) \tag{45}
$$

Finally the differential element $d\xi_1 d\xi_2$ becomes $\rho d\rho d\alpha$. This, combined with (43) and (44) allows us to recast the singular integral (42), in a completely regular form, as :

$$
I^2 = \int_E n_s(y) \hat{\Sigma}_{i,j;k}^0 (\rho, \alpha) \hat{N}_m^\alpha (\rho, \alpha) J(\xi_1, \xi_2) d\rho d\alpha \quad (46)
$$

The expression (46) takes full advantage of the CPV exclusion. The numerical evaluation of (46) can be performed with standard product Gaussian quadrature formulas, the complete procedure requiring a further coordinate change $(\rho, \alpha) \rightarrow (\nu_1, \nu_2)$ in order to recover an integral over the square $[-1,1]^2$. This involves some algebra but no major difficulty (see [1],[19]).

5.2 - Interpolation of crack-surface gradients .

To interpolate the crack-surface gradients $D_{rs}^\nu \phi_j(y)$ apply the operator $D_{rs}^\nu$ to (41)
yields :

$$
D_{rs}^{\nu} \phi_j(y,t) = \phi_{j,k}(t) D_{rs}^{\nu}(\xi_1, \xi_2) N_k(\xi_1, \xi_2) \quad (47)
$$

However, due to the discontinuity of the gradient of an interpolation such
as (41) across the element boundaries, (47) does not make sense if y is some
geometrical node $A_k$. This is the case especially when one wants to implement
the equation (31) for a collocation point x located on one mesh node, as it is
done usually. This is the reason why, extending to curved cracks the idea developed in [17] for planar cracks, we write a continuous interpolation of $D_{rs}\phi$ in terms of nodal values $(D_{rs}\phi_j)(A_k,t)$ in a form analogous to (41):

$$D_{rs}^{y}\phi_j(y) = M_k(\xi_1, \xi_2) (D_{rs}\phi_j)_k(t) \quad (48)$$

Now let us consider that the two interpolations (47) and (48) coincide in a least squares sense, thus minimizing the functional:

$$J(t) = \int_{\Gamma_c} \left[ \phi_j^k(t) D_{rs}^{y}(\xi_1, \xi_2) N_k(\xi_1, \xi_2) - (D_{rs}\phi_j)^k(t) M_k(\xi_1, \xi_2) \right]^2 dS \quad (49)$$

Taking the partial derivatives of J with respect to $(D_{rs}\phi_j)^{ik}$ leads to the matricial relationship:

$$[MM]_{pk} (D_{rs}\phi_j)^k = [MN]_{pk} \phi_j^k \quad (50a)$$

$$\text{with } [MM]_{pk} = \int_{\Gamma_c} M_p M_k \, dS \quad \text{and} \quad [MN]_{pk} = \int_{\Gamma_c} M_p D_{rs}^{y} N_k \, dS \quad (50b)$$

which gives the nodal values of the tangential gradients of $\phi_j$ as a linear function of the nodal values of $\phi_j$. The matrix equation (50) expresses a purely geometrical fact and is therefore independent of the time and of the component index j of $\phi$. We then have to compute the four matrices $[MM]$, $[MN]_{12} = -[MN]_{21}$, $[MN]_{23} = -[MN]_{32}$, $[MN]_{31} = -[MN]_{13}$ and invert one of them, namely $[MM]$ (which has a banded structure). A close examination reveals that this computational task is small compared to the solution of the complete crack problem.

**Interpolation near the crack front** The surface displacement gradients are expected to vary like $d^{-1/2}$ near the crack front, d being the distance between a point y and its orthogonal projection on the crack edge. This behaviour can be modelled using a 9-node Lagrange element with quarter-point nodes. The element jacobian is then proportional to $d^{1/2}$. This allows one to recover directly the local stress intensity factors [17].

## 6 - Concluding remarks

We have presented a unified BIE formulation entirely free of Cauchy principal value integrals and covering a very wide class of linear elasticity problems: the solids in consideration can be bounded or unbounded, cracked or uncracked, under dynamic or static loading. Our main results (21) and (31)-(35) contain and generalize all attempts which have been made in the past for the regularization of classical direct collocation BIE [7][17][18][19][20]. The numerical implementation of this regular BIE

(a) is done and successfully tested for uncracked solids in static and steady-state elasticity.

(b) is under progress for cracked solids in static and steady-state elasticity, and for time-domain elastodynamics.

The present approach can be developed the same way for two-dimensional elastodynamics and to three- or two-dimensional wave problems governed by the scalar wave equation.
# Appendix 1 - Elastodynamic kernels

## A1.1 - Some preliminary relations

$$c_L = \sqrt{\frac{\lambda + 2\mu}{\rho}} \quad c_T = \sqrt{\frac{\mu}{\rho}} \quad \text{(longitudinal and transverse wave velocities)}$$

$$\beta = \frac{c_T}{c_L} = \frac{k_L}{k_T} = \sqrt{\frac{1 - 2\nu}{2(1 - \nu)}} \quad \text{with} \quad k_{\alpha} = \frac{\omega}{c_{\alpha}} \quad (\alpha = L \text{ or } T)$$

$$r^2 = \|x-y\|^2 = (y_3-x_3)^2 + (y_2-x_2)^2 + (y_1-x_1)^2$$

$$r_i = \frac{\partial r}{\partial y_i} = \frac{y_i - x_i}{r} \quad r_{ij} = 1 \quad r_{jj} = \frac{1}{r} (\delta_{ij} - x_i r_{,j}) \quad \text{(A.1)}$$

$$r_{ijk} = \frac{1}{r^2} ( 3 r_i r_j r_k - \delta_{ij} r_k - \delta_{ik} r_j - \delta_{jk} r_i )$$

## A1.2 - Stokes displacement tensor U and associated stress tensor Σ (time-dependent case)

The Stokes tensors associated to the concentrated body force of (twice continuously differentiable) time-dependent magnitude f(t) acting on x along e_k-direction are given by:

$$U_{ijk}(x,t,y|f) = \frac{1}{4\pi\mu r} \left[ \text{Br}_{i}r_{,k} + (A + f(t - \frac{r}{c_{T}}))\delta_{ik} \right] \quad \text{(A.2)}$$

$$\Sigma_{ijk}(x,t,y|f) = \frac{1}{4\pi r^2} \left[ 2\text{Cr}_{i}r_{,k}r_{,s} + \delta_{is}r_{,k}F_{L} + (\delta_{ik}r_{,s} + \delta_{sk}r_{,i})F_{T} \right] \quad \text{(A.3)}$$

with the auxiliary quantities:

$$A = c_T^2 \int \frac{1}{1/c_T} \lambda f(t - \lambda r) \, d\lambda$$

$$B = -3A - \left[ f(t - \frac{r}{c_T}) - \beta^2 f(t - \frac{r}{c_L}) \right]$$

$$C = -5B + \left[ f(t - \frac{r}{c_T}) + \frac{r}{c_T} f(t - \frac{r}{c_T}) \right] - \beta^2 \left[ f(t - \frac{r}{c_L}) + \frac{r}{c_L} f(t - \frac{r}{c_L}) \right] \quad \text{(A.4)}$$

$$F_T = 2B - \left[ f(t - \frac{r}{c_T}) + \frac{r}{c_T} f(t - \frac{r}{c_T}) \right]$$

$$F_L = 2B - (1 - 2\beta^2) \left[ f(t - \frac{r}{c_L}) + \frac{r}{c_L} f(t - \frac{r}{c_L}) \right]$$

## A1.3 - Helmholtz displacement tensor U^w and associated stress tensor Σ^w (steady-state case)

$$U_{ijk}^{w}(x,y) = \frac{1}{4\pi\mu r} \left[ \text{Br}_{i}r_{,k} + (A + e^{ik \cdot T})\delta_{ijk} \right] \quad \text{(A.5)}$$

$$\Sigma_{ijk}^{w}(x,y) = \frac{1}{4\pi\mu r^2} \left( 2\text{Cr}_{i}r_{,k}r_{,s} + \delta_{is}r_{,k}F_{L} + (\delta_{ik}r_{,s} + \delta_{sk}r_{,i})F_{T} \right) \quad \text{(A.6)}$$
with the auxiliary quantities :

$$
\begin{align*}
A &= A_T e^{ik_T r} - \beta^2 A_L e^{ik_L r} & A_\alpha &= \frac{i}{k_\alpha r} - \frac{1}{(k_\alpha r)^2} \\
B &= B_T e^{ik_T r} - \beta^2 B_L e^{ik_L r} & B_\alpha &= -3A_\alpha - 1 \\
C &= C_T e^{ik_T r} - \beta^2 C_L e^{ik_L r} & C_\alpha &= -5B_\alpha - (k_\alpha r)^2 A_\alpha \tag{A.7} \\
F_T &= (k_T r)^2 A_T e^{ik_T r} + 2B & F_L &= (1-2\beta^2) (k_L r)^2 A_L e^{ik_L r} + 2B
\end{align*}
$$

A1.4 - Kelvin displacement tensor U$^a$ and associated stress tensor Σ$^0$ (static case)

$$
U_{ij}^{0}(x,y) = \frac{1}{4\pi\mu r} \left[ \frac{1-\beta^2}{2} r_i r_j + \frac{1+\beta^2}{2} \delta_{ik} \right] \quad (A.8)
$$

$$
\Sigma_{i,j,k}^{0}(x,y) = \frac{1}{4\pi\mu r^2} \left( 3(1-\beta^2)r_i r_j r_k + \beta^2(\delta_{ik}r_j + \delta_{jk}r_i - \delta_{ki}r_j) \right) \quad (A.9)
$$

A1.5 - Taylor expansions of the dynamic kernels for small r

(a) Stokes tensors .

$$
A = -\frac{1}{2}(1-\beta^2)f(t) + \frac{1}{3}(1-\beta^3)\frac{r}{c_T}\dot{f}(t) - \frac{1}{8}(1-\beta^4)\frac{r^2}{c_T^2}\ddot{f}(t) + o(r^2) \quad (A.10)
$$

$$
A + f\left(t - \frac{r}{c_T}\right) = \frac{1}{2}(1+\beta^2)f(t) - \frac{1}{3}(2+\beta^3)\frac{r}{c_T}\dot{f}(t) + \frac{1}{8}(3+\beta^4)\frac{r^2}{c_T^2}\ddot{f}(t) + o(r^2)
$$

$$
B = \frac{1}{2}(1 - \beta^2)f(t) - \frac{1}{8}(1 - \beta^4)\frac{r^2}{c_T^2}\ddot{f}(t) + o(r^2)
$$

$$
C = -\frac{3}{2}(1-\beta^2)f(t) + \frac{1}{8}(1-\beta^4)\frac{r^2}{c_T^2}\ddot{f}(t) + o(r^2)
$$

$$
F_T = -\beta^2 f(t) + \frac{1}{4}(1+\beta^4)\frac{r^2}{c_T^2}\ddot{f}(t) + o(r^2)
$$

$$
F_L = \beta^2 f(t) - \frac{1}{4}(1 - 2\beta^2 + 3\beta^4)\frac{r^2}{c_T^2}\ddot{f}(t) + o(r^2)
$$

The expansions (A.10) show that :

(i) the lowest-order term yields $f(t)$ times the singular static Kelvin kernel (compare with (A.9)).

(ii) the differences $U[x,t,y|f] - f(t)U^0(x,y)$ and $\Sigma[x,t,y|f] - f(t)\Sigma^0(x,y)$ remain finite for arbitrary small $r$. Thus the relations (5a-5b).

(b) Helmholtz kernels .

In this case the force function $f(t) = \exp(-i\omega t)$ is $C^\infty$. The Taylor expansion thus becomes an infinite series expansion :

$$
A = -\frac{1-\beta^2}{2} - \sum_{p \ge 1} \frac{(1-\beta^{p+2})}{p!(p+2)} (ik_T r)^p
$$

$$
A + e^{ik_T r} = \frac{1+\beta^2}{2} + \sum_{p \ge 1} \frac{(1+p+\beta^{p+2})}{p!(p+2)} (ik_T r)^p
$$
$$B = \frac{1-\beta^2}{2} + \sum_{p \ge 1} \frac{(1-\beta^{p+2})(1-p)}{p!(p+2)} (ik_T r)^p$$

$$C = -\frac{3}{2}(1-\beta^2) + (k_T r)^2 \sum_{p \ge 0} \frac{(1-\beta^{p+4})p-1}{p!(p+2)(p+4)} (ik_T r)^p$$

$$F_T = -\beta^2 - (k_T r)^2 \sum_{p \ge 0} \frac{(p+2+2\beta^{p+4})}{p!(p+2)(p+4)} (ik_T r)^p \quad (A.11)$$

$$F_L = \beta^2 + (k_T r)^2 \sum_{p \ge 0} \frac{(2 - (p+4)\beta^{p+2} + 2(p+3)\beta^{p+4})}{p!(p+2)(p+4)} (ik_T r)^p$$

A term-by-term comparison of the previous expansions with the series expansion of $exp(ik_T r)$ allows one to prove their absolute convergence for every $k_T r$. The expansions (A.11) show that:

(i) the lowest-order term yields the singular static Kelvin kernel (compare with (A.8) and (A.9)).

(ii) the differences $U^w - U^0$ and $\Sigma^w - \Sigma^0$ remain finite for arbitrary small $r$ and tend towards zero like $w$ and $w^2$ respectively for small $w$. Thus the relations (5a-5b).

# Appendix 2 - Shape functions

For brevity we restrict ourselves to the shape functions associated with the 8-node quadrilateral element but the modified shape functions $\hat{N}_m^x$ can be derived as well for any polynomial interpolation.

## A2.1 - Geometry / displacement shape functions

$$ \begin{aligned} N_1(\xi_1, \xi_2) &= 1/4 (1-\xi_1)(1-\xi_2)(-1-\xi_1-\xi_2) & N_2(\xi_1, \xi_2) &= 1/2 (1-\xi_1^2)(1-\xi_2) \\ N_3(\xi_1, \xi_2) &= 1/4 (1+\xi_1)(1-\xi_2)(-1+\xi_1-\xi_2) & N_4(\xi_1, \xi_2) &= 1/2 (1-\xi_2^2)(1+\xi_1) \end{aligned} \quad (A.12) $$

$$N_5(\xi_1, \xi_2) = 1/4 (1+\xi_1)(1+\xi_2)(-1+\xi_1+\xi_2) \quad N_6(\xi_1, \xi_2) = 1/2 (1-\xi_1^2)(1+\xi_2)$$

$$N_7(\xi_1, \xi_2) = 1/4 (1-\xi_1)(1+\xi_2)(-1-\xi_1+\xi_2) \quad N_8(\xi_1, \xi_2) = 1/2 (1-\xi_2^2)(1-\xi_1)$$

## A2.2 - Reduced shape functions $\hat{N}_m^x(\rho, \alpha)$

$$ \hat{N}_1^x(\rho, \alpha) = 1/4[\cos\alpha(b_1-a_1)a_2 + \sin\alpha(b_1-a_2)a_1 + \rho(a_3\cos^3\alpha + a_1\sin^2\alpha + (a_1+a_3-b_1)\cos\alpha\sin\alpha) - \rho^2\cos\alpha\sin\alpha(\cos\alpha+\sin\alpha)] \quad (A.13a) $$

$$ \hat{N}_2^x(\rho, \alpha) = 1/2[\cos\alpha(a_1-a_3)a_2 - a_1a_3\sin\alpha + \rho((a_3-a_1)\cos\alpha\sin\alpha - a_2\cos^2\alpha) + \rho^2\cos^2\alpha\sin\alpha] $$

with

$a_1 = 1 - \eta_1$

$a_3 = 1 + \eta_1$

$b_3 = -1 - \eta_1 - \eta_2$

$b_3 = -1 + \eta_1 + \eta_2$

$a_2 = 1 - \eta_2$

$a_4 = 1 + \eta_2$

$b_4 = -1 + \eta_1 - \eta_2$

$b_4 = -1 + \eta_1 - \eta_2$

(A.13b)

To obtain $\hat{N}_3^x$, $\hat{N}_4^x$ replace $(a_3, a_3, a_3, b_3)$ by $(a_3, a_3, a_4, b_3)$ and $(\cos\alpha, \sin\alpha)$ by $(\sin\alpha, -\cos\alpha)$ in $\hat{N}_1^x$, $\hat{N}_2^x$, and so on by circular permutation.
References

1 BONNET M. - Méthode des équations intégrales régularisées en élastodynamique tridimensionnelle. *Thèse de doctorat, publiée dans le Bulletin EDF/DER série C, n° 1/2, 1987.*

2 BREBBIA C.A., TELLES J.C.F., WROBEL L.C. - Boundary element techniques. Theory and application in engineering. Springer - Verlag, 1984.

3 BREBBIA C.A. (Editor) - Topics in boundary element research, vols 1 & 2. Springer - Verlag 1984.

4 BUI H.D. - Sur les équations intégrales en thermoélastoplasticité tridimensionnelle. *Compte Rendus de l'Académie des Sciences, tome 286, 1978.*

5 BUI H.D. - An integral equation method for solving the problem of a plane crack of arbitrary shape. *J. Mech. Phys. Solids, vol. 25, 1977, p.29-39.*

6 BUI H.D. - New trends in BIE methods. *Computational Mechanics '86, Tokyo 1986, Springer Verlag, vol 2, pp X1.3-8.*

7 BUI H.D., LORET B., BONNET M. - Régularisation des équations intégrales de l'élastodynamique et de l'élastostatique. *Compte Rendus de l'Académie des Sciences, tome 300, 1985.*

8 ERINGEN A.C., SUHUBI E.S. - Elastodynamics (vol II - linear theory). Academic Press 1975.

9 GUNTHER N.M. - La théorie du potentiel et ses applications aux problèmes fondamentaux de le physique mathématique. *Paris, Gauthier-Villars, 1934.*

10 KOBAYASHI S. - Fundamentals of boundary integral equations methods in elastodynamics. In [5] vol. 2.

11 KUPRADZE V.D. - Dynamical problems in elasticity. *Progress in solids mechanics, vol. III, North Holland, 1963.*

12 KUPRADZE V.D. (Editor) - Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity. *North Holland, 1979.*

13 LACHAT J.C., WATSON J.O. - Effective numerical treatment of boundary integral equations: a formulation for three-dimensional elastostatics. *Int. J. Num. Meth. in Eng., vol. 10, 1976, pp.991-1005.*

14 LE VAN A. - Equations intégrales pour fissures 3D. Contribution à la régularisation. *Thèse ENSM Nantes, 1988.*

15 MANOLIS G., BESKOS D.E. - Boundary element methods in elastodynamics. *Umwin Hyman, London 1988.*

16 NEDELEC J.C. - Approximation des équations intégrales en mecanique et en physique. *Cours Ecole Polytechnique, 1977.*

17 POLCH E.Z., CRUSE A.T., HUANG C.J. - Traction BIE solutions for flat cracks.

18 RIZZO F.J., SHIPPY D.J. - An advanced boundary integral equation method for three-dimensional elasticity. *Int. J. Num. Meth. in Eng., vol. 11, 1977, p.1753.*

19 RIZZO F.J., SHIPPY D.J., REZAYAT M. - A boundary integral equation method for radiation and scattering of elastic waves in three-dimensions. *Int. J. Num. Meth. in Eng., vol. 21, 1985, p.115-129.*

20 SLADEK J., SLADEK V. - Dynamic stress-intensity factors studied by boundary integro-differential equations. *Int. J. Num. Meth. in Eng., vol. 23, 1986, p.425-436.*

21 WENDLAND W.L. - Asymptotic accuracy and convergence for point collocation methods. *in [5] vol 2*
<figure><img src="image_1.png" /><figcaption>Figure 1: a multiply-connected bounded body (left) an unbounded body with a cavity (with edges) and a crack (right)</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2: an indented half-space</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: orientation convention for a crossed crack.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4: the eight-node Serendip element and the reduced polar coordinate system (α, α).</figcaption></figure>