# Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods

Marc Bonnet, Huy Duong Bui

► To cite this version:

Marc Bonnet, Huy Duong Bui. Regularization of the displacement and traction BIE for 3D elastodynamics using indirect methods. Advances in Boundary Element Techniques, pp.1-29, 1993, 10.1007/978-3-642-51027-4_1. hal-00112506

HAL Id: hal-00112506

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

Submitted on 1 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Regularization of the Displacement and Traction BIE for 3D Elastodynamics Using Indirect Methods

Marc BONNET

Laboratoire de Mécanique des Solides, CNRS, Ecole Polytechnique,
91128 Palaiseau Cedex, FRANCE.

Huy Duong BUI

Electricité de France, 92141 Clamart Cedex, FRANCE and
Laboratoire de Mécanique des Solides, CNRS, Ecole Polytechnique,
91128 Palaiseau Cedex, FRANCE.

## Abstract

This paper deals with regularization techniques developed in order to overcome the strongly singular (displacement BIE) or hypersingular (traction BIE) character of the boundary integral equations commonly used for three-dimensional elastodynamics. More specifically, we address in-direct regularization techniques, which rely upon singularity exclusion and separate evaluation of strongly singular or hypersingular integrals containing the singular kernel only. In many, but not all, cases, this evaluation can be made using an auxiliary static problem together with the fact that the dynamic and static kernels share the same singular term. In this paper we first give a bibliographical outline of the regularization problem. Then the regularized displacement and traction BIE for general 3D situations in transient elastodynamics, including crack problems, are derived and stated. We also recall related results of interest. Then numerical implementation considerations are considered. For completeness, a brief survey of Galerkin-type and direct approaches is also given. Finally some numerical examples in elastodynamics are given.

## 1 Introduction.

The boundary element method (BEM) is now well-established and applied in numerous engineering fields (structural analysis, geomechanics, acoustics, among others) [9], [10], [4], [26], [48]. However the integral equations of elastodynamics or acoustics involve strongly singular or hypersingular integrals which, in classical elastic potential theory, are interpreted as Cauchy principal value (CPV) (Kupradze [31], [32]) or Hadamard finite part (FP). As pointed out in [35], [11], [12], the usual techniques of integral calculus (differentiation under the integral sign, change of variable) applied to CPV or FP integrals may give false results. This has been a serious obstacle to a satisfactory numerical implementation of BEM because these singular integrals correspond to near-diagonal terms in the matrix of the boundary element discretized problem, and hence the accuracy of the numerical BEM solution depends strongly on an accurate evaluation of the singular integrals. Indeed this singular character is also an advantage, since the strong singularity of the kernels, if accurately evaluated, leads to a well-conditioned discrete problem.

The regularization of the singular integral equations arising in elastic or acoustic potential theory is an old problem. Giraud [19], then Mikhlin [35] considered multidimensional singular integral equations of the form

$$ (I - \lambda K)\phi = f \quad \text{with} \quad (K\phi)(x) = \int_S \phi(y)K(x, y)dS_y \quad (x, y \in S) \quad (1) $$

with unknown density $\phi(y)$ and a kernel $K(x, y)$ integrable in the CPV sense. They state the regularization problem as follows: find the singular kernel $K'(x, z, \lambda)$ such that

$$ (I + \lambda K')(I - \lambda K)\phi = (I + \lambda K')f \quad (2) $$
Combining Giraud and Mikhlin methods, Kupradze [31] exhibits a construction of the singular part of integral operator K'. He is then able to prove that the regularized form (2) of the singular integral equation is a Fredholm integral equation of the second kind, and hence to state the existence and uniqueness of the solutions of the singular integral equations associated to the basic elasticity problems. In this sense, the singular integral equations arising in elastic or scalar potential theory are well-posed.

However, this approach of regularization, though very fruitful for fundamental considerations, is not amenable to efficient numerical methods for solving singular multidimensional integral equations, the construction of operator K' being very complicated for arbitrary surfaces. For this reason other regularization approaches appeared, better suited to numerical implementation: indirect approach, variational approach and direct approach.

The present paper is devoted to the application of indirect approach to the regularization of displacement and traction boundary integral equations (BIE), with emphasis on the latter. Basically, the indirect approach consists in a singularity isolation: $\phi(\mathbf{x})$ is subtracted and added to $\phi(\mathbf{y})$ in the operator K:

$$ (K\phi)(\mathbf{x}) = \int_S [\phi(\mathbf{y}) - \phi(\mathbf{x})] K(\mathbf{x}, \mathbf{y}) dS_y + \phi(\mathbf{x}) \int_S K(\mathbf{x}, \mathbf{y}) dS_y $$

For a strongly singular integral operator K, the first integral above is weakly singular, while the second one remains to be evaluated by some means. This idea itself is extensively used in classical potential theory. In the BEM literature, indirect regularization of collocation elastostatic BIE appears in Rizzo, Shippy [41]; the isolated CPV is evaluated (or, more precisely, its actual computation is avoided) by means of a rigid-body identity. This idea has been extended for broader classes of problems governed by displacement boundary integral equations (DBIE): frequency-domain elastodynamics (Bui et al. [13], Bonnet [5], Rizzo et al. [42]) and time-domain elastodynamics [6]. During the same period, the need of numerical solutions for (hypersingular) traction BIEs became manifest, especially for crack problems. The first step was to convert these hypersingular BIEs into CPV strongly singular BIEs, the unknown density becoming the tangential gradient of the crack opening displacement (Bui [12], Weaver [47] for planar cracks, Sladek & Sladek [44] for curved cracks). The singularity isolation and evaluation in this case cannot be done using a rigid-body identity, and the final singularity isolation relies upon an analytical treatment of the residual integral followed by a limiting process (Bonnet [5], [6]). Regularization methods for hypersingular BIEs for scalar potential theory and frequency-domain elastodynamic crack problems are also developed by Krishnasamy et al. [29], [30], Nishimura & Kobayashi [37]. Leblond [33] considers the second-order regularization for 2D elastostatic crack problems.

The present paper focuses on the regularization, using the indirect approach, of the (hypersingular) traction BIE (TBIE), or, equivalently, of the gradient BIE (GBIE), for time-domain elastodynamics. The TBIE are classically used for the BIE modelling of cracks in linear brittle fracture mechanics [27], [12], [14], [16], [24], [37], [38], [45], [47], owing to the well-known degeneracy of the DBIE applied to cracks. However, the TBIE may also prove useful for the study of non-cracked solids:

* On portions of the boundary where displacements are given (Dirichlet boundary conditions), the classical DBIE is of the 'first kind' structure and shows some ill-conditioning. The use of a TBIE for the collocation points located on the Dirichlet part of the boundary restores the mathematically desirable 'second-kind' structure and may lead to a better conditioned numerical solution algorithm.

* Traction BIEs may be used to compute the complete stress tensor on the boundary, either for stress analysis purposes or for the computation of residual-type error indicators.

The purpose of the present paper is twofold:
* To establish and state the regularized TBIE for general elastodynamic problems, involving bounded as well as unbounded elastic media. The case of a crack embedded in an infinite medium is treated separately.

* To show, via the regularization process being used, that the results do not depend on the limiting process used to define the TBIE, i.e. the shape of the exclusion neighbourhood around the (singular) collocation point, in contrast with other expositions of the subject. This gives both a rigorous proof of the results and, in the authors' opinion, a better understanding of the singularities involved and of the nature of the resulting regularized TBIE.

In accordance with the latter consideration, the TBIE here are introduced and investigated without reference to a specific limiting process such as finite-part integrals. The BIE are not treated as limiting forms of interior representations when the collocation point is moved to the boundary $\partial\Omega$ of the domain $\Omega$ under consideration (as was earlier done for cracks in [5], [6]). Instead they are viewed as limiting forms of exterior representations when a vanishing exclusion neighbourhood is removed around the collocation point, the latter being kept fixed on $\partial\Omega$ during the regularization process. This presentation viewpoint follows Guiggiani et al. [21].

The regularization of time-domain (strongly singular) elastodynamic DBIE is also treated here, both for the sake of reference and in order to introduce some ideas and notations of later use.

The regularization of static BIEs plays a central role, since the static and dynamic kernels possess the same singularities.

In this paper, after some preliminaries and definitions (section 2), the indirect regularization of static DBIE is treated (section 3). Then comes the regularization of static TBIE, which constitutes the bulk of the present paper. First uncracked elastic media, either bounded or unbounded, are considered and two different regularizations strategies are developed: second-order regularization (section 5) and integration by parts followed by a first-order regularization (section 6). Then, in section 7, the former strategy is applied to the TBIE associated to a curved crack. A key preliminary for sections 6 and 7 is the introduction of tangential differential operators, integration by parts formulas for curved surfaces and integral identities involving the Kelvin static kernels, which is done in section 4. Next, the regularized DBIE and TBIEs are stated for time-domain elastodynamics and discussed in section 8. This is followed by a discussion on the actual computation of the resulting (weakly singular) element integrals using BEM discretization (section 9) and a brief overview of other approaches for dealing with hypersingular BIE, namely the variational and the direct approach (section 10). Finally two numerical examples illustrating the use of regularized DBIE and TBIE in elastodynamics are presented in section 11.

## 2 Preliminaries.

We consider the dynamic response of a 3D elastic body $\Omega$ (either bounded or unbounded (see figure 1), subjected to surface loadings. Since the regularization procedures would be unaffected by the presence of nonzero body forces and initial conditions, which appear in displacement and traction BIE as weakly singular domain integrals only, zero body forces and initial conditions are considered here. The displacement field $u(y,t)$ ($y \in \Omega$ and $t \in T^+ = [0, +\infty]$) is then governed by the homogeneous elastodynamic equilibrium equation together with Hooke's law:

$$ C_{isab}u_{a,b}s(y,t) - \rho\ddot{u}(y,t) = 0 \qquad (3) $$

$$ C_{ispq} = \lambda\delta_{is}\delta_{pq} + \mu(\delta_{ip}\delta_{sq} + \delta_{iq}\delta_{ps}) \quad (\lambda, \mu: \text{Lamé constants}) \qquad (4) $$

together with appropriate boundary conditions and, if $\Omega$ is unbounded, classical elastodynamic radiation conditions [31]. The comma, as in (3), indicates partial differentiation with respect to the components of $y$. Einstein summation convention is used, unless explicitly stated otherwise, throughout the present paper.
<figure><img src="image_2.png" /><figcaption>Figure 1: Geometrical notations and conventions.</figcaption></figure>

The boundary of $\Omega$ is denoted $\Gamma$. The surface $\Gamma$ represents the boundary of either the finite body and/or the internal cavities or inclusions, according to the geometry under consideration. 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 several disjoint surfaces (e.g. a bounded solid with cavities).

The case of a crack in an infinite elastic space will be considered separately in section 7.

## 2.1 Stokes tensors $U$ and $\Sigma$ [18].

Let $x \in \mathbb{R}^3$ and $f$ be a twice continuously differentiable function of time that vanishes for $t < 0$. Let $U_i^k[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 components $U_i^k$ define a second order tensor: the Stokes' displacement tensor $U$, which possesses the symmetry properties:

$$ \begin{gathered} U_i^k[x, t, y|f(t)] = U_i^k[y, t, x|f(t)] = U_i^k[x, t, y|f(t)] \\ \frac{\partial}{\partial x_r} U_i^k[x, t, y|f(t)] = -\frac{\partial}{\partial y_r} U_i^k[x, t, y|f(t)] \end{gathered} \quad (5) $$

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

$$ \Sigma_{is}^k[x, t, y|f(t)] = C_{ispq} U_{p,q}^k[x, t, y|f(t)] \quad (6) $$

In equation (6) and throughout the present paper, the comma used with two-point kernels denotes differentiation with respect to the field point $y$. The analytical expressions of $U[x, t, y|f(t)]$ and $\Sigma[x, t, y|f(t)]$ in the 3-D case are given by formulas (95), (96) in Appendix A.

The special cases $f(t) = e^{-i\omega t}$ and $f(t) = 1$ yield respectively the time-harmonic tensors (known as the Helmholtz tensors) $\Sigma(x, y; \omega)$, $U(x, y; \omega)$ and the static tensors (known as the Kelvin tensors) $\Sigma(x, y)$ and $U(x, y)$. The expressions of Kelvin tensors in the 3-D case are given by formulas (93), (94) respectively, in Appendix A.

The Stokes and Kelvin tensors satisfy the (dynamical) equilibrium equation:

$$ \Sigma_{is,s}^k[x, t, y|f(t)] - \rho U_i^k[x, t, y|\tilde{f}(t)] = -\delta(y-x)f(t) \quad (7) $$

$$ \Sigma_{is,s}^k(x, y) = -\delta(y - x) \quad (8) $$

Let $r = \|x - y\|$ be the euclidean distance between $x$ and $y$. The Stokes tensors (and their time-harmonic and static counterparts as well) exhibit a well-known singular behaviour for $r$ arbitrarily small:

$$ U_i^k[x, t, y|f(t)] = O(1/r) \qquad U_{i,r}^k[x, t, y|f(t)] = O(1/r^2) \qquad \Sigma_{is}^k[x, t, y|f(t)] = O(1/r^2) \quad (9) $$
## 2.2 Basic integral identities for uncracked domain $\Omega$.

Let $z$ be a fixed point, either interior or exterior to $\Omega$ (ie. not on the boundary $\Gamma$). The integral boundary representation theorem for time-domain elastodynamics reads [18]:

$$ \kappa 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 = 0 \quad (10) $$

where $\kappa = 1$ ($z$ interior to $\Omega$) or $\kappa = 0$ ($z$ exterior to $\Omega$). Equation (10) stems from an 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) and a subsequent time-convolution. The latter is performed analytically ([18]) and results in kernels such as $\Sigma_{is}^k[z, t, y|u_i(y, t)]$, where $u_i(y, t)$ is substituted to the force function $f(t)$ in (93), (94).

As $z$ is not on the boundary, the integrals in (10) are $C^\infty$ functions of $z$. In particular, they may be differentiated with respect to $z_r$, yielding:

$$ \kappa u_{k,r}(z, t) - \int_{\Gamma} n_s(y) \Sigma_{is,r}^{k}[x, t, y|u_i(y, t)] dS_y + \int_{\Gamma} U_{i,r}^{k}[x, t, y|t_i(y, t)] dS_y = 0 \quad (11) $$

In (11), use has been made of (5), which allows the exchange of differentiations with respect to the source point $z$ and the integration point $y$. The above identity (11) yields an integral representation of the displacement gradients. The interior stress tensor $\sigma(z)$ may then be obtained in terms of the boundary elastic fields using Hooke's law (4).

## 2.3 Definition of BIE via a limiting process.

In view of the singular character (9) of the fundamental tensors for $x$ arbitrarily close to $y$, a limiting process of some nature is necessary if one is to derive *boundary integral equations* from the identities (10), (11).

Let $x$ be a fixed point on $\Gamma$. Following Guiggiani et al. [21] and using the same notations, a neighbourhood $v_\epsilon(x) = v_\epsilon$ of $x$ is removed from $\Omega$, defining the domain $\Omega_\epsilon$. The neighbourhood $v_\epsilon$ has a size of the same order as $\epsilon$, and hence vanishes with $\epsilon$. The common practice is to take for $v_\epsilon$ the sphere of radius $\epsilon$ centered at $x$ and to interpret the subsequent limiting processes as Cauchy principal values (CPV) of Hadamard finite parts (FP). On the contrary, the shape of $v_\epsilon$ throughout the present paper is arbitrary, in order to show that the final (regularized) BIEs does not depend on a specific limit process.

The following notations are introduced (see figure 2): $s_\epsilon = \Omega \cap \partial v_\epsilon$, $e_\epsilon = \partial \Omega \cap v_\epsilon$ (so that $\partial \Omega_\epsilon = (\Gamma - e_\epsilon) + s_\epsilon$) and $c_\epsilon = \partial e_\epsilon$. The direct collocation displacement boundary integral equation (DBIE) and gradient boundary integral equation (GBIE) are stated as the limit for $\epsilon \to 0$ of:

$$ \int_{(\Gamma - e_{\epsilon}) + s_{\epsilon}} \{n_s(y) \Sigma_{is}^{k}[x, t, y|u_i(y, t)] - U_i^k[x, t, y|t_i(y, t)]\} dS_y = 0 \quad (12) $$

$$ \int_{(\Gamma - \epsilon_\epsilon) + s_\epsilon} \{n_s(y) \Sigma_{is,r}^k[x, t, y|u_i(y, t)] - U_{i,r}^k[x, t, y|t_i(y, t)]\} dS_y = 0 \quad (13) $$

respectively (n denoting the unit normal of $\partial \Omega_\epsilon$ exterior to $\Omega_\epsilon$. Equations (12), (13) above are identities (10) and (11) applied to the domain $\Omega_\epsilon$. Indeed the regularization process will show the very existence of the limit, which is not a priori obvious.

## 2.4 Transfer of the singularity of Stokes' tensors into Kelvin tensors.

The Taylor expansion of the Stokes' tensors for small $r$ (equation (97) in Appendix A) reveals the following crucial properties:

$$ U_i^k[x, t, y|f(t)] - f(t)U_i^k(x, y) = \dot{f}(t)O(1) \quad (14) $$

$$ U_{i,r}^k[x, t, y|f(t)] - f(t)U_{i,r}^k(x, y) = \ddot{f}(t)O(1) \quad (14) $$

$$ \Sigma_{is}^k[x, t, y|f(t)] - f(t)\Sigma_{is}^k(x, y) = \dddot{f}(t)O(1) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Exclusion neighbourhood $v_\epsilon$ and related notations used for the limiting process (12), (13).</figcaption></figure>

In other words the singularity (9) of the dynamic tensors is entirely included into their static counterparts U and $\Sigma$. This property is also valid for the Helmholtz tensors. In view of properties (14), the singularity of the dynamic kernels $\Sigma_{i\sigma}^k$, $\Sigma_{i\sigma,r}^k$, $U_{i,r}^k$ in (12), (13) can be transfered in integrals involving the static Kelvin tensors. Accordingly, the regularization of static DBIE and GBIE will be investigated first (sections 3, 5, 6, 7). Then the results for time-domain elastodynamics will be stated in section 8 and discussed.

### 3 Regularized elastostatic displacement BIE.

The static DBIE is stated as the limiting case, for $\epsilon \to 0$, of:

$$ \lim_{\epsilon \to 0} \int_{(\Gamma - e_{\epsilon}) + s_{\epsilon}} \{u_i(y)n_s(y)\Sigma_{i\sigma}^k(x, y) - t_i(y)U_i^k(x, y)\} dS_y = 0 \quad (15) $$

#### 3.1 Rigid-body identity.

Let $\epsilon > 0$ take a temporarily fixed value, and $D_\epsilon$ denote the *bounded* domain with boundary $\partial\Omega_\epsilon$ (with *outward* unit normal $n^D$, consistently with the conventions used herein). Thus:

* if $\Omega_\epsilon$ is bounded: $D_\epsilon = \Omega_\epsilon$, $n^D = n$, $x$ is exterior to $D_\epsilon$.

* if $\Omega_\epsilon$ is unbounded: $D_\epsilon = \mathbb{R}^3 - \Omega_\epsilon$, $n^D = -n$, $x$ is interior to $D_\epsilon$.

Consider, as an auxiliary solution, a rigid displacement of $D_\epsilon$ defined by:

$$ u(y) = u^0(\text{constant}), \quad t(y) = 0 \quad \text{on } \partial D_\epsilon \quad (16) $$

This rigid-body motion satisfies the elastostatic equations, thus equation (10) holds true for this case, upon substitution of $n$ by $n^D$. Using the unit normal $n$ exterior to $\Omega_\epsilon$, i.e. reversing the normal in the case $\Omega$ unbounded, it reads:

$$ -\kappa u_k^0 + u_i^0 \int_{(\Gamma - e_\epsilon) + s_\epsilon} n_s(y) \Sigma_{i\sigma}^k(x, y) dS_y = 0 \quad \text{with} \quad \begin{cases} \kappa = 1 & \Omega \text{ unbounded} \\ \kappa = 0 & \Omega \text{ bounded} \end{cases} \quad (17) $$

The meaning of $\kappa$ in (17) is consistent with the notation introduced in subsection 2.2, this significance is maintained in the sequel. The above considerations hold true for $x$ being an edge or corner point as well as a regular point of $\partial\Omega$. Since the identity (10) for an infinite domain $\Omega$ implies that the displacement and stress fields satisfy decay conditions at infinity, the rigid-body displacement had to be considered for a *bounded* domain, hence the introduction of $D_\epsilon$, and the subsequent appearance of the coefficient $\kappa$ in (17).
## 3.2 Regularization of static DBIE.

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

$$ \kappa u_k(\mathbf{x}) + \int_{\Gamma - e_{\epsilon}} \left\{ \left( (u_i(\mathbf{y}) - u_i(\mathbf{x})) n_s(\mathbf{y}) \Sigma_{is}^k(\mathbf{x}, \mathbf{y}) - t_i(\mathbf{y}) U_i^k(\mathbf{x}, \mathbf{y}) \right) \right\} dS_y = 0 \quad (18) $$

Assume $u_i \in C^{0,\alpha}$ at $\mathbf{x}$, where $C^{m,\alpha}$ denotes the set of functions $m$ times continuously differentiable such that $\exists (\alpha, C) > 0, |u_{i,m}(\mathbf{x}) - u_{i,m}(\mathbf{y})| \le \|\mathbf{x} - \mathbf{y}\|^\alpha$. Under this assumption, one has:

$$ (u_i(\mathbf{y}) - u_i(\mathbf{x})) \Sigma_{is}^k(\mathbf{x}, \mathbf{y}) \sim \|\mathbf{x} - \mathbf{y}\|^{\alpha-2} \quad (19) $$

Now the limiting process $\epsilon \to 0$ in (18) is investigated. The integrals on $s_\epsilon$ and on $\Gamma - e_\epsilon$ are considered separately.

* Because of (19) and since $\|\mathbf{x} - \mathbf{y}\| \sim \epsilon$, $dS_y \sim \epsilon^2 d\Omega$ ($\Omega$ being the solid angle from $\mathbf{x}$), the integrals over $s_\epsilon$ vanish in the limit $\epsilon \to 0$.

* On $\Gamma - e_\epsilon$, the integrands are of order $\|\mathbf{x} - \mathbf{y}\|^{\alpha-2}$ and $\|\mathbf{x} - \mathbf{y}\|^{-1}$. The limit for $\epsilon \to 0$ of both integrals over $\Gamma - e_\epsilon$ are the corresponding ordinary improper integrals over $\Gamma$.

Hence, taking the limit $\epsilon \to 0$ in (18) gives the regularized DBIE as follows:

$$ \kappa u_k(\mathbf{x}) + \int_{\Gamma} \left\{ (u_i(\mathbf{y}) - u_i(\mathbf{x})) \Sigma_{is}^k(\mathbf{x}, \mathbf{y}) - t_i(\mathbf{y}) U_i^k(\mathbf{x}, \mathbf{y}) \right\} dS_y = 0 \quad (20) $$

where $\kappa$ assumes the same meaning than in (17).

# 4 Tangential differential operators and integration by parts.

## 4.1 Tangential differential operators.

Let $S$ be a twice continuously differentiable ($C^2$) surface, closed or open, of unit normal $n$. Let $\nu$ denote the unit outward normal to $\partial S$ lying in the tangent plane to $S$ and $\tau$ the unit tangent to $\partial S$ defined as $\tau = n \wedge \nu$. Consider a scalar field $f(\mathbf{y})$, $\mathbf{y} \in S$. The function $f$ may be undefined outside $S$ (e.g. $S$ is a crack and $f$ is the crack opening displacement, or $f = n_i(\mathbf{y})$). Therefore the cartesian derivatives $f_i$ are generally meaningless. The domain of definition of $f$ is extended in a neighbourhood $V$ of $S$ by introducing a continuation $\hat{f}$ of $f$ outside $S$ defined as:

$$ (\mathbf{y} \in V) \qquad \hat{f}(\mathbf{y}) = f(P(\mathbf{y})) \quad (21) $$

where $P(\mathbf{y})$ is the orthogonal projection of $\mathbf{y}$ onto $S$. Clearly the restriction of $\hat{f}$ to $S$ is equal to $f$. Moreover the normal derivative of $\hat{f}$ is equal to zero, i.e. the vector $\text{grad} \hat{f}$ is tangent to $S$; therefore it may be used to define the tangential gradient $\text{grad}_S f$ of the function $f$;

$$ \text{grad}_S f = \text{grad} \hat{f} = \text{grad} \hat{f} \quad (22) $$

If $f$ is an arbitrary scalar function defined in $V$, one has, consistently with (22):

$$ \text{grad}_S f = \text{grad} f - n f_{,n} = e_r \hat{D}_r f = e_r (f_{,r} - n_r f_{,n}) \quad (23) $$

which defines the tangential partial derivatives $\hat{D}_r f$ (using the notation $(\cdot)_{,n} = \partial/\partial n(\cdot)$). In the following, the symbol $(\cdot)$ will be omitted, keeping in mind when necessary the extension (21).
## 4.2 Integration by parts: variants of the Stokes' identity.

The classical Stokes' identity for a vector field U defined over V reads:

$$ \int_S n \cdot \text{rot}(U)dS = \int_{\partial S} U \cdot \tau ds \quad (24) $$

Let $f$ and $v$ be respectively a scalar and a vector field on $S$. Application of Stokes' identity (24) to the vector field $U = (n \wedge e_j)f$ yields the following identity:

$$ \int_S (-n_r K f + D_r f) dS = \int_{\partial S} f v_r ds \quad (25) $$

where $K(y) = D_r n_r(y)$ is twice the mean curvature of $S$. The operator

$$ D_{rs}f = (n_r f_s - n_s f_r) \quad (26) $$

is also introduced. From (23), $D_{rs}f = n_r D_s f - n_s D_r f$: $D_{rs}f$ is a tangential differential operator. The particular choice $U = n \wedge (e_j \wedge n)f$ in (24) leads to another identity relating surface and contour integrals:

$$ \int_S D_{rs}f dS = \epsilon_{jrs} \int_{\partial S} f \tau_j ds \quad \text{for any fixed pair } r, s, \quad r, s = 1, 2, 3 \quad (27) $$

where $\epsilon_{jrs}$ denotes the permutation symbol of the indices $j, r, s$. Identity (27) is very interesting for BEM formulations. It allows integration by parts on a closed or open surface using ordinary partial derivatives (i.e. without separation of tangential and normal derivatives), thanks to eqn. (26). It will be used in section 6.

The contour integrals in the right-hand sides of (25), (27) vanish if $S$ is closed and piecewise smooth (i.e. made of several smooth open surfaces), provided $f(y)$ is continuous at the edges.

## 4.3 Integral identities involving Kelvin tensors.

The regularization approach developed in sections 6, 7 will involve the integrals:

$$ A_{is}^{k}(x, S) = \int_{S} \Sigma_{is}^{k}(x, y) dS_{y} \qquad B_{ir}^{k}(x, S) = \int_{S} U_{i,r}^{k}(x, y) dS_{y} \quad (28) $$

where $x$ is a fixed point not on the surface $S$, which ensures the validity of any analytical treatment (such as integration by

parts, for the present matters) performed on these integrals.

Examination of formulas (93) and (94) reveal that the two first integrals $A_{is}^k(x, S)$, $B_{ir}^k(x, S)$ are expressed in terms of two basic integrals $I_a(x, S)$, $J_{abc}(x)$ ($a, b, c = 1, 2, 3$), as follows:

$$ A_{is}^{k}(x, S) = -\frac{1}{4\pi} \left[ \beta^{2}(\delta_{ik}I_{s}(x, S) + \delta_{sk}I_{i}(x, S) - \delta_{is}I_{k}(x, S)) + 3(1 - \beta^{2})J_{isk}(x, S) \right] \quad (29) $$

$$ B_{ir}^{k}(x, S) = -\frac{1}{8\pi\mu} \left[ (1 + \beta^2)\delta_{ik}I_{r}(x, S) + (1 - \beta^2)(3J_{ikr}(x, S) - \delta_{rk}I_{i}(x, S) - \delta_{ir}I_{k}(x, S)) \right] \quad (29) $$

$$ I_a(x, S) = \int_S r, a \frac{dS_y}{r^2}; \qquad J_{abc}(x, S) = \int_S r, a, r, b, r, c \frac{dS_y}{r^2}; \qquad \beta^2 = \frac{\mu}{\lambda + 2\mu} \quad (30) $$

The integrals $I_a(x, S)$, $J_{abc}(x)$ can be transformed, using identities (25), (27), into the following expressions (their derivation is given in Appendix B):

$$ I_a(x, S) = \int_S \left(\frac{1}{r}, n - K(y)\right) n_a(y) \frac{dS_y}{r} - \int_{\partial S} \nu_a(y) \frac{dS_y}{r} \quad (31) $$

$$ \begin{aligned} 3J_{abc}(x, S) = & \delta_{ac}I_b(x, S) + \delta_{bc}I_a(x, S) + \int_S (r, b, r, p, D_{pa}n_c - D_a(n_b n_c) - D_b(n_a n_c)) \frac{dS_y}{r} \\ & + \int_S (2n_a n_b - r, a, r, b)n_c K \frac{dS_y}{r} + \delta_{ab} \int_S n_c r, n \frac{dS_y}{r^2} - 2 \int_S n_a n_c n_b r, n \frac{dS_y}{r^2} \\ & + \int_{\partial S} (\nu_a n_c n_b + \nu_b n_c n_a - r, a, r, b) \nu_c \frac{dS_y}{r} + \int_{\partial S} (\nu_p n_a - \nu_a n_p) r, p, r, b n_c \frac{dS_y}{r} \end{aligned} \quad (32) $$
with $r, n = n_p r, p, K = D_r n_r$ and in which the dependencies $n = n(y), \nu = \nu(y), K = K(y), r, i = \partial/\partial y_i r(x, y)$ are made implicit for brevity.

The reformulations (31), (32) of $I_a(x, S)$, $J_{abc}(x)$ involve surface integrals with integrands being of order $O(1/r)$ ($S$ being $C^2$, $1/r^2 r, n$ is of order $O(1/r)$) for $x$ close to the integration point $y$. Besides, in view of the symmetry of $J_{abc}(x, S)$ with respect to the indices $a, b, c$ in (30), formula (32) is invariant under any permutation of $a, b, c$ in either side. However, this symmetry does not hold for the surface and contour integrals taken separately.

The integral

$$K_i^k(x, S) = \int_S n_s(y) \Sigma_{is}^k(x, y) dS_y \quad (33)$$

will also be used. Although it is not directly expressible using $I_a(x, S)$, $J_{abc}(x, S)$ integrals, it can be transformed along similar lines (see Appendix B) into:

$$K_i^k(x, S) = -\frac{1}{4\pi} \left[ \delta_{ik} \int_S r, n \frac{dS_y}{r^2} + 3(1 - \beta^2) \int_S r, i r, s r, n \frac{dS_y}{r} + \beta^2 \epsilon_{ikp} \int_{\partial S} \tau_p \frac{dS_y}{r} \right] \quad (34)$$

# 5 Second-order regularized elastostatic GBIE and TBIE.

## 5.1 Definition of gradient BIE via a limiting process.

The notations are those introduced in subsection 2.3. The elastostatic gradient boundary integral equation (GBIE) results from application of identity (11) to the domain $\Omega_c$ and then taking the limit $\epsilon \to 0$ (which existence itself will result from the regularization approach). Repeating the considerations of subsection 2.3, the static GBIE is stated as the limit for $\epsilon \to 0$ of

$$\int_{(\Gamma - \epsilon_c) + s_c} \{ u_i(y) n_s(y) \Sigma_{is,r}^k(x, y) - t_i(y) U_{i,r}^k(x, y) \} dS_y = 0 \quad (35)$$

Application of Hooke's tensor (4) to (35) yields the traction boundary integral equation (TBIE).

Two approaches are considered in this section: the second-order regularization and an integration by parts followed by a first-order regularization.

## 5.2 Second-order regularization.

This approach, like the first-order regularization of DBIE, uses a simple elastostatic solution as an auxiliary problem. Owing to the $||y - x||^{-3}$ singularity of the derivatives $\Sigma_{is,r}^k[x, t, y|u_i(y, t)]$ of the second Stokes' tensor, it is then necessary to perform a second-order regularisation, i.e. subtract $u(x)$ and $u_r(x)$ from $u(y)$.

Consider, as an auxiliary solution of elastostatic equilibrium equation, the superposition of a rigid displacement and a uniform strain displacement of $D_c$, defined by:

$$v_i(y) = u_i^0 + A_{ij} r_j t_i(y) = C_{isab} n_s(y) A_{ab} \quad \text{on } \partial D_c \quad (\text{with } r_j = y_j - x_j) \quad (36)$$

$(u_i^0, A_{ij}$ constants). Equation (11) holds true for this elastic state. It reads:

$$\kappa A_{kr} + \int_{(\Gamma - \epsilon_c) + s_c} \{[u_i^0 + A_{ij} r_j] n_s(y) \Sigma_{is,r}^k(x, y) - C_{isab} n_s(y) A_{ab} U_{i,r}^k(x, y)\} dS_y = 0 \quad (37)$$

where $\kappa = 0$ ($\Omega_c$ bounded) or $\kappa = 1$ ($\Omega_c$ unbounded), and using as in (17) the unit normal exterior to $\Omega_c$. Identity (37) remains in particular true for the choice $u_i^0 = u_i(x)$ and $A_{ij} = u_{i,j}(x)$, ie the unknown displacement and displacement gradient at the fixed point $x$. With this choice, subtraction of (37) from (35) gives:

$$\begin{aligned} & -\kappa u_{k,r}(x) + \int_{(\Gamma - \epsilon_c) + s_c} [u_i(y) - u_i(x) - u_{i,j}(x)r_j] n_s(y) \Sigma_{is,r}^k(x, y) dS_y \\ & - \int_{(\Gamma - \epsilon_c) + s_c} [t_i(y) - C_{isab} n_s(y) u_{a,b}(x)] U_{i,r}^k(x, y) dS_y = 0 \end{aligned} \quad (38)$$
Assume $u_i \in C^{1,\alpha}$, ($i = 1,2,3$) at $x$. Then one has $|u_i(y) - u_i(x) - u_{i,j}(x)r_j| \le C \|x - y\|^{1+\alpha}$, and thus:

$$ 
\begin{aligned}
& [u_i(y) - u_i(x) - u_{i,j}(x)r_j] \Sigma_{is,r}^k(x, y) \sim \|x - y\|^{\alpha-2} \\
& C_{isab} n_s(y) [u_{a,b}(y - u_{a,b}(x))] U_{i,r}^k(x, y) \sim \|x - y\|^{\alpha-2}
\end{aligned} 
\quad (39) 
$$

Now the integrals on $s_\epsilon$ and on $\Gamma - e_\epsilon$ in (38) are considered separately.

* Because of (39) and since $\|x - y\| \sim \epsilon$, $dS_y \sim \epsilon^2 d\Omega$ ($\Omega$ being the solid angle from $x$), the integrals over $s_\epsilon$ vanish for $\epsilon \to 0$.

* On $\Gamma - e_\epsilon$, both integrands are of order $\|x - y\|^{\alpha-2}$ and $\|x - y\|^{-1}$, for $y$ close to $x$. The limit of both integrals over $\Gamma - e_\epsilon$ are then ordinary improper integrals over $\Gamma$.

Hence, taking the limit $\epsilon \to 0$ in (38) yields the following regularized GBIE:

$$ 
\begin{aligned}
& -\kappa u_{k,r}(x) + \int_{\Gamma} [u_i(y) - u_i(x) - u_{i,j}(x)r_j] n_s(y) \Sigma_{is,r}^k(x, y) dS_y \\
& \qquad - \int_{\Gamma} [t_i(y) - C_{isab} n_s(y) u_{a,b}(x)] U_{i,r}^k(x, y) dS_y = 0
\end{aligned} 
\quad (40) 
$$

where $\kappa$ assumes the same meaning than in (17). Since the unit normal at $(x)$ does not appear in equation (40), this result is valid for $(x)$ being a smooth point of $\Gamma$ as well as an edge or corner point, provided the regularity requirement $u_i \in C^{1,\alpha}$ at $x$ is met.

## 5.3 Reformulation of (40) for BEM discretization.

In a boundary element point of view, the regularization is achieved if the property (39) is made explicit when considering the discretization of the surface fields $u$ and $t$ using shape functions, i.e. if the appropriate cancellations occur in the discretized

$$ u_i(y) - u_i(x) - u_{i,j}(x)r_j, \quad t_i(y) - C_{isab} n_s(y) u_{a,b}(x) $$

(see section 9). As these quantities are only defined on $\Gamma$ in a BEM context, one has to rearrange them using:

$$ 
\begin{aligned}
u_{i,j}(x) &= D_r u_i(x) + n_j(x) u_{i,n}(x) \\
C_{isab} n_s(y) u_{a,b}(x) &= t_i(x) + C_{isab} u_{a,b}(x) [n_s(y) - n_s(x)]
\end{aligned} 
\quad (41) 
$$

Provided $u_i \in C^{1,\alpha}$ at $x$ and $\Gamma$ is a (piecewise) Lyapunov surface of exponent $\alpha$, one has:

$$ 
\begin{aligned}
u_i(y) - u_i(x) - r_j D_j u_i(x) &= n_j(x) u_{i,n}(x) r_j = O(\|x - y\|^{1+\alpha}) \\
t_i(y) - t_i(x) &= C_{isab} u_{a,b}(x) [n_s(y) - n_s(x)] = O(\|x - y\|^\alpha)
\end{aligned} 
\quad (42) 
$$

The regularized GBIE (40) can be accordingly recast in an alternative form, which also involves only weakly singular integrals thanks to (42):

$$ 
\begin{aligned}
& \kappa u_{k,r}(x) + \int_{\Gamma} \{[u_i(y) - u_i(x) - r_j D_j u_i(x)] n_s(y) \Sigma_{is,r}^k(x, y) - [t_i(y) - t_i(x)] U_{i,r}^k(x, y)\} dS_y \\
& \qquad + \int_{\Gamma} \{u_{i,n}(x) n_j(x) r_j n_s(y) \Sigma_{is,r}^k(x, y) - C_{isab} u_{a,b}(x) [n_s(y) - n_s(x)] U_{i,r}^k(x, y)\} dS_y = 0
\end{aligned} 
\quad (43) 
$$

The version (43) of the regularized GBIE is better suited to BEM interpolation than (40).
# 6 Integration by parts and first-order regularization of static GBIE.

In this section, an alternative regularization approach for the GBIE is developed in which, instead of a second-order regularization, a first-order regularization follows a preliminary integration by parts.

## 6.1 Integration by parts of the hypersingular static kernel.

The use of identity (27) above, together with equation (8), leads to the following result, which holds for any closed surface S and for x not located on S:

$$ 
\begin{aligned}
\int_S u_i(y)n_s(y)\Sigma_{i\epsilon,r}^k(x,y)dS_y &= \int_S u_i(y)\{D_{sr}\Sigma_{is}^k(x,y) + n_r(y)\Sigma_{is,s}^k(x,y)\}dS_y \\
&= \int_S D_{rs}u_i(y)\Sigma_{is}^k(x,y)dS_y
\end{aligned} 
\quad (44) $$

This integration by parts pattern is very useful in the elasticity BIE methods context. It appears in [44] and [32] and is used in [7] for other purposes.

Application of (44) to the first integral in static GBIE (35) (x being located outside $\Omega_\epsilon$) leads to:

$$ \int_{(\Gamma - e_\epsilon) + s_\epsilon} \{D_{rs}u_i(y)\Sigma_{is}^k(x,y) - t_i(y)U_{i,r}^k(x,y)\} dS_y = 0 \quad (45) $$

The same integration by parts is also performed on identity (37) obtained from the auxiliary solution (36), giving:

$$ 
\begin{aligned}
& \kappa u_{k,r}(\mathbf{x}) + \int_{(\Gamma - e_{\epsilon}) + s_{\epsilon}} (n_{r}(\mathbf{y})u_{i,s}(\mathbf{x}) - n_{s}(\mathbf{y})u_{i,r}(\mathbf{x}))\Sigma_{is}^{k}(\mathbf{x}, \mathbf{y})dS_y \\
& \quad - \int_{(\Gamma - e_{\epsilon}) + s_{\epsilon}} C_{isab}n_{s}(\mathbf{y})u_{a,b}(\mathbf{x})U_{i,r}^{k}(\mathbf{x}, \mathbf{y})dS_y = 0
\end{aligned} 
\quad (46) $$

## 6.2 First-order regularized static TBIE

Eqn. (46) is subtracted from (45), to get:

$$ 
\begin{aligned}
& -\kappa u_{k,r}(\mathbf{x}) + \int_{(\Gamma - e_{\epsilon}) + s_{\epsilon}} [n_{r}(\mathbf{y})(u_{i,s}(\mathbf{y}) - u_{i,s}(\mathbf{x})) - n_{s}(\mathbf{y})(u_{i,r}(\mathbf{y}) - u_{i,r}(\mathbf{x}))]\Sigma_{is}^{k}(\mathbf{x}, \mathbf{y})dS_y \\
& \quad - \int_{(\Gamma - e_{\epsilon}) + s_{\epsilon}} [t_{i}(\mathbf{y}) - C_{isab}n_{s}(\mathbf{y})u_{a,b}(\mathbf{x})]U_{i,r}^{k}(\mathbf{x}, \mathbf{y})dS_y = 0
\end{aligned} 
\quad (47) $$

Following similar lines than in subsection 5.2, identity (47) above leads to a regularized GBIE. However, since quantities $u_{i,s}(\mathbf{y}) - u_{i,s}(\mathbf{x})$ and $C_{isab}n_s(\mathbf{y})u_{a,b}(\mathbf{x})$ involve complicated combinations of tractions and (tangential gradients of) surface displacements, the result is not very convenient in a BEM point of view. The goal of the present subsection is to rearrange equation (47) to get a regularized GBIE explicitely expressed in terms of tractions and tangential gradients of surface displacements.

The substitution $n_i(\mathbf{y}) = (n_i(\mathbf{y}) - n_i(\mathbf{x})) + n_i(\mathbf{x})$ is made in the integrals over $S - e_\epsilon$ in (47). After a suitable rearrangement of the terms, one gets:

$$ 
\begin{aligned}
& -\kappa u_{k,r}(\mathbf{x}) + \int_{\Gamma - e_{\epsilon}} \{[D_{rs}u_i(\mathbf{y}) - D_{rs}u_i(\mathbf{x})]\Sigma_{is}^k(\mathbf{x}, \mathbf{y}) - [t_i(\mathbf{y}) - t_i(\mathbf{x})]U_{i,r}^k(\mathbf{x}, \mathbf{y})\} dS_y \\
& + \int_{s_{\epsilon}} [n_r(\mathbf{y})(u_{i,s}(\mathbf{y}) - u_{i,s}(\mathbf{x})) - n_s(\mathbf{y})(u_{i,r}(\mathbf{y}) - u_{i,r}(\mathbf{x}))]\Sigma_{is}^k(\mathbf{x}, \mathbf{y})dS_y \\
& - \int_{s_{\epsilon}} (t_i(\mathbf{y}) - C_{isab}n_s(\mathbf{y})u_{a,b}(\mathbf{x}))U_{i,r}^k(\mathbf{x}, \mathbf{y})dS_y + I_{kr}(\mathbf{x}, \epsilon) = 0
\end{aligned} 
\quad (48) $$
$$ \mathcal{I}_{kr}(\mathbf{x}, \epsilon) = - \int_{\Gamma - e_{\epsilon}} [(n_r(\mathbf{y}) - n_r(\mathbf{x}))u_{i,s}(\mathbf{x}) - (n_s(\mathbf{y}) - n_s(\mathbf{x}))u_{i,r}(\mathbf{x})] \Sigma_{is}^k(\mathbf{x}, \mathbf{y}) dS_y \quad (49) \\ + \int_{\Gamma - e_{\epsilon}} C_{isab}(n_s(\mathbf{y}) - n_s(\mathbf{x})) u_{a,b}(\mathbf{x}) U_{i,r}^k(\mathbf{x}, \mathbf{y}) dS_y $$

Because $u_i \in C^{1,\alpha}$ and property (39),

* The two first integrals in (48), in the limit $\epsilon \to 0$, are ordinary improper integrals over $G$.

* the integrals over $s_\epsilon$ vanish for $\epsilon \to 0$.

The limiting process thus reduces to the investigation of $\mathcal{I}_{kr}(\mathbf{x}, \epsilon)$. Inserting Hooke's law, noticing the symmetry $C_{isab} = C_{abis}$ of Hooke's tensor and using the notations of section 4 in (49) leads to:

$$ \mathcal{I}_{kr}(\mathbf{x}, \epsilon) = u_{i,r}(\mathbf{x}) K_i^k(\mathbf{x}, \Gamma - e_\epsilon) + D_{rs} u_i(\mathbf{x}) A_{is}^k(\mathbf{x}, \Gamma - e_\epsilon) - n_s(\mathbf{x}) C_{isab} u_{a,b}(\mathbf{x}) B_{ir}^k(\mathbf{x}, \Gamma - e_\epsilon) \\ - C_{isab} u_{i,s}(\mathbf{x}) \int_{\Gamma - e_\epsilon} [n_r(\mathbf{y}) U_{a,b}^k(\mathbf{x}, \mathbf{y}) - n_b(\mathbf{y}) U_{a,r}^k(\mathbf{x}, \mathbf{y})] dS_y \quad (50) $$

In (50), the first integral is integrated by parts using identity (27):

$$ \int_{\Gamma - e_{\epsilon}} [n_r(\mathbf{y}) U_{a,b}^{k}(\mathbf{x}, \mathbf{y}) - n_b(\mathbf{y}) U_{a,r}^{k}(\mathbf{x}, \mathbf{y})] dS_y = \int_{\Gamma - e_{\epsilon}} D_{rs} U_{a}^{k}(\mathbf{x}, \mathbf{y}) dS_y \\ = \int_{c_{\epsilon}} U_{a,r}^{k}(\mathbf{x}, \mathbf{y}) \epsilon_{rbi} \tau_i(\mathbf{y}) dS_y \quad (51) $$

Using (51) and the results of section 4, the integrals in (50) are split into surface integrals and contour integrals over $c_\epsilon$.

* In the limit $\epsilon \to 0$, all surface integrals in identities (31), (32) are ordinary improper integrals over $\Gamma$. Moreover, since $\mathbf{x}$ is a smooth point of $G$, it can be shown that, from (34):

$$ K_i^k(\mathbf{x}, \Gamma) = -\left(\frac{1}{2} - \kappa\right)\delta_{ik} \quad \text{with} \quad \begin{cases} \kappa = 1 & \Omega \text{ unbounded} \\ \kappa = 0 & \Omega \text{ bounded} \end{cases} \quad (52) $$

Hence the contribution of the surface integrals in (50) equals

$$ -(\frac{1}{2} - \kappa)u_{i,k}(\mathbf{x}) + D_{rs}u_i(\mathbf{x})A_{is}^k(\mathbf{x}, \Gamma) - t_i(\mathbf{x})B_{ir}^k(\mathbf{x}, \Gamma) $$

* Using (51), the total contour integral $\partial \mathcal{I}_{kr}(\mathbf{x}, \epsilon)$ over $c_\epsilon$ arising in (50) equals:

$$ \begin{aligned} \partial \mathcal{I}_{kr}(\mathbf{x}, \epsilon) &= -C_{isab} u_{i,s}(\mathbf{x}) \int_{c_\epsilon} U_{a,r}^k(\mathbf{x}, \mathbf{y}) \epsilon_{rbi} \tau_i(\mathbf{y}) dS_y + u_{i,r}(\mathbf{x}) \partial K_i^k(\mathbf{x}, \Gamma - e_\epsilon) \\ &\quad + D_{rs} u_i(\mathbf{x}) \partial A_{is}^k(\mathbf{x}, \Gamma - e_\epsilon) - n_s(\mathbf{x}) C_{isab} u_{a,b}(\mathbf{x}) \partial B_{ir}^k(\mathbf{x}, \Gamma - e_\epsilon) \\ &= -C_{isab} u_{i,s}(\mathbf{x}) \int_{c_\epsilon} U_{a,r}^k(\mathbf{x}, \mathbf{y}) \epsilon_{rbi} \tau_i(\mathbf{y}) dS_y \\ &\quad + u_{i,r}(\mathbf{x}) \partial K_i^k(\mathbf{x}, \Gamma - e_\epsilon) - u_{i,r}(\mathbf{x}) n_s(\mathbf{x}) \partial A_{is}^k(\mathbf{x}, \Gamma - e_\epsilon) \\ &\quad + C_{isab} u_{i,s}(\mathbf{x}) [n_r(\mathbf{x}) \partial B_{ab}^k(\mathbf{x}, \Gamma - e_\epsilon) - n_b(\mathbf{x}) \partial B_{ar}^k(\mathbf{x}, \Gamma - e_\epsilon)] \end{aligned} \quad (53) $$

where $\partial A_{is}^k(\mathbf{x}, S)$, $\partial B_{is}^k(\mathbf{x}, S)$, $\partial K_i^k(\mathbf{x}, S)$ symbolically collect all contour integrals over $c_\epsilon$ arising in $A_{is}^k(\mathbf{x}, S)$, $B_{is}^k(\mathbf{x}, S)$, $K_i^k(\mathbf{x}, S)$ given by identities (31), (32) and (34).
Using (29), (31) and (32), one has:

$$
\begin{equation}
\begin{aligned}
n_r(\mathbf{x}) \partial B_{ab}^k(\mathbf{x}, \Gamma - \epsilon_c) - n_b(\mathbf{x}) \partial B_{ar}^k(\mathbf{x}, \Gamma - \epsilon_c) &= \\
&= -\frac{1}{8\pi\mu} \left\{ -(1 + \beta^2)\delta_{ak} \left[ n_r(\mathbf{x}) \int_{c_\epsilon} \nu_b(\mathbf{y}) \frac{d\mathbf{sy}}{r} - n_b(\mathbf{x}) \int_{c_\epsilon} \nu_r(\mathbf{y}) \frac{d\mathbf{sy}}{r} \right] \right. \\
&\qquad \left. + (1 - \beta^2) \left[ n_r(\mathbf{x}) \int_{c_\epsilon} (n_b[\nu_a n_k + \nu_k n_a + (\nu_p n_a - \nu_a n_p)r_{,p}r_{,k}] - r_{,a}r_{,k}\nu_b) \frac{d\mathbf{sy}}{r} \right. \right. \\
&\qquad \left. - n_b(\mathbf{x}) \int_{c_\epsilon} (n_r[\nu_a n_k + \nu_k n_a + (\nu_p n_a - \nu_a n_p)r_{,p}r_{,k}] - r_{,a}r_{,k}\nu_r) \frac{d\mathbf{sy}}{r} \right] \right\}
\end{aligned}
\tag{54}
\end{equation}
$$

The point x being smooth on Γ, the curve c_ε becomes a plane curve for vanishing ε. Thus,
νᵢ(y) = O(1), nᵢ(y) = O(1), νᵢ(y) - νᵢ(x) = O(ε), nᵢ(y) - nᵢ(x) = O(ε) and dsy/r = O(1)dθ
(θ denoting the polar angle of origin x in the tangent plane to Γ at x). Thus, each integrand
in (54) behaves like O(1)dθ for small ε. This, together with nᵣ(x)n₉(y) - n₉(x)nᵣ(y) = O(ε),
gives:

$$
\int_{c_{\epsilon}} [n_{r}(\mathbf{x})n_{b} - n_{b}(\mathbf{x})n_{r}][\nu_{a}n_{k} + \nu_{k}n_{a} + (\nu_{p}n_{a} - \nu_{a}n_{p})r_{,p}r_{,k}] \frac{dsy}{r} = O(\epsilon) \quad (55)
$$

Hence:

$$
\begin{align*}
n_r(\mathbf{x}) \partial B_{ab}^k(\mathbf{x}, \Gamma - \epsilon_c) - n_b(\mathbf{x}) \partial B_{ar}^k(\mathbf{x}, \Gamma - \epsilon_c) & \\
&= \frac{1}{8\pi\mu} \int_{c_\epsilon} \{(1 + \beta^2)\delta_{ak} + (1 - \beta^2)r_{,a}r_{,k}\} (n_r(\mathbf{x})\nu_b - n_b(\mathbf{x})\nu_r) \frac{d\mathbf{sy}}{r} + O(\epsilon) \\
&= \frac{\epsilon_{rbi}}{8\pi\mu} \int_{c_\epsilon} \{(1 + \beta^2)\delta_{ak} + (1 - \beta^2)r_{,a}r_{,k}\} \tau_i \frac{d\mathbf{sy}}{r} + O(\epsilon) \\
&= \epsilon_{rbi} \int_{c_\epsilon} \tau_i U_a^k(\mathbf{x}, \mathbf{y}) \frac{d\mathbf{sy}}{r} + O(\epsilon) \tag{56}
\end{align*}
$$

using $n_r(\mathbf{x})\nu_b(\mathbf{y}) - n_b(\mathbf{x})\nu_r(\mathbf{y}) = n_r(\mathbf{y})\nu_b(\mathbf{y}) - n_b(\mathbf{y})\nu_r(\mathbf{y}) + O(\epsilon) = \epsilon_{rbi}\tau_i(\mathbf{y}) + O(\epsilon)$. On the
other hand, from (29), (31), (32):

$$
\begin{align}
n_s(\mathbf{x}) \partial A_{is}^k(\mathbf{x}, \Gamma - \epsilon_c) &= -\frac{1}{4\pi} \left[ \beta^2 \int_{c_\epsilon} [n_i(\mathbf{x})\nu_k - n_k(\mathbf{x})\nu_i - \delta_{ik}n_s(\mathbf{x})\nu_s] \frac{d\mathbf{sy}}{r} \right. \nonumber \\
&\qquad \left. + (1-\beta^2) \int_{c_\epsilon} [n_k(\mathbf{x})\nu_i(n_p(\mathbf{x})n_p - 1) + n_i n_k(n_p(\mathbf{x})\nu_p)] \frac{d\mathbf{sy}}{r} \right. \nonumber \\
&\qquad \left. + (1-\beta^2) \int_{c_\epsilon} [n_q(\mathbf{x})r_{,q}((n_i\nu_p - n_p\nu_i)n_k r_{,p} - \nu_k r_{,i})] \frac{d\mathbf{sy}}{r} \right] \nonumber \\
&= -\frac{1}{4\pi} \beta^2 \epsilon_{ikq} \int_{c_\epsilon} \tau_q \frac{d\mathbf{sy}}{r} + O(\epsilon) \tag{57}
\end{align}
$$

$$
\text{using } n_s(\mathbf{x})\nu_s(\mathbf{y}) = O(\epsilon), n_p(\mathbf{x})n_p(\mathbf{y}) - 1 = O(\epsilon), n_q(\mathbf{x})r_{,q}(\mathbf{x}, \mathbf{y}) = O(\epsilon).
$$

Collecting results (34), (56) and (57), one has:

$$
\begin{align}
\partial \mathcal{I}_{kr}(\mathbf{x}, \epsilon) &= -C_{isab} u_{i,s}(\mathbf{x}) \int_{c_\epsilon} U_{a,r}^k(\mathbf{x}, \mathbf{y}) \epsilon_{rbi} \tau_i(\mathbf{y}) d\mathbf{sy} + C_{isab} u_{i,s}(\mathbf{x}) \int_{c_\epsilon} U_{a,r}^k(\mathbf{x}, \mathbf{y}) \epsilon_{rbi} \tau_i(\mathbf{y}) d\mathbf{sy} \nonumber \\
&\quad + u_{i,r}(\mathbf{x}) \frac{\beta^2}{4\pi} \epsilon_{ikq} \int_{c_\epsilon} \tau_q \frac{d\mathbf{sy}}{r} - u_{i,r}(\mathbf{x}) \frac{\beta^2}{4\pi} \epsilon_{ikq} \int_{c_\epsilon} \tau_q \frac{d\mathbf{sy}}{r} + O(\epsilon) \nonumber \\
&= O(\epsilon) \tag{58}
\end{align}
$$

that is, the total contour integral $\partial \mathcal{I}_{kr}(\mathbf{x}, \epsilon)$ over $c_\epsilon$ arising in (50) vanishes for $\epsilon \to 0$. Equation (58) is the key step of the current regularization of the static GBIE: thanks to it, the integral $\mathcal{J}_{kr}(\mathbf{x}, \epsilon)$ reduces to a sum of weakly singular integrals over $\Gamma$, according to eqns. (31), (32):

$$
\lim_{\epsilon \to 0} \mathcal{I}_{kr}(\mathbf{x}, \epsilon) = -\left(\frac{1}{2} - \kappa\right) u_{i,k}(\mathbf{x}) + D_{rs} u_i(\mathbf{x}) A_{is}^k(\mathbf{x}, \Gamma) - t_i(\mathbf{x}) B_{ir}^k(\mathbf{x}, \Gamma) \quad (59)
$$
The result (59) completes the investigation of the limiting process (35) and finally allows the statement of the following regularized static GBIE:

$$ -\frac{1}{2}u_{k,r}(x) + D_{rs}u_i(x)A_{is}^k(x, \Gamma) - t_i(x)B_{ir}^k(x, \Gamma) + \int_{\Gamma-e_\epsilon} [D_{rs}u_i(y) - D_{rs}u_i(x)] \Sigma_{is}^k(x, y)dS_y - \int_{\Gamma-e_\epsilon} [t_i(y) - t_i(x)] U_{i,r}^k(x, y)dS_y = 0 \quad (60) $$

The densities of the integrals in GBIE (60) are the tangential gradients of u and the tractions t. This is a nice feature in a boundary element point of view, because the boundary integrals involve independently interpolated fields defined on the boundary $\Gamma$.

# 7 Regularized elastostatic first-order traction BIE for curved cracks.

## 7.1 Elastostatic traction BIE for cracks in infinite elastic solids.

Consider a perfect crack embedded in an infinite elastic space $\Omega$. The two faces $S^+$, $S^-$ of the crack are geometrically identical and have opposite unit normals $n^+$, $n^-$, such that $n^-$ is oriented from $S^+$ to $S^-$. The surface $S = S^+ = S^-$ is piecewise twice continuously differentiable. The response of the infinite elastic space to static tractions $T_i^\pm$ applied on the two crack faces is considered.

It is well-known that the DBIE applied to this crack problem is ill-posed, because the (given) tractions appear only by their sum $T_i^\pm$. Hence, the traction BIE is needed. Moreover, for $x \in S$, the kernels are singular at $S^+ \ni y^+ = x$ and $S^- \ni y^- = x$. Thus, the limiting process which defines the GBIE or TBIE in this case must be modified accordingly. Let $x$ be a regular interior point of $S$. Consider a neighbourhood $v_\epsilon$ of $x$, vanishing with $\epsilon$ and of arbitrary shape. The closed surface $s_\epsilon = \partial v_\epsilon$ is split into two parts $s_\epsilon^+$ and $s_\epsilon^-$, according to figure 3. Let $e_\epsilon^\pm = S^\pm \cap v_\epsilon$. Equation (45) can then be written for the domain $\Omega_\epsilon$ exterior to the cavity (of zero thickness, except for $v_\epsilon$) bounded by $\partial\Omega_\epsilon = ((S - e_\epsilon)^+ + s_\epsilon^+) + ((S - e_\epsilon)^- + s_\epsilon^-)$, the limit of this cavity for $\epsilon \to 0$ being the crack:

$$ \lim_{\epsilon \to 0} \left\{ \int_{S-e_\epsilon} \left\{ D_{rs}\phi_i(y,t)\Sigma_{is}^k(x,y) - S_i(y,t)U_{i,r}^k(x,y) \right\} dS_y + \int_{s_\epsilon^+ + s_\epsilon^-} \left\{ D_{rs}u_i(y,t)\Sigma_{is}^k(x,y) - t_i(y,t)U_{i,r}^k(x,y) \right\} dS_y \right\} = 0 \quad (61) $$

using $n = n^-$ and introducing the crack opening displacement (COD) $\phi_i(y,t) = [u_i(y,t)]$ and the sum of tractions $S_i(y) = T_i^+(y^+, t) + T_i^-(y^-, t) = C_{isab}n_s(y)[u_{a,b}(y,t)]$, where $[f(y)] = f(y^+) - f(y^-)$ is the jump operator.

Identity (27) have been used to integrate by parts the terms containing derivatives of the static stress kernels. Although the surface $S$ is open, the contour integrals arising from (27) vanish because $\phi_i(y,t) = 0$ on $\partial S$.

## 7.2 Regularization of the static traction BIE for cracks.

As the idealized crack has no interior, identities (37) or (46) are not applicable. The singularities must be isolated and then evaluated separately. The first step is to put $n_i(y) = (n_i(y) - n_i(x)) + n_i(x)$ in the integral over $S - e_\epsilon$ and $u_{i,j}(y) = (u_{i,j}(y) - u_{i,j}(x)) + u_{i,j}(x)$ in the integral over $s_\epsilon = s_\epsilon^+ + s_\epsilon^-$ in (61). Upon this manipulation (for a fixed value of $\epsilon$), (61) becomes:

$$ \int_{S-e_\epsilon} \left\{ [D_{rs}\phi_i(y) - D_{rs}\phi_i(x)] \Sigma_{is}^k(x, y) - [S_i(y) - S_i(x)] U_{i,r}^k(x, y) \right\} dS_y + \int_{s_\epsilon^+ + s_\epsilon^-} \left[ n_r(y)(u_{i,s}(y) - u_{i,s}(x)) - n_s(y)(u_{i,r}(y) - u_{i,r}(x)) \right] \Sigma_{is}^k(x, y)dS_y - \int_{s_\epsilon^+ + s_\epsilon^-} (t_i(y) - C_{isab}n_s(y)u_{a,b}(x)) U_{i,r}^k(x, y)dS_y + J_{kr}(x, \epsilon) = 0 \quad (62) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Geometrical notations and conventions for the crack. Exclusion neighbourhood $v_\epsilon$ and related notations used for the limiting process (61).</figcaption></figure>

$$ \mathcal{J}_{kr}(\mathbf{x}, \epsilon) = \int_{s_{\epsilon}^{+} + s_{\epsilon}^{-}} \left\{ \begin{aligned} & [n_r(\mathbf{y})u_{i,s}(\mathbf{x}) - n_s(\mathbf{y})u_{i,r}(\mathbf{x})] \sum_{is}^{k}(\mathbf{x}, \mathbf{y}) dS_y - C_{isab}n_s(\mathbf{y})u_{a,b}(\mathbf{x})U_{i,r}^{k}(\mathbf{x}, \mathbf{y}) \\ & + D_{rs}\phi_i(\mathbf{x})A_{is}^{k}(\mathbf{x}, S - e_{\epsilon}) - S_i(\mathbf{x})B_{i,r}^{k}(\mathbf{x}, S - e_{\epsilon}) \end{aligned} \right\} dS_y \quad (63) $$

Because $u_i \in C^{1,\alpha}$ and property (39),

* The integral over $S - e_\epsilon$ in (48), in the limit $\epsilon \to 0$, is the corresponding ordinary improper integrals over $S$.

* the first two integrals over $s_\epsilon^+ + s_\epsilon^-$ vanish for $\epsilon \to 0$.

The limiting process thus reduces to the investigation of the limit of $\mathcal{J}_{kr}(\mathbf{x}, \epsilon)$. The latter can be recast in a form similar to (50)-(51) and using similar considerations:

$$ \mathcal{J}_{kr}(\mathbf{x}, \epsilon) = -C_{isab} \left\{ \begin{aligned} & u_{i,s}^{+}(\mathbf{x}) \int_{c_{\epsilon}} U_{a}^{k}(\mathbf{x}, \mathbf{y}) \epsilon_{rbi} \tau_{i}^{+}(\mathbf{y}) dS_{y} + u_{i,s}^{-}(\mathbf{x}) \int_{c_{\epsilon}} U_{a}^{k}(\mathbf{x}, \mathbf{y}) \epsilon_{rbi} \tau_{i}^{-}(\mathbf{y}) dS_{y} \\ & + (n_{r}(\mathbf{x})[u_{i,s}(\mathbf{x})] - n_{s}(\mathbf{x})[u_{i,r}(\mathbf{x})]) A_{is}^{k}(\mathbf{x}, S - e_{\epsilon}) \\ & - n_{s}(\mathbf{x}) C_{isab}[u_{a,b}(\mathbf{x})] B_{ir}^{k}(\mathbf{x}, S - e_{\epsilon}) + u_{i,r}^{+}(\mathbf{x}) K_{i}^{k}(\mathbf{x}, s_{\epsilon}^{+}) + u_{i,r}^{-}(\mathbf{x}) K_{i}^{k}(\mathbf{x}, s_{\epsilon}^{-}) \end{aligned} \right\} \quad (64) $$

The integrals in (64) are split into surface integrals and contour integrals over $c_\epsilon$.

* In the limit $\epsilon \to 0$, all surface integrals in identities (31), (32), (34) are ordinary improper integrals over $\Gamma$. Moreover, since $\mathbf{x}$ is a smooth point of $G$, it can be shown from (34) that:

$$ K_i^k(\mathbf{x}, s_\epsilon^\pm) = -\frac{1}{2}\delta_{ik} - \frac{1}{4\pi} \int_{c_\epsilon} \tau_p^\pm \frac{ds_y}{r} \quad (65) $$

Hence the contribution of the surface integrals in (64) equals

$$ -\frac{1}{2}(u_{i,k}^{+}(\mathbf{x}) + u_{i,k}^{-}(\mathbf{x})) + D_{rs}\phi_{i}(\mathbf{x})A_{is}^{k}(\mathbf{x}, S) - S_{i}(\mathbf{x})B_{ir}^{k}(\mathbf{x}, S) $$

* From (64), (65) and using $\tau^+ = -\tau^-$, the total contour integral $\partial \mathcal{J}_{kr}(\mathbf{x}, \epsilon)$ over $c_\epsilon$ arising in (64) is:

$$ \begin{aligned} \partial \mathcal{J}_{kr}(\mathbf{x}, \epsilon) = & -C_{isab}[u_{i,s}(\mathbf{x})] \int_{c_{\epsilon}} U_{a,r}^{k}(\mathbf{x}, \mathbf{y}) \epsilon_{rbi} \tau_{i}^{+}(\mathbf{y}) dS_{y} \\ & + [u_{i,r}(\mathbf{x})] \partial K_{i}^{k}(\mathbf{x}, \Gamma - e_{\epsilon}) - [u_{i,r}(\mathbf{x})] n_{s}(\mathbf{x}) \partial A_{is}^{k}(\mathbf{x}, \Gamma - e_{\epsilon}) \\ & + C_{isab}[u_{i,s}(\mathbf{x})] [n_{r}(\mathbf{x}) \partial B_{ab}^{k}(\mathbf{x}, \Gamma - e_{\epsilon}) - n_{b}(\mathbf{x}) \partial B_{ar}^{k}(\mathbf{x}, \Gamma - e_{\epsilon})] \end{aligned} \quad (66) $$
where $\partial A_{is}^k(x, S)$, $\partial B_{is}^k(x, S)$, $\partial K_i^k(x, S)$ symbolically collect all contour integrals over $c_e$ arising in $A_{is}^k(x, S)$, $B_{is}^k(x, S)$, $K_i^k(x, S)$ given by identities (31), (32) and (34). In (66), the substitution $D_{rs}\phi_i(x) = n_r(x)[u_{i,s}(x)] - n_s(x)[u_{i,r}(x)]$ has been made; its validity stems from the fact that $D_{rs}\phi_i(x)$ involves only tangential derivatives, hence the (tangential) differentiation and the jump [-] operators may be exchanged.

The contour integral $\partial \mathcal{J}_{kr}(x, \epsilon)$ (equation (66)) is the same as the contour integral $\partial \mathcal{I}_{kr}(x, \epsilon)$ (equation (53)), except that the jumps of gradients $[u_{i,j}(x)]$ appear instead of the gradients $u_{i,j}(x)$ themselves. Hence, the analysis conducted in subsection 6 is applicable again to $\partial \mathcal{J}_{kr}(x, \epsilon)$ and leads to:

$$ \lim_{\epsilon \to 0} \partial \mathcal{J}_{kr}(x, \epsilon) = 0 \qquad (67) $$

that is, the total contour integral $\partial \mathcal{J}_{kr}(x, \epsilon)$ over $c_e$ arising in (50) vanishes for $\epsilon \to 0$.

Equation (67) is the key step of the current regularization of the static crack TBIE: thanks to it, the integral $\mathcal{J}_{kr}(x, \epsilon)$ reduces to a sum of weakly singular integrals over $S$ and of contour integrals over $\partial S$, according to eqns. (31), (32):

$$ \lim_{\epsilon \to 0} \mathcal{J}_{kr}(x, \epsilon) = -\frac{1}{2}(u_{i,k}^{+}(x) + u_{i,k}^{-}(x)) + D_{rs}\phi_{i}(x)A_{is}^{k}(x, S) - S_{i}(x)B_{ir}^{k}(x, S) \qquad (68) $$

The result (68) completes the investigation of the limiting process (62) and finally allows the statement of the following regularized static GBIE:

$$ \frac{1}{2}(u_{i,k}^{+}(x) + u_{i,k}^{-}(x)) + D_{rs}\phi_{i}(x)A_{is}^{k}(x, S) - S_{i}(x)B_{ir}^{k}(x, S) + \int_{S} [D_{rs}\phi_{i}(y) - D_{rs}\phi_{i}(x)] \Sigma_{is}^{k}(x, y)dS_{y} - \int_{\Gamma - e_{\epsilon}} [S_{i}(y) - S_{i}(y)] U_{i,r}^{k}(x, y)dS_{y} = 0 \qquad (69) $$

Application of the tensor $C_{lpkr} n_l^\pm(x)$ to (69) leads to the TBIEs where the given tractions $T_l^\pm(x)$ are related to unknown COD $\phi(y)$ and the known $S(y)$. In the common case where an incident stress tensor $\sigma^I(x)$ is known and the superposition principle is applied, the tractions $T_l^\pm(y)$ are of opposite sign: $S(y) = 0, y \in S$ and the static regularized TBIE reads, from (69):

$$ T_l^\pm(x) = -C_{lpkr} \left\{ \int_S [D_{rs}\phi_i(y) - D_{rs}\phi_i(x)] \Sigma_{is}^k(x, y)dS_y + D_{rs}\phi_i(x)A_{is}^k(x, S) \right\} \qquad (70) $$

## 8 Time-domain elastodynamic DBIE, GBIE and TBIE and discussion.

The properties (14) of Stokes' tensors are used to transfer the singularity of the dynamical kernels in integrals involving the static Kelvin tensors. For example, the limiting process in (12) which defines the elastodynamic DBIE may be rewritten:

$$ \int_{\Gamma} n_s(y) \{ \Sigma_{is}^k[x, t, y | u_i(y, t)] - u_i(y, t) \Sigma_{is}^k(x, y) \} dS_y + \lim_{\epsilon \to 0} \int_{(\Gamma - e_{\epsilon}) + s_{\epsilon}} \{ u_i(y, t) \Sigma_{is}^k(x, y) - U_i^k[x, t, y | t_i(y, t)] \} dS_y = 0 \qquad (71) $$

Indeed, because of properties (14) of the Stokes tensors, such considerations can be applied as well to the other BIE considered in previous sections. Hence all results obtained for elastostatics in sections 3, 5, 6 and 7 can be extended to time-domain elastodynamics.
## 8.1 Regularization of the elastodynamic DBIE (12).

Identity (17) is written for $u^0 = u(x, t)$ and subtracted from (71). As a result, the regularized elastodynamic DBIE reads:

$$ \kappa u_k(x, t) + \int_{\Gamma} n_s(y) \left[ \Sigma_{is}^{k}[x, t, y|u_i(y, t)] - \Sigma_{is}^{k}(x, y)u_i(y, t) \right] dS_y \\ + \int_{\Gamma} n_s(y) \Sigma_{is}^{k}(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 \quad (72) $$

## 8.2 Second-order regularization of the elastodynamic GBIE (13) and TBIE.

Identity (37) is written for the static displacement field $U(y)$ defined by:

$$ U_i(y) = u_i(y, t) + u_{i,j}(y, t)(y_j - x_j) \quad (73) $$

and subtracted from (13). This manipulation leads to the second-order regularization of the elastodynamic BIE:

$$ \begin{aligned} -\kappa u_{k,r}(\mathbf{x}) + \int_{\Gamma} \{ &n_s(\mathbf{y}) [\Sigma_{is,r}^{k}[\mathbf{x}, t, y|u_i(\mathbf{y}, t)] - u_i(\mathbf{y}, t)\Sigma_{is,r}^{k}(\mathbf{x}, y)] \} dS_y \\ &+ \int_{\Gamma} [u_i(\mathbf{y}, t) - u_i(\mathbf{x}, t) - u_{i,j}(\mathbf{x}, t)r_j] n_s(\mathbf{y})\Sigma_{is,r}^{k}(\mathbf{x}, y)dS_y \\ &- \int_{\Gamma} [U_{i,r}^{k}[\mathbf{x}, t, y|t_i(\mathbf{y}, t)] - t_i(\mathbf{y}, t)U_{i,r}^{k}(\mathbf{x}, y)] dS_y \\ &- \int_{\Gamma} [t_i(\mathbf{y}, t) - C_{isab}n_s(\mathbf{y})u_{a,b}(\mathbf{x}, t)] U_{i,r}^{k}(\mathbf{x}, y)dS_y = 0 \end{aligned} \quad (74) $$

## 8.3 First-order regularization of the elastodynamic GBIE (13) and TBIE.

The elastodynamic counterpart of identity (44) is obtained using identity (27) above, together with equation (8). It holds any closed surface $S$ and for $x$ not located on $S$:

$$ \begin{aligned} & \int_S n_s(y) \Sigma_{is,r}^k [\mathbf{x}, t, y|u_i(\mathbf{y}, t)] dS_y \\ &= \int_S D_{sr} \Sigma_{is}^k [\mathbf{x}, t, y|u_i(\mathbf{y}, t)] + \int_S n_r(y) \Sigma_{is,s}^k [\mathbf{x}, t, y|u_i(\mathbf{y}, t)] dS_y \\ &= \int_S \Sigma_{is}^k [\mathbf{x}, t, y|D_{rs} u_i(\mathbf{y}, t)] dS_y + \int_S U_i^k [\mathbf{x}, t, y|\ddot{u}_i(\mathbf{y}, t)] dS_y \end{aligned} \quad (75) $$

The elastodynamic version of the first-order regularized GBIE (60) results from application of identities (75)-(44) to the integral

$$ \int_{\Gamma} n_s(y) \left[ \Sigma_{is,r}^{k} [\mathbf{x}, t, y | u_i(\mathbf{y}, t)] - u_i(\mathbf{y}, t) \Sigma_{is,r}^{k}(\mathbf{x}, y) \right] dS_y $$

This leads to state the first-order regularization of the elastodynamic GBIE as follows:

$$ \begin{aligned} & -\frac{1}{2} u_{k,r}(\mathbf{x}, t) + D_{rs} u_i(\mathbf{x}, t) A_{is}^k(\mathbf{x}, \Gamma) - t_i(\mathbf{x}, t) B_{ir}^k(\mathbf{x}, \Gamma) \\ & + \int_{\Gamma} \left[ \Sigma_{is}^k[\mathbf{x}, t, y|D_{rs} u_i(\mathbf{y}, t)] - D_{rs} u_i(\mathbf{y}, t) \Sigma_{is}^k(\mathbf{x}, y) \right] dS_y \\ & - \int_{\Gamma} \left[ U_{i,r}^k[\mathbf{x}, t, y|t_i(\mathbf{y}, t)] - t_i(\mathbf{y}, t) U_{i,r}^k(\mathbf{x}, y) \right] dS_y + \rho \int_{\Gamma} U_i^k[\mathbf{x}, t, y|\ddot{u}_i(\mathbf{y}, t)] n_r(\mathbf{y}) dS_y \\ & + \int_{\Gamma} \left[ \{D_{rs} u_i(\mathbf{y}, t) - D_{rs} u_i(\mathbf{x}, t)\} \Sigma_{is}^k(\mathbf{x}, y) - [t_i(\mathbf{y}, t) - t_i(\mathbf{x}, t)] U_{i,r}^k(\mathbf{x}, y) \right] dS_y = 0 \end{aligned} \quad (76) $$
## 8.4 Regularized elastodynamic GBIE and TBIE for cracks.

Using the same arguments than in the previous subsections, the elastodynamic counterpart of the regularized GBIE for cracks is:

$$ 
\begin{aligned}
& \frac{1}{2}(u_{i,k}^{+}(\mathbf{x}, t) + u_{i,k}^{-}(\mathbf{x}, t)) + D_{rs}\phi_{i}(\mathbf{x}, t)A_{is}^{k}(\mathbf{x}, S) - S_{i}(\mathbf{x}, t)B_{ir}^{k}(\mathbf{x}, S) \\
& \qquad + \int_{S} [\Sigma_{is}^{k}[\mathbf{x}, t, y|D_{rs}\phi_{i}(\mathbf{y}, t)] - D_{rs}\phi_{i}(\mathbf{y}, t)\Sigma_{is}^{k}(\mathbf{x}, y)] dS_{y} \\
& \qquad - \int_{S} [U_{i,r}^{k}[\mathbf{x}, t, y|S_{i}(\mathbf{y}, t)] - S_{i}(\mathbf{y}, t)U_{i,r}^{k}(\mathbf{x}, y)] dS_{y} + \rho \int_{S} U_{i}^{k}[\mathbf{x}, t, y|\tilde{\phi}_{i}(\mathbf{y}, t)]n_{r}(\mathbf{y})dS_{y} \\
& \qquad + \int_{S} [D_{rs}\phi_{i}(\mathbf{y}, t) - D_{rs}\phi_{i}(\mathbf{x}, t)]\Sigma_{is}^{k}(\mathbf{x}, y)dS_{y} - \int_{S} [S_{i}(\mathbf{y}, t) - S_{i}(\mathbf{x}, t)]U_{i,r}^{k}(\mathbf{x}, y)dS_{y} = 0
\end{aligned}
\quad (77) $$

## 8.5 Discussion of the results.

The above BIE (72), (74), (76), (77) are regularized BIE expressed with only weakly singular integrals, provided the surface fields satisfy the appropriate regularity requirements. In particular, as shown in section 4, $A_{is}^k(\mathbf{x}, S)$, $B_{ir}^k(\mathbf{x}, S)$ which appear in (76), (77) are made of contour integrals (in the case of a crack) and weakly singular surface integrals. These regularized BIEs hold for bounded solids as well as for infinite elastic media. The boundedness or unboundedness of $\Omega$ in all regularized BIEs is entirely taken into account, given the coefficient $\kappa$ (in (72) and (74)) and the orientation convention chosen for $n$.

**Comments about the limiting process.** Throughout the analysis conducted in sections 3,5, 6, 7, the BIE are defined as limiting cases of exterior representations as the size $\epsilon$ of an exclusion neighbourhood $v_\epsilon$ vanishes. During the limit process, the collocation point $\mathbf{x}$ remains fixed and is located on the boundary of $\Omega$. Moreover, the shape of $v_\epsilon$ is arbitrary throughout the limit process and may vary as $\epsilon \to 0$. These features are in contrast with other presentations of hypersingular BIEs and their regularization [6], [30], where the hypersingular BIEs:

* are defined as limiting cases of internal representations for an internal point $\mathbf{x}$.

* are formulated using FP integrals (before any regularization) or CPV integrals (where an integration by parts yields a first regularization).

In the latter the final BIE seems to depend on a particular limit process (CPV, FP), in which exclusion neighbourhoods of specific shapes have to be considered. On the contrary, the present exposition shows that the resulting regularized BIEs are truly expressed in terms of ordinary improper integrals, as the value of such integrals does not depend on the choice of exclusion neighbourhood. Such considerations were previously put forward by Guiggiani et al. [21] in their direct treatment of hypersingular integrals, see subsection 10.2.

The present application of indirect approach to DBIE and GBIE results in the derivation of regularized BIE, with weakly singular integrals; moreover, the regularization process shows the very existence of the limits which define the BIEs and their independance with respect to $v_\epsilon$. Thus the introduction of specific limit processes, such as FP integrals, for the study of GBIE is by no means essential. In a sense, the hypersingularity of the original GBIE, or the strong singularity of the original DBIE, is 'apparent' rather than actual: the regularization process consists essentially in subtracting and adding the same appropriate quantity to the original BIE (either directly or by means of an auxiliary solution), thus the final (weak) singularity may be viewed as more 'essential' than the original strong- or hypersingularity.

Establishment of the first-order regularized GBIE (76) and (77) make an essential use of identities (31), (32) and (34) given in section 4:

* The surface integrals in (31), (32) and (34) are weakly singular, as a result of integrations by parts.
* In the same time, although the contour integrals over the curve $c_\epsilon$ do not vanish individually with $\epsilon$, the nonvanishing quantities cancel each other in the very combinations $\partial J_{kr}(\mathbf{x}, \epsilon)$ (eqn. (53)), $\partial J_{kr}(\mathbf{x}, \epsilon)$ (eqn. (66)) of such integrals involved in the GBIE. This highlights the fact that the GBIE, i.e. the limiting process (35), is independent on the shape of the exclusion neighbourhood $v_\epsilon$. It is even unnecessary to select a certain shape for $v_\epsilon$, e.g. $v_\epsilon = \epsilon V$, and maintain the selected shape $V$ while taking the limit $\epsilon \to 0$. This may be view as another manifestation of the fact that the GBIE is not 'essentially' hypersingular.

**Regularity requirements for the density function and implications for BEM interpolation.** Another result of the regularization process is the regularity requirements on **u** and **t** under which the resulting BIE are actually weakly singular.

* In view of equation (72), $u(\mathbf{y}, t) \in C^{0,\alpha}$ at $\mathbf{y} = \mathbf{x}$ is necessary. This requirement is always fulfilled by the usual conformal BEM interpolations using e.g. polynomial shape functions, then regularized DBIE can be used in the usual BEM framework. This has been done in [5], [40], [42], among others.

* On the contrary, the GBIEs (74), (76), (77) require $u(\mathbf{y}, t) \in C^{1,\alpha}$ (or, equivalently, a $C^{0,\alpha}$ regularity of the displacement gradients) at $\mathbf{y} = \mathbf{x}$, and consequently $t(\mathbf{y}, t) \in C^{0,\alpha}$. If the latter condition is easy to fulfill with usual BEM discretizations, this is not the case of the former, because conformal $C^{1,\alpha}$ boundary elements are very difficult to develop for general surfaces in $\mathbb{R}^3$. This problem is by no means negligible: Krishnasamy et al. [30] point out that neglecting to ensure the $C^{1,\alpha}$ regularity at $\mathbf{y} = \mathbf{x}$ leads to scale-dependent results. The other choices available by now include:

- The use of nonconformal interpolations for **u**, i.e. using boundary elements with displacement nodes away from the element edges. Then the $C^{1,\alpha}$ requirement is met. It is generally difficult, and maybe undesirable, to make the number of collocation points match exactly the number of unknowns. This results in an expected increase of storage and CPU time, because of the additional equations and of the subsequent use of least-squares solvers which are computationally more expensive than Gaussian elimination.

- The parallel use of a $C^{0,\alpha}$ interpolation for **u** and another $C^{0,\alpha}$ for the gradients $u_{,j}$. The necessary (linear) relation between the displacement and gradient nodal values is obtained by stating that the two interpolations are nearest in the least-squares sense. This has been done by Polch, Cruse & Huang ([38]), which studied elastostatic planar crack problems using regularized TBIE.

Moreover, in the first-order regularized GBIEs (76), (77), the continuity of the displacement gradient at $\mathbf{x}$ is necessary for the free-term to make sense.

The regularity requirements for the density functions which appear along the regularization process are consistent with known theoretical results and other approaches:

- The Lyapunov-Tauber theorem in elastic potential theory [32] states that the derivative of the double-layer elastic potential ($K\phi)(\mathbf{x})$, with $K_i^k(\mathbf{x}, \mathbf{y}) = n_s(\mathbf{y})\Sigma_{is}^k(\mathbf{x}, \mathbf{y})$ at $\mathbf{x} \in \Gamma$ exists only if $\phi_i \in C^{1,\alpha}$ at $\mathbf{x}$.

- The direct approach for the evaluation of hypersingular integrals (Guiggiani et al., [21]) leads to the same regularity requirements.

**Collocation at an edge or corner point.** The regularized DBIE (72) and the second-order regularization (74) of GBIE are valid if collocated at an edge or corner point $\mathbf{x}$, provided, in the latter case, the total displacement gradients $u_i \in C^{1,\alpha}$ at $\mathbf{x}$. In both cases, no apparent free term arise due to the geometrical singularity. Reformulation (43) of (40), and its elastodynamic counterpart as well, remains valid at edge or corner points, the integrals being taken separately
on each regular component of $\Gamma$ and the quantities $D_j u_i(x)$, $n_s(x)$ being attached to each regular component of $\Gamma$.

**Frequency-domain elastostatics.** The regularized BIE for frequency-domain elastodynamics may be obtained simply by removing the time argument in $u_i$ and $t_i$ and doing the substitutions $u_i(y)\Sigma(x, y, \omega)$ to $\Sigma[x, t, y|u_i(y, t)]$ and the like, in equations (72), (74), (76), (77). The transfer of singularity from Helmholtz to Kelvin kernels uses series expansions (in infinite series of $k_\alpha r$) of the Helmholtz kernels [5], [6].

# 9 Numerical implementation of the regularization approach.

A detailed discussion of the whole discretization procedure of the time-domain BEM is beyond the scope of this paper. This section focuses on the treatment, at the spatial interpolation of geometry and unknowns level, which allows full use of the regularization approach and ensures a proper numerical evaluation of the singular integrals. Hence this discussion can be restricted without harm to elastostatic DBIE and GBIE.

## 9.1 Discretization of $\Gamma$ and $u$.

The surface $\Gamma$ is divided into boundary elements. They are mapped, in a standard fashion, on a reference element $E_0$, which is generally the square $\xi = (\xi_1, \xi_2) \in [-1, 1]^2$ or the triangle $0 \le \xi_1 + \xi_2 \le 1$. Let the discussion be restricted to the consideration of a single element $E$.

The location of a point $y$ on $E$ is expressed in terms of $n$ shape functions $N^k$ and $n$ geometrical nodes $A^k$ ($k = 1...n$):

$$ Oy = N^k(\xi)OA^k \quad (\xi \in E_0) \quad (78) $$

The nodes $A_k$ are located on the boundary $\partial E$ of $E$, so that the geometry interpolation is conformal. Then the natural basis ($a_\alpha$), metric tensor ($g_{\alpha\beta}$) and unit normal $n$ on $E$ are given by:

$$ \begin{gathered} a_{\alpha}(\xi) = N_{,\alpha}^{k}(\xi)OA_{k} \qquad g_{\alpha\beta}(\xi) = a_{\alpha}(\xi) \cdot a_{\beta}(\xi) \qquad (\xi \in E_0) \\ \sqrt{\hat{g}(\xi)}n(\xi) = a_1 \wedge a_2 \qquad \hat{g}(\xi) = (g_{11}g_{22} - g_{12}^2)(\xi) \end{gathered} \quad (79) $$

The displacement $u(y)$ on $E$ is interpolated using $m$ shape functions $M^q$ and $m$ nodal values $u^q$ ($q = 1...m$). According to the discussion of subsection 8.5, the displacement nodes may be either boundary nodes or internal nodes (nonconformal interpolation) and $m \neq q$ is possible.

The shape functions discussed here are polynomials of $(\xi_1, \xi_2)$.

## 9.2 Numerical evaluation of singular integrals in DBIE.

Singular integrals occur if $E$ contains the collocation point $x$, which may be neither a geometrical node nor a displacement node. The regularization procedure leads to singular integrals of the form:

$$ I^s = \int_E n_s(y) \Sigma_{is}^k(x, y) (u_i(y) - u_i(x)) dS_y \quad (80) $$

Let $\eta = (\eta_1, \eta_2)$ denote the antecedent of $x$ on $E_0$. The following definition is introduced: for any polynomial $P(\xi)$, the polynomials $P_\alpha^I$, $P_{\alpha\beta}^{II}$ are constructed as:

$$ \begin{aligned} P(\xi) &= P(\eta) + (\xi_{\alpha} - \eta_{\alpha})P_{\alpha}^{I}(\xi; \eta) \\ P(\xi) &= P(\eta) + (\xi_{\alpha} - \eta_{\alpha})P_{,\alpha}(\eta) + 1/2(\xi_{\alpha} - \eta_{\alpha})(\xi_{\beta} - \eta_{\beta})P_{\alpha\beta}^{II}(\xi; \eta) \end{aligned} \quad (81) $$
Following a common practice in BEM (see e.g. [42]), set $ξ₁ = ρ \cos α$, $ξ₂ = ρ \sin α$. Then $dξ₁dξ₂ = ρdρdα$ and, from (81):

$$ 
\begin{aligned}
M^q(\xi) - M^q(\eta) &= (\xi_\alpha - \eta_\alpha) M_\alpha^{I,q}(\xi; \eta) = \hat{M}^q(\rho, \alpha; \eta) \\
r(x, y) &= \| (\xi_\alpha - \eta_\alpha) N_\alpha^{I,k}(\xi; \eta) OA^k \| = \rho \hat{r}(\rho, \alpha; \eta) \\
\Sigma_{is}^k(x, y) &= \frac{1}{\rho^2} \hat{\Sigma}_{is}^k(\rho, \alpha; \eta)
\end{aligned} 
\quad (82) $$

where $\hat{r}(\rho, \alpha; \eta) \neq 0$ and $(\hat{\Sigma}(\rho, \alpha; \eta)$ is regular at $\rho = 0$. Hence integral (80) is recast in a completely regular form, as:

$$ I^s = u_i^q \int_E n_s(y) \hat{\Sigma}_{is}^k(\rho, \alpha; \eta) \hat{M}^q(\rho, \alpha; \eta) \sqrt{\hat{g}(\xi)} d\rho d\alpha \quad (83) $$

Expression (83) takes full advantage of the regularization. The numerical evaluation of (83) can be performed with standard product Gaussian quadrature formulas, the complete procedure requiring a further coordinate change $(\rho, \alpha) \to (v_1, v_2)$ in order to recover an integral over the square $[-1, 1]^2$ [5], [42].

## 9.3 Numerical evaluation of singular integrals in GBIE and TBIE.

A typical integral occurring in second-order regularized GBIE is:

$$ I^s = \int_E [u_i(y) - u_i(x) - r_j D_j u_i(x)] n_s(y) \Sigma_{is,r}^k(x, y) dS_y \quad (84) $$

Using the intrinsic expression of grad_S, definition (81) and $a_\alpha(\eta) \cdot a^\gamma(\eta) = \delta_\gamma^\alpha$, one has on E:

$$ 
\begin{aligned}
r_j D_j u_i(x) &= (\xi_\alpha - \eta_\alpha) \{a_\alpha(\eta) + 1/2(\xi_\beta - \eta_\beta) N_{\alpha\beta}^{II,k}(\xi; \eta) OA^k\} \cdot u_i^q M_{\gamma,\alpha}^q a^\gamma(\eta) \\
&= (\xi_\alpha - \eta_\alpha) u_i^q \{M_{\alpha,\gamma}^q + 1/2(\xi_\beta - \eta_\beta) N_{\alpha\beta}^{II,k}(\xi; \eta) M_{\gamma,\alpha}^q (OA^k \cdot a^\gamma(\eta))\} \\
u_i(y) - u_i(x) &= u_i^q (\xi_\alpha - \eta_\alpha) \{M_{\alpha,\gamma}^q + 1/2(\xi_\beta - \eta_\beta) N_{\alpha\beta}^{II,q}(\xi; \eta) M_{\gamma,\alpha}^q\}
\end{aligned} 
\quad (85) $$

$$ 
\begin{aligned}
u_i(y) - u_i(x) - r_j D_j u_i(x) \\
&= 1/2 u_i^q (\xi_\alpha - \eta_\alpha)(\xi_\beta - \eta_\beta) \{ M_{\alpha\beta}^{II,q}(\xi;\eta) - N_{\alpha\beta}^{II,k}(\xi;\eta) M_{\gamma,\alpha}^q (OA^k \cdot a^\gamma(\eta)) \}
\end{aligned} 
\quad (86) $$

Then, because of (82) and (86), $\Sigma_{is,r}^k(x, y) \sim \rho^{-3}$, $u_i(y) - u_i(x) - r_j D_j u_i(x) \sim \rho^2$ and $dS_y \sim \rho d\rho d\alpha$, integral (84) may be recast in a completely regular form in the system $(\rho, \alpha)$. This requires the analytical derivation of $M_{\alpha\beta}^{II,q}(\xi; \eta)$ and $N_{\alpha\beta}^{II,k}(\xi; \eta)$ for given shape functions $M^q$ and $N^k$.

In first-order regularized GBIEs (76), (77), the following type of singular integral occur:

$$ I^s \doteq \int_E [D_{rs} u_i(y) - D_{rs} u_i(x)] \Sigma_{is}^k(x, y) dS_y \quad (87) $$

Using the intrinsic expression of grad_S, definition (81), $a_\alpha(\eta) \cdot a^\gamma(\eta) = \delta_\gamma^\alpha$ and $dS_y = \sqrt{\hat{g}(\xi)} d\xi_1 d\xi_2$, one has on E:

$$ D_{rs} u_i(y) dS_y = \epsilon_{prs} [a_2(\xi) M_{,1}^q(\xi) - a_1(\xi) M_{,2}^q(\xi)] \cdot e_p d\xi_1 d\xi_2 \quad (88) $$

and $(D_{rs}u_i(y) - D_{rs}u_i(x))dS_y$ is rewritten:

$$ 
\begin{aligned}
(D_{rs}u_i(y) - D_{rs}u_i(x))dS_y \\
&= \left\{ (1 - \sqrt{\hat{g}(\eta)/\hat{g}(\xi)}) [a_2(\eta)M_{,1}^q(\eta) - a_1(\eta)M_{,2}^q(\eta)] \right. \\
&\qquad \left. + [(a_2(\xi) - a_2(\eta))M_{,1}^q(\xi) - (a_1(\xi) - a_1(\eta))M_{,2}^q(\xi)] \right. \\
&\qquad \left. + a_2(\eta)(M_{,1}^q(\xi) - M_{,1}^q(\eta)) - a_1(\eta)(M_{,2}^q(\xi) - M_{,1}^q(\eta)) \right\} \cdot e_p \epsilon_{prs} d\xi_1 d\xi_2
\end{aligned} 
\quad (89) $$
Then, using (79) and definition (81), one can put:

$$ 
\begin{align*}
a_{\alpha}(\xi) - a_{\alpha}(\eta) &= (\xi_{\beta} - \eta_{\beta})(N_{,\alpha}^{k})_{\beta}^{I}(\xi; \eta)OA_{k} \\
M_{,\alpha}^{q}(\xi) - M_{,\alpha}^{q}(\eta) &= (\xi_{\beta} - \eta_{\beta})(M_{,\alpha}^{q})_{\beta}^{I}(\xi; \eta) \\
\hat{g}(\xi) - \hat{g}(\eta) &= (\xi_{\beta} - \eta_{\beta})\hat{g}_{\beta}^{I}(\xi; \eta) \\
1 - \sqrt{\hat{g}(\eta)/\hat{g}(\xi)} &= -(\xi_{\beta} - \eta_{\beta})\frac{\hat{g}_{\beta}^{I}(\xi; \eta)}{1 + \sqrt{\hat{g}(\eta)/\hat{g}(\xi)}} \tag{90}
\end{align*}
$$

Then, because of (82) and (86), $(D_{rs}u_i(y) - D_{rs}u_i(x))dS_y \sim \rho^2 d\rho d\alpha$. Integral (87) may be recast in a completely regular form in the system $(\rho, \alpha)$. This requires the analytical derivation of $(N_{,\alpha}^k)_{\beta}^I(\xi; \eta)$, $(M_{,\alpha}^q)_{\beta}^I(\xi; \eta)$ and $\hat{g}_{\beta}^I(\xi; \eta)$ for given shape functions $M^q$ and $N^k$.

# 10 Overview of variational and direct approaches for hypersingular BIEs.

For completeness, a brief account of two other approaches to formulate and deal with hypersingular BIEs, namely the variational BIE approach and the direct approach for evaluation of hypersingular integrals.

## 10.1 Variational approach for static TBIE.

Roughly speaking, equation (1) is multiplied by a test function $\psi(x)$ belonging to some appropriate function space and integrated on $S$ with respect to $x$. The regularization here stems from the fact that the singular kernel $K(x, y)$, which appears in the bilinear form of the variational equation, is integrable over $(x, y) \in S \times S$, even in the case $K(x, y) = \sum_{i,r}^k (x, y)$. However, in the actual numerical computation of the (discrete counterpart of) the bilinear form, the two integrations over $S$ have to be done sequentially. As a consequence, the bilinear form has to be rewritten in terms of weakly singular kernels (i.e. integrable over $S$). This is done using two integrations by parts (one for each variable $x$ and $y$), e.g. in [36], [10], [23]. Actual statements of variational BIE formulations are generally derived from the variational theorems of elasticity.

This approach is investigated by Nedelec and co-workers (see e.g. Nedelec [36], Bamberger & Ha Duong [1]-[2], Becache [3]), Hamdi [23], Polizzotto, Maier and co-workers (see e.g. [39], [34], [43]), and others. This approach is conceptually better than the usual collocation method: it allows convergence study and lead to symmetric BEM matrices. Moreover, the variational TBIE requires $C^{0,\alpha}$ interpolation of the densities instead of $C^{1,\alpha}$ (at $x$) for the collocation TBIE. The usual shape functions may then be used. On the other hand, derivation of the appropriate expressions for the bilinear forms as well as the implementation of this approach for general engineering problems is more involved than using collocation. Numerical applications of this approach may be found e.g. in [23] (exterior acoustics), [10], [15] (elastodynamic crack problems), [43].

## 10.2 Direct computation of hypersingular integrals.

The CPV and FP integrals are defined in terms of specific limit processes, using exclusion neighbourhood of specified shape and vanishing size. In the direct approach, no prior regularization is performed but the limiting process is carefully preserved in the mapping between the physical element and the reference element used for the numerical evaluation of element integrals, which involve shape functions. This approach allows the direct computation of any strongly singular or hypersingular integral arising in BIE methods. The development of this approach (Guiggiani and co-workers [20], [21]) is recent. In [21], the GBIE (35) for elastostatic problems is formulated (using the notations of section 5) becomes:

$$ c_{pqr} u_{p,q}(\mathbf{x}) + \lim_{\epsilon \to 0} \int_{\Gamma - \epsilon} \left\{ u_i(\mathbf{y}) n_s(\mathbf{y}) \sum_{is,r}^k (\mathbf{x}, \mathbf{y}) - t_i(\mathbf{y}) U_{i,r}^k(\mathbf{x}, \mathbf{y}) + u_i(\mathbf{x}) \frac{b_{kri}(\mathbf{x})}{\epsilon} \right\} dS_y = 0 \quad (91) $$
(where $c_{pqr}u_{p,q}(\mathbf{x})$ and $b_{kr_i}(\mathbf{x})$ are known) and they show that formulation (91) does not depend on the shape chosen for $v_\epsilon$, and hence on the shape of $e_\epsilon$. Numerical examples on element integrals are very good.

# 11 Numerical examples.

## 11.1 Incident plane wave on a spherical cavity. [5].

In this example the regularized DBIE is applied to the problem of scattering of a time-harmonic (pulsation $\omega$) incident plane longitudinal wave by a spherical cavity of radius $R$ and center $O$. The wave propagates along the $x_3$-direction and has amplitude $u^I$. In view of the axisymmetry of this problem, the only nonzero displacement components are $u_r, u_\theta$, where $(O, r, \theta, \phi)$ denotes a spherical coordinate system, the points $\phi = 0, \phi = \pi$ lying on $Ox_3$. This problem has an exact solution [18].

The moduli of the surface ($r/R = 1$) and far-field ($r/R = 100$) displacements are presented, in the form of a polar diagram, on figure 4 for the case $\omega R\sqrt{\rho/(\lambda + 2\mu)} = 3$. Using the symmetry capabilities of our code [8], only one-eighth of the spherical surface is meshed, with 12 eight-noded quadrilateral elements and 49 nodes in the present case. The results show very good agreement between the numerical values (symbols) and the analytical solution (curves).

## 11.2 Dynamical propagation of a crack in antiplane strain [27].

This example illustrates the use of the regularized TBIE in time-domain elastodynamics.

The spontaneous propagation, in antiplane strain, of a straight crack $C$ extending along the $x_1$-axis in an infinite elastic space is considered. The normal direction of $C$ is the $x_2$-direction. Via the superposition principle, $C$ is loaded by a shear traction: $T^\pm(\mathbf{y}) = \pm\tau H(t)e_3$, where $\tau$ is a constant and $H(t)$ is the Heaviside step function. The propagation of $C$ is governed by $K_{III} = K_{III}^c$, where $K_{III}$ and $K_{III}^c$ respectively denote the mode III dynamical stress intensity factor and the toughness of the elastic material. The left end of the crack is kept fixed, so that only the right end propagates, according to the above criterion. Let $\ell(t)$ denote the length of $C$ at time $t$, the initial length being $\ell(0)$.

The regularized TBIE for cracks in antiplane strain is obtained by integrating equation (77) from $y_3 = -\infty$ to $y_3 = +\infty$.

The crack is discretized into $J$ elements of equal length $\Delta x$, while the time interval $[0,T]$ of interest is split into $I$ equal time intervals $\Delta t$. The ratio $c\Delta t/\Delta x$ has been set to the value $1/2$, which allows to perform most of the integrations analytically. This has been done in order to circumvent inaccuracy problems arising when dealing with numerical integration of time-domain kernels and related to causality considerations [25].

The only nonzero component $\phi_3(\mathbf{y},t) = \phi(\mathbf{y},t)$ is interpolated linearly in space and in time. The regularized TBIE is collocated at the midpoint of each boundary element and at $t = i\Delta t$, $i = 1,2,\dots,I$. The linear system of equations which is to be solved at each time step has $J-1$ unknowns and $J$ equations, hence it is considered in a least squares sense, using the Householder factorization of the matrix (LINPACK software library, [17]).

The propagation is simulated by adding a new element if the propagation criterion is met, so that $I$ increases during the time-stepping scheme. Details about the discratization procedure and the numerical treatment of the propagation criterion may be found in [27].

The numerical results for the spontaneously propagating crack compare very favourably with the analytical solution for the semi-infinite spontaneously propagating crack in antiplane strain given by Kostrov [28]. Moreover their accuracy is better than those obtained for the same problem in [46] using a finite-difference method, especially for the initiation phase, which is delayed due to the poor stress resolution of the FD technique. Figure 5 shows our numerical results for $\ell(t) - \ell(0)$ compared with those of [46] and with the analytical value given in [28].
<figure><img src="image_5.png" /><figcaption>Figure 4: Surface and far-field (for $r/R = 100$) displacements for the problem of subsection 11.1. The ratios $u_r/u^l$ and $u_\theta/u^l$ are depicted in polar diagrams; $\theta$ ranges in $[0, \pi]$ for each component, because of the symmetry of the problem under consideration.</figcaption></figure>

## 12 Concluding remarks.

The regularization of the strongly singular and hypersingular collocation boundary integral equations arising for 3D general elastodynamic problems has been derived and stated. The main results are equations (72), (74), (76) and (77), in which all integrals are ordinary improper integrals, which can be computed accurately using standard numerical quadrature methods. The regularized TBIEs allow the modelling of cracks and the computation of the entire stress tensor on the boundary.

The limiting process used to establish the regularized BIEs shows that they do not depend on a specific limit process such as Cauchy Principal Value or Finite Part. Moreover, the regularity requirements on the densities are natural consequences of the regularization process.

Integration by parts identities and transformation of integrals over open surfaces of the Kelvin tensors play a key role in the derivation of the first-order regularized GBIE and TBIE.

As a global conclusion, all usual BIE are regularizable using indirect approach, and the occurrence of highly singular kernels in the BIE associated to elasticity problems (and scalar potential problems as well) should not be a concern when implementing a BE method.
<figure><img src="image_6.png" /><figcaption>Figure 5: Crack-tip location during the spontaneous rupture propagation of a semi-infinite crack. Our results are compared to Kostrov's [28] exact solution and Virieux & Madariaga [46] numerical results (obtained using FDM). The results are normalized, and $K'$ denotes the nondimensional quantity $K' = K_{III}^c / (\mu\sqrt{\pi}\Delta x)$.</figcaption></figure>

Numerical implementation of these regular BIE is done and successfully tested for uncracked solids in static and steady-state elasticity, and for propagating cracks in 2-D time-domain elastodynamics. Numerical implementation of regularized TBIEs for 3-D situations is under investigation. In the authors' opinion, the study of interpolation strategies for GBIE/TBIE (in connexion with the $C^{1,\alpha}$ requirement at the collocation point) and of accurate spatial numerical integration algorithms for time-domain BIEs (accuracy problems arise from causality considerations) deserve attention.

The present regularization approach can be developed the same way for 2-D elastodynamics and 3-D or 2-D scalar potential or wave problems.

Références

[1] BAMBERGER A., HA DUONG T. - Formulation variationnelle espace-temps pour le calcul par potentiel retardé de la diffraction d'une onde acoustique (I). Math. Meth. in the Appl. Sci. vol. 8 pp 405-435, 1986.

[2] BAMBERGER A., HA DUONG T. - Formulation variationnelle pour le calcul de la diffraction d'une onde acoustique par une surface rigide. Math. Meth. in the Appl. Sci. vol. 8 pp 598-608, 1986. 1987.

[3] BECACHE E. - Résolution par une méthode d'équations intégrales d'un problème de diffraction d'ondes élastiques transitoires par une fissure. Thèse de doctorat, Université Paris 6, 1991.

[4] BESKOS D.E. (Editor) - Boundary Element Methods in Mechanics. Computational Methods in Mechanics, vol. 3, North Holland 1987.

[5] 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, nr. 1/2, 1987.

[6] BONNET M. - Regularized Boundary Integral Equations for Three-dimensional Bounded or Unbounded Elastic Bodies Containing Curved Cracks of Arbitrary Shape Under Dynamic Loading. In "Boundary Element Techniques: Applications in Engineering" (ed. C.A. Brebbia & N.G. Zamani), Computational Mechanics Publications (Southampton), 1989.

[7] BONNET M. - Shape differentiation of regularized BIE: application to 3-D crack analysis by the virtual crack extension approach. In "Boundary Elements in Mechanical and Electrical Engineering" (C.A. Brebbia and A. Chaudouet, eds.), Computational Mechanics Publications, Springer-Verlag, 1990.
[8] BONNET M. - On the use of geometrical symmetry in the boundary element methods for 3D elasticity. *In Boundary Element Technology*, C.A. Brebbia editor, Computational Mechanics Publications, Elsevier, 1991.

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

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

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

[12] 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.

[13] 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, Paris, tome 300, 1985.

[14] BUI H.D. - A boundary integral equation approach to fracture mechanics in three-dimensional and related problems. 6th National Japonese Conf. on Boundary Elements Methods., Tokyo, December 1989).

[15] CORTEY-DUMONT Ph - Simulation numérique de problèmes de diffraction d'ondes par une fissure. Thèse d'Etat, Université Paris VI, 1985.

[16] T.A. CRUSE - Boundary Element Analysis in Computational Fracture Mechanics. Kluwer Academic Publishers, 1988.

[17] DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.B. - LINPACK users'guide. SIAM, Philadelphia, 1979.

[18] ERINGEN A.C., SUHUBI E.S. - Elastodynamics (vol II - linear theory. Academic Press 1975.

[19] GIRAUD G. Equations á intégrales principales, étude suivie d'une application. Ann. Ec. Normale, (3), LI, fasc. 4, 1936, Paris.

[20] GUIGGIANI M., GIGANTE A. - A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method. ASME J. Appl. Mech. (to appear).

[21] GUIGGIANI M., KRISHNASAMY G., RUDOLPHI T.J., RIZZO F.J. - A general algorithm for the numerical solution of hypersingular boundary integral equations. ASME j. Appl. Mech., to appear.

[22] HA DUONG T. - Equations intégrales pour la résolution numérique de problèmes de diffraction d'ondes acoustiques dans $R^3$. Thèse d'Etat, Université Paris VI, 1987.

[23] HAMDI M.A. - Formulation variationnelle par équations intégrales pour le calcul de champs acoustiques linéaires proches et lointains. Thèse d'Etat, UTC Compiègne, 1982.

[24] HIROSE S. & ACHENBACH J.D. - Time-domain boundary element analysis of elastic wave interaction with a crack. Int. J. Num. Meth. in Eng. vol. 28, pp 629-644, 1989. foundations pp

[25] KARABALIS D.L. - A simplified 3-D time-domain BEM for dynamic soil-structure interaction problems. Eng. Anal. with Boundary Elements, vol.8 no 3, pp.139-145, 1991.

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

[27] KOLLER M.G., BONNET M., MADARIAGA R. - Modelling of dynamical crack propagation using time-domain boundary integral equations. Submitted to Wave Motion, August 1991.

[28] KOSTROV B.V. - On the crack propagation with variable velocity. Intern. J. Fract., vol.11 n°1, pp.47-55, 1975.

[29] KRISHNASAMY G., SCHMERR L.W., RUDOLPHI T.J. & RIZZO F.J. - Hypersingular boundary integral equations: some applications in acoustic and elastic wave scattering. ASME J. Appl. Mech., vol. 57, pp. 404-414, 1990.

[30] KRISHNASAMY G., RIZZO F.J., RUDOLPHI T.J. - Hypersingular boundary integral equations: their occurrence, interpretation, regularization and computation. To appear in "Developments in Boundary Element Methods" - vol. 7: Advanced Dynamic Analysis, chap. 7 (P.K. Banerjee and S. Kobayashi, eds.).

[31] KUPRADZE V.D. - Dynamical problems in elasticity. Progress in solids mechanics, vol. III, North Holland, 1963.
[32] KUPRADZE V.D. (Editor) - Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity. North Holland, 1979.

[33] LEBLOND J.B. - Equations intégrales régularisées pour un corps bidimensionnel contenant une fissure de forme quelconque. C.R.A.S.

[34] MAIER G., DILIGENTI M. & CARINI A. - A variational formulation of the boundary integral equation method in elastodynamics. In Boundary element methods in engineering, pp.268-275, Springer, 1990. equation

[35] MIKHLIN S.G. - Multi-dimensional singular integrals and singular integral equations. Pergamon Press, 1965.

[36] NEDELEC J.C. - Integral equations with non integrable kernels. Integral equations and operator theory, vol. 5, pp 562-572, 1982.

[37] NISHIMURA N. & KOBAYASHI S. - A regularized boundary integral equation method for elastodynamic crack problems. Comp. Mech. vol. 4 pp 319-328, 1989.

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

[39] POLIZZOTTO C. - An energy approach to the boundary element method. Part I: elastic solids. Comp. Meth. Appl. Mech. Eng. vol. 69, 1988, pp 167-184.

[40] REZAYAT M., SHIPPY D.J. & RIZZO F.J. - On time-harmonic elastic wave analysis by the boundary element method for moderate to high frequencies. Comp. Meth. Appl. Mech. Eng. vol. 55, 1986, pp 349-367.

[41] 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.

[42] 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.

[43] SIRTORI S., MAIER G., NOVATI G., MICCOLI S. - A Galerkin symmetric boundary element method in elasticity: formulation and implementation. 1990, to appear.

[44] SLADEK J., SLADEK V. - Three-dimensional curved crack in an elastic body. Int. J. Solids Struct. vol. 19, pp 425-436, 1983.

[45] 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.

[46] VIRIEUX J., MADARIAGA R. - Dynamic faulting studied by a finite difference method. Bull. Seism. Soc. Am., vol.72 n°1, pp.345-369, 1982. by

[47] WEAVER J. - Three-dimensional crack analysis. Int. J. Solids Struct. vol. 13, pp 321-330, 1977.

[48] WENDLAND W.L. - Asymptotic accuracy and convergence for point collocation methods. in [12] vol 2.

## A Elastodynamic kernels.

$$r^2 = \|x - y\|^2 = (y_1 - x_1)^2 + (y_2 - x_2)^2 + (y_3 - x_3)^2$$
$$r_{,i} = \frac{y_i - x_i}{r} \quad r_{,i}r_{,i} = 1 \quad r_{,ij} = \frac{1}{r}(\delta_{ij} - r_{,i}r_{,j}) \quad r_{ijk} = \frac{1}{r^2}(3r_{,i}r_{,j}r_{,k} - \delta_{ij}r_{,k} - \delta_{ik}r_{,j} - \delta_{jk}r_{,i})$$
$$c_L = \sqrt{(\lambda + 2\mu)/\rho} \quad c_T = \sqrt{\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)}} = \sqrt{\mu/(\lambda+2\mu)} \quad \text{with} \quad k_\alpha = \frac{\omega}{c_\alpha} \quad (\alpha = \text{L or T}) \qquad (92)$$

### A.1 Kelvin displacement tensor U and stress tensor Σ (static case).

$$U_i^k(x, y) = \frac{1}{8\pi\mu r} [(1 - \beta^2)r_{,i}r_{,k} + (1 + \beta^2)\delta_{ik}] \qquad (93)$$

$$\Sigma_{is}^k(x, y) = \frac{1}{4\pi r^2} [3(1 - \beta^2)r_{,i}r_{,k}r_{,s} + \beta^2(\delta_{ik}r_{,s} + \delta_{sk}r_{,i} - \delta_{is}r_{,k})] \qquad (94)$$
## A.2 Stokes displacement tensor $U$ and stress tensor $\Sigma$ (time-dependent case).

$$ U_i^k(x, t, y|f) = \frac{1}{4\pi\mu r} [Br_i, r, k + (A + f(t - r/c_T))\delta_{ik}] \quad (95) $$

$$ \Sigma_{is}^k(x, t, y|f) = \frac{1}{4\pi r^2} [2Cr_i, r, k, r, s + \delta_{is}r, kF_L + (\delta_{ik}r, s + \delta_{sk}r, i)F_T] $$

$$ \begin{aligned} A &= c_T^2 \int_{1/c_T}^{1/c_L} \lambda f(t - \lambda r) d\lambda \\ B &= -3A - (f(t - v_T) - \beta^2 f(t - v_L)) \\ C &= -5B + (f(t - v_T) + v_T \dot{f}(t - v_T)) - \beta^2 (f(t - v_L) + v_L \dot{f}(t - v_L)) \end{aligned} \quad (96) $$

$$ \begin{aligned} F_T &= 2B - (f(t - v_T) + v_T \dot{f}(t - v_T)) \\ F_L &= 2B - (1 - 2\beta^2) (f(t - v_L) + v_L \dot{f}(t - v_L)) \\ v_T &= r/c_T \qquad v_L = r/c_L \end{aligned} $$

The following Taylor expansions for small $r$ hold:

$$ \begin{aligned} A &= -\frac{1}{2}(1 - \beta^2)f(t) + \frac{1}{3}(1 - \beta^3)v_T \dot{f}(t) - \frac{1}{8}(1 - \beta^4)v_T^2 \ddot{f}(t) + o(r^2) \\ A + f(t - r/c_T) &= \frac{1}{2}(1 + \beta^2)f(t) - \frac{1}{3}(2 + \beta^3)v_T \dot{f}(t) + \frac{1}{8}(3 + \beta^4)v_T^2 \ddot{f}(t) + o(r^2) \\ B &= \frac{1}{2}(1 - \beta^2)f(t) - \frac{1}{8}(1 - \beta^4)v_T^2 \ddot{f}(t) + o(r^2) \\ C &= -\frac{3}{2}(1 - \beta^2)f(t) + \frac{1}{8}(1 - \beta^4)v_T^2 \ddot{f}(t) + o(r^2) \\ F_T &= -\beta^2 f(t) + \frac{1}{4}(1 + \beta^4)v_T^2 \ddot{f}(t) + o(r^2) \\ F_L &= \beta^2 f(t) - \frac{1}{4}(1 - 2\beta^2 + 3\beta^4)v_T^2 \ddot{f}(t) + o(r^2) \end{aligned} \quad (97) $$

They show that:

* the lowest-order terms yields $f(t)$ times the singular static Kelvin kernels (see (93), (94)).

* the differences $U[x, t, y|f] - f(t)U(x, y)$ and $\Sigma[x, t, y|f] - f(t)\Sigma(x, y)$ remain finite for arbitrary small $r$. (14).

The properties (14) of the Stokes tensors are a consequence of the expansions (97) and the above remarks. The similar expansions (in infinite series of $k_\alpha r$) of the Helmholtz kernels can be found e.g. in [5], [6].

## B Proof of the identities (31) (32) (34).

* $I_a(x, S)$ is easily established by noticing that:

$$ \frac{r_{,a}}{r^2} = -\left(\frac{1}{r}\right)_{,a} = -\left[D_a \frac{1}{r} + n_a \left(\frac{1}{r}\right)_{,n}\right] \quad (98) $$

and applying identity (25) with $f = n_a 1/r$.

* From (92), one has:

$$ 3J_{abc}(x, S) = -\delta_{ab}I_c(x, S) - \delta_{ac}I_b(x, S) - \delta_{bc}I_a(x, S) + \int_S r_{,abc}dS_y \quad (99) $$
Then, putting $r_{,abc} = D_c r_{,ab} + n_c n_p r_{,abp}$ and applying identity (23) with $f = r_{,ab}$ leads to

$$
\begin{align*}
\int_S r_{,abc} \frac{dS_y}{r^2} &= \int_S K n_c r_{,ab} dS_y + \int_S n_c n_p r_{,abp} dS_y + \int_{\partial S} r_{,ab} \nu_c dS_y \\
&= \int_S K n_c r_{,ab} dS_y + \int_S n_c (D_{pa} r_{,bp} + n_a r_{,bpp}) dS_y + \int_{\partial S} r_{,ab} \nu_c dS_y \\
&= \int_S K n_c r_{,ab} dS_y + \int_S n_c \left( D_{ba} \frac{1}{r} + n_a \left( \frac{2}{r} \right)_{,b} - D_{pa} (r_{,b} r_{,p}) \right) dS_y + \int_{\partial S} r_{,ab} \nu_c dS_y \\
&= \int_S K n_c r_{,ab} dS_y + \int_S n_c \left( n_a \left( \frac{1}{r} \right)_{,b} + n_b \left( \frac{1}{r} \right)_{,a} - D_{pa} (r_{,b} r_{,p}) \right) dS_y + \int_{\partial S} r_{,ab} \nu_c dS_y \quad (100)
\end{align*}
$$

because of $r_{,bpp} = (2/r)_{,b}$. Then, upon application of identity (23) to integrands $n_a n_c (1/r)$ and $n_a n_b (1/r)$ and of identity (26) to integrand $D_{pa}(r_b r_p)$, one gets:

$$
\begin{align}
\int_S r_{,abc} \frac{dS_y}{r^2} &= \int_S K n_c r_{,ab} dS_y + \int_S 2n_c n_a n_b K \frac{ds_y}{r} - \int_S 2n_c n_a n_b \frac{ds_y}{r^2} \nonumber \\
&\quad - \int_S (D_b(n_c n_a) + D_b(n_c n_a)) \frac{ds_y}{r^2} + \int_S r_{,b} r_{,p} D_{pa} n_c \nonumber \\
&\quad + \int_{\partial S} (\nu_b n_c n_a + \nu_b n_c n_a) \frac{ds_y}{r} - \epsilon_{paq} \int_{\partial S} n_c r_{,b} r_{,p} \tau_q dS_y + \int_{\partial S} r_{,ab} \nu_c dS_y \tag{101}
\end{align}
$$

Finally, the terms in (100) are rearranged using (31) and noticing that $\epsilon_{paq}\tau_q = n_p\nu_a - n_a\nu_p$.
Expression (32) of $J_{abc}(x, S)$ is then readily get from (99).

• Expression (34) of integral $K_i^k(x, S)$ results from the following manipulation:

$$
\int_S (n_k r_i - n_i r_k) \frac{dS_y}{r^2} = \int_S D_{ik} \frac{1}{r} dS_y = \epsilon_{ikp} \int_{\partial S} \tau_p \frac{dS_y}{r} \quad (102)
$$