PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 96 (110) (2014), 5–22 DOI: 10.2298/PIM1410005A

A NUMERICAL STUDY OF ENERGETIC BEM-FEM APPLIED TO WAVE PROPAGATION IN 2D MULTIDOMAINS

A. Aimi, L. Desiderio, M. Diligenti, and C. Guardasoni

Abstract. Starting from a recently developed energetic space-time weak formulation of boundary integral equations related to wave propagation problems defined on single and multidomains, a coupling algorithm is presented, which allows a flexible use of finite and boundary element methods as local discretization techniques, in order to efficiently treat unbounded multilayered media. Partial differential equations associated to boundary integral equations will be weakly reformulated by the energetic approach and a particular emphasis will be given to theoretical and experimental analysis of the stability of the proposed method.

1. Introduction

The coupling of finite and boundary element methods (FEM and BEM) for the solution of time-dependent problems is attractive because it allows for an optimal exploitation of the respective advantages of both methods. A mathematical survey of the coupling of FEM and BEM is given in [16, 22, 28]. The BEM, also when formulated directly in the space-time domain, has attracted particular interest for its accuracy, the simplicity of imposing the interface conditions in problems defined on multidomains, the implicit fulfillment of the infinity radiation conditions and the low cost of discretization when problems are defined over unbounded domains and the classical numerical methods as finite difference or finite element cannot efficiently determine the solution.

When one deals with regions having different material properties (e.g., layered soils [11, 25]) or even different physics (e.g., in solid-fluid coupling [17] or wavesoil-structure interaction [27]) domain decomposition is needed.

The idea of subdividing the computational domain into subregions where the most appropriate solution technique is applied, is computationally very actractive. Therefore this approach has been addressed in many publications, mainly in the context of BEM-FEM coupling. In fact the use of boundary integral equations (BIEs) and BEMs is complex and not particularly efficient in presence of nonlinearities localized in small parts of the domain. In this case, the classical differential models and numerical techniques, such as the Finite Difference Method (FDM) and FEM help efficiently deal with the nonlinear part of the problem, but require, in general a fine discretization of the entire domain with a significant increase in computational cost, even if, when using fully unstructured grids, local mesh refinement is in principle feasible, at least in absence of strong inhomogeneities. Anyway, in this context, BEM and FEM methods for the approximation of boundary integral equations systems and systems of partial differential equations are complementary and a suitable coupling of these two techniques can take advantage of what both offer.

2010 Mathematics Subject Classification: 65M38.

Key words and phrases: wave propagation, boundary element and finite element methods, 2D multidomain.

Within engineering calculations, the BEM and the FEM are well established tools for the numerical approximation of problems for which analytical solutions are mostly unknown or available under unrealistic modeling only. Applications occur in elasticity problems, elastodynamics, electromagnetics, wave scattering. In these problems we have to solve a given differential equation in two adjacent domains with specified interface conditions. In the last decades, contributions to BEM-FEM coupling, in the context of hyperbolic problems, started to appear [1, 21, 23, 29], especially analyzing stability issues. In this work, taking advantage of a recently developed energetic space-time weak formulation of BIEs related to wave propagation problems defined on single and multidomains (see in particular [3,4] and references therein), a coupling algorithm is presented, which allows a flexible use of FEM and BEM as local discretization techniques, in order to efficiently treat unbounded multilayered media. In principle, both the frequency-domain and time-domain BEM can be used for hyperbolic boundary value problems. Most earlier contributions concerned direct formulations of BEM in the frequency domain, often using the Laplace or Fourier transforms and addressing wave propagation problems. Conversely, time-domain BEM yields directly the unknown time-dependent quantities. In this last approach, the construction of the BIEs uses the fundamental solution of the hyperbolic partial differential equation and jump relations. In this direction the most interesting results are given by the weak formulation due to Bamberger and Ha Duong [13, 14].

Partial differential equations associated to BIEs will be weakly reformulated by the energetic approach and a particular emphasis will be given to theoretical and experimental analysis of the stability of the proposed method.

The paper is structured as follow: in Sections 2 and 3 the model problem and its energetic weak formulation are introduced; Section 4 is dedicated to the Galerkin BEM-FEM discretization phase; in Section 5 we obtain theoretical stability results using energy arguments. Al last, in Section 6, several numerical results are presented and discussed.

2. Model problem

Let Ω ⊂ ℝ^{2} be an open bounded domain, with a sufficiently smooth boundary ∂Ω with inward normal n. Let Ω_{1}∪Ω_{2} = ℝ^{2}∖Ω be a decomposition of ℝ^{2} ∖ Ω , with Ω_{1} unbounded and Ω_{2} bounded nonoverlapping subdomains such that Ω̄_{1} ∩ Ω̄_{2} = Γ, Ω̄ ∩ Ω̄_{1} = ∅ as depicted in Figure 1. Note that the boundary of the unbounded subdomain Ω_{1} is just the interface Γ.

Figure 1. Spatial domain for the model problem.

Having denoted with u_{i}(x, t) the unknown function in the i-th subdomain, which depends on space and time, and with p_{i}(x, t) := μ_{i}∂u_{i}(x, t)/∂n_{x}, which depends on a unitary (outward) normal vector and on μ_{i}, a typical constant related to the material constituting Ω_{i}, we want to solve the following wave propagation model problem

(2.1)  Δu_{i}(x, t) − (1/c_{i}^{2})ü_{i}(x, t) = f_{i}(x, t), x ∈ Ω_{i}, t ∈ [0, T], i = 1, 2
(2.2)  u_{i}(x, t) = 0,  x ∈ Ω_{i}, i = 1, 2
(2.3)  u̇_{i}(x, 0) = 0,  x ∈ Ω_{i}, i = 1, 2
(2.4)  p_{2}(x, t) = p̄(x, t), x ∈ Γ_{N} := ∂Ω, t ∈ [0, T],

where overhead dots indicate derivatives with respect to time, c_{i} is the propagation velocity of a perturbation in the i-th subdomain, p̄(x, t) is a given function, the assigned PDE right-hand sides f_{1}(x, t) ≡ 0 and f_{2}(x, t) are suitably connected. Moreover on the common boundary or interface Γ the matching conditions read

(2.5)  u_{1}(x, t) = u_{2}(x, t), p_{1}(x, t) = −p_{2}(x, t), x ∈ Γ, t ∈ [0, T].

where u_{i} ∈ H^{1}([0, T]; H^{1}(Ω)), i.e., the unknown functions u_{i} are understood in a weak sense.

Since the goal of this work is to approximate u_{1} using a BEM approach and u_{2} using a FEM technique, we have to obtain a boundary integral reformulation of the problem (2.1) in Ω_{1}.

As it is well known, following [19, 24], problem (2.1)–(2.3) in the subdomain Ω_{1} can be rewritten as a strong system of two BIEs in the boundary unknowns the functions p_{1}(x, t) and u_{1}(x, t) on Γ

\begin{aligned}\frac{1}{2}u_{1}(x, t) & = \frac{1}{μ_{1}} (Vp_{1})(x, t) − (Ku_{1})(x, t), \\ \frac{1}{2}p_{1}(x, t) & = (K^{∗}p_{1})(x, t) − μ_{1}(Du_{1})(x, t). \end{aligned} ; (2.6)

where

(Vp_{1})(x, t)= ∫_{Γ} ∫_{0}^{t} G(r; t − τ)p_{1}(ξ, τ)dτdγ_{ξ},
(Ku_{1})(x, t) = ∫_{Γ} ∫_{0}^{t} \frac{∂G}{∂n_{ξ}} (r; t − τ)u_{1}(ξ, τ)dτdγ_{ξ},
(K^{∗}p_{1})(x, t) = ∫_{Γ} ∫_{0}^{t} \frac{∂G}{∂n_{x}} (r; t − τ)p_{1}(ξ, τ)dτdγ_{ξ}
(Du_{1})(x, t) = ∫_{Γ} ∫_{0}^{t} \frac{∂^{2}G}{∂n_{x}∂n_{ξ}} (r; t − τ)u_{1}(ξ, τ)dτdγ_{ξ}

and

G(r; t − τ) = \frac{c_{1}}{2π} \frac{H[c_{1}(t − τ) − r]}{\sqrt{c_{1}^{2}(t − τ)^{2} − r^{2}}} ,

is the fundamental solution of the two dimensional wave operator, with H[·] the Heaviside function and r = ∥r∥_{2} = ∥x − ξ∥_{2}. Note that in (2.6) the operator K^{∗} is the adjoint of the Cauchy singular operator K.

Of course, problem (2.6) has to be coupled with the differential one specified for Ω_{2}, under coupling conditions (2.5) at the interface. In particular, we are here interested in a direct space-time weak formulation for the coupling of the integrodifferential problem on Ω_{1} ∪ Ω_{2}, and this will be done in the next Section.

3. Energetic weak formulation for the coupling

We start remarking that the solution of (2.1) in Ω_{1} satisfies the following energy identity

(3.1)  E_{Ω_{1}} (u_{1}, T) := \frac{1}{2} ∫_{Ω_{1}} [ \frac{1}{c_{1}^{2}} u̇_{1}^{2}(x, T) + |∇u_{1}(x, T)|^{2}]dxdt

= \frac{1}{μ_{1}} ∫_{Γ} ∫_{0}^{T} u̇_{1}(x, t)p_{1}(x, t)dtdγ_{x}

which can be obtained multiplying by u̇_{1} equation (2.1) specified for i = 1 and integrating by parts over Ω_{1} × [0, T]. Then, the energetic weak formulation of system (2.6) consists in:

finding u_{1} ∈ H^{1}([0, T]; H^{\frac{1}{2}} (Γ)) and p_{1} ∈ H^{0}([0, T]; H^{− \frac{1}{2}} (Γ)) such that

\begin{aligned}\frac{1}{2}⟨u̇_{1}, q_{1}⟩ & = \frac{1}{μ_{1}} ⟨ ˙ (Vp_{1}), q_{1}⟩ − ⟨ ˙ (Ku_{1}), q_{1}⟩, \\ −\frac{1}{2}⟨p_{1}, v̇_{1}⟩ & = −⟨K^{∗}p_{1}, v̇_{1}⟩ + μ_{1}⟨Du_{1}, v̇_{1}⟩, \end{aligned} ; (3.2)

where ⟨·, ·⟩ = ⟨·, ·⟩_{L^{2} ([0,T]×Γ)} and q_{1}(x, t), v_{1}(x, t) are suitable test functions, belonging to the same functional space of p_{1}(x, t) and u_{1}(x, t), respectively.

In particular, the first equation in (2.6) has been differentiated with respect to time and projected with the L^{2}([0, T] × Γ) scalar product by means of functions belonging to H^{0}([0, T]; H^{− \frac{1}{2}} (Γ)), while the second equation in (2.6) has been projected with the L^{2}([0, T] × Γ) scalar product by means of functions belonging to H^{1}([0, T]; H^{\frac{1}{2}} (Γ)), derived with respect to time.

For the energetic weak formulation in Ω_{2}, let us multiply differential equation (2.1) for the time derivative of the test function v_{2}(x, t) ∈ H^{1}([0, T]; H^{1}(Ω_{2})) and integrate by parts in space obtaining, after a multiplication by the coefficient μ_{2}

(3.3)  − μ_{2}A(u_{2}, v_{2}) + ⟨v̇_{2|_{Γ}} , p_{2|_{Γ}} ⟩ = μ_{2}F(v_{2}) − ⟨v̇_{2|_{Γ_{N}}} , p̄⟩_{L^{2} ([0,T]×Γ_{N} )} ,

where

(3.4)  A(u_{2}, v_{2}) := ∫_{0}^{T} ∫_{Ω_{2}} [∇v̇_{2}(x, t) · ∇u_{2}(x, t) + \frac{1}{c_{2}^{2}} v̇_{2}(x, t)ü_{2}(x, t)]dxdt

(3.5)  F(v_{2}) := ∫_{0}^{T} ∫_{Ω_{2}} v̇_{2}(x, t)f_{2}(x, t)dxdt .

Now, remembering interface condition (2.5) and using the further coupling condition at interface for test functions v_{1}(x, t) = v_{2|_{Γ}} (x, t), combining (3.3) with the second weak BIE in (3.2), we finally obtain the following energetic weak formulation of the coupled problem

\frac{1}{μ_{1}} ⟨ ˙ (Vp_{1}), q_{1}⟩ − ⟨ ˙ (Ku_{2|_{Γ}} ), q_{1}⟩− \frac{1}{2}⟨u̇_{2|_{Γ}} , q_{1}⟩ = 0

−\frac{1}{2}⟨p_{1}, v̇_{2|_{Γ}} ⟩ − ⟨K^{∗}p_{1}, v̇_{2|_{Γ}} ⟩ + μ_{1}⟨Du_{2|_{Γ}} , v̇_{2|_{Γ}} ⟩ − μ_{2}A(u_{2}, v_{2})
= μ_{2}F(v_{2}) − ⟨v̇_{2|_{Γ_{N}}} , p̄⟩_{L^{2} ([0,T]×Γ_{N} )} .

(3.6)

At every time instant, the unknowns are p_{1} over the interface Γ and u_{2} in Ω_{2}. Note that, integrating in the sense of distributions the scalar products in the first weak equation of (3.6), we can write also:

−\frac{1}{μ_{1}} ⟨Vp_{1}, q̇_{1}⟩ + ⟨Ku_{2|_{Γ}} , q̇_{1}⟩ + \frac{1}{2}⟨u_{2|_{Γ}} , q̇_{1}⟩ = 0
−\frac{1}{2}⟨p_{1}, v̇_{2|_{Γ}} ⟩ − ⟨K^{∗}p_{1}, v̇_{2|_{Γ}} ⟩ + μ_{1}⟨Du_{2|_{Γ}} , v̇_{2|_{Γ}} ⟩ − μ_{2}A(u_{2}, v_{2})
= μ_{2}F(v_{2}) − ⟨v̇_{2|_{Γ_{N}}} , p̄⟩_{L^{2} ([0,T]×Γ_{N} )} .

(3.7)

Let us conclude this Section with some energy considerations. At first, let us consider system (3.2) with q_{1} = p_{1} and v_{1} = u_{1} and change the sign in the second equation; then, summing up the two equation and remembering (3.1) one obtains

⟨u̇_{1}, p_{1}⟩ = \frac{1}{μ_{1}} ⟨ ˙ (Vp_{1}), p_{1}⟩ − μ_{1}⟨Du_{1}, u_{1}⟩ = μ_{1}E_{Ω_{1}} (u_{1}, T).

On the other side, considering v_{2} = u_{2} in (3.4), one gets

A(u_{2}, u_{2}) = E_{Ω_{2}} (u_{2}, T).

Following similar arguments, starting from (3.6) the following relation appears

μ_{1}E_{Ω_{1}} (u_{1}, T) + μ_{2}E_{Ω_{2}} (u_{2}, T) = −μ_{2}F(u_{2}) + ⟨u̇_{2|_{Γ_{N}}} , p̄⟩_{L^{2} ([0,T]×Γ_{N} )} ,

from which one can deduce a-priori stability estimates for regular solutions u_{1} and u_{2} bounding from above the related energies by means of the problem data [8].

4. Space-time Galerkin discretization

For time discretization we consider a uniform decomposition of the time interval [0, T] with time step Δt = T/N_{Δt} , N_{Δt} ∈ ℕ^{+} , generated by the N_{Δt} + 1 time-knots: t_{k} = kΔt, k = 0, …, N_{Δt} , and we choose piecewise constant shape function for the time approximation of p_{1} and piecewise linear shape functions for the time approximation of u_{2}, although higher degree shape functions can be used. In particular, time shape functions, for k = 0, …, N_{Δt} − 1, will be defined by

ψ̄_{k} (t) = H[t − t_{k}] − H[t − t_{k+1} ], for the approximation of p_{1},
ψ̂_{k} (t) = R(t − t_{k}) − R(t − t_{k+1} ), for the approximation of u_{2},

where R(t − t_{k}) = \frac{t−t_{k}}{Δt} H[t − t_{k}] is the ramp function.

For the space discretization we consider the bounded subdomain Ω_{2} suitably approximated by a polygonal domain and a mesh T_{Δx} = {e_{1}, …, e_{M_{Δx}} } constituted by M_{Δx} triangles, with diam(e_{i}) ⩽ Δx, e_{i} ∩ e_{j} = ∅ if i ≠ j and such that ⋃ _{i=1}^{M_{Δx}} ē_{i} = Ω̄_{2} . The restriction of the mesh T_{Δx} defines on the polygonal approximation of Γ a mesh T_{Γ,Δx} formed by M_{1} nonoverlapping segments.

The functional background compels one to choose spatially shape functions belonging to L^{2}(Γ) for the approximation of p_{1} and to C^{0}(Ω_{2}) for the approximation of u_{2}. Hence, we will choose piece-wise constant basis functions φ̄_{j} (x), j = 1, …, M_{1} related to T_{Γ,Δx} for the approximation of p_{1} over the interface and piecewise linear continuous functions φ̂_{j} (x), j = 1, …, M_{2} related to T_{Δx} for the approximation of u_{2} in Ω_{2}. The approximate solutions of the problem at hand will be expressed as

(4.1)  p_{1}(x, t) ≃ ∑_{k=0}^{N_{Δt} −1} ψ̄_{k} (t)Φ̄_{k} (x) := ∑_{k=0}^{N_{Δt} −1} ψ̄_{k} (t) ∑_{j=1}^{M_{1}} ᾱ_{kj} φ̄_{j} (x), x ∈ Γ,

(4.2)  u_{2}(x, t) ≃ ∑_{k=0}^{N_{Δt} −1} ψ̂_{k} (t)Φ̂_{k} (x) := ∑_{k=0}^{N_{Δt} −1} ψ̂_{k} (t) ∑_{j=1}^{M_{2}} α̂_{kj} φ̂_{j} (x), x ∈ Ω_{2}.

The Galerkin BEM-FEM discretization coming from energetic weak formulation (3.7) produces the linear system

𝔼α = b,

where matrix 𝔼 has a block lower triangular Toeplitz structure, since its elements depend on the difference t_{h} − t_{k} and in particular they vanish if t_{h} < t_{k}. Each block has dimension M := M_{1} + M_{2}. If we denote by 𝔼^{(ℓ)} the block obtained when t_{h} − t_{k} = ℓΔt, ℓ = 0, …, N_{Δt} − 1, the linear system can be written as

(4.3) [ \begin{matrix}𝔼^{(0)} & 0 & 0 & … & 0 \\ 𝔼^{(1)} & 𝔼^{(0)} & 0 & … & 0 \\ 𝔼^{(2)} & 𝔼^{(1)} & 𝔼^{(0)} & … & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 𝔼^{(N_{Δt} −1)} & 𝔼^{(N_{Δt} −2)} & 𝔼^{(N_{Δt} −3)} & … & 𝔼^{(0)} \end{matrix} ] [ α ^{(0)} ; α ^{(1)} ; α ^{(2)} ; ⋮ ; α ^{(N_{Δt} −1)} ] = [ b^{(0)} ; b^{(1)} ; b^{(2)} ; ⋮ ; b^{(N_{Δt} −1)} ]

where for ℓ = 0, …, N_{Δt} − 1 and j = 1, …, M:

α ^{(ℓ)} = ( α_{j}^{(ℓ)}) ,  α ^{(ℓ)} = (ᾱ_{ℓ1} , …, ᾱ_{ℓM_{1}} , α̂_{ℓ1} , …, α̂_{ℓM_{2}} )^{⊤} ,  b^{(ℓ)} = ( b_{j}^{(ℓ)}) ,

Note that each block has a symmetric 2 × 2 block substructure of the type

(4.4)  𝔼^{(ℓ)} = [ \begin{matrix}𝔼_{Γ}^{(ℓ)} & 𝔼_{Γ,F EM}^{(ℓ)} \\ 𝔼_{F EM,Γ}^{(ℓ)} & 𝔼_{F EM}^{(ℓ)} \end{matrix} ]

where we can recognize the contribution of the coupling on the interface Γ and of the pure energetic FEM inside Ω_{2}, and it generally presents a highly sparse structure due to the FEM contribution.

The solution of (4.3) is obtained with a block forward substitution, i.e., at every time instant t_{ℓ} = ℓΔt, ℓ = 0, …, N_{Δt} − 1, one computes:

z^{(ℓ)} = b^{(ℓ)} − ∑_{j=1}^{ℓ} 𝔼^{(j)} α ^{(ℓ−j)}

and then solves the reduced linear system:

(4.5)  𝔼^{(0)}α ^{(l)} = z^{(l)}.

Procedure (4.5) is a time-marching technique, where the only matrix to be inverted is the nonsingular symmetric block 𝔼^{(0)}, therefore LU factorization needs to be performed only once and saved. At each time step the solution of (4.5) requires only a forward and backward substitution phases, while all the other blocks are used to update at every time step the right-hand side. Owing to this procedure, we can construct and store only the blocks 𝔼^{(0)}, …, 𝔼^{(N_{Δt} −1)} with a considerable reduction of computational cost and memory requirement.

Now, referring to the left-hand sides of weak problem (3.7), having set Δ_{hk} = t_{h}−t_{k}, after a double analytic integration in the time variable, we have the following results

• matrix elements coming from the discretization of ⟨Vp_{1}, q̇_{1}⟩ are of the form

(4.6) ∑_{α,β=0}^{1} (−1)^{α+β} ∫_{Γ} φ̄_{i} (x) ∫_{Γ} H[c_{1}Δ_{h+α,k+β} − r]S_{V}(r; Δ_{h+α,k+β} )φ̄_{j} (ξ)dγ_{ξ}dγ_{x}

where

S_{V}(r; Δ_{h+α,k+β} ) = \frac{1}{2π} [ log( c_{1}Δ_{h+α,k+β} + \sqrt{c_{1}^{2}Δ_{h+α,k+β}^{2} − r^{2}} ) − logr]

• matrix elements coming from the discretization of ⟨Ku_{2|_{Γ}} , q̇_{1}⟩ are of the form

(4.7) ∑_{α,β=0}^{1} (−1)^{α+β} ∫_{Γ} φ̄_{i} (x)∫_{Γ}H[c_{1}Δ_{h+α,k+β} − r]S_{K}(r; Δ_{h+α,k+β} )φ̂_{j|_{Γ}} (ξ)dγ_{ξ}dγ_{x}

where

S_{K}(r; Δ_{h+α,k+β} ) = \frac{1}{2πc_{1}Δt} \frac{r · n_{ξ}}{r^{2}} \sqrt{c_{1}^{2}Δ_{h+α,k+β}^{2} − r^{2}}

• matrix elements coming from the discretization of ⟨Du_{2|_{Γ}} , v̇_{2|_{Γ}} ⟩ are of the form

(4.8) ∑_{α,β=0}^{1} (−1)^{α+β} ∫_{Γ} φ̂_{i|_{Γ}} (x)∫_{Γ}H[c_{1}Δ_{h+α,k+β} − r]S_{D}(r; Δ_{h+α,k+β} )φ̂_{j|_{Γ}} (ξ)dγ_{ξ}dγ_{x}

where

S_{D}(r; Δ_{h+α,k+β} ) := \frac{1}{2πΔt^{2}} \{ \frac{(r · n_{x})(r · n_{ξ})}{r^{2}}\frac{Δ_{h+α,k+β} \sqrt{c_{1}^{2}Δ_{h+α,k+β}^{2} − r^{2}}}{c_{1}r^{2}}

+\frac{n_{x} · n_{ξ}}{2c_{1}^{2}}[ log( c_{1}Δ_{h+α,k+β} + \sqrt{c_{1}^{2}Δ_{h+α,k+β}^{2} − r^{2}} )

− logr − \frac{c_{1}Δ_{h+α,k+β} \sqrt{c_{1}^{2}Δ_{h+α,k+β}^{2} − r^{2}}}{r^{2}} ]\} .

We observe in (4.6)–(4.8) a space weak singularity of type O(logr), a space strong singularity of type O(1/r) and a space hypersingularity of type O(1/r^{2}) as r → 0, which are typical of integral kernels related to 2D elliptic problems, and also the presence of the Heaviside function which analytically represents the wave front propagation. Hence, the numerical treatment of space strong singularity and hypersingularity have been operated through quadrature schemes widely used in the context of Galerkin BEM coming from elliptic problems (see [9]), coupled with a suitable regularization technique, after a careful subdivision of the integration domain due to the presence of the Heaviside function. We refer the interested reader to [10] for the description of the accurate and fast evaluation of such type of double integrals on the interface Γ.

For what concerns matrix elements coming from < u_{2|_{Γ}} , q̇_{1} > after an analytic integration in time, one has

δ_{hk} ∫_{Γ} φ̂_{j|_{Γ}} (x)φ̄_{i} (x)dγ_{x}

where δ_{hk} is the Kronecker symbol. An analogous result holds for the discretization of ⟨p_{1}, v̇_{2|_{Γ}} ⟩.

Finally, the matrix elements coming from A(u_{2}, v_{2}) will be of the type

ρ_{hk} ∫_{Ω_{2}} ∇φ̂_{j} (x) · ∇φ̂_{i} (x)dx + \frac{χ_{hk}}{(c_{2}Δt)^{2}} ∫_{Ω_{2}} φ̂_{j} (x)φ̂_{i} (x)dx,

where

ρ_{hk} = \{ \begin{matrix}1/2, & h = k \\ 1, & h > k \\ 0, & h < k \end{matrix} and χ_{hk} = \{ \begin{matrix}1, & h = k \\ −1, & h = k + 1 \\ 0, & elsewhere \end{matrix}

hence they involve the evaluation of standard stiffness and mass matrices elements related to the bounded subdomain Ω_{2}.

5. Energy estimates for the numerical scheme

Let us consider the second energetic weak equation in (3.7), which involves all the problem unknowns, and, to simplify the following notation, p̄ = 0. Now, using (4.1) and (4.2) as approximate solutions and v_{2}(x, t) = ψ̂_{h} (t)Φ̂_{h} (x) as a test function, remembering the second equation in (2.6), we obtain

μ_{2} ∑_{k=0}^{N_{Δt} −1} [(∇Φ̂_{h} , ∇Φ̂_{k} )_{L^{2} (Ω_{2} )} ∫_{0}^{T} ˙ ψ̂_{h}(t)ψ̂_{k} (t)dt + \frac{(Φ̂_{h} , Φ̂_{k} )_{L^{2} (Ω_{2} )}}{c_{2}^{2}} ∫_{0}^{T} ˙ ψ̂_{h}(t) ¨ψ̂_{k}(t)dt]

+ ∑_{k=0}^{N_{Δt} −1} (Φ̂_{h|_{Γ}} , Φ̄_{k} )_{L^{2} (Γ)} ∫_{0}^{T} ˙ ψ̂_{h}(t)ψ̄_{k} (t)dt = μ_{2}(Φ̂_{h} , F_{h})_{L^{2} (Ω_{2} )} , (5.1)

where

F_{h}(x) = − ∫_{0}^{T} ˙ ψ̂_{h}(t)f_{2}(x, t)dt .

Performing analytically time integrals in the left-hand side of (5.1), we get

μ_{2} ∑_{k=0}^{h−1} (∇Φ̂_{h} , ∇Φ̂_{k} )_{L^{2} (Ω_{2} )} + \frac{μ_{2}}{2} (∇Φ̂_{h} , ∇Φ̂_{h} )_{L^{2} (Ω_{2} )} + \frac{μ_{2}}{c_{2}^{2}} ( Φ̂_{h} , \frac{Φ̂_{h} − Φ̂_{h−1}}{(Δt)^{2}} )_{L^{2} (Ω_{2} )}

(5.2)  + (Φ̂_{h|_{Γ}} , Φ̄_{h} )_{L^{2} (Γ)} = μ_{2}(Φ̂_{h} , F_{h})_{L^{2} (Ω_{2} )} .

At this stage we need to observe that, from (4.1) and (4.2),

p_{1}(x, t_{h}) ≃ p_{1}^{h} := Φ̄_{h} (x), u_{2}(x, t_{h}) ≃ u_{2}^{h} := ∑_{k=0}^{h−1} Φ̂_{k} (x),

hence we can rewrite the numerical scheme (5.2) as

μ_{2}( ∇(u_{2}^{h+1} − u_{2}^{h}), ∇\frac{u_{2}^{h+1} + u_{2}^{h}}{2} )_{L^{2} (Ω_{2} )} + \frac{μ_{2}}{c_{2}^{2}} ( u_{2}^{h+1} − u_{2}^{h}, \frac{u_{2}^{h+1} − 2u_{2}^{h} + u_{2}^{h−1}}{(Δt)^{2}} )_{L^{2} (Ω_{2} )}

(5.3)  + (u_{2|_{Γ}}^{h+1} − u_{2|_{Γ}}^{h} , p_{1}^{h})_{L^{2} (Γ)} = μ_{2}(u_{2}^{h+1} − u_{2}^{h}, F_{h})_{L^{2} (Ω_{2} )} .

Now, let us introduce the discrete energy at time instant t_{h+1} in the subdomain Ω_{i}, i = 1, 2, defined as

(5.4)  E_{i}^{h+1} := \frac{1}{2}[ \frac{1}{c_{i}^{2}} ∥ \frac{u_{i}^{h+1} − u_{i}^{h}}{Δt} ∥_{L^{2} (Ω_{i} )}^{2} + ∥∇u_{i}^{h+1} ∥_{L^{2} (Ω_{i} )}^{2} ];

further we observe that

(∇(u_{2}^{h+1} − u_{2}^{h}), ∇(u_{2}^{h+1} + u_{2}^{h}))_{L^{2} (Ω_{2} )} = ∥∇u_{2}^{h+1} ∥_{L^{2} (Ω_{2} )}^{2} − ∥∇u_{2}^{h}∥_{L^{2} (Ω_{2} )}^{2}

and

(u_{2}^{h+1} − u_{2}^{h}, u_{2}^{h+1} − 2u_{2}^{h} + u_{2}^{h−1} )_{L^{2} (Ω_{2} )}

= ∥u_{2}^{h+1} − u_{2}^{h}∥_{L^{2} (Ω_{2} )}^{2} − (u_{2}^{h+1} − u_{2}^{h}, u_{2}^{h} − u_{2}^{h−1} )_{L^{2} (Ω_{2} )}

⩾ ∥u_{2}^{h+1} − u_{2}^{h}∥_{L^{2} (Ω_{2} )}^{2} − \frac{1}{2}∥u_{2}^{h+1} − u_{2}^{h}∥_{L^{2} (Ω_{2} )}^{2} − \frac{1}{2}∥u_{2}^{h} − u_{2}^{h−1} ∥_{L^{2} (Ω_{2} )}^{2}

= \frac{1}{2}[∥u_{2}^{h+1} − u_{2}^{h}∥_{L^{2} (Ω_{2} )}^{2} − ∥u_{2}^{h} − u_{2}^{h−1} ∥_{L^{2} (Ω_{2} )}^{2} ].

Combining these results with (5.3) and (5.4), we finally obtain

(5.5)  μ_{2}[E_{2}^{h+1} − E_{2}^{h}] + (u_{2|_{Γ}}^{h+1} − u_{2|_{Γ}}^{h} , p_{1}^{h})_{L^{2} (Γ)} ⩽ μ_{2}(u_{2}^{h+1} − u_{2}^{h}, F_{h})_{L^{2} (Ω_{2} )} .

Applying the Cauchy–Schwarz inequality to the right-hand side of (5.5), we deduce

μ_{2}[E_{2}^{h+1} − E_{2}^{h}] + (u_{2|_{Γ}}^{h+1} − u_{2|_{Γ}}^{h} , p_{1}^{h})_{L^{2} (Γ)} ⩽ μ_{2}Δt∥ \frac{u_{2}^{h+1} − u_{2}^{h}}{c_{2}Δt} ∥_{L^{2} (Ω_{2} )} c_{2}∥F_{h}∥_{L^{2} (Ω_{2} )} ,

wherefrom, applying the Young inequality and remembering (5.4), we obtain

(5.6) μ_{2}[E_{2}^{h+1} − E_{2}^{h}] + (u_{2|_{Γ}}^{h+1} − u_{2|_{Γ}}^{h} , p_{1}^{h})_{L^{2} (Γ)} ⩽ μ_{2}ΔtE_{2}^{h+1} + μ_{2}c_{2}^{2} \frac{Δt}{2} ∥F_{h}∥_{L^{2} (Ω_{2} )}^{2} .

Let us now sum inequality (5.6), for h = 0, …, n − 1, with 1 ⩽ n ⩽ N_{Δt} , observing that E^{0} = 0

(5.7) μ_{2}E_{2}^{n}+Δt ∑_{h=0}^{n−1} ( \frac{u_{2|_{Γ}}^{h+1} − u_{2|_{Γ}}^{h}}{Δt} , p_{1}^{h}) _{L^{2} (Γ)} ⩽ μ_{2}Δt[E_{2}^{n}+∑_{h=0}^{n−1} E_{2}^{h}+ \frac{c_{2}^{2}}{2} ∥F_{h}∥_{L^{2} (Ω_{2} )}^{2} ].

Let us note that, remembering the coupling conditions on the interface, for the second term in the left-hand side of (5.7) it holds, for sufficiently small Δt,

Δt ∑_{h=0}^{n−1} ( \frac{u_{1}^{h+1} − u_{1}^{h}}{Δt} , p_{1}^{h})_{L^{2} (Γ)} ≃ ∫_{Γ} ∫_{0}^{t_{n}} u̇_{1}(x, t)p_{1}(x, t) dt dγ_{x} = μ_{1}E_{Ω_{1}} (u_{1}, t_{n}),

hence we can write

(5.8)  Δt ∑_{h=0}^{n−1} ( \frac{u_{1}^{h+1} − u_{1}^{h}}{Δt} , p_{1}^{h})_{L^{2} (Γ)} = μ_{1}E_{1}^{n} ⩾ 0

and finally from (5.7) and (5.8)

(5.9)  μ_{2}(1 − Δt)E_{2}^{n} + μ_{1}E_{1}^{n} ⩽ μ_{2}Δt[ ∑_{h=0}^{n−1} E_{2}^{h} + \frac{c_{2}^{2}}{2} ∑_{h=0}^{n−1} ∥F_{h}∥_{L^{2} (Ω_{2} )}^{2} ].

For the discrete energy in the subdomain Ω_{2} we obtain, under the natural assumption 0 < Δt < 1

E_{2}^{n} ⩽ \frac{Δt}{1 − Δt} [ ∑_{h=0}^{n−1} E_{2}^{h} + \frac{c_{2}^{2}}{2} ∑_{h=0}^{n−1} ∥F_{h}∥_{L^{2} (Ω_{2} )}^{2} ]

and applying at last Gronwall’s discrete lemma we deduce, for every time instant t_{n}, the following upper bound

(5.10)  E_{2}^{n} ⩽ \frac{c_{2}^{2}Δt}{2(1 − Δt)} exp( \frac{t_{n}}{1 − Δt} ) ∑_{h=0}^{n−1} ∥F_{h}∥_{L^{2} (Ω_{2} )}^{2} .

If we finally consider the discrete energy in the subdomain Ω_{1}, from (5.9) we can write

(5.11)  E_{1}^{n} ⩽ \frac{μ_{2}}{μ_{1}} Δt[ ∑_{h=0}^{n−1} E_{2}^{h} + \frac{c_{2}^{2}}{2} ∑_{h=0}^{n−1} ∥F_{h}∥_{L^{2} (Ω_{2} )}^{2} ]

and using (5.10) to bound from above every term in the first sum on the right-hand side of (5.11), we succeed in proving a complete stability result for our numerical scheme.

6. Numerical results

We illustrate first the performance of the proposed method with a classical benchmark taken from [2]. We consider a circular hole of radius R_{C} = 1 in an infinite linear and homogeneous membrane. The Ω_{2} region is defined by R_{C} < r < R_{I}, with R_{I} = 2, the Ω_{1} region is defined by r > R_{I} and the interface Γ is the circumference r = R_{I} (see Figure 2). Here μ_{1} = μ_{2} = 1 and c_{1} = c_{2} = 1. Meshes in Ω_{2} and over Γ are conformal, i.e., segments for Γ coincide with sides of triangles for Ω_{2} (see Figure 3). The hole is subjected to an uniform traction p̄(x, t) = e^{−t} . The observation time interval is [0, 20]. For the discretization, we have used different temporal and spatial steps, Δt and Δx respectively.

Figure 2. The problem domain

The model describes an explosive phenomenon in which the solution is the highest in the first instants of time and then rapidly decreases almost to exhaustion.

Table 1. Solution u(A, t) with c = 1 and \overline{p}(x, t) = e^{−t}
Δx Δt  u(A, 0.4)  u(A, 0.8)  u(A, 1.2)  u(A, 1.6)  u(A, 2.0)
0.4 0.4 3.880 · 10^{−1} 4.338 · 10^{−1} 5.214 · 10^{−1} 4.922 · 10^{−1} 4.885 · 10^{−1}
0.2 3.197 · 10^{−1} 4.502 · 10^{−1} 5.049 · 10^{−1} 5.171 · 10^{−1} 4.915 · 10^{−1}
0.1 3.099 · 10^{−1} 4.460 · 10^{−1} 5.029 · 10^{−1} 5.267 · 10^{−1} 4.996 · 10^{−1}
0.2 0.2 3.064 · 10^{−1} 4.481 · 10^{−1} 5.054 · 10^{−1} 5.129 · 10^{−1} 4.887 · 10^{−1}
0.1 3.088 · 10^{−1} 4.495 · 10^{−1} 5.085 · 10^{−1} 5.207 · 10^{−1} 4.974 · 10^{−1}
0.1 0.1 3.055 · 10^{−1} 4.498 · 10^{−1} 5.092 · 10^{−1} 5.190 · 10^{−1} 4.967 · 10^{−1}
0.05 0.05 3.027 · 10^{−1} 4.492 · 10^{−1} 5.100 · 10^{−1} 5.213 · 10^{−1} 5.013 · 10^{−1}
Table 2. Solution u(A, t) with c = 1 and \overline{p}(x, t) = e^{−t}
Δx Δt  u(A, 4)  u(A, 8)  u(A, 12)  u(A, 16)  u(A, 20)
0.4 0.4 3.365 · 10^{−1} 1.552 · 10^{−1} 9.438 · 10^{−1} 6.798 · 10^{−1} 5.322 · 10^{−1}
0.2 3.385 · 10^{−1} 1.542 · 10^{−1} 9.433 · 10^{−1} 6.798 · 10^{−1} 5.321 · 10^{−1}
0.1 3.376 · 10^{−1} 1.536 · 10^{−1} 9.430 · 10^{−1} 6.799 · 10^{−1} 5.322 · 10^{−1}
0.2 0.2 3.396 · 10^{−1} 1.549 · 10^{−1} 9.483 · 10^{−2} 6.834 · 10^{−2} 5.348 · 10^{−2}
0.1 3.387 · 10^{−1} 1.544 · 10^{−1} 9.480 · 10^{−1} 6.834 · 10^{−1} 5.349 · 10^{−2}
0.1 0.1 3.389 · 10^{−1} 1.546 · 10^{−1} 9.492 · 10^{−2} 6.843 · 10^{−2} 5.356 · 10^{−2}
Table 3. Solution u(B, t) with c = 1 and \overline{p}(x, t) = e^{−t}
Δx Δt  u(B, 0.4)  u(B, 0.8)  u(B, 1.2)  u(B, 1.6)  u(B, 2.0)
0.4 0.4 6.984 · 10^{−3} 4.676 · 10^{−2} 1.389 · 10^{−1} 2.500 · 10^{−1} 3.252 · 10^{−1}
0.2 1.128 · 10^{−3} 2.664 · 10^{−2} 1.355 · 10^{−1} 2.595 · 10^{−1} 3.363 · 10^{−1}
0.1 9.125 · 10^{−4} 9.831 · 10^{−3} 1.323 · 10^{−1} 2.658 · 10^{−1} 3.424 · 10^{−1}
0.2 0.2 1.296 · 10^{−3} 2.818 · 10^{−2} 1.327 · 10^{−1} 2.600 · 10^{−1} 3.388 · 10^{−1}
0.1 3.320 · 10^{−5} 1.244 · 10^{−2} 1.277 · 10^{−1} 2.673 · 10^{−1} 3.460 · 10^{−1}
0.1 0.1 7.421 · 10^{−5} 1.363 · 10^{−2} 1.258 · 10^{−1} 2.679 · 10^{−1} 3.472 · 10^{−1}
0.05 0.05 4.748 · 10^{−7} 4.794 · 10^{−3} 1.214 · 10^{−1} 2.732 · 10^{−1} 3.515 · 10^{−1}
Table 4. Solution u(B, t) with c = 1 and \overline{p}(x, t) = e^{−t}
Δx Δt  u(B, 4)  u(B, 8)  u(B, 12)  u(B, 16)  u(B, 20)
0.4 0.4 3.232 · 10^{−1} 1.584 · 10^{−1} 9.525 · 10^{−2} 6.828 · 10^{−2} 5.335 · 10^{−2}
0.2 3.261 · 10^{−1} 1.569 · 10^{−1} 9.510 · 10^{−2} 6.827 · 10^{−2} 5.335 · 10^{−2}
0.1 3.273 · 10^{−1} 1.562 · 10^{−1} 9.506 · 10^{−2} 6.827 · 10^{−2} 5.335 · 10^{−2}
0.2 0.2 3.270 · 10^{−1} 1.576 · 10^{−1} 9.559 · 10^{−2} 6.862 · 10^{−2} 5.362 · 10^{−2}
0.1 3.270 · 10^{−1} 1.576 · 10^{−1} 9.559 · 10^{−2} 6.862 · 10^{−2} 5.362 · 10^{−2}
0.1 0.1 3.287 · 10^{−1} 1.571 · 10^{−1} 9.566 · 10^{−2} 6.870 · 10^{−2} 5.369 · 10^{−2}

In Tables 1–4, the values of the numerical solution u(x, t) at points A and B obtained by varying the mesh (see Figure 3) of the domain and of the time step are Figure 3. Triangulation of the domain Ω_{2} with Δx = 0.4 (left) and its refinement (right)

reported. The results show both the efficiency of the numerical scheme proposed and the stability of the numerical solution.

In Figure 4 the solution u(x, t) on the segment \overline{AB} for different time instants is represented. The curves, shown in Figure 5, reproduce the numerical solution of the problem in the points A and B. They tend to coincide with increasing time.

Figure 4. Solution u(x, t) on the segment \overline{AB} for different time instants (c = 1 and \overline{p}(x, t) = e^{−t} )

When the domain Ω_{2}∪Ω_{1} is not homogeneous, the interface Γ is no longer fictitious and separates two regions with different characteristics, and different speeds.

Now, we consider again a circular crown Ω_{2} and an unbounded region Ω_{1} complementary to it (see Figure 2). Suppose then that the characteristic speeds in Figure 5. Graph of u(x, t) at points A and B (c = 1 and \overline{p}(x, t) = e^{−t} )

Figure 6. Graph of u(x, t) at the point A for different values of speed c_{i} and \overline{p}(x, t) = H[t]

Ω_{1} and Ω_{2} are one double of the other, in particular c_{1} = 1, c_{2} = 2 and subsequently c_{1} = 2, c_{2} = 1. We will refer to numerical solution obtained with a triangulation built with Δx = 0.2 and a time step Δt = 0.2. In this example we take on Γ the Neumann datum \overline{p}(x, t) = H[t].

From Figure 6 we see that, when the perturbation is localized in Ω_{2}, the solution u(A, t) behaves as in the case of a monodomain with the same physical characteristics of Ω_{2}. When the wave crosses the interface Γ diffraction and reflection phenomena occur, as a result of which the value of the solution is no longer Figure 7. Graph of u(x, t) at the point A for different values of speed c_{i} and \overline{p}(x, t) = H[t]

Figure 8. Unbounded domain with two holes

assimilated to one of the two monodomain. For growing time, the solution u(A, t) tends to take on the behavior that would have in a monodomain with the same physical characteristics Ω_{1}.

Remark. When the two subdomains have the same physical properties, the BIEs used on the interface in the coupling system can be seen as an exact NRBC [18] assigned on an arbitrary truncation boundary chosen in an unbounded domain in order to deal with standard domain methods (in this case finite elements). Hence, the interface acts as an absorbing layer preventing waves from being reflected backwards. In this framework, let us consider the unbounded plain domain depicted in Figure 9. Solution in the domain Ω_{2} in different time instants

Figure 10. Solution in the domain Ω_{2} in different time instants

figure 8, exterior to two holes. The truncation line is chosen as the circumference of radius 3 (dotted line). In figures 9 and 10 we report the computed solution inside the domain Ω_{2}.

References

1. T. Abboud, P. Joly, J. Rodriguez, I. Terrasse, Coupling discontinuous Galerkin methods and retarded potentials for transient wave propagation on unbounded domains, J. Comput. Physics 230(15) (2011), 5877–5907.

2. A. I. Abreu, J. A. M. Carrer, W. J. Mansur, Scalar wave propagation in 2D: a BEM formulation based on the operational quadrature method, Eng. Anal. Bound. Elem. 27 (2003), 101–105.

3. A. Aimi, M. Diligenti, A new spacetime energetic formulation for wave propagation analysis in layered media by BEMs, Int. J. Numer. Methods Eng. 75 (2008), 1102–1132.

4. A. Aimi, S. Gazzola, C. Guardasoni, Energetic Boundary Element Method analysis of wave propagation in 2D multilayered media, Math. Methods Appl. Sci. 35 (2012), 1140–1160.

5. A. Aimi, M. Diligenti, A. Frangi, C. Guardasoni, A stable 3D energetic Galerkin BEM approach for wave propagation interior problems, Eng. Anal. Bound. Elem. 36 (2012), 1756– 1765.

6. , Neumann exterior wave propagation problems: computational aspects of 3D energetic Galerkin BEM, Comp. Mech. 51 (2013), 475–493.

7. A. Aimi, M. Diligenti, C. Guardasoni, S. Panizzi, Energetic BEM-FEM coupling for wave propagation in layered media, Commun. Appl. Ind. Math. 3(2) (2012), 1–21.

8. A. Aimi, S. Panizzi, BEM-FEM coupling for the 1D Klein-Gordon equation, Numer. Methods Partial Differ. Equations, to appear; DOI:10:1002/num.21888 (2014).

9. A. Aimi, M. Diligenti, G. Monegato, New numerical integration schemes for applications of Galerkin BEM to 2D problems, Int. J. Numer. Methods Eng. 40 (1997), 1977–1999.

10. A. Aimi, M. Diligenti, C. Guardasoni, Numerical schemes for space-time hypersingular integrals in energetic Galerkin BEM, Numer. Alg. 55 (2010), 145–170.

11. V. S. Almeida, J. B. Paiva, Static analysis of soil/pile interaction in layered soil by BEM/BEM coupling, Adv. Eng. Softw. 38 (2007), 835–845.

12. H. Antes, G. Beer, W. Moser, Soil-structure interaction and wave propagation problems in 2D by a Duhamel integral based approach and the convolution quadrature method, Comput. Mech. 36(6) (2005), 431–443.

13. A. Bamberger, T. Ha Duong, Formulation variationnelle espace-temps pour le calcul par potentiel retardé de la diffraction d’une onde acoustique, I, Math. Methods Appl. Sci. 8 (1986), 405–435.

14. , Formulation variationnelle pour le calcul de la diffraction d’une onde acoustique par une surface rigide, Math. Methods Appl. Sci. 8 (1986), 598–608.

15. A. Buffa, M. Costabel, C. Schwab, Boundary element methods for Maxwell’s equations on non-smooth domains, Numer. Math. 92(4) (2002), 679–710.

16. M. Costabel, E. Stephan, Coupling of finite and boundary element methods for an elastoplastic interface problem, SIAM J. Numer. Anal. 27(5) (1990), 1212–1226.

17. O. Czygan, O. von Estorff, Fluid-structure interaction by coupling BEM and nonlinear FEM, Eng. Anal. Bound. Elem. 26 (2002), 773–779.

18. S. Falletta, G. Monegato, An exact non reflecting boundary condition for 2D time-dependent wave equations problems, Wave Motion 51 (2014), 168–192.

19. A. Frangi, Casual shape functions in the time domain boundary element method, Comp. Mech. 25 (2000), 533–541.

20. T. Ha Duong, B. Ludwig, I. Terrasse, A Galerkin BEM for transient acoustic scattering by an absorbing obstacle, Int. J. Num. Methods Eng. 57 (2003), 1845–1882.

21. S. T. Lie, G. Y. Yu, Stability improvement to BEM/FEM coupling scheme for 2D scalar wave problems, Adv. Eng. Softw. 33 (2002), 17–26.

22. C. Johnson, J. C. Nedelec, On the coupling of boundary integral and finite element methods, Math. Comp. 35(152) (1980), 1063–1079.

23. W. J. Mansur, G. Yu, J. A. M. Carrer, S. T. Lie, E. F. N. Siqueira, The θ scheme for timedomain BEM/FEM coupling applied to the 2D scalar wave equation, Comm. Numer. Methods Eng. 16(6) (2000), 439–448.

24. G. Maier, M. Diligenti, A. Carini, A variational approach to boundary element elastodynamic analysis and extension to multidomain problems, Comput. Methods Appl. Mech. Eng. 92 (1991), 193–213.

25. M. Sari, I. Demir, Wave modelling through layered media using the BEM, J. Appl. Sci. 6(8) (2006), 1703–1711.

26. E. Stephan, Coupling of boundary element methods and finite element methods; in: E. Stein, R. de Borst, T. Hughes (Eds.), Encyclopedia of Computational Mechanics 1, Wiley, (2004), 375–412.

27. J. L. Wegner, M. M. Yao, X. Zhang, Dynamic wave-soil-structure interaction analysis in the time domain, Computer and Structures 83 (2005), 2206–2214.

28. W. L. Wendland, On the coupling of finite elements and boundary elements; in: G. Kuhn, H. Mang (eds.), Discretization Methods in Structural Mechanics, IUTAM-IACM Symposium 1989, Springer-Verlag, Berlin, (1990), 405–414.

29. G. Y. Yu, A symmetric BEM/FEM coupling procedure for 2D elastodynamic problems, J. Appl. Mech. 70 (2003), 451–454.

Department of Mathematics and Computer Science University of Parma Parma Italy

alessandra.aimi@unipr.it
mauro.diligenti@unipr.it
chiara.guardasoni@unipr.it

POEMS (UMR 7231, CNRS-ENSTA-INRIA), ENSTA Paris France

Luca.Desiderio@ensta-paristech.fr