A Dirichlet-to-Neumann finite element method for axisymmetric elastostatics in a semi-infinite domain

Accepted Manuscript A Dirichlet-to-Neumann finite element method for axisymmetric elastostatics in a semi-infinite domain

Eduardo Godoy, Valeria Boccardo, Mario Durán

PII: DOI: Reference:

S0021-9991(16)30493-4 http://dx.doi.org/10.1016/j.jcp.2016.09.066 YJCPH 6877

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

11 February 2016 14 September 2016 30 September 2016

Please cite this article in press as: E. Godoy et al., A Dirichlet-to-Neumann finite element method for axisymmetric elastostatics in a semi-infinite domain, J. Comput. Phys. (2016), http://dx.doi.org/10.1016/j.jcp.2016.09.066

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights • • • • •

A numerical method for elasticity in axisymmetric semi-infinite domains is proposed The method couples finite elements with Dirichlet-to-Neumann boundary conditions The lack of an explicit closed-form expression for the DtN map needs to be overcome This is done by using a procedure that combines analytical and numerical techniques The numerical experiments confirm the effectiveness and accuracy of the method

A Dirichlet-to-Neumann finite element method for axisymmetric elastostatics in a semi-infinite domain Eduardo Godoya,∗, Valeria Boccardob,c , Mario Dur´ ana a

b

INGMAT R&D Centre, Jos´e Miguel de la Barra 412, 4to piso, Santiago, Chile. Facultad de Ingenier´ıa, Pontificia Universidad Cat´ olica de Chile, Av. Vicu˜ na Mackenna 4860, Macul, Santiago, Chile. c Facultad de Ciencias, Universidad Mayor, Av. Manuel Montt 367, Providencia, Santiago, Chile.

Abstract The Dirichlet-to-Neumann finite element method (DtN FEM) has proven to be a powerful numerical approach to solve boundary-value problems formulated in exterior domains. However, its application to elastic semi-infinite domains, which frequently arise in geophysical applications, has been rather limited, mainly due to the lack of explicit closed-form expressions for the DtN map. In this paper, we present a DtN FEM procedure for boundary-value problems of elastostatics in semi-infinite domains with axisymmetry about the vertical axis. A semi-spherical artificial boundary is used to truncate the semi-infinite domain and to obtain a bounded computational domain, where a FEM scheme is employed. By using a semi-analytical procedure of solution in the unbounded residual domain lying outside the artificial boundary, the exact nonlocal boundary conditions provided by the DtN map are numerically approximated and efficiently coupled with the FEM scheme. Numerical results are provided to demonstrate the effectiveness and accuracy of the proposed method. Keywords: Dirichlet-to-Neumann map, Semi-infinite domain, Finite elements, Elasticity

1. Introduction Boundary-value problems (BVPs) formulated in unbounded domains arise in many fields of application, such as acoustics, geophysics, seismology, oceanography, and meteorology, among others. When solving numerically such problems, the unboundedness of the domain represents a major difficulty, since most standard numerical methods require a bounded domain. Different mathematical and numerical approaches have been devised to solve BVPs in unbounded domains. According to Givoli [1, 2], they are classified into four main categories: boundary integral methods, infinite element methods, absorbing layer methods, and artificial boundary condition (ABC) methods. The advantages and limitations of each one of them are discussed in [1]. The present work is concerned with the latter category, also referred to as artificial boundary method [3, 4]. In the standard ABC method, the original unbounded domain is truncated by introducing an artificial boundary enclosing the particular area of interest, thus defining a bounded computational domain, where standard numerical methods may be used to solve the BVP. This is possible on the condition that suitable boundary conditions are set on the artificial boundary (the ABCs), which are supposed to properly represent the unbounded residual domain that was eliminated. ∗

Corresponding Author Email addresses: [email protected] (Eduardo Godoy), [email protected] (Valeria Boccardo), [email protected] (Mario Dur´ an) Preprint submitted to Journal of Computational Physics

October 3, 2016

In general, many different choices of ABCs are possible. We focus herein on ABCs based on the Dirichlet-to-Neumann (DtN) map and their use in combination with finite elements, which results in the Dirichlet-to-Neumann finite element method (DtN FEM) [2, 5]. The main advantage of this approach is that the DtN map provides exact ABCs, in such a way that the resulting BVP in the computational domain is mathematically equivalent to the original unbounded BVP. Therefore, the use of the FEM to solve the former results in highly accurate and robust numerical schemes. The DtN FEM has been mainly developed for BVPs formulated in exterior domains, i.e. the complement of a compact set. It is customary to assume a circular or spherical artificial boundary, in order to be able to apply the method of separation of variables to solve the resulting BVP in the residual domain. Thus, an analytical solution in series form is obtained, with coefficients that are computed exactly. From it, an explicit closed-form expression for the DtN map is deduced, given as an exact nonlocal relation between the solution and its normal derivative on the artificial boundary, which is used as a boundary condition that completes the mathematical formulation of the BVP in the computational domain, making it available for numerical solution by finite elements. This procedure has been extensively used to solve wave BVPs, namely the Helmholtz equation [6, 7, 8], the time-dependent wave equation [9, 10, 11], and the elastodynamic equation, both in frequency domain [12, 13] and time domain [14, 15, 16]. In this context, ABCs are also called non-reflecting boundary conditions, since they are aimed at preventing any spurious reflection of waves from the artificial boundary. In addition, some important efforts have been made in order to extend the class of allowable artificial boundaries to more general shapes than a circle or a sphere [17, 18, 19]. The DtN FEM has also been successfully applied to solve elliptic BVPs in exterior domains. Han and Wu provided nonlocal and local approximations of the exact ABCs for the 2D Laplace equation [20] and 2D linear elastostatics [21], including FEM approximations and their error estimates. Explicit closed-form expressions for the DtN map associated with 2D and 3D Laplace equation and 2D linear elastostatics were given by Givoli and Keller [22]. More recently, Han and Zheng provided the exact ABCs for 3D linear elastostatics, as well as a FEM approximation and its error estimate [23]. Further error estimates were given by Han and Bao for general 2D scalar elliptic problems [24, 25] and 2D elastostatics [26]. A particularly important field of application of unbounded BVPs is geophysics, where the ground is usually modelled as an unbounded elastic domain. A typical example is the problem of determining deformations and stresses around structures of interest in the ground, such as an excavation in mining or civil engineering. Even though an exterior domain may be used as a model of the ground in a first approximation, a more realistic model needs to take into account the ground surface, which is customarily assumed, for simplicity, to be a plane boundary that extends to infinity, where a traction-free (Neumann) boundary condition holds. The resulting unbounded domain is called semi-infinite. To apply a DtN FEM approach to such a domain, a natural approach would be to consider as artificial boundary a half circle or half sphere surrounding the structure of interest, in analogy to the exterior case. This procedure works properly for some standard 2D scalar BVPs, such as the Laplace equation [1] and the Helmholtz equation [27] with Dirichlet or Neumann boundary conditions holding on the infinite plane boundary, since in these cases the BVP in the residual semi-infinite domain can be solved by separation of variables. However, the same procedure cannot be directly applied to linear elastostatics in semi-infinite domains, since in this case the method of separation of variables fails in solving the BVP in the residual domain in a fully analytical way, and thus it is not possible to get an explicit closed-form expression for the DtN map. Givoli and Vigdergauz [27] proposed an alternative to deal with this drawback in the 2D case. They considered a semi-circular artificial boundary and solved the BVP in the semi-infinite residual domain using complex analysis techniques, resulting in a solution in series form, whose coefficients can only be computed in an approximated way. A different approach was considered by Han et 2

al. [28], who employed the direct method of lines with semi-discretisation to solve the BVP in the residual domain lying outside a half circular artificial boundary, obtaining an approximation of the DtN map. Bao and Han [29] proposed a variation of this approach, where the artificial boundary is a polygonal line instead of a half circle, and additionally an error estimate for the method was provided. No previous works on DtN FEM procedures for 3D linear elastostatics in semi-infinite domains were found in the literature. In this paper, we present a DtN FEM for a semi-infinite elastic domain in 3D. A geometrical perturbation on the plane surface representing the structure of interest is considered, which is assumed to be axisymmetric about the vertical axis. Assuming further that the elastic domain is only under axisymmetric loading, the whole problem is treated as a BVP of axisymmetric elastostatics. The semi-infinite domain is truncated by introducing a semi-spherical artificial boundary surrounding the perturbation. Then a standard FEM discretisation in the resulting computational domain with exact ABCs in the artificial boundary is established at a theoretical level. As it is not possible to obtain a fully explicit closed-form expression for the DtN map, we proceed in a similar way to [27], solving the BVP in the residual domain just for particular Dirichlet data on the artificial boundary, in order to approximate those boundary integral terms involving the DtN map that occur in the finite element formulation. The BVP in the residual domain for each required case is solved by a semi-analytical method, based to a large extent on the procedure reported by the authors in a recent paper [30]. By applying separation of variables, a general analytical solution is calculated, expressed as a series with unknown coefficients, which are approximated by minimising a quadratic energy functional appropriately chosen. The minimisation yields a symmetric and positive definite linear system of equations for a finite number of coefficients, which is efficiently solved by exploiting its particular block structure, in such a way that the coefficients are computed by mere forward and backward substitutions. This procedure allows an approximate but effective coupling of the exact nonlocal ABCs provided by the DtN map with the FEM scheme. The method is validated by solving a model problem whose exact solution is available. The relative error between the numerical and the exact solution is analysed for different locations of the artificial boundary, number of terms in the series of the DtN map and mesh sizes. An application to a simplified geophysical case is also presented, in order to show the potentiality of the method. 2. Mathematical formulation 2.1. Generalities Let us consider the lower half-space with a geometrical perturbation on its plane surface which is local, i.e. it is restricted to a bounded region. This type of domain is oftentimes referred to as a locally perturbed half-space. In addition, the perturbation is assumed to be rotationally symmetric about a vertical axis, in such a way that the whole semi-infinite domain is actually an axisymmetric solid (or solid of revolution). Therefore, its geometry is completely described by its generating cross section, which consists in a 2D domain that we denote by Ω. The boundary of Ω consists of three parts: A vertical boundary of axisymmetry coinciding with the axis of revolution, denoted by Γs , a horizontal unperturbed boundary coinciding with the infinite plane surface of the half-space, denoted by Γ∞ , and a perturbed bounded boundary which is assumed piecewise smooth, denoted by Γp (see Figure 1). The domain Ω will be described with the aid of axisymmetric cylindrical coordinates (ρ, z) or spherical coordinates (r, φ) as appropriate, with the origin placed at the point of intersection between the axis of revolution and the plane boundary of the unperturbed half-space. Variables ρ, z, r and φ are linked by the relations r2 = ρ2 + z 2 ,

z φ = arctan , ρ 3

ρ = r sin φ,

z = r cos φ.

z

φ

r

Γ∞

ρ

Γp

Ω Γs

Figure 1: Axisymmetric semi-infinite domain.

ˆ ˆ, z ˆ , rˆ and φ. The unit vectors associated with variables ρ, z, r and φ are denoted respectively by ρ They are linked by the relations ˆ cos φ, ˆ = rˆ sin φ + φ ρ

ˆ sin φ, ˆ = rˆ cos φ − φ z ˆ =ρ ˆ cos φ − z ˆ sin φ. φ

ˆ sin φ + z ˆ cos φ, rˆ = ρ

(1a) (1b)

Thus an arbitrary position in Ω is expressed either as (ρ, z) or (r, φ), depending on the chosen system of coordinates. Moreover, we denote by θ the azimuthal angle, which is not necessary to describe geometrically the 2D domain Ω, but it is required in the axisymmetric elastostatic model, as we shall see below. 2.2. Axisymmetric elastostatic model We assume that Ω is occupied by an isotropic, homogeneous, linear elastic medium. Under the condition that the 3D axisymmetric domain is only subject to axisymmetric loading, the resulting elastic deformations will keep the axisymmetric nature of the problem. A generic displacement field defined in Ω is denoted by u and its associated stress tensor is denoted by σ = σ(u). By virtue of the axisymmetry, u has only components uρ and uz (resp. ur and uφ ), whereas σ has normal components σρ and σz (resp. σr and σφ ), a shear component σρz (resp. σrφ ), and an additional non-vanishing normal component σθ (see [31] for details). It is assumed that σ is given in terms of u by the isotropic Hooke’s law, that is, σ(u) = λ (∇ · u) I + μ(∇u + ∇uT ),

(2)

where λ, μ > 0 are the Lam´e constants of the elastic solid and I is the identity matrix. The Hooke’s law (2) is written by components as uρ ∂uρ ∂uz +λ +λ , ∂ρ ρ ∂z uρ ∂uρ ∂uz σθ (u) = λ + (λ + 2μ) + λ , ∂ρ ρ ∂z

σρ (u) = (λ + 2μ)

4

(3a) (3b)

uρ ∂uρ ∂uz + λ + (λ + 2μ) , ∂ρ ρ ∂z  ∂u ∂uz  ρ σρz (u) = μ + , ∂z ∂ρ σz (u) = λ

(3c) (3d)

in cylindrical coordinates and as  ur λ  ∂uφ ∂ur + 2λ + + cot φ uφ , ∂r r r ∂φ  ∂uφ ur 1 ∂ur σφ (u) = λ + 2(λ + μ) + (λ + 2μ) + λ cot φ uφ , ∂r r r ∂φ  ur 1  ∂uφ ∂ur σθ (u) = λ + 2(λ + μ) + λ + (λ + 2μ) cot φ uφ , ∂r r r ∂φ   1 ∂u ∂u u φ φ r σrφ (u) = μ + − , r ∂φ ∂r r σr (u) = (λ + 2μ)

(4a) (4b) (4c) (4d)

in spherical coordinates (see [31] for details). We consider the following elastostatic BVP in Ω: Find u : Ω → R2 such that ∇ · σ(u) = 0

in Ω,

(5a)

σ(u)ˆ z=0

on Γ∞ ,

(5b)

ˆ =f σ(u)n

on Γp ,

(5c)

on Γs ,

(5d)

as r → ∞,

(5e)

ˆ =u·ρ ˆ=0 σ(u)ˆ ρ·z



|u| = O r

−1



ˆ denotes the unit outward normal vector on Γp and f : Γp → R2 is a given piecewise where n continuous vector function. The equation of linear elastostatics (also known as the Navier equation) is given by (5a) in generic form, but its axisymmetric form is to be used, which will be made explicit in the next sections. Concerning the boundary conditions, Γ∞ is assumed traction-free as expressed by (5b), whereas a prescribed traction acting on Γp is considered, represented by the nonzero righthand side f in (5c). To allow for the axisymmetry, it is customary to assume Γs to be free of shear traction and constrained against normal displacement, as expressed by (5d). The BVP is completed with a decaying condition at infinity (5e). Remark 1. A particular application of (5) is obtained when Ω is subject to gravity. Denoting by  the mass density of the elastic medium and by g the acceleration of gravity, the lithostatic displacement ug is defined as gz 2 ˆ, z (6) ug (z) = − 2(λ + 2μ) which corresponds to the displacement of the unperturbed elastic half-space due to gravity. Substitution of (6) in (2) yields the lithostatic stress tensor σ(ug ). The displacement of the perturbed ˆ in half-space Ω due to gravity is obtained by adding to ug the solution u of (5) with f = −σ(ug )n (5c) (See [30] for more details). Remark 2. Although we have assumed as axisymmetric loading a prescribed traction on Γp , expressed by a non-homogeneous Neumann boundary condition (5c), the DtN FEM procedure to be employed for solving (5) also works if we assume instead as axisymmetric loading some prescribed displacement on Γp (i.e. a non-homogeneous Dirichlet boundary condition), or even a combination of both (e.g. mixed or Robin boundary conditions). 5

3. FEM formulation in the computational domain 3.1. Equivalent bounded boundary-value problem In order to apply the DtN FEM to solve (5), an artificial boundary is introduced to truncate Ω, which is chosen to be a quarter circle of radius R, denoted by ΓR (cf. Figure 2). Notice that if Ω is rotated about the axis of revolution to generate the 3D semi-infinite solid, ΓR gives rise to a hemispherical surface of radius R. The bounded domain lying inside ΓR is denoted by Ωi , which corresponds to the computational domain. The resulting bounded parts of boundaries Γ∞ and Γs are respectively denoted by Γi∞ and Γis (cf. Figure 2). In order to solve (5), we formulate the following mathematically equivalent BVP in Ωi : Find u : Ωi → R2 such that in Ωi ,

(7a)

σ(u)ˆ z=0

on Γi∞ ,

(7b)

ˆ =f σ(u)n

on Γp ,

(7c)

ˆ =u·ρ ˆ=0 σ(u)ˆ ρ·z

Γis ,

(7d)

on ΓR ,

(7e)

∇ · σ(u) = 0

on

σ(u)ˆ r = −Mu

where (7e) is the exact ABC on ΓR , expressed in terms of the DtN map M, which we assume for the moment to be known. The definition of M, as well as its numerical approximation, will be discussed in detail in the next section.

z

Γi∞

ρ

Γp

R

Ωi ΓR

Γis

Figure 2: Axisymmetric computational domain.

Remark 3. We have assumed that the medium occupying Ω is isotropic, homogeneous, and linear elastic. Nevertheless, the DtN FEM also applies if more general hypotheses in this regard are assumed, but restricted to the the computational domain, e.g. the medium in Ωi could be inhomogeneous, anisotropic, nonlinear, and even inelastic. Moreover, nonzero right-hand sides in (7a) or (7b) are allowed as well.

6

3.2. Weak formulation Throughout this section, we work in axisymmetric cylindrical coordinates (ρ, z). To state a weak (or variational) formulation of (7), we consider a Hilbert space containing physically admissible displacement fields in Ωi , defined as   V = v = (vρ , vz ) ∈ [H 1 (Ωi )]2 : vρ = 0 on Γis , where H 1 (Ωi ) is the usual Sobolev space of L2 functions with L2 first derivatives. Notice that the boundary condition of zero normal displacement in Γs is included in V. The weak formulation of (7) is expressed as: Find u ∈ V such that a(u, v) + b(u, v) = (f , v)Γp

∀ v ∈ V,

(8)

where a corresponds to the bilinear form associated with the weak formulation of the elastostatic equation, expressed in axisymmetric cylindrical coordinates (cf. [32]) as 

 ∂vρ vρ ∂vρ ∂vz ∂vz + σθ (u) + σz (u) + σρz (u) + ρ dρ dz, (9) σρ (u) a(u, v) = ∂ρ ρ ∂z ∂z ∂ρ Ωi whereas b is a bilinear form involving the DtN map, defined as  Mu · v dΓR , b(u, v) =

(10)

ΓR

and the right-hand side of (8) is

 (f , v)Γp =

Γp

f · v dΓp .

By substituting (3) in (9), expanding terms and rearranging appropriately, a is reexpressed as   



∂uρ ∂vρ ∂uz ∂vz ∂uρ ∂vz ∂uz ∂vρ + ρ dρ dz + λ + ρ dρ dz a(u, v) = (λ + 2μ) ∂ρ ∂ρ ∂z ∂z ∂ρ ∂z ∂z ∂ρ Ωi Ωi 

 

∂vρ ∂vz ∂uρ ∂uz +μ + + ρ dρ dz (11) ∂z ∂ρ ∂z ∂ρ Ωi  

  uρ v ρ ∂uρ ∂uz ∂vρ ∂vz +λ + vρ + uρ + dρ dz + (λ + 2μ) dρ dz. ∂ρ ∂z ∂ρ ∂z Ωi Ωi ρ The first three integrals in (11) are terms analogous to the 2D elastostatic case (except for the differential element), whereas the last two integrals are additional terms that arise due to the axisymmetry. 3.3. FEM discretisation In what follows, the weak formulation (8) is discretised by using standard conforming P 1 finite elements. Let us assume for simplicity that Γp is a polygonal line. We denote by Th a family of

regular triangular meshes of Ωi such that Ωi = T ∈Th T . The indexing parameter h measures the size of  each triangular element T and is defined as h = max{diam T : T ∈ Th }, with diam T = sup{ (ρ2 − ρ1 )2 + (z2 − z1 )2 : (ρ1 , z1 ), (ρ2 , z2 ) ∈ T }. For each triangular mesh T h , let V h be a finite-dimensional vector subspace of V consisting of continuous piecewise linear vector functions, defined as    V h = v h ∈ V : v h ∈ [C 0 (Ωi )]2 , v h T ∈ [P1 (T )]2 ∀ T ∈ T h , 7

where C 0 (Ωi ) denotes the space of continuous functions in Ωi and P1 (T ) stands for the space of polynomials of degree at most one defined in T . The discrete version of (8) is expressed as: Find uh ∈ V h such that ∀ vh ∈ V h. (12) a(uh , v h ) + b(uh , v h ) = (f , v h )Γp Let I be the set of all nodes of the triangular mesh T h . For each i ∈ I, we introduce a nodal shape (or basis) function ψi ∈ V h which has unit value at node i and zero value at all other nodes. Let also Is ⊂ I be the set of nodes lying on Γis . A nodal basis of V h is thus built by considering ˆ for nodes i ∈ I \ Is , to allow for the condition of zero normal displacement vector functions ψi ρ ˆ for all nodes i ∈ I. Consequently, the solution of (12) is given by (7d), and vector functions ψi z expressed as a linear combination of the nodal basis functions as   uh (ρ, z) = dρi ψi (ρ, z)ˆ ρ+ dzi ψi (ρ, z)ˆ z, (13) i∈I

i∈I\Is

where dρi and dzi are respectively the unknown nodal values of components uhρ and uhz of the solution uh . Replacing (13) in (12) and substituting v h by each one of the nodal basis functions ˆ and ψi z ˆ , we arrive at the finite element matrix form of the problem, expressed as ψi ρ Kd = F , 

with a

b

K =K +K ,

a

K = 

and d=

K aρρ K aρz

(14)

K azρ K azz dρ dz





 F =

K =

, 

,

b

Fρ

K bρρ K bρz K bzρ K bzz

 ,



Fz

,

where ˆ , ψj ρ ˆ ), [K aρρ ]ij = a(ψi ρ

ˆ , ψj z ˆ ), [K aρz ]ij = a(ψi ρ

ˆ , ψj ρ ˆ ), [K azρ ]ij = a(ψi z

ˆ , ψj z ˆ ), [K azz ]ij = a(ψi z

ˆ , ψj ρ ˆ ), [K bρρ ]ij = b(ψi ρ

ˆ , ψj z ˆ ), [K bρz ]ij = b(ψi ρ

ˆ , ψj ρ ˆ ), [K bzρ ]ij = b(ψi z

ˆ , ψj z ˆ ), [K bzz ]ij = b(ψi z

(16)

[dz ]j = dzj ,

(17)

ˆ ) Γp . [F z ]i = (f , ψi z

(18)

[dρ ]j = dρj , ˆ ) Γp , [F ρ ]i = (f , ψi ρ

(15)

a

The entries of matrix K and vector F are computed by standard numerical integration techniques. However, this is not possible in the case of matrix K b . Substituting (10) in (16), we reexpress the entries of K b as   b b ˆ · Mψj ρ ˆ dΓR , ˆ · Mψj z ˆ dΓR , ψi ρ [K ρz ]ij = ψi ρ (19a) [K ρρ ]ij = Γ R Γ R ˆ · Mψj ρ ˆ dΓR , ˆ · Mψj z ˆ dΓR . ψi z [K bzz ]ij = ψi z (19b) [K bzρ ]ij = ΓR

ΓR

Computing these terms is far from being straightforward, since they involve the DtN map M, which has to be calculated in an accurate way. Notice that these terms are different from zero only if i, j ∈ IR , where IR ⊂ I denotes the set of mesh nodes lying on ΓR . A correct calculation of the nonzero entries of matrix K b is critical for obtaining an accurate solution of (7), since K b accounts for the contribution from the exact ABCs in the standard finite element scheme used for numerical discretisation. This crucial issue is thoroughly treated in the next section. 8

Remark 4. As the domain Ωi is meshed with triangular elements, the artificial boundary ΓR is approximated by a polygonal line. This will limit the accuracy that the method can achieve, despite the choice of the basis functions. Provided that we are assuming piecewise linear basis functions, the associated error of the method will not be larger than that of finite elements. However, higher-order basis functions (e.g. piecewise quadratic) will not provide higher order accuracy. 4. Approximation of the exact artificial boundary conditions 4.1. Definition of the DtN map Next, we provide the mathematical definition of the DtN map M introduced in the previous section. Let us consider the residual semi-infinite domain lying outside the artificial boundary ΓR , denoted by Ωe . The resulting unbounded parts of boundaries Γ∞ and Γs are respectively denoted by Γe∞ and Γes (cf. Figure 3). Let H s (ΓR ) be the standard Sobolev space on ΓR of real order s. Given any v ∈ [H 1/2 (ΓR )]2 , we obtain Mv ∈ [H −1/2 (ΓR )]2 as  Mv = −σ(u)ˆ r Γ , (20) R

where u is the solution of the following BVP in

Ωe :

Find u : Ωe → R2 such that in Ωe ,

∇ · σ(u) = 0

(21a)

Γe∞ ,

(21b)

u=v

on ΓR ,

(21c)

ˆ =u·ρ ˆ=0 σ(u)ˆ ρ·z

Γes ,

(21d)

as r → ∞.

(21e)

σ(u)ˆ z=0

on



|u| = O r

−1



on

As already mentioned, unlike most exterior BVPs, in this case the method of separation of variables fails in obtaining a full analytical closed-form expression for the solution u of (21) in function of a generic Dirichlet datum v on ΓR . This drawback is overcome by solving (21) approximately, just for those v required to compute the nonzero entries of matrix K b . By virtue of (19), the required ˆ and v = ψj z ˆ , for every j ∈ IR . right-hand sides in (21c) are v = ψj ρ 4.2. General analytical solution In what follows, the method of separation of variables is applied to calculate a general analytical solution in series form satisfying (21a), (21b), (21d) and (21e). For this, we proceed analogously as in [30], where the solution originally proposed by Eubanks [33] was completed and improved. Throughout this section, we work in axisymmetric spherical coordinates (r, φ). The residual domain of Figure 3 is defined as Ωe = {(r, φ) : R < r < ∞, π/2 < φ < π}. First, a solution of (21a) is sought under the form 2μu = ∇(Φ + zΨ) − 4(1 − ν)Ψˆ z,

(22)

where ν = λ/2/(λ + μ) denotes the Poisson’s ratio of the elastic solid and Φ, Ψ are the so-called Boussinesq potentials, which are scalar harmonic functions. Decomposition (22) is reexpressed by components in (r, φ) as ∂Ψ ∂Φ (r, φ) + r cos φ (r, φ) − (3 − 4ν) cos φ Ψ(r, φ), ∂r ∂r ∂Ψ 1 ∂Φ 2μ uφ (r, φ) = (r, φ) + cos φ (r, φ) + (3 − 4ν) sin φ Ψ(r, φ). r ∂φ ∂φ 2μ ur (r, φ) =

9

(23a) (23b)

φ

Γe∞

r

R ΓR

Ωe Γes

Figure 3: Axisymmetric semi-infinite residual domain.

The Boussinesq potentials Φ, Ψ are then sought as general solutions to the Laplace equation in Ωe and with decay rates at infinity such that the solution u defined in (22) or (23) fulfils (21e). For this, it suffices to look for Ψ satisfying ΔΨ = 0



Ψ=O r

−1



in Ωe ,

(24a)

as r → ∞,

(24b)

in Ωe ,

(25a)

as r → ∞,

(25b)

and Φ satisfying ΔΦ = 0



|∇Φ| = O r

 −1

where Δ denotes the Laplacian. Applying standard separation of variables in (r, φ) to (24a) and (25a), and ruling out unbounded solutions at infinity, we arrive at the sets of potentials Ψn (r, φ) = r−(n+1) Pn (cos φ),

Φn (r, φ) = r−(n+1) Pn (cos φ),

n = 0, 1, 2, . . . ,

where Pn (·) denotes the Legendre polynomial of order n. Notice that when r tends to infinity, the asymptotically dominant term corresponds to the case n = 0. Then, as it holds that   Ψ0 (r, φ) = O r−1 as r → ∞, the set of potentials Ψn satisfies (24b). However, we have that   |∇Φ0 (r, φ)| = O r−2 as r → ∞, which means that the set of potentials Φn satisfies a more restrictive condition than (25b), so it is not a sufficiently general set of solutions of (25). We remedy this by incorporating into the set a logarithmic potential, defined as Φ−1 (r, φ) = ln(r − r cos φ), 10

(26)

which is continuous in Ωe , since the points where the logarithm is not well defined are not included in Ωe , and satisfies (25a). Furthermore, it holds that   |∇Φ−1 (r, φ)| = O r−1 as r → ∞, so the set of potentials Φn with the case n = −1 included fulfils (25b) exactly. The logarithmic potential (26) arises when solving the problem of a concentrated force acting normal to the free surface of an elastic half-space (known as the Boussinesq problem, cf. [31]). It is also known as the Boussinesq’s elementary solution of the second kind (cf. [34]). To get a solution of (21a), we build two sets of displacements by setting particular combinations of potentials Φn and Ψn in (22) or (23). The first set is obtained by setting Φ = Φn and Ψ = 0 for n = −1, 0, 1, . . ., i.e. 2μu(1) n = ∇Φn ,

(27)

whereas the second set is obtained by setting Φ = −(n − 4 + 4ν)Φn−1 and Ψ = (2n + 1)Ψn for n = 0, 1, 2, . . ., i.e.   ˆ . 2μu(2) (28) n = −(n − 4 + 4ν)∇Φn−1 + (2n + 1) ∇(zΨn ) − 4(1 − ν)Ψn z In fact, as stated in [33], this particular combination of potentials Φn and Ψn considerably simplifies the forthcoming calculations. Expanding expressions in (27) and (28), and rearranging terms, both sets of displacements are conveniently reexpressed as −(n+2) (1) u(1) wn (φ), n (r, φ) = r

n = −1, 0, 1, . . .

(29a)

−(n+1) (2) wn (φ), u(2) n (r, φ) = r

n = 0, 1, 2, . . .

(29b)

(1)

(2)

for certain vector functions wn and wn depending only on φ. Substituting components of (29a) and (29b) in (4) yields their associated stress tensors, which are similarly expressed as −(n+3) (1) σ (1) τ n (φ), n (r, φ) = r

n = −1, 0, 1, . . .

(30a)

−(n+2) (2) τ n (φ), σ (2) n (r, φ) = r

n = 0, 1, 2, . . .

(30b)

(1)

(2)

(1)

for certain tensor functions τ n and τ n depending only on φ. Full expressions of functions wn , (2) (1) (2) wn , τ n and τ n are provided in [30]. The general solution of (21a) then corresponds to an infinite linear combination of displacements defined in (29a) and (29b), which we conveniently express as u(r, φ) =

∞ 

  (2) (2) (1) r−(n+2) a(1) n w n (φ) + an+1 w n+1 (φ) ,

(31)

n=−1 (1)

(2)

for arbitrary real coefficients an and an . This solution satisfies the decaying condition at infinity (21e) as well. The associated stress tensor is written in analogous way as σ(r, φ) =

∞ 

  (2) (2) (1) r−(n+3) a(1) n τ n (φ) + an+1 τ n+1 (φ) .

(32)

n=−1

ˆ Conseˆ = −φ. In addition, the solution u has to satisfy (21b). On Γe∞ (φ = π/2) it holds that z quently, ˆ =0 σ(u)ˆ z = −σrφ (u)ˆ r − σφ (u)φ

⇒ 11

σφ (u) = σrφ (u) = 0

on Γe∞ .

Imposing this condition on (32) yields     (2) (2) (1) a(1) n [τn ]φ π/2 + an+1 [τn+1 ]φ π/2 = 0,     (2) (2) (1) a(1) n [τn ]rφ π/2 + an+1 [τn+1 ]rφ π/2 = 0. (1)

(33a) (33b)

(2)

To evaluate functions τ n , τ n at φ = π/2, it is necessary to distinguish between the cases n even and n odd. By doing so in (33) and expanding, we obtain the following relations between (1) (2) coefficients an and an : (1)

(2)

a−1 = (3 − 2ν)a0 , (1)

(34a)

(2)

(2n + 1)a2n = α2n a2n+1 (1)

(2)

(2n + 2)a2n+1 = α2n+2 a2n+2 where

n = 0, 1, 2, . . . ,

(34b)

n = 0, 1, 2, . . . ,

(34c)

α2n = (2n + 1)2 − 2(1 − ν). (1)

Incorporating relations (34) into (31)-(32) and using the fact that coefficients an arbitrary, allows us to reexpress u and σ as u(r, φ) =

∞ 

An

 R 2n+2

n=0

1 σ(r, φ) = R

 ∞ n=0

r An

∞ 

w(A) n (φ) +

Bn

 R 2n+3

n=−1

 R 2n+3 r

τ (A) n (φ)

+

∞  n=−1

r Bn

w(B) n (φ),

 R 2n+4 r (A)

τ (B) n (φ) (B)

where An and Bn are arbitrary real coefficients. Vector functions wn , wn (A) (B) τ n , τ n are defined as (1)

(2)

w(A) n (φ) = α2n w 2n (φ) + (2n + 1)w 2n+1 (φ) (B) w−1 (φ)

=

w(B) n (φ) = τ (A) n (φ) = (B)

τ −1 (φ) = τ (B) n (φ) =

(1) (2) (3 − 2ν)w−1 (φ) + w0 (φ), (1) (2) α2n+2 w2n+1 (φ) + (2n + 2)w2n+2 (φ) (1) (2) α2n τ 2n (φ) + (2n + 1)τ 2n+1 (φ) (1) (2) (3 − 2ν)τ −1 (φ) + τ 0 (φ), (1) (2) α2n+2 τ 2n+1 (φ) + (2n + 2)τ 2n+2 (φ)

are

(35) ,

(36)

and tensor functions

n = 0, 1, 2, . . . ,

(37a) (37b)

n = 0, 1, 2, . . . ,

(37c)

n = 0, 1, 2, . . . ,

(37d) (37e)

n = 0, 1, 2, . . . .

Explicit expressions for all these functions are provided next:   2μ[wn(A) ]r (φ) = −(2n + 1) α2n P2n (cos φ) + γ2n P2n+2 (cos φ) ,     (cos φ) + 2n P2n+2 (cos φ) , 2μ[wn(A) ]φ (φ) = − sin φ α2n P2n   (B) 2μ[w−1 ]r (φ) = − 1 − 2ν + 4(1 − ν) cos φ ,   (B) 2μ[w−1 ]φ (φ) = sin φ 3 − 4ν − (1 − 2ν)q(φ) ,   2μ[wn(B) ]r (φ) = −(2n + 2) α2n+2 P2n+1 (cos φ) + γ2n+1 P2n+3 (cos φ) ,     (cos φ) + 2n+1 P2n+3 (cos φ) , 2μ[wn(B) ]φ (φ) = − sin φ α2n+2 P2n+1 12

(2)

and an

(37f)

(38a) (38b) (39a) (39b) (39c) (39d)

  [τn(A) ]r (φ) = (2n + 1)(2n + 2) α2n P2n (cos φ) + β2n P2n+2 (cos φ) ,  (cos φ) − (2n + 1)(2n + 2) [τn(A) ]φ (φ) = (α2n + 2n )P2n+1   × α2n P2n (cos φ) + (α2n − 2n + 2 − 4ν)P2n+2 (cos φ) ,  [τn(A) ]θ (φ) = −(4n + 3) (2n + 1)(2n + 2)(1 − 2ν)P2n+2 (cos φ)   (cos φ) , + (2n − 1 + 2ν)P2n+1     (cos φ) + (2n + 1)α2n+1 P2n+2 (cos φ) , [τn(A) ]rφ (φ) = sin φ (2n + 2)α2n P2n (B)

[τ−1 ]r (φ) = 1 − 2ν + 2(2 − ν) cos φ,   (B) [τ−1 ]φ (φ) = −(1 − 2ν) 1 + cos φ − q(φ) ,   (B) [τ−1 ]θ (φ) = −(1 − 2ν) cos φ + q(φ) ,   (B) [τ−1 ]rφ (φ) = −(1 − 2ν) sin φ 1 − q(φ) ,   [τn(B) ]r (φ) = (2n + 2)(2n + 3) α2n+2 P2n+1 (cos φ) + β2n+1 P2n+3 (cos φ) ,    (cos φ) − (2n + 2)(2n + 3) [τn(B) ]φ (φ) = α2n+2 + 2n+1 P2n+2   × α2n+2 P2n+1 (cos φ) + (α2n+1 − 2n + 1 − 4ν)P2n+3 (cos φ) ,  [τn(B) ]θ (φ) = −(4n + 5) (2n + 2)(2n + 3)(1 − 2ν)P2n+3 (cos φ)   (cos φ) , + (2n + 1 + 2ν)P2n+2     (cos φ) + (2n + 2)P2n+3 (cos φ) , [τn(B) ]rφ (φ) = sin φ α2n+2 (2n + 3)P2n+1

(40a) (40b) (40c) (40d) (41a) (41b) (41c) (41d) (41e) (41f) (41g) (41h)

where n = 0, 1, 2, . . ., q(φ) = 1/(1 − cos φ) and β2n = (2n + 2)(2n + 5) − 2ν, γ2n = (2n + 2)(2n + 5 − 4ν), 2n = (2n + 1)(2n − 2 + 4ν). It can be easily verified that by virtue of (38b), (39b), (39d), (40d), (41d) and (41h), u and σ satisfy the condition r > R, uφ (r, π) = σrφ (r, π) = 0 which is equivalent to satisfy (21d). Consequently, the general solution in series form given in (35) and (36) satisfies (21a), (21b), (21d) and (21e). Notice that it is a fully analytical solution, since no numerical approximation has been introduced up to now. Full details about the procedure to obtain this solution are found in [30]. 4.3. Numerical enforcement of exact boundary conditions on the artificial boundary In what follows, the general analytical solution calculated in the previous subsection is numerically forced to satisfy (21c). At a theoretical level, there exist particular values of coefficients An and Bn such that (21c) holds, which of course depend on the particular Dirichlet datum v. However, these values are determined by an infinite system of simultaneous linear equations, so in general there is no way to calculate them analytically. To overcome this drawback, we compute numerically a finite number of coefficients, which, when substituted in (35) and (36), give rise to a semi-analytical solution of (21). The numerical procedure to compute coefficients An and Bn is described next for a generic right-hand side v in (21c). In the next subsection, this procedure will ˆ and v = ψj z ˆ . First, we truncate the series be applied to the particular right-hand sides v = ψj ρ

13

in (35) and (36) at a finite order N and we call uN and σ N the resulting displacement vector and stress tensor, that is, uN (r, φ) =

N 

An

 R 2n+2

n=0

r

w(A) n (φ) +

N 

Bn

 R 2n+3

n=−1

r

w(B) n (φ),

N N  R 2n+4  1   R 2n+3 (A) (B) σ N (r, φ) = An τ n (φ) + Bn τ n (φ) , R r r n=0

(42)

(43)

n=−1

and we consider the quadratic energy functional J defined as   1 1 ˆ · uN dΓR − ˆ · v dΓR , σN n σN n J (uN ) = 2R ΓR R ΓR

(44)

ˆ = −ˆ where n r is the unit normal vector on ΓR , pointing outwards Ωe . The first term in (44) is quadratic in uN and represents the surface elastic potential energy on ΓR , whereas the second term is linear in uN and is related to the Dirichlet boundary condition assumed on ΓR . The functional J has been defined in such a way that its minimisation provides numerical values for coefficients A0 , A1 , A2 , . . . , AN and B−1 , B0 , B1 , . . . , BN so that (21c) holds approximately. By expanding terms ˆ = −ˆ in (44), replacing n r and making explicit the integrals (notice that dΓR = R2 sin φ dφ), we arrive at   R π  [σN ]r (R, φ) [uN ]r (R, φ) + [σN ]rφ (R, φ) [uN ]φ (R, φ) sin φ dφ J (uN ) = − 2 π/2  π (45)   +R [σN ]r (R, φ) vr (R, φ) + [σN ]rφ (R, φ) vφ (R, φ) sin φ dφ. π/2

Let us define the column vectors A ∈ RN +1 and B ∈ RN +2 as those whose entries are the coefficients A0 , A1 , A2 , . . . , AN and B−1 , B0 , B1 , . . . , BN , respectively. By evaluating the corresponding components of uN and σ N at r = R in (42) and (43), substituting in (45) and expanding terms appropriately, J is reexpressed as an explicit function of A and B as J (A, B) =

N N N N N N 1   (AB) 1   (BA) 1   (AA) Qnk An Ak + Qnk An Bk + Qnk Bn Ak 2 2 2 n=0 k=0

1 + 2 where (αβ)

Qnk

 =−

N 

n=−1 k=−1



π π/2



and yn(α)

N 

=−

n=0 k=−1

(BB)

Qnk

Bn Bk −

N 

yn(A) An −

n=0

n=−1 k=0

N 

yn(B) Bn ,

n=−1

 (β) (β) [wn(α) ]r (φ)[τk ]r (φ) + [wn(α) ]φ (φ)[τk ]rφ (φ) sin φ dφ π

π/2



(46)

 [τn(α) ]r (φ)vr (R, φ) + [τn(α) ]rφ (φ)vφ (R, φ) sin φ dφ

α, β = A, B,

(47)

α = A, B.

(48)

Substitution of (38a), (38b), (39a), (39b), (39c), (39d), (40a), (40d), (41a), (41d), (41e) and (41h) (αβ) in (47), yields expressions for coefficients Qnk in terms of explicit integrals. These expressions are too cumbersome to be presented here, however, all the involved integrals are computed exactly with the aid of formulae provided in the Appendix A of [30]. The resulting explicit formulae for 14

(αβ)

(AA)

(AA)

(BA)

coefficients Qnk are provided in Appendix A. In particular, it holds that Qnk = Qkn , Qkn = (AB) (BB) (BB) (α) Qnk and Qnk = Qkn . Concerning coefficients yn , they depend on the particular right-hand side v assumed in (21c), so their calculation is discussed in the next subsection. From coefficients (AA) (AB) (BB) (A) (B) Qnk , Qnk , Qnk , yn and yn used as entries, we assemble the matrices Q(AA) , Q(AB) , Q(BB) and the vectors y (A) , y (B) , respectively. By virtue of (A.1)-(A.3), Q(AA) is a tridiagonal symmetric matrix of size N + 1, Q(AB) is a non-square full matrix of size (N + 1) × (N + 2), and Q(BB) is a symmetric matrix of size N + 2, which is almost tridiagonal, except for its first row and column (A) (B) that are full. Figure 4 shows schematically the structure of these matrices. Vectors yn and yn are of length N + 1 and N + 2, respectively. Additionally, we assemble the square matrix Q of size 2N + 3 and the vectors x, y of length 2N + 3 as (A) A y Q(AA) Q(AB) , x= , y= , Q= B y (B) [Q(AB) ]T Q(BB) which allow us to reexpress J in (46) as a quadratic form in x 1 J (x) = xT Qx − xT y. 2

(49)

It is straightforward to verify that Q is a symmetric matrix. It is in addition a positive definite matrix, due to the elliptic nature of the equation of elastostatics. Therefore, the quadratic form J has a unique global minimum that is reached when ∇J (x) = Qx − y = 0, or equivalently, Qx = y,

(50)

which is a linear system of equations for the coefficients An and Bn (stored within the vector x). To solve (50), we take advantage of the symmetry and positive definiteness of matrix Q. Expressing once again Q, x and y in terms of their block components, (50) is rewritten as

Q(AA)

Q(AA) [Q(AB) ]T

Q(AB) Q(BB)



Q(AB)

A B



=

y (A) y (B)

,

(51)

Q(BB)

N +1

N +2

N +1

N +2

N +2

Figure 4: Structure of matrices Q(AA) ,Q(AB) and Q(BB) .

15

and the solution of (51) is expressed explicitly with the aid of the Schur-Banachiewicz blockwise inversion formula (cf. [35]):    (BB) ]−1 [Q(AB) ]T [Q(AA) ]−1 y (A) A = [Q(AA) ]−1 + [Q(AA) ]−1 Q(AB) [Q (52a)  (BB) ]−1 y (B) , − [Q(AA) ]−1 Q(AB) [Q  (BB) ]−1 y (B) ,  (BB) ]−1 [Q(AB) ]T [Q(AA) ]−1 y (A) + [Q B = −[Q

(52b)

 (BB) denotes the Schur complement of Q(BB) in Q, defined as where Q  (BB) = Q(BB) − [Q(AB) ]T [Q(AA) ]−1 Q(AB) . Q

(53)

Therefore, we observe from (52) that, in order to get A and B, it suffices to invert matrices Q(AA)  (BB) , being both positive definite thanks to the positive definiteness of Q. As Q(AA) is a and Q tridiagonal matrix, it is inverted by using the Thomas algorithm for tridiagonal systems (cf. [36],  (BB) , we use its Cholesky factorisation, which is Algorithm 4.3.6), whereas in order to invert Q computed using a standard algorithm for that purpose (cf. [36], Algorithm 4.2.1). Hence, both inversions are very fast and efficient, since in practice they are performed by applying successive forward and backward substitutions. The resulting algorithm to compute A and B for a given Dirichlet datum v consists of the following steps: 1. 2. 3. 4. 5. 6.

Compute the coefficients for tridiagonal inversion of Q(AA) . Use them to evaluate [Q(AA) ]−1 y (A) and [Q(AA) ]−1 Q(AB) .  (BB) using (53). Calculate the Schur complement Q  (BB) . Compute the Cholesky factorisation of Q  (BB) ]−1 y (B) and [Q  (BB) ]−1 [Q(AB) ]T . Use it to evaluate [Q Evaluate A and B as indicated in (52a) and (52b).

This algorithm has to be applied repeatedly for each pair of vectors y (A) and y (B) necessary in the computations. Nevertheless, as matrices Q(αβ) remain unchanged, each one of them, together with  (BB) , are computed and stored only once. the coefficients for inversion of Q(AA) and Q 4.4. Numerical approximation of integral terms involving the DtN map in the FEM formulation In what follows, the numerical procedure described in the previous subsection is applied to (α) ˆ and approximate the terms in (19). First, coefficients yn in (48) are to be calculated for v = ψj ρ ˆ for every j ∈ IR . Combining with (1a), we obtain that these two cases are respectively v = ψj z equivalent to consider y (A) = C jρ , y (B) = D jρ , and

y (A) = C jz ,

y (B) = D jz ,

where C jρ , C jz ∈ RN +1 and D jρ , D jz ∈ RN +2 are the vectors of components  π   (A) j Cρ,n =− n = 0, . . . , N, [τn ]r (φ) sin φ + [τn(A) ]rφ (φ) cos φ ψj (R, φ) sin φ dφ π/2  π   (A) j Cz,n = − n = 0, . . . , N, [τn ]r (φ) cos φ − [τn(A) ]rφ (φ) sin φ ψj (R, φ) sin φ dφ π/2  π   (B) j Dρ,n = − n = −1, . . . , N, [τn ]r (φ) sin φ + [τn(B) ]rφ (φ) cos φ ψj (R, φ) sin φ dφ π/2  π   (B) j Dz,n = − n = −1, . . . , N, [τn ]r (φ) cos φ − [τn(B) ]rφ (φ) sin φ ψj (R, φ) sin φ dφ π/2

16

(54a) (54b) (54c) (54d)

respectively. By substituting (40a), (40d), (41a), (41d), (41e) and (41h) in (54) and expanding, components of vectors C jρ , C jz , D jρ , D jz are rewritten as   j j j = −(2n + 1)(2n + 2) α2n ls,2n + β2n ls,2n+2 Cs,n − (2n + 2)α2n mjs,2n − (2n + 1)α2n+1 mjs,2n+2 j Ds,−1 j Ds,n

=

j −3 ls,1

− (1 −

n = 0, . . . , N,

2ν)qsj ,

  j j = −(2n + 2)(2n + 3) α2n+2 ls,2n+1 + β2n+1 ls,2n+3   − α2n+2 (2n + 3)mjs,2n+1 + (2n + 2)mjs,2n+3

(55a) (55b)

n = 0, . . . , N,

(55c)

where s = ρ, z and  π  π j j = sin2 φPn (cos φ)ψj (R, φ)dφ, lz,n = sin φ cos φPn (cos φ)ψj (R, φ)dφ, (56a) lρ,n π/2 π/2  π  π j 2  j mρ,n = sin φ cos φPn (cos φ)ψj (R, φ)dφ, mz,n = − sin3 φPn (cos φ)ψj (R, φ)dφ, (56b) π/2 π/2  π j qρ = (1 + cos φ)ψj (R, φ)dφ, qzj = 0. (56c) π/2

To compute these integrals, we observe that the mesh Th restricted to ΓR gives rise to a partition of the quarter circle into segments (usually equispaced). Let J be the number of these segments. Then the set IR necessarily consists of J + 1 mesh nodes, which we assume for simplicity to be correlatively numbered, that is, IR = {1, 2, . . . , J, J + 1}. We associate with each node j ∈ IR an angle φj , in such a way that π/2 = φ1 < φ2 < . . . < φJ < φJ+1 = π. The P 1-nodal shape functions ψj , restricted to ΓR , do not depend on the radius R, and have the following explicit expressions

 φ2 − φ ψ1 (φ) = 1[φ1 ,φ2 ] (φ) , (57a) φ2 − φ1



  φ − φj−1 φj+1 − φ ψj (φ) = 1[φj−1 ,φj ] (φ) + 1[φj ,φj+1 ] (φ) j = 2, . . . , J, (57b) φj − φj−1 φj+1 − φj 

φ − φJ ψJ+1 (φ) = 1[φJ ,φJ+1 ] (φ) , (57c) φJ+1 − φJ where 1[a,b] (·) stands for the indicator function of the interval [a, b]. Substituting (57) in (56), s , ms and q ρ are decomposed as coefficients ln,j n,j j 1 1,+ ls,n = ls,n ,

m1s,n = m1,+ s,n ,

1 1,+ qρ,n = qρ,n ,

j j,− j,+ ls,n = ls,n + ls,n ,

j,+ mjs,n = mj,− s,n + ms,n ,

j j,− j,+ qρ,n = qρ,n + qρ,n ,

J+1 ls,n

=

J+1,− ls,n ,

mJ+1 s,n

=

mJ+1,− , s,n

J+1 qρ,n

=

j = 2, . . . , J,

J+1,− qρ,n ,

where for j = 2, . . . , J + 1, j,− lρ,n

mj,− ρ,n qρj,−

 φ − φj−1 dφ, = sin φPn (cos φ) φj − φj−1 φj−1 

 φj φ − φj−1 dφ, = sin2 φ cos φPn (cos φ) φj − φj−1 φj−1 

 φj φ − φj−1 dφ, = (1 + cos φ) φj − φj−1 φj−1 

φj



2

17

(58a) (58b) (58c)

j,− lz,n

mj,− z,n



 φ − φj−1 = sin φ cos φPn (cos φ) dφ, φj − φj−1 φj−1 

 φj φ − φj−1 3  dφ, =− sin φPn (cos φ) φj − φj−1 φj−1 

φj

(58d) (58e)

and for j = 1, . . . , J, j,+ lρ,n

mj,+ ρ,n qρj,+ j,+ lz,n

mj,+ z,n

 φj+1 − φ dφ, = sin φPn (cos φ) φj+1 − φj φj 

 φj+1 φj+1 − φ 2  dφ, = sin φ cos φPn (cos φ) φj+1 − φj φj 

 φj+1 φj+1 − φ dφ, = (1 + cos φ) φj+1 − φj φj

  φj+1 φj+1 − φ = sin φ cos φPn (cos φ) dφ, φj+1 − φj φj 

 φj+1 φj+1 − φ 3  dφ. =− sin φPn (cos φ) φj+1 − φj φj 

φj+1



2

(59a) (59b) (59c) (59d) (59e)

Integrals in (58) and (59) are computed analytically with the aid of symbolic computation techniques. The resulting expressions become more and more unwieldy as n increases, so only a few of j the coeffithem are provided in Appendix B by the way of example. We denote by Ajs,n and Bs,n cients An and Bn obtained by applying the algorithm of the previous subsection with y (A) = C js and y (B) = D js , respectively. Notice that as for each j = 1, 2, . . . , J, J + 1 two cases need to be considered (s = ρ and s = z), the algorithm is to be applied 2(J + 1) times in total. Additionally, let us define the vectors Ajs ∈ RN +1 and B js ∈ RN +2 as those whose components are the coefficients j , respectively, which by virtue of (52) are given by Ajs,n and Bs,n    (BB) ]−1 [Q(AB) ]T [Q(AA) ]−1 C j Ajs = [Q(AA) ]−1 + [Q(AA) ]−1 Q(AB) [Q s (60a) (AA) −1 (AB)  (BB) −1 j − [Q ] Q [Q ] D , s

 (BB) ]−1 [Q(AB) ]T [Q(AA) ]−1 C j + [Q  (BB) ]−1 D j . B js = −[Q s s

(60b)

Thus, using the definition of the DtN map (20) and combining with (43) yields the approximations N N     1  j  (A) (A) j (B) (B) ˆ ˆ ˆ (φ) ≈ − Aρ,n [τn ]r (φ)ˆ r + [τn ]rφ (φ)φ + Bρ,n [τn ]r (φ)ˆ r + [τn ]rφ (φ)φ , Mψj ρ R n=0 n=−1 N N     1  j  (A) (A) j (B) (B) ˆ ˆ ˆ (φ) ≈ − Mψj z Az,n [τn ]r (φ)ˆ r + [τn ]rφ (φ)φ + Bz,n [τn ]r (φ)ˆ r + [τn ]rφ (φ)φ . R n=0

n=−1

By replacing these expressions in (19), combining with (1a) and expanding, we arrive at [K bρρ ]ij

≈−R −R

N 

Ajρ,n

n=0 N  n=−1



j Bρ,n

π/2



 [τn(A) ]r (φ) sin φ + [τn(A) ]rφ (φ) cos φ ψi (R, φ) sin φ dφ



π

π π/2

 [τn(B) ]r (φ) sin φ + [τn(B) ]rφ (φ) cos φ ψi (R, φ) sin φ dφ,



18

[K bρz ]ij

≈−R −R

N 

Ajz,n

n=0 N 



π/2



j Bz,n

−R

N 

Ajρ,n

n=0 N 



≈−R −R

N 

n=0 N 





 [τn(B) ]r (φ) cos φ − [τn(B) ]rφ (φ) sin φ ψi (R, φ) sin φ dφ,



π π/2



π π/2

j Bz,n

n=−1

 [τn(A) ]r (φ) cos φ − [τn(A) ]rφ (φ) sin φ ψi (R, φ) sin φ dφ

π/2

j Bρ,n

Ajz,n

 [τn(B) ]r (φ) sin φ + [τn(B) ]rφ (φ) cos φ ψi (R, φ) sin φ dφ,



π

n=−1

[K bzz ]ij



π π/2

n=−1

[K bzρ ]ij ≈ − R

 [τn(A) ]r (φ) sin φ + [τn(A) ]rφ (φ) cos φ ψi (R, φ) sin φ dφ



π



π π/2

 [τn(A) ]r (φ) cos φ − [τn(A) ]rφ (φ) sin φ ψi (R, φ) sin φ dφ 

 [τn(B) ]r (φ) cos φ − [τn(B) ]rφ (φ) sin φ ψi (R, φ) sin φ dφ,

and combining with (54), these terms are conveniently reexpressed as [K bρρ ]ij

≈R

 N

i Cρ,n Ajρ,n

+

n=0

[K bρz ]ij

≈R

 N

≈R

 N

i Cρ,n Ajz,n

+

≈R

 N

n=0

N 

i j Dρ,n Bz,n

n=−1 i Cz,n Ajρ,n

+

n=0

[K bzz ]ij

i j Dρ,n Bρ,n

n=−1

n=0

[K bzρ ]ij

N 

N 

i j Dz,n Bρ,n

n=−1 i Cz,n Ajz,n

+

N 

i j Dz,n Bz,n



  = R C iρ · Ajρ + D iρ · B jρ ,

(61a)

  = R C iρ · Ajz + D iρ · B jz ,

(61b)

  = R C iz · Ajρ + D iz · B jρ ,

(61c)

  = R C iz · Ajz + D iz · B jz .

(61d)

n=−1

It should be observed that vectors C iρ , C iz , D iρ , D iz actually do not depend on the radius R, and so do vectors Aiρ , Aiz , B iρ , B iz . Hence, by virtue of (61), the entries of matrix K b are linear in R. 5. Numerical experiments 5.1. Model problem We introduce next a model problem whose exact solution is available. Let us assume the perturbed boundary Γp to be a quarter circle of radius a > 0. In this case, it holds in axisymmetric spherical coordinates that Ω = {(r, φ) : a < r < ∞, π/2 < φ < π}, Γp = {(r, φ) : r = a, π/2 < φ < π}. ˆ : Γp → R2 as We define the function f = fr rˆ + fφ φ   fr (φ) = 2μc 1 − 2ν + 2(2 − ν) cos φ , 19

(62a)



fφ (φ) = 2μc (1 − 2ν)

 cos φ sin φ , 1 − cos φ

(62b)

ˆ : Ω → R2 , defined as where c is a scale factor. Then the displacement u = ur rˆ + uφ φ  c a2  1 − 2ν + 4(1 − ν) cos φ , r

 c a2 1 − 2ν uφ (r, φ) = 3 − 4ν − sin φ, r 1 − cos φ ur (r, φ) = −

(63a) (63b)

is an exact solution of (5) with f given by (62). This solution is obtained from (35) with R replaced by a and setting An = Bn = 0 for all n = 0, 1, 2, . . . and B−1 = 2μc a. The stress components associated with (63) are obtained by setting these same values of coefficients An and Bn in (36) with R replaced by a, or simply by substituting (63) in (4), yielding  2μc a2  1 − 2ν + 2(2 − ν) cos φ , 2 r 

cos2 φ 2μc a2 (1 − 2ν) , σφ (r, φ) = r2 1 − cos φ 

2μc a2 (1 − 2ν) 1 + cos φ − cos2 φ , σθ (r, φ) = − r2 1 − cos φ 

2μc a2 (1 − 2ν) cos φ sin φ . σrφ (r, φ) = r2 1 − cos φ σr (r, φ) =

ˆ + uz z ˆ , with In axisymmetric cylindrical coordinates, (63) reads u = uρ ρ

 1 − 2ν c a2 ρ z , uρ (ρ, z) = − 2 + (ρ + z 2 )1/2 (ρ2 + z 2 )1/2 − z ρ2 + z 2

 z2 c a2 2(1 − ν) + , uz (ρ, z) = − 2 ρ2 + z 2 (ρ + z 2 )1/2

(64a) (64b) (64c) (64d)

(65a) (65b)

and associated stress components given by

 1 − 2ν 2μc a2 3ρ2 z σρ (ρ, z) = 2 , + (ρ + z 2 )1/2 (ρ2 + z 2 )1/2 − z (ρ2 + z 2 )2

 1 2μc a2 (1 − 2ν) z , + σθ (ρ, z) = − (ρ2 + z 2 )1/2 (ρ2 + z 2 )1/2 − z ρ2 + z 2 6μc a2 z 3 , σz (ρ, z) = 2 (ρ + z 2 )5/2 6μc a2 ρz 2 , σρz (ρ, z) = 2 (ρ + z 2 )5/2

(66a) (66b) (66c) (66d)

which follow from substituting (65) in (3). This exact solution will be used as a benchmark to test the performance of the DtN FEM approach. 5.2. Implementation aspects To apply the DtN FEM to solve the model problem, the radius of the artificial boundary R must be, of course, greater than a. The resulting computational domain Ωi is described in axisymmetric spherical coordinates (r, φ) as Ωi = {(r, φ) : a < r < R, π/2 < φ < π}. 20

In general, in any DtN FEM approach, the radius R is to be chosen large enough so that any possible irregularity (anisotropy, inhomogeneity) is enclosed by ΓR , but at the same time small enough to minimise the size of Ωi and thus the number of finite elements considered (cf. [27]). In our case, for the purpose of solving the model problem in order to test the accuracy of the DtN FEM, we shall consider values lying between R = 1.2a and R = 2a. In order to maintain a constant mesh quality, we employ structured triangular meshes of Ωi generated as follows: We assume equispaced partitions of intervals [a, R] and [π/2, π] consisting respectively of I and J segments, delimited by nodes satisfying a = r1 < r2 < . . . < rI < rI+1 = R and π/2 = φ1 < φ2 < . . . < φJ < φJ+1 = π (the latter partition already introduced in Section 4), in such a way that ri = a + (i − 1)(R − a)/I and φj = (π/2)(1 + (j − 1)/J). The cartesian product of both partitions gives rise to a structured quadrilateral mesh of Ωi , which consists of (I + 1)(J + 1) nodes and IJ elements. If each quadrilateral element is split by its diagonal, two triangles are created. Then the resulting triangular mesh consists of (I + 1)(J + 1) nodes and 2IJ elements in total. In order to ensure good quality triangles, we set the ratio I/J to be proportional to the quotient between the length of interval [a, R] and the length of quarter circle of radius (R + a)/2, which yields the relation   4 R−a J , (67) I= π R+a where · stands for the nearest integer function. Figure 5 shows two examples of such structured meshes for R = 1.5a. The mesh size h is calculated using relatively simple geometrical considerations, arriving at the formula

 π 1/2   1 2 h= (R − a) + 2IR (I − 1)R + a 1 − cos . I 2J

(a)

(b)

Figure 5: Two structured triangular meshes for R = 1.5a. (a) I = 5, J = 20. (b) I = 10, J = 40.

21

(68)

To choose an optimal value of the truncation order N for fixed R and h, we use the simple crude formula first proposed in [6] for the Laplace equation, and subsequently modified and used in [22] and [27] for the 2D elasticity equation. According to this formula, and given that we are using P 1 finite elements, the optimal N is estimated as   ln (h/a) Nopt = − 2 , (69) ln (R/a) where · stands for the ceiling function. This formula provides a good estimation of N , as we shall see below. 5.3. Accuracy and effectiveness Next, the model problem is solved by the DtN FEM procedure in the structured meshes described above. Numerical results are presented and the accuracy of the method is tested by comparing the numerical solution to the exact one given in (65) and (66), evaluated at the same mesh. The effect of varying parameters R, N and h is investigated. Throughout this subsection, we assume a perturbation radius a = 600 m and an elastic solid with Young’s modulus E = 70 GPa and Poisson’s ratio ν = 0.3. The Lam´e constants are immediate from E and ν through the standard formulae λ = Eν/(1 + ν)/(1 − 2ν) and μ = E/2/(1 + ν). In order to get displacements and stresses of physically realistic magnitudes, we set c = 10−4 . Given a fixed radius of artificial boundary R and a fixed mesh size h, the truncation order N is calculated from formula (69) unless otherwise is indicated. Figure 6 shows the components of the computed displacement vector uh for R = 1.5a and a mesh with I = 61 and J = 240. The associated stress tensor σ h is numerically computed from uh by substituting the displacement components in (3), with the derivatives calculated numerically in the mesh. The components of σ h are presented in Figure 7.

(a)

(b)

Figure 6: Computed displacement field. (a) Component uhρ (ρ, z). (b) Component uhz (ρ, z).

22

(a)

(b)

(c)

(d)

Figure 7: Computed stress tensor. (a) Component σρh (ρ, z). (b) Component σθh (ρ, z). (c) Component σzh (ρ, z). (d) h Component σρz (ρ, z).

To test preliminarily the effectiveness of our DtN FEM approach, we set R = 1.5a and the numerical solutions calculated in meshes of decreasing size are compared against the exact solution at some specific points of the domain. Table 1 presents numerical values of some displacement and stress components evaluated at the extreme points of ΓR . From this table, we see that, as the mesh size decreases, the values of the numerical solution approximate the values of the exact one at the considered points, just as expected. Moreover, in the case of the displacement, this approximation is considerably faster. 23

Table 1: Some components of the solution evaluated at points (0, R) and (−R, 0).

Solution Num., I = 3, J = 12 Num., I = 15, J = 58 Num., I = 26, J = 103 Num., I = 38, J = 149 Num., I = 49, J = 194 Num., I = 61, J = 240 Exact

uρ (0, R) [m] -0.01709 -0.01608 -0.01603 -0.01602 -0.01601 -0.01601 -0.01600

uz (−R, 0) [m] -0.09664 -0.09607 -0.09603 -0.09601 -0.09601 -0.09601 -0.09600

σρ (0, R) [MPa] 0.29365 0.83454 0.88596 0.90597 0.91680 0.92379 0.95726

σz (−R, 0) [MPa] -7.58848 -7.28673 -7.24649 -7.22431 -7.21383 -7.20685 -7.17949

In order to study more rigorously the accuracy of the DtN FEM procedure and to test the effect of varying R, N and h, we employ the relative errors between the numerical and the exact solution, in L2 -norm and H 1 -seminorm, defined respectively as E0,Ωi =

uh − Πh u 0,Ωi

Πh u 0,Ωi

,

E1,Ωi =

|uh − Πh u|1,Ωi |Πh u|1,Ωi

,

where · 0,Ωi and | · |1,Ωi stand for the usual L2 -norm and H 1 -seminorm in Ωi , respectively, and Πh u denotes the Lagrange interpolation of the exact solution u over the mesh of size h. Firstly, we test how the location of the artificial boundary affects the accuracy of the method. For this purpose, the radius R is varied from 1.2a to 2a in increments of 0.1a. In each one of the respective meshes, numbers I and J are set in order to keep the size h relatively constant. The parameters of the resulting meshes are given in Table 2, including the optimal value of N in terms of R and h given by (69). Figure 8 presents the relative errors in L2 -norm and H 1 -seminorm as functions of R, where we observe that both of them decrease as the artificial boundary moves away from the perturbation, and the error in L2 -norm is about 2 orders of magnitude lower than the error in H 1 -seminorm. Nevertheless, both errors are very small, even for a relatively close artificial boundary (when R = 1.2a, the respective values are 0.00782% and 0.78831%). Therefore, the location of the artificial boundary affects to some extent the accuracy of the DtN FEM procedure in solving the model problem. A more distant location (i.e. a larger computational domain) results in a better accuracy, however, the accuracy achieved for a relatively close location (i.e. a small computational domain) is more than acceptable. Table 2: Parameters of structured meshes for testing the effect of R.

R 1.2a 1.3a 1.4a 1.5a 1.6a 1.7a 1.8a 1.9a 2.0a

I 15 23 31 39 48 56 65 74 83

J 128 136 145 153 162 170 179 188 196

No. nodes 2064 3288 4672 6160 7987 9747 11880 14175 16548

No. elements 3840 6256 8990 11934 15552 19040 23270 27824 32536

24

h [m] 11.88283 11.89933 11.91552 11.99243 11.92557 12.01760 11.98914 11.97479 12.00785

Nopt 44 30 24 20 17 15 14 13 12

−3

x 10

1

9

0.9

8

0.8

Relative error [%]

Relative error [%]

10

7 6 5 4 3

0.7 0.6 0.5 0.4

800

900

1000 R [m]

1100

0.3

1200

(a)

800

900

1000 R [m]

1100

1200

(b)

Figure 8: Relative errors between the numerical and the exact solution for increasing R. (a) Error in L2 -norm. (b) Error in H 1 -seminorm.

Secondly, we test the effect of the truncation order on the accuracy of the DtN FEM procedure. We assume R = 1.5a and the numerical simulations are performed in a mesh with I = 61 and J = 240, which yields a mesh size h = 7.66128 m. The truncation order N is varied from 2 to 50 in increments of 2. The relative errors in L2 -norm and H 1 -seminorm are computed as functions of N and presented in Figure 9. We observe that both of them are decreasing functions in N , just as expected. In addition, both relative errors exhibit very small values, even for a low truncation order N , with the error in L2 -norm being far smaller than the error in H 1 -seminorm. When N is increased from 2 to 50, the reduction in both relative errors is minor (from 0.00287% to 0.00281% and from 0.43285% to 0.43240%, respectively). For the considered mesh, the optimal value of N given by (69) is Nopt = 22, which has associated relative errors of 0.00282% in L2 -norm and 0.43246% in H 1 -seminorm. If N is increased from Nopt to 50, the reduction in both relative errors is even more minor, so the estimation for optimal truncation order provided by (69) is by far acceptable. Consequently, the value of the truncation order N does not have an important effect on the accuracy of the DtN FEM approach in solving the model problem. Nevertheless, it should be noted that the model problem corresponds to a case where the geometry of the perturbation is quite simple and the solution corresponds to a smooth and regular function. A perturbation with a more complex geometry, together with a solution exhibiting stronger spatial variations, is likely to require larger values of N in order to achieve a good accuracy. Finally, we test the effect of the mesh size h on the accuracy of the method and we investigate numerically its convergence as h tends to zero. For this purpose, we assume R = 1.5a and nine structured meshes of decreasing size are considered, whose parameters are given in Table 3, including the optimal N provided by (69). Numbers I and J have been set in order to obtain approximately equispaced decreasing values of h in logarithmic scale. The relative errors between the numerical and the exact solution, in L2 -norm and H 1 -seminorm, are presented in a single log-log plot in Figure 10, where we observe that both of them decrease as the mesh size h is diminished, 25

−3

2.88

x 10

0.433

2.86

Relative error [%]

Relative error [%]

0.4329

2.84

2.82

0.4328 0.4327 0.4326 0.4325 0.4324

0

10

20

30

40

0.4323

50

0

10

20

N

30

40

50

N

(a)

(b)

Figure 9: Relative errors between the numerical and the exact solution for increasing N . (a) Error in L2 -norm. (b) Error in H 1 -seminorm.

Table 3: Parameters of structured meshes for testing the effect of h.

I 3 4 6 8 12 19 28 41 61

J 10 15 22 33 49 73 108 161 240

No. nodes 44 80 161 306 650 1480 3161 6804 14942

No. elements 60 120 264 528 1176 2774 6048 13202 29280

h [m] 166.51964 117.30296 79.98915 56.25572 37.87096 24.85467 16.85528 11.40243 7.66128

Nopt 7 9 10 12 14 16 18 20 22

with the error in L2 -norm being lower than the error in H 1 -seminorm. Therefore, the numerical solution computed by the DtN FEM converges to the exact solution as the mesh is refined, which confirms the effectiveness of the method. In fact, the relative error in L2 -norm associated with the finest mesh is less than 0.003%. This very good agreement between the numerical and the exact solution is achieved with a relatively low number of terms in the series of the DtN map (N = 22). In addition, the curves of both relative errors are nearly straight lines in logarithmic scale, which suggests that the DtN FEM has a constant rate of convergence. The computed slopes are nearly 2 and 1, respectively, which predicts that the rate of convergence of the method is 2. The numerical evidence thus indicates that the DtN FEM approach presented in this work possesses second-order accuracy.

26

2

10

L2−norm H1−seminorm

1

Relative error [%]

10

0

10

−1

10

−2

10

−3

10

0

10

1

2

10

10

3

10

h [m] Figure 10: Relative errors between the numerical and the exact solution for decreasing h.

5.4. An application to a geophysical problem In what follows, we apply the DtN FEM procedure introduced throughout this work to solve numerically a simplified model of a geophysical problem. An important area in mining and geotechnical engineering is slope stability analysis, particularly in designing open-pit mines. The FEM has been commonly used for engineers to study numerically the stability of slopes, both in 2D and 3D. However, the unboundedness of the surrounding domain remains as an important drawback, which is often not appropriately dealt with. In order to show a more suitable treatment of the unboundedness in this type of problems, we propose a simplified mathematical model to analyse the static stability of a single slope. The mathematical domain is assumed to be semi-infinite and axisymmetric, which will allow us to approach to some extent the roughly axisymmetric shape of an open-pit mine. The surrounding solid medium is supposed elastic and piecewise homogeneous isotropic. Figure 11 depicts schematically the domain Ω to be considered, which is divided into two regions with different physical parameters: A layer of constant thickness around the slope (Ω1 ) and the rest of the medium (Ω0 ) which corresponds to the unbounded region. We denote H the height of the slope, d the thickness of the layer, b the horizontal distance between the axis of symmetry and the bottom of the slope, and α the inclination angle of the slope (cf. Figure 11). Let E0 , ν0 and E1 , ν1 be the Young’s moduli and Poisson’s ratios of the elastic media occupying Ω0 and Ω1 , respectively. The corresponding Lam´e constants λ0 , μ0 and λ1 , μ1 are obtained from the aforementioned standard formulae. Let also 0 and 1 be the mass densities in Ω0 and Ω1 . We assume that the whole domain is subject to the gravity force, which acts everywhere in the direction of −ˆ z . Therefore, in order to solve (5), we consider the lithostatic displacement ug (associated with the solid medium in the unbounded ˆ in (5c), exactly as indicated in Remark 1. The physical region Ω0 ) and we set f = −σ(ug )n solution is then obtained by adding ug to the numerical solution uh of (5) (see also [30]). In the numerical simulations to be presented below, we set H = 800 m, b = 400 m and d = 100 m. The assumed values for the physical constants are E0 = 20 GPa, ν0 = 0.2, 0 = 2500 kg/m3 , E1 = 10 GPa, ν1 = 0.3, 1 = 2000 kg/m3 , and an acceleration of gravity g = 9.81 m/s2 . Three values of slope inclination are considered, namely slopes of 0.75 (α = tan−1 0.75 ≈ 0.20483π), 1.25 27

d H

Ω1 b

g α Ω0

d

Figure 11: Mathematical domain in the simplified model to analyse the stability of a slope.

(α = tan−1 1.25 ≈ 0.28522π) and 1.75 (α = tan−1 1.75 ≈ 0.33475π). The artificial boundary is located at R = 3000 m. For each one of the three inclinations, unstructured triangular meshes of Ω = Ω0 ∪ Ω1 of similar size are employed, whose parameters are given in Table 4. The optimal N for each mesh is obtained from (69) with the length a set to the horizontal distance between the axis of symmetry and the top of the slope (a = b + H cot α). The resulting values of N (also indicated in Table 4) are assumed in each simulation. The displacement fields for the slopes of 0.75, 1.25 and 0.75 are respectively presented in Figures 12, 13 and 14. We observe in all considered Table 4: Parameters of unstructured meshes for each slope inclination.

Slope 0.75 1.25 1.75

No. nodes 14018 14435 14648

No. elements 27551 28393 28820

(a)

h [m] 38.02701 37.52022 38.69456

Nopt 11 7 5

(b)

Figure 12: Computed displacement fields for a slope of 0.75. (a) Component uhρ (ρ, z). (b) Component uhz (ρ, z).

28

(a)

(b)

Figure 13: Computed displacement fields for a slope of 1.25. (a) Component uhρ (ρ, z). (b) Component uhz (ρ, z).

(a)

(b)

Figure 14: Computed displacement fields for a slope of 1.75. (a) Component uhρ (ρ, z). (b) Component uhz (ρ, z).

cases that the domain deforms downwards due to the gravity force, and this deformation is larger for a lower slope. Furthermore, the upper part of the slope tends to deform horizontally towards the axis of symmetry, also with a larger deformation for a lower slope. The slope stability can be estimated by using the von Mises stress, which in axisymmetric cylindrical coordinates is computed as  2 . σV M = σρ2 + σθ2 + σz2 − σρ σθ − σθ σz − σρ σz + 3σρz (70) The von Misses stress was evaluated along the slope surface for each inclination studied. Figure 15 presents a plot of the three curves of Von Mises stress on the slope in function of the height, which remains constant in all the cases considered. We observe that for the lowest slope (0.75), a maximum value of von Mises stress arises at a height of about 135 m with respect to the bottom of the pit, whereas for the other two slopes, both of them steeper than the first one, the maximum of von Mises stress occurs at the bottom of the slope.

29

12 Slope of 0.75 Slope of 1.25 Slope of 1.75

Stress [MPa]

10 8 6 4 2 0

0

100

200

300

400 Height [m]

500

600

700

800

Figure 15: Von Mises stress on the slope surface in function of the height for different slope inclinations.

6. Concluding remarks In this paper, we have presented a DtN FEM to solve boundary-value problems of linear elastostatics in an axisymmetric semi-infinite domain. The method proceeds in a similar way to the standard version for exterior domains, however, a special semi-analytical procedure was devised to overcome the lack of an explicit closed-form expression for the DtN map. This constitutes a drawback that is inherent to elasticity problems formulated in semi-infinite domains, which commonly arise in geophysical applications. First, a semi-spherical artificial boundary was used to truncate the original semi-infinite domain, dividing it into a bounded computational domain and a semi-infinite residual domain. The basic idea of the method was to approximate numerically only those integral terms occurring in the FEM formulation, established in the computational domain, that involve the DtN map. The method of separation of variables conveniently applied in the residual domain yields a general analytical solution in series form, depending on unknown coefficients, which are to be determined in order to satisfy the Dirichlet boundary condition on the artificial boundary. As this is not possible in an analytical way, we proceeded numerically, truncating the series and establishing a suitable quadratic functional depending on the solution, whose minimisation provides a linear system of equations for a finite number of coefficients. The particular block structure of the associated matrix allows us to solve the linear system very efficiently, using in practice mere forward and backward substitutions to calculate the coefficients of the series. This procedure was applied repeatedly to approximate all the terms involving the DtN map in the FEM formulation, resulting in an effective coupling between the exact boundary conditions provided by the DtN map and the FEM scheme. The performance of the DtN FEM was studied by solving a model problem whose exact solution is available. The numerical solution was compared against the exact one in structured triangular meshes by analysing the relative errors in L2 -norm and H 1 -seminorm. The respective effects of the location of the artificial boundary, number of terms in the series of the DtN map and mesh size on the accuracy of the method were tested. Improvements of accuracy were observed as the artificial boundary was located farer away from the perturbation, however, a good accuracy was achieved

30

even for a close artificial boundary. The number of terms in the series of the DtN map did not show to have an important effect on the accuracy of the method; a good agreement between the numerical and the exact solution was achieved using a relatively low number of terms (less than 25). Decreases in both relative errors were observed as the mesh size was reduced, which demonstrates the effectiveness and accuracy of the DtN FEM procedure. Furthermore, the numerical evidence showed that the method possesses second-order accuracy. Finally, the DtN FEM was applied to solve a simplified geophysical problem, in order to illustrate the potentiality of the method. Acknowledgements This research received funding from Proyecto CYTED-RED Tem´ atica MAVS P711RT0278. The authors also thank Professor J.-C. N´ed´elec for his valuable support of this work. Appendix A. The entries of matrix Q(AA) are given by (2n + 1)(2n + 2)(4n + 3) (4n + 5)   × (2n + 3)(16n3 + 32n2 + 22n + 8ν 2 − 12ν + 9) + 4(1 − ν 2 )

= 2μQ(AA) nn

(AA)

(AA)

2μQn,n+1 = 2μQn+1,n =

(AA)

2μQnk

(2n + 1)(2n + 2)(2n + 3)(2n + 4)(4n + 3) (4n + 5)   × (2n + 2)(2n + 4) − 1 + 2ν

=0

n = 0, . . . , N, (A.1a)

n = 0, . . . , N − 1, (A.1b)

n, k = 0, . . . , N,

|n − k| > 1. (A.1c)

The entries of matrices Q(AB) and Q(BA) are given by (AB)

4(4n + 3)P2n (0) (2n − 1)(2n + 2)(2n + 4) × (2n(2n + 3)ν(ν − 2) − 3 + 5ν − 4ν 2 )

(BA)

2μQn,−1 = 2μQ−1,n = −

(AB)

2μQnk

(BA)

= 2μQkn

=

4(2k + 1)(2k + 3)(4k + 5) (2k + 6 + 2n)(2k + 3 − 2n)(2k + 4 + 2n) (2n + 1)(4n + 3) ηnk P2n (0)P2k (0) × (2k − 1 − 2n)(2k + 1 − 2n)

n = 0, . . . , N, (A.2a)

n, k = 0, . . . , N, (A.2b)

with ηnk = 51 + 58k + 32k 2 n2 + 188kn + 138n + 56k 2 n + 72n2 + 16k 2 + 104kn2     − (2k − 1 − 2n) 4n2 (2k + 4) − 4k((2k + 5)n + 2k + 6) − 21 − 8k(2k + 2) ν − (2k + 3 − 2n)(2k + 6 + 2n)(2k − 1 − 2n)ν 2 . The entries of matrix Q(BB) are given by 2 (BB) 2μQ−1,−1 = 2(1 − 2ν)2 ln 2 + (2 + 5ν − 6ν 2 ), 3

(A.3a)

31

(BB)

(BB)

2μQn,−1 = 2μQ−1,n = 2(7 + 2ν)δn0

  2(1 − 2ν)(4n + 5)(2n + 1) (2n + 4)ν − 1 P2n (0) − (2n + 2)(2n + 4)

n = 0, . . . , N, (A.3b)

(4n + 5)(2n + 2)(2n + 3) (4n + 7) 4 3 × (32n + 208n + 492n2 + 4(125 − 2ν + 4ν 2 )n + 187 − 20ν + 28ν 2 )

2μQ(BB) = nn

(BB)

(BB)

2μQn+1,n = 2μQn,n+1 =

(BB)

2μQnk

(4n + 5)(2n + 2)(2n + 3)(2n + 4)(2n + 5) (4n + 7) × ((2n + 4)(2n + 6) − 1 + 2ν)

=0

n = 0, . . . , N, (A.3c)

n = 0, . . . , N − 1, (A.3d)

n, k = 0, . . . , N,

|n − k| > 1, (A.3e)

where δnk denotes the Kronecker delta. Appendix B. j,± j,± , mj,± are given next Explicit expressions for a few coefficients ls,n s,n and qs

φj − φj−1 1 cos2 φj − cos2 φj−1 − cos φj sin φj − , 4 2 4(φj − φj−1 ) φj+1 − φj cos2 φj+1 − cos2 φj 1 + cos φj sin φj + , = 4 2 4(φj+1 − φj )

j,− = lρ,0 j,+ lρ,0

j,− = mj,− lρ,1 ρ,1 =

(2 + sin2 φj ) cos φj − (2 + sin2 φj−1 ) cos φj−1 1 sin3 φj + , 3 9(φj − φj−1 )

(2 + sin2 φj+1 ) cos φj+1 − (2 + sin2 φj ) cos φj 1 j,+ 3 sin , = mj,+ = − φ − lρ,1 j ρ,1 3 9(φj+1 − φj ) φj − φj−1 cos φj sin φj (6 cos2 φj − 7) j,− lρ,2 − =− 32 16   (cos2 φj − cos2 φj−1 ) 3(cos2 φj + cos2 φj−1 ) − 7 − , 32(φj − φj−1 ) cos φj sin φj (6 cos2 φj − 7) φj+1 − φj j,+ lρ,2 + =− 32 16   (cos2 φj+1 − cos2 φj ) 3(cos2 φj + cos2 φj+1 ) − 7 + , 32(φj+1 − φj ) φj − φj−1 cos φj − cos φj−1 + sin φj + , qρj,− = 2 φj − φj−1 φj+1 − φj cos φj+1 − cos φj − sin φj − , qρj,+ = 2 φj+1 − φj j,− =− lz,0

sin φj cos φj − sin φj−1 cos φj−1 cos2 φj − sin2 φj + , 4 4(φj − φj−1 )

sin φj+1 cos φj+1 − sin φj cos φj cos2 φj − sin2 φj − , 4 4(φj+1 − φj ) (2 + cos2 φj ) sin φj − (2 + cos2 φj−1 ) sin φj−1 1 , = − cos3 φj + 3 9(φj − φj−1 )

j,+ = lz,0 j,− lz,1

32

j,+ lz,1 =

mj,− z,1 =

(2 + cos2 φj+1 ) sin φj+1 − (2 + cos2 φj ) sin φj 1 cos3 φj − , 3 9(φj+1 − φj ) (6 + sin2 φj ) sin φj − (6 + sin2 φj−1 ) sin φj−1 (2 + sin2 φj ) cos φj − , 3 9(φj − φj−1 )

(6 + sin2 φj+1 ) sin φj+1 − (6 + sin2 φj ) sin φj (2 + sin2 φj ) cos φj + , 3 9(φj+1 − φj ) 8(cos2 φj − 1)(3 cos2 φj + 1) + 7 =− 64 cos φj sin φj (6 cos2 φj + 1) − cos φj−1 sin φj−1 (6 cos2 φj−1 + 1) + , 64(φj − φj−1 ) 8(cos2 φj − 1)(3 cos2 φj + 1) + 7 = 64 cos φj+1 sin φj+1 (6 cos2 φj+1 + 1) − cos φj sin φj (6 cos2 φj + 1) − . 64(φj+1 − φj )

mj,+ z,1 = − j,− lz,2

j,+ lz,2

References [1] D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier, Amsterdam, 1992. [2] D. Givoli, Recent advances in the DtN FE method, Arch. Comput. Method E. 6 (2) (1999) 71–116. [3] H. Han, X. N. Wu, A survey on artificial boundary method, Sci. China Math. 56 (2013) 2439–2488. [4] H. Han, X. N. Wu, Artificial Boundary Method, Tsinghua University Press and SpringerVerlag, Beijing-Berlin-Heidelberg, 2013. [5] D. Givoli, Exact representations on artificial interfaces and applications in mechanics, Appl. Mech. Rev. 52 (11) (1999) 333–349. [6] J. B. Keller, D. Givoli, Exact non-reflecting boundary conditions, J. Comput. Phys. 82 (1989) 172–192. [7] M. J. Grote, J. B. Keller, On nonreflecting boundary conditions, J. Comput. Phys. 122 (1995) 231–243. [8] A. S. Deakin, H. Rasmussen, Nonreflecting boundary condition for the Helmholtz equation, Comput. Math. Appl. 41 (2001) 307–318. [9] M. J. Grote, J. B. Keller, Exact nonreflecting boundary condition for the time dependent wave equation, SIAM J. Appl. Math. 55 (1995) 280–297. [10] M. J. Grote, J. B. Keller, Nonreflecting boundary conditions for time-dependent scattering, J. Comput. Phys. 127 (1996) 52–65. [11] U. E. Aladl, A. S. Deakin, H. Rasmussen, Nonreflecting boundary condition for the wave equation, J. Comput. Appl. Math. 138 (2002) 309–323. [12] D. Givoli, J. B. Keller, Non-reflecting boundary conditions for elastic waves, Wave Motion 12 (1990) 261–279. 33

[13] I. Harari, Z. Shohet, On non-reflecting boundary conditions in unbounded elastic solids, Comput. Methods Appl. Mech. Engrg. 163 (1998) 123—-139. [14] M. J. Grote, J. B. Keller, Exact nonreflecting boundary condition for elastic waves, SIAM J. Appl. Math. 60 (2000) 803–819. [15] M. J. Grote, Nonreflecting boundary conditions for elastodynamic scattering, J. Comput. Phys. 161 (2000) 331–353. [16] G. G¨ achter, M. J. Grote, Dirichlet-to-Neumann map for three-dimensional elastic waves, Wave Motion 37 (2003) 293–311. [17] D. P. Nicholls, N. Nigam, Exact non-reflecting boundary conditions on general domains, J. Comput. Phys. 194 (2004) 278–303. [18] D. P. Nicholls, N. Nigam, Error analysis of an enhanced DtN-FE method for exterior scattering problems, Numer. Math. 105 (2006) 267–298. [19] L. Chindelevitch, D. P. Nicholls, N. Nigam, Error analysis and preconditioning for an enhanced DtN-FE algorithm for exterior scattering problems, J. Comput. Appl. Math. 204 (2007) 493– 504. [20] H. Han, X. N. Wu, Approximation of infinite boundary condition and its application to finite elements methods, J. Comput. Math. 3 (1985) 179–192. [21] H. Han, X. N. Wu, The approximation of the exact boundary conditions at an artificial boundary for linear elastic equations and its application, Math. Comput. 59 (199) (1992) 21–37. [22] D. Givoli, J. B. Keller, A finite element method for large domains, Comput. Methods Appl. Mech. Engrg. 76 (1989) 41—-66. [23] H. Han, C. X. Zheng, Artificial boundary method for the three-dimensional exterior problem of elasticity, J. Comput. Math. 23 (2005) 603—-618. [24] H. Han, W. Bao, Error estimates for the finite element approximation of problems in unbounded domains, SIAM J. Numer. Anal. 37 (2000) 1101–1119. [25] W. Bao, H. Han, High-order local artifcial boundary conditions for problems in unbounded domains, Comput. Methods Appl. Mech. Engrg. 188 (2000) 455–471. [26] H. Han, W. Bao, Error estimates for the finite element approximation of linear elastic equations in an unbounded domain, Math. Comput. 70 (2001) 1437–1459. [27] D. Givoli, S. Vigdergauz, Artificial boundary conditions for 2D problems in geophysics, Comput. Methods Appl. Mech. Engrg. 110 (1993) 87—-101. [28] H. Han, W. Bao, T. Wang, Numerical simulation for the problem of infinite elastic foundation, Comput. Methods Appl. Mech. Engrg. 147 (1997) 369–385. [29] W. Bao, H. Han, The direct method of lines for the problem of infinite elastic foundation, Comput. Methods Appl. Mech. Engrg. 175 (1999) 157–173. [30] V. Boccardo, E. Godoy, M. Dur´an, An efficient semi-analytical method to compute displacements and stresses in an elastic half-space with a hemispherical pit, Adv. Appl. Math. Mech. 7 (3) (2015) 295–322. 34

[31] M. H. Sadd, Elasticity: Theory, Applications, and Numerics, Elsevier Butterworth-Heinemann, Boston, 2005. [32] B. Szabo, I. Babuˇska, Finite Element Analysis, John Wiley & Sons, New York, 1991. [33] R. A. Eubanks, Stress concentration due to a hemispherical pit at a free surface, J. Appl. Mech. 21 (1954) 57–62. [34] Y. A. Amenzade, Theory of Elasticity, Mir Publishers, Moscow, 1979. [35] A. Bj¨ orck, Numerical Methods for Least Square Problems, SIAM, Philadelphia, 1996. [36] G. H. Golub, C. F. V. Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, 1996.

35