Boundary integral equations for isotropic linear elasticity

TL;DR

This paper develops boundary integral equations for 3D isotropic linear elasticity in biphasic models, analyzes spectra for spherical inclusions, and proposes a Galerkin discretization for multiple inclusions with numerical validation.

Contribution

It introduces a new integral equation approach and Galerkin discretization for multiple spherical inclusions with different Lamé parameters in isotropic elasticity.

Findings

Eigenfunctions of boundary operators for spherical inclusions identified

Spectra of boundary operators explicitly computed for spherical cases

Numerical tests demonstrate effectiveness of the proposed discretization

Abstract

This articles first investigates boundary integral operators for the three-dimensional isotropic linear elasticity of a biphasic model with piecewise constant Lam\'e coefficients in the form of a bounded domain of arbitrary shape surrounded by a background material. In the simple case of a spherical inclusion, the vector spherical harmonics consist of eigenfunctions of the single and double layer boundary operators and we provide their spectra. Further, in the case of many spherical inclusions with isotropic materials, each with its own set of Lam\'e parameters, we propose an integral equation and a subsequent Galerkin discretization using the vector spherical harmonics and apply the discretization to several numerical test cases.

Tables

Table 1. Table 1: Degree N g subscript 𝑁 𝑔 N_{g} and number of points T g subscript 𝑇 𝑔 T_{g} of Lebedev quadrature rules such that spherical haromonics up to degree N g subscript 𝑁 𝑔 N_{g} are integrated exactly with T g subscript 𝑇 𝑔 T_{g} points.

$N_g$	3	5	7	9	11	13	15	17	19	21	23	25	27	29
$T_g$	6	14	26	38	50	74	86	110	146	170	194	230	266	302
$N_g$	31	35	41	47	53	59	65	71	77	83	89	95	101	107
$T_g$	350	434	590	770	974	1202	1454	1730	2030	2354	2702	3074	3470	3890

Table 2. Table 2: Relative error 𝚁𝚎 𝚁𝚎 {\tt Re} of the approximation to the one-sphere model.

$N$	Case 1		Case 2		Case 3		Case 4
2	0		0		1.985e-01		5.375e-01
5	0		0		4.020e-03		6.797e-02
8	0		0		4.796e-09		1.370e-04
11	0		0		0		7.892e-13
14	0		0		0		7.097e-13

Table 3. Table 3: Geometric configuration of the case study for convergence with respect to the degree of vector spherical harmonics N 𝑁 N involving three spheres.

Set	Sphere	Center	Radius	Lamé constants	Stress tensor	Transmission
$J_1$	${\Gamma }_1$	$(1,0,0)$	0.1	${\mu }_1=10,{\lambda }_1=10$	———	$⟦𝒯u⟧=0$
$J_2$	${\Gamma }_2$	$(-1,0,0)$	0.1	———	${𝒯}_{{𝐧}_2}^{+}u=-{\sigma }_2$	———
$J_2$	${\Gamma }_3$	$(0,0,0)$	2	${\mu }_0=1,{\lambda }_0=1$	${𝒯}_{{𝐧}_3}^{+}u={\sigma }_3$	———

Table 4. Table 4: Number of spheres w.r.t the radius R 𝑅 R

Radius of the big sphere	1	2	3	4	5
Number of total spheres	2	28	94	252	486

Equations

$${\bf L}u:=-\mathop{\rm div}\Big{(}2\mu e(u)+\lambda\mathop{\rm Tr}\big{(}e(u)\big{)}\mathop{\rm Id}\Big{)},$$

$$\lambda =K-\frac{2}{3}\mu ,$$

$$\nu =\frac{\lambda }{2(\mu +\lambda )},$$

$$H_{\mathrm{loc}}^s(\Omega )$$

$$H_{\mathrm{comp}}^s(\Omega )=\bigcup _K\left\{u\in H_{\mathrm{loc}}^s(\Omega )|\mathrm{supp}(u)\subset K\right\},$$

$$V_0({\Omega }^{-})=\left\{v\in H^1({\Omega }^{-})|{\int }_{{\Omega }^{-}}v=0\right\}.$$

$$\Vert {\lambda }_k{\Vert }_{H^{\frac{1}{2}}(\Gamma )}^2:=\Vert {\lambda }_k{\Vert }_{L^2(\Gamma )}^2+{\int }_{\Gamma }{\int }_{\Gamma }\frac{|{\lambda }_k(x)-{\lambda }_k(y){|}^2}{|x-y{|}^3}dxdy.$$

$${\gamma }^{\mp }:H_{\mathrm{loc}}^1({\Omega }^{\mp }{)}^3\to H^{\frac{1}{2}}(\Gamma {)}^3$$

$$[[\phi ]]={\gamma }^{-}\phi -{\gamma }^{+}\phi .$$

$$\mathcal{T}\phi :=2\mu e(\phi )+\lambda \mathrm{Tr}e(\phi )\mathrm{Id},\phi \in H_{\mathrm{loc}}^1({\Omega }^{\mp }{)}^3.$$

$${\mathcal{T}}_{\mathbf{n}}^{\mp }\phi :={\gamma }^{\mp }(\mathcal{T}\phi \mathbf{n}),$$

$$[[\mathcal{T}\phi ]]={\mathcal{T}}_{\mathbf{n}}^{-}\phi -{\mathcal{T}}_{\mathbf{n}}^{+}\phi .$$

$$\mathbf{L}{G}_i(x)=\delta (x){\mathbf{e}}_i,$$

$$G_{ij}(x):=\frac{1}{8\pi \mu |x|}\left(\frac{\lambda +3\mu }{\lambda +2\mu }{\delta }_{ij}+\frac{\lambda +\mu }{\lambda +2\mu }\frac{x_i{x}_j}{|x{|}^2}\right),$$

$$\mathbf{L}u=0,\text{mbox}in{H}^{-1}(\Omega {)}^3$$

$$(\mathcal{S}\varphi )(x):={\int }_{\Gamma }G(x-y)\varphi (y)dy,x\in {\mathbb{R}}^3\text{textbackslash}\Gamma .$$

$$(\mathcal{V}\varphi )(x):={\int }_{\Gamma }G(x-y)\varphi (y)dy,x\in \Gamma ,$$

$$\mathbf{L}\mathcal{S}\varphi =0\text{mbox}in{\mathbb{R}}^3\text{textbackslash}\Gamma .$$

$$\mathcal{D}\phi (x)={\int }_{\Gamma }{\mathcal{T}}_{\mathbf{n},y}(G)(x-y)\phi (y)dy,x\in {\mathbb{R}}^3\backslash \Gamma ,$$

$$(\mathcal{K}\phi )(x)={\int }_{\Gamma }{\mathcal{T}}_{\mathbf{n},y}(G)(x-y)\phi (y)dy,x\in \Gamma ,$$

$$({\mathcal{K}}^{\ast }\varphi )(x)={\int }_{\Gamma }{\left({\mathcal{T}}_{\mathbf{n},x}(G)\right)}^{\top }(x-y)\varphi (y)dy,x\in \Gamma ,$$

$$\mathbf{L}\mathcal{D}\phi =0\text{mbox}in{\mathbb{R}}^3\text{textbackslash}\Gamma ,$$

$$\mathcal{N}\psi (x)={\int }_{{\mathbb{R}}^3}G(x-y)\psi (y)dy,x\in {\mathbb{R}}^3,$$

$$\psi =\mathbf{L}\mathcal{N}\psi =\mathcal{N}\mathbf{L}\psi ,\text{mbox}in{\mathcal{D}}^{\prime }({\mathbb{R}}^3{)}^3.$$

$$\mathcal{S}=\mathcal{N}{\gamma }^{\ast },\mathcal{D}=\mathcal{N}{\mathcal{T}}_{\mathbf{n}}^{\ast }.$$

$${\mathcal{T}}_{\mathbf{n}}^{-}\mathcal{S}\varphi =\frac{1}{2}\varphi +{\mathcal{K}}^{\ast }\varphi ,{\mathcal{T}}_{\mathbf{n}}^{+}\mathcal{S}\varphi =-\frac{1}{2}\varphi +{\mathcal{K}}^{\ast }\varphi ,\text{mbox}on{H}^{-\frac{1}{2}}(\Gamma {)}^3.$$

$$\begin{array}[]{l}\llbracket\mathcal{S}\phi\rrbracket=0,\quad\llbracket\mathcal{D}\varphi\rrbracket=-\varphi,\quad\mbox{on }H^{\frac{1}{2}}(\Gamma)^{3},\\ \llbracket\mathcal{T}\mathcal{S}\phi\rrbracket=\phi,\quad\llbracket\mathcal{T}\mathcal{D}\varphi\rrbracket=0,\quad\mbox{on }H^{-\frac{1}{2}}(\Gamma)^{3}\end{array}$$

$$⟨\mathcal{V}\varphi ,\varphi {⟩}_{H^{\frac{1}{2}}(\Gamma )\times H^{-\frac{1}{2}}(\Gamma )}>c\Vert \varphi {\Vert }_{H^{-\frac{1}{2}}(\Gamma )}^2,\forall \varphi \in H^{-\frac{1}{2}}(\Gamma {)}^3.$$

$$\nabla f(x)\in {\mathbb{R}}^3,\text{mbox}with(\nabla f{)}_i=\frac{\partial f}{\partial x_i},\nabla F(x)\in {\mathbb{R}}^{3\times 3},\text{mbox}with(\nabla F{)}_{ij}=\frac{\partial F_i}{\partial x_j}.$$

$${\nabla }_{\mathrm{s}}=\hat{\theta }\frac{\partial }{\partial \theta }+\hat{\varphi }\frac{1}{\sin \theta }\frac{\partial }{\partial \varphi }$$

Full text

1 Introduction

We consider three-dimensional boundary value or interface problems of the isotropic elasticity equation related to the following operator:

[TABLE]

where the strain tensor reads $e(u)={\frac{1}{2}}(\nabla u+\nabla u^{\top})$. It is obvious to see that the operator ${\bf L}$ is self-adjoint on $L^{2}(\mathbb{R}^{3})^{3}$.

In the definition of the operator (1.1), $\mu,\lambda\in\mathbb{R},\mu>0,2\mu+3\lambda>0$ are the so-called (constant) Lamé parameters. The parameter $\mu$ denotes the shear modulus which describes the tendency of the object to deform at a constant volume when being imposed with opposing forces. The other Lamé parameter $\lambda$ has no physical meanings but is introduced to simplify the definition of the operator (1.1). Indeed, it is related to the bulk modulus $K$ through the relation

[TABLE]

where the bulk modulus $K$ represents the object’s tendency to deform in all directions when acted on by opposing force from all directions. We refer to [12] for more detailed descriptions of the Lamé parameters. It is sometimes useful to introduce Poisson’s ratio $\nu$ which is defined by

[TABLE]

and whose admissible range is $(-1,1/2)$. The material is extremely compressible in the limit $\nu\to-1$ while extremely incompressible in the other limit $\nu\to 1/2$ [18].

A model of linear elasticity with appropriate boundary conditions can be approximated by the classic finite element method, see for example [15, 19] just to name a few contributions from an abundant body of literature, for the general case with non-homogeneous source term. On the other hand, displacement fields $u$ being homogeneous solutions, i.e., ${\bf L}u=0$ within a given domain, can also be represented by isotropic elastic potentials [3, 13] and elasticity in piecewise constant isotropic media can then be treated as integral equations for specified interface conditions. At the origin of the integral formulation lies the definitions of layer potentials and their corresponding integral operators [3, 20] based on the Green’s function [1] in the context of the isotropic linear elasticity.

In particular, on a unit sphere, one can introduce the vector spherical harmonics forming an orthonormal basis of $[L^{2}(\mathbb{S}^{2})]^{3}$ and which are eigenfunctions of the corresponding double and single layer boundary operators based on the Green’s function [1] of isotropic linear elasticity. The vector spherical harmonics were introduced in [9, 10] as an extension of the scalar spherical harmonics [16, 23] to the vectorial case. They were further used in the discretization of different physical models such as the Navier-Stokes equations [8] or Maxwell’s equations [6, 2]. However, they are not widely used and only sparely reported in literature, in particular in the context of isotropic elasticity. We demonstrate in this article that the corresponding integral operators have interesting spectral properties which can be made explicit by employing the vector spherical harmonics.

Our main motivation for this work is the derivation of an integral equation to model elastic materials represented by piecewise constant Lamé constants with spherical inclusions following similar principles that were presented in [14, 4, 5] in the case of scalar diffusion. The particular choice of the vector spherical harmonics as basis functions for a Galerkin discretization thereof leads then to an efficient and stable numerical scheme by exploiting the spectral properties of the involved integral operators. A similar physical model was introduced in  [22] with an algebraic formula of the approximate solution. However, with the spectral properties of the layer potentials and integral operators at hand, our approach first introduces an integral formulation for the exact solution and thus a rigorous mathematical framework. In a second step, we then propose the Galerkin discretization. The mathematical framework lays out the basis to derive a rigorous error analysis which we plan in the future.

We summarize the main contributions and organization of this work as follows:

• •

In Section 2 and 3, we give an introduction and overview of the layer potentials and corresponding boundary integral operators of the isotropic linear elasticity operator (1.1) on an arbitrary bounded domain with Lipschitz boundary which are sparely reported in the literature.

• •

Analytical properties of layer potentials and boundary integral operators are presented and proven in Section 3.4.

• •

On the unit sphere, we introduce the vector spherical harmonics in Section 4 and prove spectral properties of the boundary operators and layer potentials of this particular basis.

• •

As an application, we consider a piecewise constant elastic model with spherical inclusions and derive a integral equation in Section 6 that is then discretized by means of the vector spherical harmonics and tested numerically in Section 7.

2 Preliminaries

Denote $\mathbb{S}^{2}$ the unit sphere and $B$ the unit ball in $\mathbb{R}^{3}$. Let throughout this paper $\Omega^{-}\subset\mathbb{R}^{3}$ denote a bounded domain with Lipschitz boundary $\Gamma=\partial\Omega^{-}$ and outward pointing normal vector field ${\bf n}:\Gamma\to\mathbb{S}^{2}$. Further, we denote by $\Omega^{+}$ the unbounded set $\mathbb{R}^{3}\backslash\overline{\Omega^{-}}$.

2.1 Notations

We will first introduce some standard notions in the context of integral equations which can be found in standard textbooks (see, for example, [17, 20, 21]).

Let $\Omega$ be a domain with Lipschitz boundary, e.g., $\Omega=\Omega^{-}$ or $\Omega=\Omega^{+}$ (unbounded). Following the conventions and notation of [20], we define for $s\in\mathbb{R}$

[TABLE]

see Definition 2.6.1 in [20], and note that this consist of a slightly unconventional definition of $H_{\mathop{\rm loc}}^{\ell}(\Omega)$, see also Remark 2.6.2. We further define, see Definition 2.6.5 in [20], for $s\in\mathbb{R}$

[TABLE]

where the union is taken over all relatively compact subsets $K\subset\Omega$, and introduce

[TABLE]

Next, we denote by $H^{\frac{1}{2}}(\Gamma)^{3}$ the Sobolev space of order $\frac{1}{2}$ with the usual Sobolev-Slobodeckij norm $\|\lambda\|^{2}_{H^{\frac{1}{2}}(\Gamma)}:=\sum_{k=1}^{3}\|\lambda_{k}\|^{2}_{H^{\frac{1}{2}}(\Gamma)}$ for $\lambda=(\lambda_{1},\lambda_{2},\lambda_{3})^{\top}$ and with

[TABLE]

Moreover, we define $H^{-\frac{1}{2}}(\Gamma)^{3}:=\left(H^{\frac{1}{2}}(\Gamma)^{3}\right)^{*}$ and we equip this Sobolev space with the canonical dual norm $\|\cdot\|_{H^{-\frac{1}{2}}(\Gamma)}$. We introduce

[TABLE]

as the continuous, linear and surjective interior and exterior Dirichlet trace operators respectively, see Theorem 2.6.8 [20], and define the jump operator by

[TABLE]

Further, let $\gamma\colon H_{\mathop{\rm loc}}^{1}(\mathbb{R}^{3})^{3}\rightarrow H^{\frac{1}{2}}(\Gamma)^{3}$ be given by $\gamma\varphi=\gamma^{-}\varphi=\gamma^{+}\varphi$ almost everywhere.

Consider now the stress tensor $\mathcal{T}$ associate with ${\bf L}$, as is defined by (1.1), reading

[TABLE]

For the domains $\Omega^{\mp}$, the classical normal derivative operator, satisfying

[TABLE]

for regular $\varphi$, can be extended to an operator $\mathcal{T}^{\mp}_{\bf n}:H^{1}_{\bf L}(\Omega^{\mp})^{3}\to H^{-\frac{1}{2}}(\Gamma)^{3}$, with $H^{1}_{\bf L}(\Omega)^{3}=\{u\in H^{1}(\Omega)^{3}\;|\;\mathbf{L}u\in L^{2}_{\mathop{\rm loc}}(\Omega)^{3}\}$, based on Green’s first identity. We then define the corresponding jump operator by

[TABLE]

Further, define $\mathcal{T}_{\bf n}:H^{1}_{\bf L}(\mathbb{R}^{3})^{3}\to H^{-\frac{1}{2}}(\Gamma)^{3}$ the global normal derivative operator given by $\mathcal{T}_{\bf n}\varphi=\mathcal{T}^{-}_{\bf n}\varphi=\mathcal{T}^{+}_{\bf n}\varphi$.

2.2 Fundamental solutions

Consider the matrix-valued fundamental solution $G=(G_{ij})_{ij}$ to the linear isotropic elasticity equation such that $G_{i}$, the $i$-th column of the matrix $G$ satisfies the following identity:

[TABLE]

with ${\bf L}$ defined by (1.1), $\delta$ being the Dirac distribution at the origin and ${\bf e}_{i}$ the canonical basis in $\mathbb{R}^{3}$. The Green’s function $G$ is given by [17, 21]:

[TABLE]

where we recall that $\mu,\lambda$ are the Lamé constants and $\delta_{ij}$ is the Kronecker symbol.

2.3 Rigid displacement

For a given domain $\Omega\subset\mathbb{R}^{3}$, we consider the following problem: find $u\in H^{1}_{\mathop{\rm loc}}(\Omega)^{3}$ such that

[TABLE]

with appropriate boundary conditions. Equation (2.11) holds obviously if $e(u)=0$. Indeed, we call the displacement $u\in H^{1}_{\mathop{\rm loc}}(\Omega)^{3}$ a rigid displacement if $e(u)=0$. It is well-known that the displacement $u$ is a rigid displacement if and only if it has the form $u=Ax+b$ where $A\in\mathbb{R}^{3\times 3}$ is a constant skew matrix and $b\in\mathbb{R}^{3}$ a constant vector (see, for example, [11, 17]).

3 Layer potentials

In this section, we introduce the layer potentials and associated boundary operators which only have been sparsely reported in the literature for the operator ${\bf L}$. We therefore provide a complete overview.

3.1 Single layer potentials

Using the fundamental solution (2.10), we can now define the single layer potential $\mathcal{S}:H^{-\frac{1}{2}}(\Gamma)^{3}\to H^{1}(\mathbb{R}^{3}\textbackslash\Gamma)^{3}$ associated to the isotropic elasticity operator ${\bf L}$:

[TABLE]

Further, see e.g., [3], such function $\mathcal{S}\phi$ defined on $\mathbb{R}^{3}\textbackslash\Gamma$ is continuous across the interface $\Gamma$, i.e. $\llbracket\mathcal{S}\phi\rrbracket=0$ and a single layer boundary operator $\mathcal{V}:H^{-\frac{1}{2}}(\Gamma)^{3}\to H^{\frac{1}{2}}(\Gamma)^{3}$ can be defined by restricting the single layer potential to the boundary $\Gamma$:

[TABLE]

so that $\gamma\mathcal{S}\phi=\mathcal{V}\phi$. The following result is obvious:

Lemma 3.1**.**

For $\phi\in H^{-\frac{1}{2}}(\Gamma)^{3}$ and $\mathcal{S}\phi$ defined by (3.1), let ${\bf L}$ be the isotropic elasticity operator (1.1) and we have

[TABLE]

3.2 Double layer potential

We introduce the double layer potential $\mathcal{D}:H^{\frac{1}{2}}(\Gamma)^{3}\to H^{1}(\mathbb{R}^{3}\backslash\Gamma)^{3}$, by

[TABLE]

where the subscript $y$ means that the normal derivative operator $\mathcal{T}_{\bf n}$, defined in Section 2.1, is taken with respect to the $y$-variable. We define the double layer boundary operator $\mathcal{K}:H^{\frac{1}{2}}(\Gamma)^{3}\to H^{\frac{1}{2}}(\Gamma)^{3}$ by

[TABLE]

in the sense of principal value. Further, the adjoint double layer boundary operator $\mathcal{K}^{*}:H^{-\frac{1}{2}}(\Gamma)^{3}\to H^{-\frac{1}{2}}(\Gamma)^{3}$ is given as

[TABLE]

Similar to Lemma 3.1, the following result is obvious:

Lemma 3.2**.**

For $\varphi\in H^{\frac{1}{2}}(\Gamma)^{3}$ and $\mathcal{D}\varphi$ defined by (3.3), we have

[TABLE]

where ${\bf L}$ is the isotropic elasticity operator (1.1).

3.3 Newton potential

Finally, for sake of completeness, we also give the Newton potential associated to the isotropic elasticity operator (1.1). Define $\mathcal{N}:H_{\rm comp}^{s}(\mathbb{R}^{3})^{3}\to H_{\rm loc}^{s+2}(\mathbb{R}^{3})^{3}$ for $s\in\mathbb{R}$:

[TABLE]

where $G$ is the Green’s function defined by (2.10). Following the definition for the elasticity operator ${\bf L}$ and all $\psi\in\mathcal{D}^{\prime}(\mathbb{R}^{3})^{3}$, we have

[TABLE]

Let $\gamma^{*}:H^{-\frac{1}{2}}(\Gamma)^{3}\to H^{-1}_{\mathop{\rm comp}}(\mathbb{R}^{3})^{3}$, $\mathcal{T_{\bf n}}^{*}:H^{\frac{1}{2}}(\Gamma)^{3}\to H^{-1}_{\mathop{\rm comp}}(\mathbb{R}^{3}\textbackslash\Gamma)^{3}$ be the adjoint of the trace operator $\gamma$ and the adjoint of the normal derivative operator $\mathcal{T_{\bf n}}$ respectively, defined in Section 2.1. We then give an equivalent definition of the single and double layer potential:

[TABLE]

3.4 Properties of layer potentials

We are now listing a selection of known results of layer potentials that will be used in the following. Let us first recall the following theorem given in [3] (see also [21, Section 6.7]):

Theorem 3.3**.**

Let $\phi\in H^{-\frac{1}{2}}(\Gamma)^{3}$ and the single layer potential $\mathcal{S}$, the adjoint double layer boundary operator $\mathcal{K}^{*}$ be defined by (3.1), (3.5) respectively. Then the interior and exterior normal traces of the stress tensor satisfy

[TABLE]

We now show several jump conditions relating to the boundary layer potentials above which can be found, for example, in [17, Theorem 6.10].

Theorem 3.4**.**

Let $\Omega^{-}\in\mathbb{R}^{3}$ be a bounded Lipschitz domain with boundary $\Gamma$. Consider the single and double layer potentials defined by (3.1) and (3.3) respectively. Then it holds

[TABLE]

for all $\phi\in H^{-\frac{1}{2}}(\Gamma)^{3},\varphi\in H^{\frac{1}{2}}(\Gamma)^{3}$.

We now consider the invertibility of the single layer boundary operator (3.2) (see [17, Theorem 10.7] or [21, Theorem 6.36].

Lemma 3.5**.**

Let $\Omega^{-}\subset\mathbb{R}^{3}$ be a bounded domain with Lipschitz boundary $\Gamma$. If $\mu>0$ and $\lambda\geq 0$, the single layer boundary operator $\mathcal{V}$ defined by (3.2) is coercive, i.e.

[TABLE]

Corollary 3.6** (Invertibility of the single layer boundary operator).**

Let $\Omega^{-}\subset\mathbb{R}^{3}$ be a bounded domain with Lipschitz boundary $\Gamma$ and $\mu>0$ and $\lambda\geq 0$. Then, the single layer boundary operator $\mathcal{V}:H^{-\frac{1}{2}}(\Gamma)^{3}\to H^{\frac{1}{2}}(\Gamma)^{3}$ is invertible.

4 Real vector spherical harmonics

4.1 Surface gradient

In the following, we introduce the real vector spherical harmonics. We begin with some conventions of the gradient. On a given domain $\Omega$, consider a scalar valued function $f:\Omega\to\mathbb{R}^{3}$ and a column-vector valued function $F:\Omega\to\mathbb{R}^{3}$, we define their gradients by

[TABLE]

Note in particular that $\nabla f$ is a column-vector while $\nabla F$ are row-wise gradients for each component $F_{i}$.

Restricting the considerations to the unit ball $\Omega=B$ and it surface $\partial\Omega=\mathbb{S}^{2}$, we denote by

[TABLE]

the surface gradient operator and $\hat{r},\hat{\theta},\hat{\phi}$ are radial, polar and azimuthal unit vectors which are supposed to be row vectors. Let $f$ denote a scalar function and $F$ a vector-valued function, i.e. $f\in\mathbb{R}$ and $F\in\mathbb{R}^{3}$ and with the convention of the gradient field, the surface gradient (4.2) can alternatively be written as

[TABLE]

where $\nabla$ is the gradient in $\mathbb{R}^{3}$ based on the convention (4.1). It is immediate to verify that

[TABLE]

4.2 Definition of vector spherical harmonics

The construction of the vector spherical harmonics is based on the scalar real spherical harmonics defined on the unit sphere $\mathbb{S}^{2}$ denoted by $(Y_{\ell m})_{l\geq 0}^{|m|\leq\ell}$ which are normalized such that

[TABLE]

The vector spherical harmonics $V_{\ell m},W_{\ell m},X_{\ell m}:\mathbb{S}^{2}\to\mathbb{R}^{3}$ of degree $\ell\geq 0$ and order $m$, $|m|\leq\ell$ are given by

[TABLE]

The symbol $\times$ represents the cross product in $\mathbb{R}^{3}$. We refer to Appendix A for some explicit expressions of the vector spherical harmonics for the first few degrees. The vector spherical harmonics satisfy the following orthogonal properties:

[TABLE]

The scalar spherical harmonics (and thus the vector spherical harmonics) can be extended to any sphere $\Gamma_{r}(x_{0})=\partial B_{r}(x_{0})$ by translation and scaling. We will introduce the following scaled scalar product on $\Gamma_{r}(x_{0})$ given by

[TABLE]

In practice, the exact value of the scalar product (4.6) cannot be computed explicitly in general. With a set $\{s_{t},w_{t}\}_{t=1}^{T_{g}}$ of integration points and weights on the unit sphere, the scalar product is approximated by the quadrature rule

[TABLE]

In the numerical tests below in Section 7, we will use the Lebedev quadrature points  [7], which have the property that scalar spherical harmonics up to a certain degree $N_{g}$ are integrated exactly. This relationship is displayed in Table 1. It can be noticed that the number of points increases quadratically with $N_{g}$.

Further, the family of vector spherical harmonics gives a complete basis of $L^{2}(\Gamma_{r}(x_{0}))^{3}$ and any real function $f\in L^{2}(\Gamma_{r}(x_{0}))^{3}$ can be represented as

[TABLE]

where $[v]_{\ell m},[w]_{\ell m},[x]_{\ell m}\in\mathbb{R}$.

4.3 Properties of the derivatives

We give some derivative properties of the surface gradient (4.3) that shall be useful in the upcoming analysis. In the following, let $u$ be the a scalar-valued function, $F,G$ be vector-valued functions and ${\bf A}$ be a matrix-valued function. We have the following product rule:

[TABLE]

We also have the property for the cross product:

[TABLE]

Later proof also requires the triple product

[TABLE]

Now let $h=h(r)$ be a scalar function which does not depend on the polar angles and $u=u(\theta,\phi),H=H(\theta,\phi)$ a scalar and a vector valued function respectively depending only on the polar angles. Then there holds

[TABLE]

For the scalar function $h=h(r)$, denote by $h_{r},h_{rr}$ the first and second derivative. Then we have

[TABLE]

and

[TABLE]

Equations (4.13), (4.14) are given in [9]. Finally, the following identities hold, resulting directly from the definition of the vector spherical harmonics:

[TABLE]

5 Spectral properties of the layer potentials

We will give the main results in Section 5.1, prepare some preliminary results in Section 5.2 and finally provide the proofs in Section 5.3.

5.1 Main results

Consider the single layer potential $\mathcal{S}$ and the single layer boundary operator $\mathcal{V}$ defined by (3.1) and (3.2), we have the following result.

Theorem 5.1**.**

Let $\underline{Y_{\ell m}}$ be the matrix such that

[TABLE]

Then we have:

On the unit sphere $\mathbb{S}^{2}$,

[TABLE]

where $\mathcal{V}$ is the single layer boundary operator defined by (3.2) and $A_{\mathcal{V},\ell}$ is a constant matrix given by

[TABLE] 2. 2.

When $|x|<1$, we have

[TABLE]

where $\mathcal{S}$ is the single layer potential given by (3.1) and the matrix $A_{\mathcal{S},\ell}^{in}(x)$ has the form

[TABLE] 3. 3.

When $|x|>1$, we have

[TABLE]

where the matrix $A_{\mathcal{S},\ell}^{out}(x)$ is given by

[TABLE]

The following result is a corollary of Theorem 5.1.

Corollary 5.2**.**

Let $\underline{Y_{\ell m}}$ be the a matrix defined by (5.1). Then, on the unit sphere $\mathbb{S}^{2}$, there holds

[TABLE]

where $\mathcal{K}$ is the double layer boundary operator (3.5) with its adjoint $\mathcal{K}^{*}$ and $A_{\mathcal{K}^{*},\ell}$ is a constant and diagonal matrix:

[TABLE]

Remark 5.3**.**

Recall that the Lamé constats $\mu,\lambda$ satisfy $\mu>0,2\mu+3\lambda>0$, we can verify that the eigenvalue $\tau_{\mathcal{K}^{*},\ell}^{k}$ of the adjoint double layer boundary operator is $-{1\over 2}$ if and only if $\ell=1,k=2$. And the eigenvectors associated with the eigenvalue $-{1\over 2}$ are $W_{1m}$, $m=\pm 1,0$.

The following theorem gives explicit expressions of the double layer potential.

Theorem 5.4**.**

Let $\underline{Y_{\ell m}}$ be given by (5.1) and $\mathcal{D}$ the double layer potential on the unit sphere introduced in (3.3), we have:

For $|x|<1$,

[TABLE]

where

[TABLE] 2. 2.

For $|x|>1$,

[TABLE]

where

[TABLE]

The constants $a^{in,\mathcal{D},\ell}_{ij},a^{out,\mathcal{D},\ell}_{ij}$ are listed in Appendix B.

The results of Theorems 5.1, 5.4 and Corollary 5.2 can be extended to any sphere $\Gamma_{r}(x_{0})=\partial B_{r}(x_{0})$, which is, a sphere centered at $x_{0}$ with radius $r_{0}$, by setting \underline{Y_{\ell m}}\Big{(}{x-x_{0}\over|x-x_{0}|}\Big{)} and using the following scaling in $r$:

[TABLE]

5.2 Preliminary lemmas

To prove the results of Section 5.1 in the upcoming Section 5.3, we derive first several preliminary lemma.

Lemma 5.5**.**

For a scalar function $h=h(r)$, we have the following identities

[TABLE]

where ${\bf n}$ is the outward pointing unit normal vector of the unit sphere $\mathbb{S}^{2}$.

Proof.

We consider $h(r)V_{\ell m}$ first. Following (4.9) and (4.12), we have,

[TABLE]

where we also use the fact that ${\bf n}$ is equal to the radial basis $\hat{r}$ on the unit sphere $\mathbb{S}^{2}$. Consider now the first term $h_{r}(1)(V_{\ell m}+\hat{r}V_{\ell m}^{\top})\hat{r}$. To compute $\hat{r}V_{\ell m}^{\top}\hat{r}$, we need first the relation $\nabla_{\!\rm s}Y_{\ell m}^{\top}\hat{r}=0$ following from (4.3). Then according to the definition of the vector spherical harmonics $V_{\ell m}$ given by (4.4), we have

[TABLE]

Therefore, the first term in (5.10) yields

[TABLE]

Now consider the second term $(\nabla_{\!\rm s}V_{\ell m}+\nabla_{\!\rm s}V_{\ell m}^{\top})\hat{r}$. According to (4.3), we have $\nabla_{\!\rm s}V_{\ell m}\hat{r}=0$. The definition of $V_{\ell m}$ (4.4) gives

[TABLE]

We compute the two terms separately. Using the product rule (4.9), we have

[TABLE]

For $\nabla_{\!\rm s}(Y_{\ell m}\hat{r})^{\top}\hat{r}$, we use (4.9) and have

[TABLE]

where we use the fact that $\nabla_{\!\rm s}\hat{r}$ is a symmetric matrix and $\nabla_{\!\rm s}\hat{r}\hat{r}=0$. Therefore, there holds

[TABLE]

The computation of e\big{(}h(r)V_{\ell m}\big{)}{\bf n} is thus completed in view of (5.11) and (5.12). A similar computation gives $h(r)W_{\ell m}$. Consider now $h(r)X_{\ell m}$. Similar to (5.10), we have to compute the sum:

[TABLE]

Notice that $X_{\ell m}=\hat{r}\times\nabla_{\!\rm s}Y_{\ell m}$ is orthogonal to $\hat{r}$ and it follows immediately that $X_{\ell m}^{\top}\hat{r}=0$. Further, by (4.3), there holds $\nabla_{\!\rm s}X_{\ell m}\hat{r}=0$. Then, it remains to consider $X_{\ell m}\hat{r}$ and $\nabla_{\!\rm s}X_{\ell m}^{\top}\hat{r}$. For the term $\nabla_{\!\rm s}X_{\ell m}^{\top}\hat{r}$, we use the relation (4.10) and have

[TABLE]

Both terms can be computed by (4.11):

[TABLE]

and

[TABLE]

This gives

[TABLE]

Then we get the result for $h(r)X_{\ell m}$. ∎

By (4.13), (5.9), we get, for a displacement $h(r)V_{\ell m}+g(r)W_{\ell m}+h(r)X_{\ell m}$, it holds that

[TABLE]

The flowing lemma concerns the double layer boundary operator and its adjoint. In particular, on a sphere, we have the following lemma.

Lemma 5.6**.**

Let $\mathcal{K}$ be the double layer boundary operators defined by (3.4) on a sphere and $\mathcal{K}^{*}$ its adjoint operators. Then for $v\in L^{2}(\mathbb{S}^{2})^{3}$, we have $\mathcal{K}v=\mathcal{K}^{*}v$.

Proof.

Indeed, we have

[TABLE]

and

[TABLE]

Hence, we have

[TABLE]

The same result holds for \Big{(}2\mu e_{y}(G_{i})+\lambda\mathop{\rm Tr}e_{y}(G_{i})\mathop{\rm Id}\Big{)}_{jk}{y_{k}\over|y|} by replacing $x$ by $y$. Further, for $x,y\in\mathbb{S}^{2}$, the following relation holds:

[TABLE]

Therefore, we have

[TABLE]

∎

5.3 Proof of the principal results

We are now ready to prove Theorem 5.1.

Proof of Theorem 5.1.

Consider the single layer potential $\mathcal{S}$ defined by (3.1). Note that for any $\phi\in H^{-\frac{1}{2}}(\mathbb{S}^{2})^{3}$, as announced in Lemma 3.1, $u=\mathcal{S}\phi$ satisfies the following linear isotropic elasticity system

[TABLE]

Now we determine $u=\mathcal{S}V_{\ell m}$ by means of separation of variables in spherical coordinates. That is, we propose the Ansatz displacement field $\mathcal{S}V_{\ell m}$ as a function of the spherical coordinates of form $u=\mathcal{S}V_{\ell m}=h(r)V_{\ell m}+g(r)W_{\ell m}+h(r)X_{\ell m}$ where $f,g,h$ are three scalar functions of $r$ to be determined. Using the relation

[TABLE]

and plugging the Ansatz into (5.14), together with (4.13)-(4.15), we have the following equation:

[TABLE]

Since $V_{\ell m},W_{\ell m},X_{\ell m}$ is an orthogonal basis, all the coefficients of $V_{\ell m},W_{\ell m},X_{\ell m}$ must be zero in (5.15). Let first $\ell\geq 1$ and we have six sets of analytical solutions to (5.15) reading

[TABLE]

in which $(i),(ii),(iii)$ are admissible only for $|r|>0$ while $(iv),(v),(vi)$ are unbounded when $|r|\to\infty$. Now, we consider the exterior and the interior of the unit sphere separately in which we aim to get $\mathcal{S}V_{\ell m}$. For the sake of simplicity, write

[TABLE]

Write $\mathcal{S}V_{\ell m}$ as a linear combination of the three solutions in the exterior and three in the interior of the unit sphere:

[TABLE]

where $a^{in/out}_{\ell m},b^{in/out}_{\ell m},c^{in/out}_{\ell m}\in\mathbb{R}$ are constants to be determined. In order to determine these six unknown constants, we use the jump relation given by Theorems 3.4:

[TABLE]

Since $V_{\ell m},W_{\ell m},X_{\ell m}$ is an orthogonal basis, it follows immediately from the first equality in (5.17) that

[TABLE]

Now take the inner product of the second equation (5.17) with $V_{\ell m},W_{\ell m},X_{\ell m}$. With the computations in (5.13), we have

[TABLE]

Hence, we conclude

[TABLE]

Similar computations following the same logic give the results for $W_{\ell m},X_{\ell m}$ for $\ell\geq 1$.

In the case where $\ell=0$, we only have to treat $V_{00}$ since $W_{00}=X_{00}=0$. Using (5.16)–(5.17), we have

[TABLE]

Notice that the result (5.18) is indeed consistent with the cases where $\ell\geq 1$. We have proved therefore the theorem. ∎

Corollary 5.2 can now be deduced.

Proof of Corollary 5.2.

As is shown in (3.9), we have

[TABLE]

Again, apply the computation in (5.13), we have

[TABLE]

Hence, we have

[TABLE]

Similar computations give the results for the other two components $W_{\ell m},X_{\ell m}$. Further, Lemma 5.6 provides the result for the double layer operator $\mathcal{K}$. ∎

Finally, the proof of Theorem 5.4 follows the same structure to that of Theorem 5.1 in employing the jump relations

[TABLE]

for the choices $\phi=V_{\ell m},W_{\ell m},X_{\ell m}$ respectively.

6 Application

We study here a case of an elasticity problem involving several spherical inclusions as an application of the results in the above sections, derive an integral equation formulation and propose a Galerkin formulation thereof based on the vectorial spherical harmonics.

6.1 Problem setting

Set the sets of indices $J_{1},J_{2},J$ such that $M+1\in J_{2}$, $J_{1}\cap J_{2}=\emptyset$ and

[TABLE]

and let $\Omega_{i}\subset\mathbb{R}^{3}$, $i\in J$ be non-overlapping balls, centered at $x_{i}\in\mathbb{R}^{3}$ with radius $r_{i}$, all contained in an additional ball $B_{R}$ centred at the origin with the radius $R$. Moreover, define the domains

[TABLE]

Denote the boundaries $\Gamma_{i}=\partial\Omega_{i}$, $i\in\{0\}\cup J_{1}\cup J_{2}$. Then, it holds that $\Gamma_{0}=\bigcup\limits_{i\in J_{1}\cup J_{2}}\Gamma_{i}$. Set further ${\bf n}_{0}$ as the outward pointing unit normal vector with respect to the domain $\Omega_{0}$ and ${\bf n}_{i}$ the outward pointing normal vector with respect to each domain $\Omega_{i}$, $i\in J_{1}\cup J_{2}$. Then it holds that ${\bf n}_{i}=-{\bf n}_{0}$. We refer to Figure 1 for an illustration of geometry configuration.

In the numerical example presented below, we assume that each inclusion $\Omega_{i}$, $i\in J_{1}$ is filled with an isotropic elastic medium associated with Lamé parameters $\mu_{i},\lambda_{i}$. The remaining background domain $\Omega_{0}$ is filled with medium of Lamé constants $\mu_{0},\lambda_{0}$. Further, we denote ${\mathcal{T}^{i}_{{\bf n}_{j}}}^{\pm}$ the normal derivative operator acting on the boundary $\Gamma_{j}$ with Lamé constants $\mu_{j},\lambda_{j}$ and the normal vector ${\bf n}_{j}$:

[TABLE]

where $\gamma^{\pm}_{j}$ are the exterior and interior the trace operators on $\Gamma_{j}$, following the notations given by (2.4) by taking $\Omega^{-}=\Omega_{j}$. Define the parameter $s_{i}$ by

[TABLE]

In particular, write

[TABLE]

and it is obvious that

[TABLE]

For given $f_{i}\in H^{-\frac{1}{2}}(\Gamma_{i})^{3}$, we impose the transmission condition:

[TABLE]

Further, we let the domain $\Omega$ be subjected of a given stress tensor on its boundary. Indeed, we let each part of the boundary $\Gamma_{i}$, $i\in J_{2}$ be subjected to a given stress tensor $\sigma_{i}\in H^{1/2}(\Gamma_{i})^{3}$:

[TABLE]

Then, we consider the solution $u\in V_{0}(\Omega)^{3}$ to the following interface problem:

[TABLE]

Standard arguments involving the Lax-Milgram theorem yields the well-posedness of the problem.

Remark 6.1**.**

In our case, we pre-defined the index set $J_{2}$ with the condition $M+1\in J_{2}$. However, by equipping the domain $\Omega_{M+1}$ (which is indeed the exterior domain of the sphere $B_{R}$) with Lamé parameters $\mu_{M+1},\lambda_{M+1}$, one can also relax the setting and consider the case where $M+1\in J_{1}$.

6.2 Integral equation

Let $u\in V_{0}(\Omega)^{3}$ denote the solution to (6.6) satisfying $\llbracket u\rrbracket=0$ on $\Gamma_{i}$ for $i\in J_{1}$ and define $\nu$ on $\Gamma_{0}$ by

[TABLE]

where $\gamma_{0}^{-}$ is the trace operator $\gamma_{0}^{-}\colon H^{1}(\Omega_{0})^{3}\rightarrow H^{\frac{1}{2}}(\Gamma_{0})^{3}$ defined by (2.4) for $\Omega^{-}=\Omega_{0}$. To deduce an integral equation for $\nu$, we first introduce an auxiliary problem: find a solution $v\in H^{1}(\Omega)^{3}$ to

[TABLE]

The auxiliary problem (6.7) admits a unique solution $v$ in $\Omega$ and observe that

[TABLE]

but is different in all $\Omega_{i}$, $i\in J_{1}$. Further, there exists a global density $\phi$ supported on $\Gamma_{0}$ such that

[TABLE]

where $\mathcal{S}_{G}^{0}$ is the global layer potential (3.1) with Lamé parameters $\mu_{0},\lambda_{0}$ defined on the whole boundary $\Gamma_{0}$ while $\mathcal{S}_{i}^{0}$ is the local single layer potential with Lamé parameters $\mu_{0},\lambda_{0}$ defined locally on the sphere $\Gamma_{i}$ and $\mathcal{V}_{G}^{0},\mathcal{V}^{0}_{i}$ their corresponding single layer boundary operators (3.2). Further, according to Theorem 3.4, the density $\phi_{i}$ is given by the jump relation

[TABLE]

The last equivalence in (6.9) is obtained because $u=v$ on $\Omega_{0}$. Further, both solutions $u,v$ can be represented by some local densities in each domain $\Omega_{i}$:

[TABLE]

where $\mathcal{S}_{i}$ is the local layer potential with Lamé parameters $\mu_{i},\lambda_{i}$ while $\mathcal{S}_{i}^{0}$ with $\mu_{0},\lambda_{0}$, both of which are defined locally on the sphere $\Gamma_{i}$. For the corresponding single layer boundary operators $\mathcal{V}_{i}^{0},\mathcal{V}_{i}$ defined by (3.2), we have

[TABLE]

According to Corollary 3.6, the single layer boundary operators $\mathcal{V}_{i}^{0},\mathcal{V}_{i}$ are invertible, so we have

[TABLE]

Now, according to Theorem 3.3, we have

[TABLE]

where ${\mathcal{K}_{i}^{0}}^{*}$, $\mathcal{K}_{i}^{*}$ are adjoint double layer boundary operators (3.5) defined locally on $\Gamma_{i}$ with Lamé constants $\mu_{0},\lambda_{0}$ and $\mu_{i},\lambda_{i}$ respectively. In problem (6.6), we have

[TABLE]

where $f_{i}\in H^{-\frac{1}{2}}(\Gamma_{i})^{3},\sigma_{i}\in H^{\frac{1}{2}}(\Gamma_{i})^{3}$ are given transmission conditions and Neuman boundary conditions respectively. Providing (6.9)–(6.11), we have

[TABLE]

According to (6.8), we have

[TABLE]

where $\mathcal{L}$ is an operator defined on each $\Gamma_{i}$:

[TABLE]

and the vector $\Sigma$ satisfies

[TABLE]

6.3 Galerkin approximation

Introduce $V_{N,i},i\in J_{1}\cup J_{2}$ the set spanned by vector spherical harmonics (4.4) on the sphere $\Gamma_{i}$ with a maximum degree $N$:

[TABLE]

where we write $Y_{\ell m}^{1}=V_{\ell m},Y_{\ell m}^{2}=W_{\ell m},Y_{\ell m}^{3}=X_{\ell m}$ and $Y_{\ell m}^{ki}(x)=Y_{\ell m}^{k}\left({x-x_{i}\over r_{i}}\right).$ Define also the global set

[TABLE]

We look for the approximation of $\nu_{N}\in V_{N}$ to (6.6) with

[TABLE]

In practice, we use the quadrature (4.7) to approximate the inner product and the approximate solution $\nu_{N}$ on each $\Gamma_{i}$ thus satisfies

[TABLE]

Denote by $\mathcal{M}=3(N+1)^{2}(|J_{1}|+|J_{2}|)$ the number of degrees of freedom. The $\mathbb{R}^{\mathcal{M}}$-vector collecting all the coefficients $[y]_{\ell m}^{k}$ denoted by ${\bf\Lambda}$ yields the linear system

[TABLE]

where by (4.5), (4.6), the $\mathcal{M}\times\mathcal{M}$ diagonal matrix ${\bf D}$ is given by

[TABLE]

and where $\bf{N}$ is a $\mathcal{M}\times\mathcal{M}$ matrix with coefficients

[TABLE]

The right hand side ${\bf F}\in\mathbb{R}^{\mathcal{M}}$ is given by its coefficients:

[TABLE]

To derive the entries of the matrix (6.17), recall the spectral results in Lemma 5.1 so that we have on $\Gamma_{j}$

[TABLE]

where $\tau^{k^{\prime}j}_{\mathcal{V},\ell^{\prime}},\tau^{k^{\prime}j}_{\mathcal{K}^{*},\ell^{\prime}}$ are eigenvalues of the single and double layer boundary operator with Lamé constants $\mu_{j},\lambda_{j}$ given by (5.2), (5.5) respectively. Hence, we have

[TABLE]

where the constant $C_{j\ell^{\prime}k^{\prime}}$ reads

[TABLE]

Recall that $s_{i}$ denotes the parameter defined by (6.3). Further, by (5.3), (5.4),

[TABLE]

where $[\cdot]_{k^{\prime}}$ denotes the $k^{\prime}$-th column of the obtained matrix. Hence, the coefficient $[N_{ij}]_{\ell m,\ell^{\prime}m^{\prime}}^{kk^{\prime}}$ of the matrix ${\bf N}$ (6.17) reads

[TABLE]

where $y_{ij}^{t}=x_{i}+r_{i}s_{t}-x_{j}$ and $f(j)$ takes the value

[TABLE]

In particular, when $i=j$, we use (4.5) for the exact value of the inner product and obtain:

[TABLE]

Finally, the right-hand side vector ${\bf F}$ with coefficient $[F_{i}]_{\ell m}^{k}$ is given by:

[TABLE]

where $[\Sigma_{i}]_{\ell m}^{k}\in\mathbb{R}$ is determined by the right-hand side vector $\Sigma$ in (6.12) :

[TABLE]

Remark 6.2**.**

A similar physical model called “Finite Cluster Model” was considered in [22] in where an algebraic formulation is derived through the use of M2L-operators (using the fast multipole method terminology). However, with the jump relations given in (3.10), the algebraic formulation of the “Finite Cluster Model” can be proven to be equivalent to the discrete integral formulation (6.15) presented above.

It shall be noted that the mathematical framework introduced here through the use of layer potentials and boundary operators in order to derive an integral equation (6.12) defining the exact solution and and the subsequent introduction of the Galerkin discretization (6.15) allows a mathematical analysis which will be subject of an upcoming work.

7 Numerical tests

For all following computations, we chose the number of Lebedev integration points $T_{g}$ such that, for given $N$, products of two scalar spherical harmonics of maximal degree $N$, thus spherical harmonics of degree $2N$, are integrated exactly. The number of points can then be extracted from Table 1.

7.1 One sphere model

We start with a simple model involving only one single sphere whose solution can be computed analytically in order to assess the convergence of the method in this simple setting. For simplicity, let $\mathbb{S}^{2}\subset\mathbb{R}^{3}$ be the unit sphere on which a stress tensor $\sigma\in H^{\frac{1}{2}}(\mathbb{S}^{2})^{3}$ is imposed and let the Lamé constants be $\mu_{0}=\lambda_{0}=1$. Let ${\bf n}$ be the outward pointing normal vector with respect to the unit sphere $\mathbb{S}^{2}$. The solution $u\in V_{0}(B)^{3}$ to the problem

[TABLE]

reads

[TABLE]

with $\mathcal{S}$ being the single layer potential (3.1) and $\mathcal{K}^{*}$ the adjoint of the double layer boundary operator (3.5). For the given tensor $\sigma$, if there exists an integer $\ell_{\mathop{\rm ex}}$ such that we can expand $\sigma$ by means of vector spherical harmonics up to order $\ell_{\mathop{\rm ex}}$:

[TABLE]

then the exact solution restricted to the sphere $\Lambda_{\mathop{\rm ex}}=u|_{\mathbb{S}^{2}}$ is given explicitly by

[TABLE]

where $\tau_{\mathcal{V},\ell}^{k},\tau_{\mathcal{K}^{*},\ell}^{k}$ are the eigenvalues of the single layer boundary operator and the adjoint double layer boundary operator given by (5.2) and (5.5) resp. They only concern in computing (7.2) is that the denominator tends to zero if $\tau_{\mathcal{K}^{*},\ell}^{k}$ approaches $-1/2$. Recall that according to Remark 5.3, the only possible eigenvectors of $\mathcal{K}^{*}$ with the eigenvalue $-1/2$ are $W_{1,-1},W_{1,0},W_{1,1}$. According to Appendix A, we see that they are constant and parallel to the cartesian basis ${\bf e}_{i},i=1,2,3$. To ensure that (7.2) is well defined and these modes avoided, we simply impose that

[TABLE]

We consider the following four cases:

• •

Case 1. $\sigma=-(x,y,z)^{\top}$.

• •

Case 2. $\sigma=-(x,0,0)^{\top}$.

• •

Case 3. $\sigma=-(x^{7},y^{7},z^{7})^{\top}$.

• •

Case 4. \sigma=-\big{(}\sin(2\pi x),\sin(2\pi y),\sin(2\pi z)\big{)}^{\top}.

Table 2 lists the $L^{2}$ norm of the numerical solution on the unit sphere $||\Lambda_{s}||_{L^{2}}$ in each case with different degrees of vector spherical harmonics and the relative error is defined by

[TABLE]

The exact solutions $\Lambda_{{\rm ex}}$ in the first three cases are exactly computed by the (7.1), (7.2) while the in the last cases, the “exact” solution is obtained by taking a large enough $\ell_{\mathop{\rm ex}}$ (in this case $\ell_{\mathop{\rm ex}}$ = 50).

7.2 Convergence with respect to the degree $N$

We study now the convergence of the error measured in the $L^{2}$ norm with respect to the degree $N$ of the vector spherical harmonics. Using the notation introduced in Section 6, we test a model with $M=2$ and

[TABLE]

Let $\Gamma_{1}$ be the sphere centered at $(1,0,0)$ respectively with radius $0.1$ and $\Gamma_{2}$ centered at $(-1,0,0)$ with radius $0.1$ while $\Gamma_{3}$ is centered at $(0,0,0)$ with radius $2$. The inclusion $\Omega_{1}$ is filled with a medium represented by the Lamé constants $\mu_{1}=10,\lambda_{1}=10$ while the background domain $\Omega_{0}$ uses $\mu_{0}=1,\lambda_{0}=1$ as Lamé parameters.

The interface condition $\llbracket\mathcal{T}u\rrbracket=0$ is imposed on $\Gamma_{1}$ while the spheres $\Gamma_{2},\Gamma_{3}$ are subjected to a stress tensor $\mathcal{T}_{{\bf n}_{2}}^{+}u=-\sigma_{2}$ respectively $\mathcal{T}_{{\bf n}_{3}}^{+}u=\sigma_{3}$. Table 3 illustrates the parameters of the above geometry configuration.

We now test two cases to see the relation between the relative error and the degree of the spherical harmonics for different kinds of imposed stress tensors $\sigma_{2},\sigma_{3}$:

The two stress tensors are set to be smooth functions such that

[TABLE] 2. 2.

The stress tensors are set to be piecewise smooth such that

[TABLE]

and

[TABLE]

We compute the “exact” solution to the problem with a large degree of vector spherical harmonics ($N_{\rm ex}=50$) for both cases. In the Figure 2, we illustrate the log of the relative error (7.3) with respect to the degree $N$ of spherical harmonics of the two tests above. We observe exponential convergence in the first case where the given stress tensor is regular. In the second case, the situation is less clear as an initial pre-asymptotic is followed by a very fast convergence and the asymptotic regime is not yet reached, but the absolute error is already very small.

7.3 Computational cost

Next, we study the computational cost of our numerical method by considering an “embedded model” with $M$ inclusions by increasing the value of $M$. We do the following test with Matlab on an iMac with a 2,7 GHz Intel Core i5 processor.

We consider a case where a stress tensor ${\mathcal{T}^{0}{{\bf n}_{M+1}}}^{+}u=-{1\over R}(x,y,z)^{\top}$ is imposed on a origin-centered sphere with radius $R$, denoted by $\mathbb{S}^{2}_{R}$. Inclusions are taken to be all the spheres with radii $0.1$, centered on a cubic lattice $\mathbb{Z}^{3}$ and which are contained in $\mathbb{S}^{2}_{R}$. We increase the number of inclusions $M$ by increasing the value of the radius $R$ of the big sphere. Table 4 lists the number of spheres with respect to the radius $R$ that grows of course cubically.

We fill each small inclusion with a medium associated with the Lamé constants $\mu_{i}=10,\lambda_{i}=10$, $i=1,...,M$ and take the transmission condition $\llbracket\mathcal{T}u\rrbracket=0$ on each embedded sphere. Further, the Lamé constants of the background domain are fixed to be $\mu_{0}=1,\lambda_{0}=1$. The degree of the vector spherical harmonics is chosen to be $N=3$. Further, we stop the iterative solver of the linear system when the residual is smaller than $10^{-6}$. Figure 3 illustrates the computed elastic deformation of the model computed when $R=3$. The colorcode represents the modulus of the displacement.

We report the result of the computational time in Figure 4 which illustrates that the computational cost with respect to the number of spheres grows as $O(M^{2})$. This is the normal scaling for an integral equation involving $M$ spheres, whose iterative solver requires a number of iterations that is independent of $M$ which we observe.

7.4 The effect of an inclusion

We now consider the unit ball $B_{1}$ which contains an additional inclusion $\Omega_{1}$ in form of a sphere centered at the origin with radius $0.5$. We study how the displacement on the unit sphere $\mathbb{S}^{2}$ is influenced by the compressibility of the small inclusions $\Omega_{1}$. We will use the Poisson’s ratio as the parameter defined by (1.2) describing the compressibility of a substance.

Let the stress tensor $-(x,y,z)^{\top}$ be imposed on the unit sphere $\mathbb{S}^{2}$ and fix the shear modulus of the exterior shell $\Omega_{0}=B_{1}\backslash\overline{\Omega}_{1}$ to be $\mu_{0}=1$ and the shear modulus $\mu_{1}=1$ for the inclusion $\Omega_{1}$. We test several cases where the the exterior shell and the inclusion are associated to different Poisson’s ratio $\nu_{0},\nu_{1}$. Recall that $\nu_{0},\nu_{1}\in(-1,1/2)$ according to the definition, we have the limit values of the first Lamé parameter $\lambda_{1}$:

[TABLE]

In Figure 5, we plot the $L^{2}$ norm of the displacement on the unit sphere by letting the Poisson’s ratio $\nu_{1}$ vary in $[-1,0.4998]$ with different given values of Poisson’s ratio $\nu_{0}$ of the background domain $\Omega_{0}$. In Figure 6, we give two solutions with different Poisson’s ratios: the left solution is obtained by setting $\nu_{1}=0.4995$ while the other is obtained by setting $\nu=-1$, both embedded into a background domain with $\nu_{0}=1/6$.

8 Conclusion

In this article, we have discussed the layer potentials and their corresponding integral operators on arbitrary bounded domains with Lipschitz boundary in the context of isotropic elasticity. We proved jump relations of layer potentials and the invertibility of the single layer boundary operator. In the particular case where the body is a unit ball, we present spectral properties of the boundary operators on the base of the vector spherical harmonics. We then derived a second-kind integral equation for isotropic elastic materials with spherical inclusions that was then discretized by employing the vector spherical harmonics as basis functions and exploiting the spectral properties to enhance efficiency of the discretization. In the last part, we effect some numerical tests to asses the properties of the method: the accuracy with respect to the degree of the vector spherical harmonics and the complexity of the computational cost with respect to the number of spherical inclusions. We also used the method to explore how the deformation of the elastic material is effected by the value of the Poisson’s ratio.

9 Aknowledgement

Benjamin Stamm acknowledges the funding from the German Academic Exchange Service (DAAD) from funds of the â Bundesministeriums für Bildung und Forschungâ (BMBF) for the project Aa-Par-T (Project-ID 57317909). Shuyang Xiang acknowledges the funding from the PICS-CNRS as well as the PHC PROCOPE 2017 (Project N37855ZK).

Appendix A: computation of the first few vector spherical harmonics

We first start considering the table of vector spherical harmonics up to the second order as listed below:

[TABLE]

This gives first the obvious result that

[TABLE]

Using the definition of the surface gradient (4.3), we obtain

[TABLE]

The spherical harmonics $V_{\ell m}$ up to order 2 are then given as follows

[TABLE]

The spherical harmonics $W_{\ell m}$ up to order 2 are given by

[TABLE]

And finally, the spherical harmonics $X_{\ell m}$ up to order 2 are given by

[TABLE]

Appendix B: Entries of matrices $A_{\mathcal{D},\ell}^{in}$ and $A_{\mathcal{D},\ell}^{out}$

The coefficients in $A_{\mathcal{D},\ell}^{in}$ and $A_{\mathcal{D},\ell}^{out}$ are given as follows:

[TABLE]

and

[TABLE]