Skip to main content
  • Research Article
  • Open access
  • Published:

The double absorbing boundary method for elastodynamics in homogeneous and layered media

Abstract

Background

Recently the Double Absorbing Boundary (DAB) method was introduced as a new approach for solving wave problems in unbounded domains. It has common features to each of two types of existing techniques: local high-order Absorbing Boundary Conditions (ABC) and Perfectly Matched Layers (PML). However, it is different from both and enjoys relative advantages with respect to both.

Methods

The DAB method is based on truncating the unbounded domain to produce a finite computational domain, and on applying a local high-order ABC on two parallel artificial boundaries, which are a small distance apart, and thus form a thin non-reflecting layer. Auxiliary variables are defined on the two boundaries and within the layer, and participate in the numerical scheme. In previous studies DAB was developed for acoustic waves which are solutions to the scalar wave equation. Here the approach is extended to time-dependent elastic waves in homogeneous and layered media. The equations are written in second-order form in space and time. Standard Finite Elements (FE) are used for space discretization and the damped Newmark scheme is used for time discretization.

Results

The performance of the scheme is demonstrated via numerical examples. The DAB was applied to elastodynamics problems in conjunction with the FE method to demonstrate the performance of the method.

Conclusions

DAB is a viable method for solving wave problems in unbounded domains.

Background

Even though many computational schemes for treating wave problems in unbounded domains have been proposed during the last four decades, the quest is still going on for a high fidelity stable scheme that will allow the numerical solution in a finite region of interest. It turns out that in some cases it is very difficult to find an absorbing boundary scheme, as it is often called, that is at the same time stable, sufficiently accurate, computationally efficient, robust, and can be employed in conjunction with standard interior computational schemes. This is especially true for some specific types of problems that remain a challenge in this context. Elastodynamics, which is dealt with in this paper, is one such problem.

In recent years, the two most prominent absorbing-boundary type schemes have been those based on the use of a high-order Absorbing Boundary Condition (ABC) and those based on the use of a Perfectly Matched Layer (PML). Research on these methods remains very active. A search in the article archive ISI shows that during the last 5 years, more than 200 papers were published with the words ABC or PML in the title, and almost 1500 published papers indicated ABC or PML as keywords. See the review papers [1-4].

PML was originally devised by Bérenger [5] in 1994 for electromagnetic waves, and since then it has been further developed, analyzed and used in various applications by many authors. See, e.g., references in the review paper [4]. The PML is a layer adjacent to the boundary that truncates the unbounded domain, in which the governing equations are artificially modified. It possesses two properties at the continuous level: (a) there is a perfect match between the layer and the interior domain, namely any outgoing plane wave produces zero reflection; and (b) the solution decays exponentially when it travels inside the layer. These two properties theoretically guarantee excellent performance of the PML. What may sometimes hamper this theoretical performance is the sensitivity of the PML to the discretization and the need to introduce ad-hoc damping and stretching profiles.

The first high-order ABC was devised by Collino [6] in 1993, and a few other formulations followed by other authors. See, e.g., references in the review paper [2]. High-order ABCs are local in space and time, like the classical ABCs of Engquist and Majda [7] and Bayliss and Turkel [8], but unlike those, they do not involve high-order derivatives. Therefore they can be implemented in practice up to any desired order, as opposed to the classical ones that have been implemented up to second order only. In the high-order ABC scheme, the order of the ABC is simply an input parameter. The high derivatives that initially appear when designing a high-order ABC are eliminated by introducing auxiliary variables ϕ j on the boundary.

Recently the two approaches — high-order ABC and PML — have been compared, theoretically and numerically [3,9], in the frequency domain. They were found to be equally effective, with some relative advantages for both. In fact, although usually derived by very different approaches, recent work has shown that, at the discrete level, the two methods are quite closely related. In particular, it is shown in [10] how to design a nonstandard PML with a purely imaginary mesh continuation to exactly annihilate propagating waves at any incidence angle. This nonstandard PML is formally equivalent to the high-order ABC proposed by Hagstrom and Warburton [11].

As mentioned above, each of the two classes of techniques has relative advantages. One major disadvantage of high-order ABCs is that they require special treatment at corners formed by the intersection of two flat segments of the artificial boundary, and in some cases also at corners between an artificial and a physical boundary. Such special treatment is sometimes cumbersome or even difficult to devise. In contrast, handling corners with PMLs is usually straight forward. Another disadvantage of high-order ABCs is that they are constrained not to include any normal derivative of an auxiliary variable ϕ j on the boundary, since the ϕ j are discretized in practice only on the boundary. Thus, the ABC is allowed to involve only tangential and temporal derivatives of the ϕ j . Eliminating the normal derivatives from the ABC operators is sometimes difficult and may require a lot of algebra; a case in point is elastodynamics [12]. PML is also usually easier to incorporate in an existing numerical code.

On the other hand, an important disadvantage of PML is that it is not associated with a clear notion of convergence, except under the expensive scenario of widening a layer where all physical and auxiliary fields are well-resolved. By contrast, in the case of high-order ABCs, with a fixed location of the boundary, one can approach the exact solution arbitrarily closely (up to discretization error) by increasing the order P of the ABC (with cost that increases only linearly with P). More efficient underresolved PMLs seem to be more sensitive to discretization and to the computational parameters than ABCs. A good design of an ABC at the continuous level usually guarantees good performance at the discrete level. This does not seem to be the general case for PML, where the matching between the solutions in the interior and in the layer at the discrete level is sometimes far from perfect. In addition, the theoretical analysis of a PML is usually more difficult than that for a high-order ABC for the same application. Additional discussion on the comparison of the two types of methods can be found in [1,3,9].

In [13] we presented a new method, which shares some features of both the PML and the high-order ABC, but enjoys some of the advantages that each of them lacks. In the new method, called the Double Absorbing Boundary (DAB) method, a high-order ABC is applied on two parallel artificial boundaries, which are a small distance apart. Auxiliary variables are defined on the two boundaries and in the thin layer beteween them. Like the PML, the DAB does not require special treatment of corners. The algebra involved is relatively simple, since no elimination of normal derivatives is needed. As in the method of high-order ABCs on a single boundary, DAB is clearly associated with the notion of convergence; one can approach the exact solution arbitrarily closely (up to discretization error) by increasing the order P, with only linearly-increasing cost. The numerical properties of DAB, like accuracy, stability and sensitivity to discretization, are similar to those of a high-order ABC on a single boundary.

In [13], the new method was applied to the scalar wave equation. We incorporated the DAB in a fully explicit finite difference scheme in 1D, and in a Finite Element (FE) scheme in 2D. In [14], a well-posedness proof was provided for the DAB scheme for the acoustics problem written in second-order form. The energy method was employed to obtain uniform-in-time estimates of the norm of the solution and the auxiliary functions, thus establishing the well-posedness and asymptotic stability of the DAB formulation. In addition, in [14] the DAB was applied to problems in 2D isotropic elastodynamics, written in first-order conservation form. The problem was discretized using the Lax-Wendroff finite difference scheme.

Although DAB is a general approach, and in principle can be used with any high-order ABC applied on the double boundary, in [13,14] the ABCs of the form proposed by Hagstrom and Warburton (H-W) in [11,15] were considered. The ABC formulation in [15], called the Complete Radiation Boundary Condition (CRBC), generalizes that in [11], and leads to an almost uniform-in-time error estimate for both propagating and decaying waves. The H-W ABCs have been incorporated in both FE and finite difference schemes, and have been shown to be extremely effective in a variety of situations, including those associated with dispersive, stratified, anisotropic and convective media [16-18], and where exterior sources are present (nesting) [19].

Recently, an H-W type ABC was applied to problems in elastodynamics [20,21]. This is the first known high-order ABC for elastodynamics that is long-time stable. The key to obtaining stability turns out to be the use of the Lysmer-Kuhlemeyer (LK) ABC as the termination condition of the recursive relations between the auxiliary variables ϕ j . The LK ABC is classical [22,23] and is commonly employed in solid earth geophysics computations, either as originally proposed or with some improvements; see, e.g., [24-26]. In the stable high-order ABC for elastodynamics, all the recursive relations except the last one are scalar in nature (as in the H-W ABC), while the last one is the vectorial LK condition. In [20], we have proved this combination to be stable and converging, at the continuous level.

In this paper, the DAB is incorporated in a standard FE scheme for 2D elastodynamics in homogeneous and layered media. This is the first FE implementation of DAB for elastodynamics. The design of ABCs for heterogeneous media poses an additional challenge. PML has been used for heterogeneous media in [27], for the scalar wave equation, in [28], for the MHD equations, and in [29], for elastodynamics. High order ABCs have been applied to continuously-stratified and layered media, in [16,30]. We show that the DAB formulation, when combined with the FE method, fits quite naturally to layered media, and is as simple as DAB for the homogeneous case.

Our formulation and numerical examples assume periodic boundary conditions along the boundaries perpendicular to the artificial DAB boundaries. This choice is made since some stability issues arise when the periodic conditions are replaced by some physical boundary conditions (e.g., traction free conditions). Attempts to resolve these issues are underway.

Following is the outline of the remaining sections. In Section “The problem in the computational domain and the double absorbing boundary (DAB) method” we describe the DAB method for elastodynamics. We consider the elastic equations of motion, with a possibly nonzero damping term. We first describe the DAB formulation for a homogeneous medium, and then we discuss the case of a layered medium. In Section “Finite elements: variational formulation” we describe the FE formulation of the elastodynamics problem, using the DAB. We also discuss some important computational aspects. In Section “Results and discussion” we present some numerical experiments to demonstrate the performance of the method. We conclude with some remarks in Section “Conclusions”.

Methods

The original unbounded-domain problem of elastodynamics

We consider a two-dimensional semi-infinite elastic wave guide of width b, as shown in Figure 1(a). A Cartesian coordinate system (x,y) is introduced with the origin at the southwest corner, so that the waveguide is parallel to the x direction, and y[0,b]. The south, north and west boundaries of the waveguide are denoted Γ S , Γ N and Γ W , respectively.

Figure 1
figure 1

Setup: (a) the original unbounded domain, (b) the computation domain, with the DAB layer at the right end.

In the waveguide we consider the equation of elastodynamics with mass- and stiffness-proportional damping, i.e.,

$$\begin{array}{*{20}l} \rho \ddot u_{i} + A_{M} \rho\dot u_{i} =\sigma_{ij,j}(\boldsymbol{u})+f_{i} \, \end{array} $$
((1))
$$\begin{array}{*{20}l} \sigma_{ij} = C_{ijkl} (\epsilon_{kl} + A_{K} \dot\epsilon_{kl}), \end{array} $$
((2))
$$ {\varepsilon}_{kl}=\left({u}_{k,l}+{u}_{l,k}\right)/2\kern1em . $$
((3))

Here and elsewhere, the summation convention is understood, and a comma subscript denotes partial differentiation with respect to the variable following it. We shall assign the values x and y to the indices i,j,k,l. A superposed dot indicates differentiation with respect to time. In (1)–(3), ρ is the density of the medium, u= {u i } is the unknown displacement vector field, σ= [σ ij ] is the unknown stress tensor, f= {f i } is a given body-force vector, C= [ C ijkl ] is the elastic moduli tensor, and ε= [ε ij ] is the unknown strain tensor. In (1) and (2), A M and A K are the mass-proportional and stiffness-proportional damping coefficients, which together define the Rayleigh damping that exists in the system.

We define the wave speeds c L and c T as

$$ c_{L} = \sqrt{\frac{\lambda+2\mu}{\rho}} \quad,\quad c_{T}=\sqrt{\frac{\mu}{\rho}}, $$
((4))

where λ and μ are the Lamé constants.

Boundary conditions are specified on the three boundaries: on Γ N and Γ S we impose periodic boundary conditions:

$$\begin{array}{*{20}l} & u_{i}(x,0,t) = u_{i}(x,b,t) \end{array} $$
((5))
$$\begin{array}{*{20}l} & u_{i,y}(x,0,t) = u_{i,y}(x,b,t), \end{array} $$
((6))

Periodic boundary conditions allow us to detach the study of the stability and accuracy of the ABC, which shall be introduced in the following, from the effects of corners and interactions between the absorbing boundary and other boundaries.

On Γ W any boundary condition may be imposed. We choose to prescribe zero displacement boundary conditions

$$\begin{array}{*{20}l} & u_{x}=0 \hspace{5mm} \text{on} \hspace{5mm} \Gamma_{W}\, \end{array} $$
((7))
$$ \begin{array}{ll}& {u}_y=0\kern14.22636pt \mathrm{on}\kern14.22636pt {\varGamma}_W\kern1em .\kern2em \end{array} $$
((8))

Initial conditions are also given, i.e.,

$$\begin{array}{*{20}l} & u_{i}(x,y,t=0)=u_{i0}(x,y), \end{array} $$
((9))
$$\begin{array}{*{20}l} & \dot u_{i}(x,y,t=0)=\dot u_{i0}(x,y) \end{array} $$
((10))

where u 0={u i0} and \(\dot {\boldsymbol {u}}_{0}=\{\dot {u}_{i0}\}\) are known functions.

We assume that outside a compact region, denoted Ω 0, in which C ijkl , f i , u i0 and \(\dot u_{i0}\) may, in principle, be general, the following simplified conditions hold: (a) the medium is homogeneous, namely C ijkl is constant; (b) there are no sources, namely f i =0; and (c) the initial values vanish, namely u i0=0 and \(\dot u_{i0}=0\).

The problem in the computational domain and the double absorbing boundary (DAB) method

We now truncate the semi-infinite domain by introducing the artificial boundary Γ E , located at x=x E and spanning 0≤yb. Slightly to the west of Γ E we set an interface denoted Γ I , located at x=x I , with 0≤yb. See Figure 1(b). The entire computational domain bounded by Γ N Γ W Γ S Γ E is denoted as Ω. As Figure 1(b) shows, this domain is divided by the interface Γ I into two sub-domains: the interior domain Ω I and a thin layer Ω L . We choose the location of Γ E and Γ I such that Ω 0 (the “irregular region” defined above) is strictly contained within Ω I . Thus, in the layer Ω L we have that C ijkl is constant, f i =0, and the initial conditions are zero.

The function u satisfies the equation of elastodynamics (1) in Ω, the periodic boundary conditions (5), (6) and the conditions (7), (8), and the initial conditions (9) and (10) in Ω. In the layer Ω L we shall apply a special treatment, with the goal of rendering the solution in the interior domain Ω I as close as possible to the solution of the original semi-infinite problem in that domain. Thus, Ω L will act as an absorbing or non-reflecting layer.

To this end we define a sequence of auxiliary variables \({\phi ^{0}_{i}}, \ldots, {\phi ^{P}_{i}}\), i=x,y in the layer only. Here P is a chosen parameter that will determine the order of accuracy of the absorbing layer. The first auxiliary variable is defined to be \({\phi ^{0}_{i}}=u_{i}\) in Ω L . The problem for the \({\phi ^{m}_{i}}\) is given as follows. In the layer, we require the \({\phi ^{m}_{i}}\) for a particular value of m to satisfy the same elastic wave equation as for u i , i.e. (1), with f i =0:

$$ \rho \ddot {\phi^{m}_{i}} + A_{M} \rho\dot {\phi^{m}_{i}} = \sigma_{ij,j}(\phi^{m}),\quad m=1,\ldots,P $$
((11))

We also denote by \({T^{m}_{i}}\) the traction vector on Γ I and Γ E that corresponds to the variables \({\phi ^{m}_{i}}\). Also, \(T_{i}={T^{0}_{i}}\).

All the auxiliary variables satisfy a zero initial condition:

$$ {\phi^{m}_{i}}(x,y,0)=0,\quad \dot{\phi}^{m}_{i}(x,y,0)=0, \quad m=1,\ldots,P,\quad \text{in }\Omega_{L}\ $$
((12))

On Γ S and Γ N (or more precisely the parts of these boundaries that are in the layer, i.e., for x I xx E and 0≤yb) we apply the same (periodic) boundary conditions as for u i , i.e.,

$$\begin{array}{*{20}l} & {\phi^{m}_{i}}(x,0,t) = {\phi^{m}_{i}}(x,b,t) \end{array} $$
((13))
$$ \begin{array}{ll}& {\phi}_{i,y}^m\left(x,0,t\right)={\phi}_{i,y}^m\left(x,b,t\right),\kern1em m=1,\dots, P,\kern1em \mathrm{on}{\varGamma}_S,\kern1em {\varGamma}_N\kern1em .\kern2em \end{array} $$
((14))

Now we need to define boundary conditions for the \({\phi ^{m}_{i}}\) on Γ I and on Γ E . We define the boundary conditions on both boundaries recursively. The recursive boundary conditions are defined in a manner that is in many respects similar to the ABC for elasticity defined by equation (28) in [20], with added low-order terms:

$$\begin{array}{*{20}l} & a_{m}{\dot\phi^{m}_{x}} + c_{x}\phi^{m}_{x,x} + b_{m}{\phi^{m}_{x}} = \tilde a_{m}\dot\phi^{m+1}_{x} - c_{x}\phi^{m+1}_{x,x} + \tilde b_{m}\phi^{m+1}_{x} \end{array} $$
((15))
$$\begin{array}{*{20}l} & a_{m}{\dot\phi^{m}_{y}} + c_{y}\phi^{m}_{y,x} + b_{m}{\phi^{m}_{x}} = \tilde a_{m}\dot\phi^{m+1}_{y} - c_{y}\phi^{m+1}_{y,x} + \tilde b_{m}\phi^{m+1}_{x} \\ & \quad \text{with } m=0,\ldots,P-1,\quad \text{on } \Gamma_{I},\ \Gamma_{E},\ \notag \end{array} $$
((16))

Here, c x and c y are some chosen wave speeds. In general, various combinations of c x and c y are possible in terms of c L and c T ; however, coefficients that do not satisfy c x =c y lead to an unstable formulation. In all the numerical experiments (section “3”) we took c x =c y =c L .

Also, in (15) and (16), a m , b m , \(\tilde a_{m}\) and \(\tilde b_{m}\) are the parameters of the boundary conditions as defined in the “Computational aspects” section below.

The use of the auxiliary functions \({\phi _{i}^{m}}\) and the corresponding equations (11)–(16) is motivated by the work of Hagstrom and Warburton [11], who introduced the auxiliary functions at the truncation boundary for the high-order ABCs. These ABCs have been shown to have excellent absorbing properties; see, e.g., [21]. In this work, we use (and later discretize) the \({\phi _{i}^{m}}\) within a small layer, and impose the H-W ABCs on both boundaries of this layer, employing the DAB construction devised in [13]. This frees us from the need to eliminate the normal derivatives of the auxiliary variables, which is necessary when the ABC is used on a single truncation boundary.

To complete the recursive definition, we require the following Lysmer-Kuhlemeyer (LK) termination condition on Γ E :

$$\begin{array}{*{20}l} & {T^{P}_{x}} + \rho c_{L} \dot {\phi^{P}_{x}} = \left(2\mu+\lambda\right)\phi^{P}_{x,x}+\lambda\phi^{P}_{y,y}+\rho c_{L}{\dot\phi^{P}_{x}} = 0 \end{array} $$
((17))
$$\begin{array}{*{20}l} & {T^{P}_{y}} + \rho c_{T} \dot {\phi^{P}_{y}} = \mu\left(\phi^{P}_{x,y} + \phi^{P}_{y,x}\right) + \rho c_{T} \dot {\phi^{P}_{y}} = 0. \end{array} $$
((18))

Figure 2 illustrates the “ladder” structure of the DAB scheme, namely the flow of information on the two boundaries bounding the DAB layer.

Figure 2
figure 2

The “ladder” structure of the DAB, showing the flow of information on the two boundaries bounding the layer.

Layered media

In addition to the homogeneous scheme described in the previous subsections, we also consider the case of layered media with different material properties in the different layers. The layers are assumed to be perpendicular to the absorbing boundary. The jump conditions at the interfaces between the different layers for the case of the wave equation were given in [16]. For the elastic case, a procedure similar to the one described in that paper leads to the following jump conditions:

((19))
((20))

Here denoted a jump across the interface between two layers.

By dividing the expression for the recursive relation (15) by \(a_{m}^{(n)}\), where the superscript (n) indicates the n-th layer, and denoting c=c x =c y , we obtain for each layer n:

$$\begin{array}{*{20}l} & {\dot\phi^{m}_{i}} + \frac{c^{(n)}}{a_{m}^{(n)}}\phi^{m}_{i,x} + \frac{b_{m}^{(n)}}{a_{m}^{(n)}}{\phi^{m}_{i}} = \frac{\tilde a_{m}^{(n)}}{a_{m}^{(n)}}\dot\phi^{m+1}_{i} - \frac{c^{(n)}}{a_{m}^{(n)}}\phi^{m+1}_{i,x} +\frac{\tilde b_{m}^{(n)}}{a_{m}^{(n)}}\phi^{m+1}_{i} \\ & \quad \text{with } m=0,\ldots,P-1; \quad i=x,y, \quad \text{on } \Gamma_{I},\ \Gamma_{E}.\ \notag \end{array} $$
((21))

Since the auxiliary variables \({\phi ^{m}_{i}}\) and their derivatives in the x direction and in time are continuous across the layers, we will require that the ratios \(\frac {c_{x}^{(n)}}{a_{m}^{(n)}}\), \(\frac {b_{m}^{(n)}}{a_{m}^{(n)}}\), \(\frac {\tilde a_{m}^{(n)}}{a_{m}^{(n)}}\) and \(\frac {\tilde b_{m}^{(n)}}{a_{m}^{(n)}}\) also be continuous across the layers, namely that they do not depend on n. Thus if we denote by some reference value of the corresponding parameter, by requiring that

$$ \frac{c^{(n)}}{a_{m}^{(n)}} = \frac{c^{*}}{a_{m}^{*}}; \quad \frac{b_{m}^{(n)}}{a_{m}^{(n)}} = \frac{b_{m}^{*}}{a_{m}^{*}}; \quad \frac{\tilde a_{m}^{(n)}}{a_{m}^{(n)}}=\frac{\tilde a_{m}^{*}}{a_{m}^{*}}; \quad \frac{\tilde b_{m}^{(n)}}{a_{m}^{(n)}}=\frac{\tilde b_{m}^{*}}{a_{m}^{*}}, $$
((22))

all the recursive relations (21) for the different layers can be replaced by the following recursive relation written out for a reference layer:

$$\begin{array}{*{20}l} & {\dot\phi^{m}_{i}} + \frac{c^{*}}{a_{m}^{*}}\phi^{m}_{i,x} + \frac{b_{m}^{*}}{a_{m}^{*}}{\phi^{m}_{i}} = \frac{\tilde a_{m}^{*}}{a_{m}^{*}}\dot\phi^{m+1}_{i} - \frac{c^{*}}{a_{m}^{*}}\phi^{m+1}_{i,x} +\frac{\tilde b_{m}^{*}}{a_{m}^{*}}\phi^{m+1}_{i} \\ & \quad \text{with } m=0,\ldots,P-1; \quad i=x,\quad \text{on } \Gamma_{I},\ \Gamma_{E}.\ \notag \end{array} $$
((23))

Now, when we consider the termination conditions for the different layers, we take into account that

$$ {\sigma}_{xy}^m={T}_y^m\left|{}_{\varGamma_E}={T}_x^m\right|{}_{\mathrm{interface}\ \mathrm{between}\ \mathrm{the}\ \mathrm{layers}}\kern1em . $$
((24))

Therefore, on the interface between the layers we obtain from the second termination condition (18) that:

$$ {T^{P}_{x}} + \rho c_{T}^{(n)} \dot {\phi^{P}_{y}} = 0, $$
((25))

which cannot hold since \({T^{P}_{x}}\), \({\dot \phi ^{P}_{y}}\) are continuous, while \(\rho c_{T}^{(n)}\) is not. However, in practice, any positive wave speed can be substituted for c T in (18) without loss of stability. The slightly reduced accuracy that is caused by this modification will be compensated for by the recursive relations. Thus, in the case of layered media we take the following termination conditions:

$$\begin{array}{*{20}l} & {T^{P}_{x}} + \rho c_{L}^{*} \dot {\phi^{P}_{x}} = 0 \end{array} $$
((26))
$$\begin{array}{*{20}l} & {T^{P}_{y}} + \rho c_{T}^{*} \dot {\phi^{P}_{y}} = 0. \end{array} $$
((27))

Eqs. (23), (26) and (27) constitute the modified DAB formulation for the case of layered media. The spatial discretization in this case is carried out in the same manner as for the homogeneous medium, using the starred parameters \(a_{m}^{*}\), \(\tilde a_{m}^{*}\), \(b_{m}^{*}\), \(\tilde b_{m}^{*}\), c , \(c_{L}^{*}\), \(c_{T}^{*}\) which pertain to a chosen reference layer.

Finite elements: variational formulation

The strong form of the DAB scheme consists of the elastic Eqs. (1) in Ω and (11) in Ω L , the initial conditions (9), (10) in Ω and (12) in Ω L , the periodic boundary conditions (5), (6) and (13), (14) on Γ S and Γ N , the boundary conditions (7) and (8), and the recursive boundary relations (15) and (16) on Γ I and Γ E , with the termination conditions (17) and (18) on Γ E .

To derive the weak form, we first define the following function spaces:

$$\begin{array}{*{20}l} & \mathcal{S}=\Bigl\lbrace\boldsymbol{u}\equiv\{u_{i}\}\,\Big\vert\, \{u_{i}\}\in H^{1}(\Omega), u_{i}(\Gamma_{S})=u_{i}(\Gamma_{N}), u_{i}(\Gamma_{W})=0\Bigr\rbrace \end{array} $$
((28))
$$\begin{array}{*{20}l} & \mathcal{S}_{L}=\Bigl\lbrace\phi^{m}\equiv\{{\phi^{m}_{i}}\}\,\Big\vert\, \{{\phi^{m}_{i}}\}\in H^{1}(\Omega_{L}), {\phi^{m}_{i}}(\Gamma_{S})={\phi^{m}_{i}}(\Gamma_{N}), m=1,\dotsc P\Bigr\rbrace \end{array} $$
((29))

We multiply Eq. 1 by a weight function \(\boldsymbol {w}^{0}\in \mathcal {S}\) and Eqs. (11) by the weight functions \(\boldsymbol {w}^{m}\in \mathcal {S}_{L}\), m=1,…,P, and integrate the terms that include C ijkl by parts to obtain

$$\begin{array}{*{20}l} & \int_{\Omega} {w^{0}_{i}}\rho u_{i,tt} \,d\Omega +A_{M} \int_{\Omega} {w^{0}_{i}} \rho u_{i,t} \,d\Omega\notag \\ & \quad +\int_{\Omega} w^{0}_{i,j}C_{ijkl} u_{k,l}\,d\Omega +A_{K}\int_{\Omega} w^{0}_{i,j}C_{ijkl} \dot u_{k,l}\,d\Omega + {B^{0}_{E}} =\int_{\Omega} {w^{0}_{i}} f_{i} \,d\Omega \end{array} $$
((30))
$$\begin{array}{*{20}l} & \int_{\Omega_{L}} {w^{m}_{i}}\rho \phi^{m}_{i,tt} \,d\Omega +A_{M} \int_{\Omega_{L}} {w^{m}_{i}} \rho \phi^{m}_{i,t} \,d\Omega\notag \\ & \quad +\int_{\Omega_{L}} w^{m}_{i,j}C_{ijkl} \phi^{m}_{k,l}\,d\Omega +A_{K} \int_{\Omega_{L}} w^{m}_{i,j}C_{ijkl}\dot\phi^{m}_{k,l}\,d\Omega+ {B^{m}_{I}} + {B^{m}_{E}} = 0, \notag \\ & \qquad m=1,\dotsc,P, \end{array} $$
((31))

where we have

$$\begin{array}{*{20}l} & {B^{m}_{I}}=-\int_{\Gamma_{I}}\left[{w^{m}_{x}} {T^{m}_{x}} + {w^{m}_{y}} {T^{m}_{y}}\right]d\Gamma \end{array} $$
((32))
$$\begin{array}{*{20}l} & {B^{m}_{E}}=-\int_{\Gamma_{E}}\left[{w^{m}_{x}} {T^{m}_{x}} + {w^{m}_{y}} {T^{m}_{y}}\right]d\Gamma. \end{array} $$
((33))

We rewrite the recursive relations (15) and (16) in the following form:

On the boundary Γ I :

$$\begin{array}{*{20}l} & \phi^{m}_{x,x}=\frac{1}{c_{x}}\left(\tilde a_{m-1}{\dot\phi^{m}_{x}}-a_{m-1}\dot\phi^{m-1}_{x}\right) + \frac{1}{c_{x}}\left(\tilde b_{m-1} {\phi^{m}_{x}}-b_{m-1} \phi^{m-1}_{x}\right) - \phi^{m-1}_{x,x} \end{array} $$
((34))
$$\begin{array}{*{20}l} & \phi^{m}_{y,x}=\frac{1}{c_{y}}\left(\tilde a_{m-1}{\dot\phi^{m}_{y}}-a_{m-1}\dot\phi^{m-1}_{y}\right) + \frac{1}{c_{y}}\left(\tilde b_{m-1} {\phi^{m}_{y}}-b_{m-1} \phi^{m-1}_{y}\right) - \phi^{m-1}_{y,x} \\ & \quad m=1,\ldots,P, \quad \text{on } \Gamma_{I} \notag \end{array} $$
((35))

On the boundary Γ E :

$$\begin{array}{*{20}l} & \phi^{m}_{x,x}= \frac{1}{c_{x}}\left(\tilde a_{m}\dot\phi^{m+1}_{x} - a_{m}{\dot\phi^{m}_{x}}\right)+ \frac{1}{c_{x}}\left(\tilde b_{m} \phi^{m+1}_{x} - b_{m} {\phi^{m}_{x}}\right)- \phi^{m+1}_{x,x} \end{array} $$
((36))
$$\begin{array}{*{20}l} & \phi^{m}_{y,x}= \frac{1}{c_{y}}\left(\tilde a_{m}\dot\phi^{m+1}_{y} - a_{m}{\dot\phi^{m}_{y}}\right)+ \frac{1}{c_{y}}\left(\tilde b_{m} \phi^{m+1}_{y} - b_{m} {\phi^{m}_{y}}\right)- \phi^{m+1}_{y,x} \\ & \quad m=0,\ldots,P-1, \quad \text{on } \Gamma_{E} \notag \end{array} $$
((37))

and (17) and (18) serve as the termination conditions:

$$\begin{array}{*{20}l} & {T^{P}_{x}} =-\rho c_{L} \dot {\phi^{P}_{x}} \end{array} $$
((38))
$$\begin{array}{*{20}l} & {T^{P}_{y}} =-\rho c_{T} \dot {\phi^{P}_{y}} , \quad \text{on} \Gamma_{E}, \end{array} $$
((39))

in which case, by substituting (34)–(37) in the terms \({B^{m}_{E}}\), \({B^{m}_{I}}\) in (32) and (33), these terms become:

$${} {\fontsize{9.2pt}{9.6pt}\selectfont{\begin{aligned} {B_{I}^{m}} &= \int_{\Gamma_{I}} {w^{m}_{x}}\left[(2\mu+\lambda)\phi^{m}_{x,x}+\lambda\phi^{m}_{y,y}\right]d\Gamma+ \int_{\Gamma_{I}} {w^{m}_{y}}\left[\mu\phi^{m}_{x,y}+\mu\phi^{m}_{y,x}\right]d\Gamma \\ & =\int_{\Gamma_{I}} {w^{m}_{x}}\left[\frac{2\mu+\lambda}{c_{x}}\tilde a_{m-1}{\dot\phi^{m}_{x}} -\frac{2\mu+\lambda}{c_{x}} a_{m-1}\dot\phi^{m-1}_{x} +\frac{2\mu+\lambda}{c_{x}}\tilde b_{m-1}{\phi^{m}_{x}} -\frac{2\mu+\lambda}{c_{x}} b_{m-1}\phi^{m-1}_{x} \right.\\ &\left. \qquad -(2\mu+\lambda)\phi^{m-1}_{x,x} +\lambda\phi^{m}_{y,y}\right]\,d\Gamma + \\ & \int_{\Gamma_{I}} {w^{m}_{y}}\left[\frac{\mu}{c_{y}}\tilde a_{m-1}{\dot\phi^{m}_{y}} -\frac{\mu}{c_{y}} a_{m-1}\dot\phi^{m-1}_{y} +\frac{\mu}{c_{y}}\tilde b_{m-1}{\phi^{m}_{y}} -\frac{\mu}{c_{y}} b_{m-1}\phi^{m-1}_{y} -\mu\phi^{m-1}_{y,x} +\mu\phi^{m}_{x,y}\right]\,d\Gamma, \\ & \quad m=1,\ldots,P \quad \text{on}\quad \Gamma_{I} \end{aligned}}} $$
((40))
$${} {\fontsize{9.2pt}{9.6pt}\selectfont{\begin{aligned} {B_{E}^{m}} = &-\int_{\Gamma_{E}} {w^{m}_{x}}\left[(2\mu+\lambda)\phi^{m}_{x,x}+\lambda\phi^{m}_{y,y}\right]d\Gamma -\int_{\Gamma_{E}} {w^{m}_{y}}\left[\mu\phi^{m}_{x,y}+\mu\phi^{m}_{y,x}\right]d\Gamma= \\ &-\int_{\Gamma_{E}} {w^{m}_{x}}\left[\frac{2\mu+\lambda}{c_{x}}\tilde a_{m}\dot\phi^{m+1}_{x} -\frac{2\mu+\lambda}{c_{x}} a_{m}{\dot\phi^{m}_{x}} +\frac{2\mu+\lambda}{c_{x}}\tilde b_{m}\phi^{m+1}_{x} -\frac{2\mu+\lambda}{c_{x}} b_{m}{\phi^{m}_{x}} \right.\\ & \qquad \left.-(2\mu+\lambda)\phi^{m+1}_{x,x} +\lambda\phi^{m}_{y,y}\right]\,d\Gamma \\ &-\int_{\Gamma_{E}} {w^{m}_{y}}\left[\frac{\mu}{c_{y}}\tilde a_{m}\dot\phi^{m+1}_{y} -\frac{\mu}{c_{y}} a_{m}{\dot\phi^{m}_{y}} +\frac{\mu}{c_{y}}\tilde b_{m}\phi^{m+1}_{y} -\frac{\mu}{c_{y}} b_{m}{\phi^{m}_{y}} -\mu\phi^{m+1}_{y,x} +\mu\phi^{m}_{x,y}\right]\,d\Gamma, \\ & \quad m=0,\ldots,P-1 \quad \text{on}\quad \Gamma_{E} \end{aligned}}} $$
((41))

and (38) and (39) are substituted directly into (33), with m=P:

$$ {B_{E}^{P}} = \int_{\Gamma_{E}} {w^{P}_{x}} \rho c_{L} {\dot\phi^{P}_{x}} d\Gamma + \int_{\Gamma_{E}} {w^{P}_{y}} \rho c_{T} {\dot\phi^{P}_{y}} d\Gamma \, \text{ on } \Gamma_{E}, $$
((42))

Here and elsewhere, the x-derivatives on the boundaries Γ I ,Γ E are calculated in a one-sided manner everywhere, except for the derivatives \(\phi ^{0}_{x,x}\) and \(\phi ^{0}_{y,x}\) on Γ I , which are calculated as the average of the derivative in the elements adjacent to this boundary from the left and from the right.

Thus the weak form is: find \(\boldsymbol {u}\equiv \{u_{i}\}\in \mathcal {S}(\Omega)\) and \(\phi ^{m}\equiv \{{\phi ^{m}_{i}}\}\in \mathcal {S}_{L}(\Omega _{L})\), m=1,…,P, which satisfy the initial conditions (9), (10) in Ω and (12) in Ω L , and satisfy (30) and (31) for all \(\boldsymbol {w}^{0}\equiv \{w_{i}\}\in \mathcal {S}(\Omega)\) and all \(\boldsymbol {w}^{m}\equiv \{{w^{m}_{i}}\}\in \mathcal {S}_{L}(\Omega _{L})\), m=1,…,P.

Semi-discrete form

We discretize the weak form described in the previous subsection in space using the standard Galerkin FE method. At the global level, the variables u i in Ω and \({\phi ^{m}_{i}}\) in Ω L are replaced by their finite-dimensional approximations

$$\begin{array}{*{20}l} & {u^{h}_{i}}(\boldsymbol{x},t)=\sum_{A\in \eta_{i}}d^{0h}_{Ai}(t)N_{A}(\boldsymbol{x}), \quad \boldsymbol{x}\in\Omega, i=x,y \end{array} $$
((43))
$$\begin{array}{*{20}l} & \phi^{mh}_{i}(\boldsymbol{x},t)=\sum_{A\in {\eta_{i}^{L}}}d^{mh}_{Ai}(t)N_{A}(\boldsymbol{x}), \quad \boldsymbol{x}\in\Omega_{L}, i=x,y, \quad m=1,\dotsc,P \end{array} $$
((44))

Here h is the mesh parameter, A stands for the global node number, η i is the set of nodes in Ω on which there is no essential boundary condition prescribed on the displacement in the i direction, \({\eta ^{L}_{i}}\) is the set of nodes in Ω L on which there are no essential boundary conditions prescribed on the displacement in the i direction, N A is the global-level shape function associated with the variables u i in Ω and node A, or with the variables \({\phi ^{m}_{i}}\) in Ω L and node A. We use identical bilinear shape functions for the discretization of u i and for the discretization of all the \({\phi ^{m}_{i}}\).

At the element level the analogous expansion is:

$$\begin{array}{*{20}l} & {u^{e}_{i}}\left(\boldsymbol{x},t\right)=\sum_{a=1}^{N_{en}}d^{0e}_{ai}(t)N_{a}\left(\boldsymbol{x}\right), \quad \boldsymbol{x}\in\Omega^{e}, i=x,y \end{array} $$
((45))
$$\begin{array}{*{20}l} & \phi^{me}_{i}\left(\boldsymbol{x},t\right)=\sum_{a=1}^{N_{en}}d^{me}_{ai}(t)N_{a}(\boldsymbol{x}), \quad \boldsymbol{x}\in\Omega^{e}, i=x,y, \quad m=1,\dotsc,P \end{array} $$
((46))

Here e stands for an element number, N en is the number of element nodes, a stands for the local element node number, Ω e is the domain of element e, N a is the element-level displacement shape function associated with the node a, the quantities \(d^{me}_{\textit {ai}}\), m=0,…,P are the values of \({u^{e}_{i}}\) and \(\phi ^{me}_{i}\) at node a of element e. We also write global expressions similar to (43) and (44) for the test functions \(w^{mh}_{i}\), m=0,…,P and element level expressions similar to (45) and (46) for the element-level test functions \(w^{me}_{i}\), m=0,…,P.

Substitution of the approximations (43) and (44) into the weak Eqs. 34–(42) yields a system of ordinary differential equations in time, of the form:

$$\begin{array}{*{20}l} & \boldsymbol{M}^{11}_{0} \ddot{\boldsymbol{d}}^{x} + \boldsymbol{C}^{11}_{0} \dot{\boldsymbol{d}}^{x} + {\boldsymbol{C}}^{12}_{0} \dot{\boldsymbol{d}}^{y} +\boldsymbol{k}^{11}_{0} \boldsymbol{d}^{x} +\boldsymbol{k}^{12}_{0} \boldsymbol{d}^{y} +\boldsymbol{G}^{11}_{0} \dot{\phi}^{x}_{1} + \boldsymbol{H}^{11}_{0} {\phi^{x}_{1}} = \boldsymbol{f}^{x} \end{array} $$
((47))
$$\begin{array}{*{20}l} & \boldsymbol{M}^{22}_{0} \ddot{\boldsymbol{d}}^{y} + \boldsymbol{C}^{12}_{0} \dot{\boldsymbol{d}}^{x} + \boldsymbol{C}^{22}_{0} \dot{\boldsymbol{d}}^{y} +\boldsymbol{k}^{21}_{0} \boldsymbol{d}^{x} +\boldsymbol{k}^{22}_{0} \boldsymbol{d}^{y} +\boldsymbol{G}^{22}_{0} {\dot\phi^{y}_{1}} + \boldsymbol{H}^{22}_{0} {\phi^{y}_{1}} = \boldsymbol{f}^{y} \end{array} $$
((48))

For m=1,…,P−1

$$\begin{array}{*{20}l} & \boldsymbol{M}^{11}_{m} {\ddot\phi^{x}_{m}} + \boldsymbol{C}^{11}_{m} {\dot\phi^{x}_{m}} + \boldsymbol{C}^{12}_{m} {\dot\phi^{y}_{m}} +\boldsymbol{k}^{11}_{m} {\phi^{x}_{m}} +\boldsymbol{k}^{12}_{m} {\phi^{y}_{m}} +\boldsymbol{A}^{11}_{m} \dot\phi^{x}_{m-1} + \boldsymbol{B}^{11}_{m} \phi^{x}_{m-1} \notag\\ &\qquad +\boldsymbol{G}^{11}_{m} \dot\phi^{x}_{m+1} + \boldsymbol{H}^{11}_{m} \phi^{x}_{m+1} = 0 \end{array} $$
((49))
$$\begin{array}{*{20}l} & \boldsymbol{M}^{22}_{m} {\ddot\phi^{y}_{m}} + \boldsymbol{C}^{12}_{m} {\dot\phi^{x}_{m}} + \boldsymbol{C}^{22}_{m} {\dot\phi^{y}_{m}} +\boldsymbol{k}^{21}_{m} {\phi^{x}_{m}} +\boldsymbol{k}^{22}_{m} {\phi^{y}_{m}} +\boldsymbol{A}^{22}_{m} \dot\phi^{y}_{m-1} + \boldsymbol{B}^{22}_{m} \phi^{y}_{m-1} \notag\\ &\qquad +\boldsymbol{G}^{22}_{m} \dot\phi^{y}_{m+1} + \boldsymbol{H}^{22}_{m} \phi^{y}_{m+1} = 0 \end{array} $$
((50))

and also:

$$\begin{array}{*{20}l} & \boldsymbol{M}^{11}_{P} {\ddot\phi^{x}_{P}} + \boldsymbol{C}^{11}_{P} {\dot\phi^{x}_{P}} + \boldsymbol{C}^{12}_{P} {\dot\phi^{y}_{P}} +\boldsymbol{k}^{11}_{P} {\phi^{x}_{P}} +\boldsymbol{k}^{12}_{P} {\phi^{y}_{P}} +\boldsymbol{A}^{11}_{P} \dot{\phi}^{x}_{P-1} + \boldsymbol{B}^{11}_{P} \phi^{x}_{P-1} = 0 \end{array} $$
((51))
$$\begin{array}{*{20}l} & \boldsymbol{M}^{22}_{P} {\ddot\phi^{y}_{P}} + \boldsymbol{C}^{12}_{P} {\dot\phi^{x}_{P}} + \boldsymbol{C}^{22}_{P} {\dot\phi^{y}_{P}} +\boldsymbol{k}^{21}_{P} {\phi^{x}_{P}} +\boldsymbol{k}^{22}_{P} {\phi^{y}_{P}} +\boldsymbol{A}^{22}_{P} \dot{\phi}^{y}_{P-1} + \boldsymbol{B}^{22}_{P} \phi^{y}_{P-1} = 0 \end{array} $$
((52))

with the initial conditions

$$\begin{array}{*{20}l} & \boldsymbol{d}^{i}(t=0)={\boldsymbol{d}^{i}_{0}}; \quad \dot{\boldsymbol{d}}^{i}(t=0)={\boldsymbol{v}_{0}^{i}}; \quad i=x,y \end{array} $$
((53))
$$\begin{array}{*{20}l} & {\phi^{i}_{m}}(t=0)=\mathbf{0}; \quad {\dot\phi^{i}_{m}}(t=0)=\mathbf{0}; \quad i=x,y, \quad m=1,\dotsc,P \end{array} $$
((54))

Here d i and \({\phi _{m}^{i}}\) are the vectors whose entries are the unknown nodal values of u i in Ω and of \({\phi _{m}^{i}}\) in Ω L , respectively; a dot indicates differentiation with respect to time.

The element level expressions may be extracted from (34)–(42). The first (mass) and second (stiffness) terms in each of (30) and (31) contribute to the global mass and stiffness matrices respectively via the following arrays at the element level (we denote these arrays by the superscript OM):

$$\begin{array}{*{20}l} &{} M_{0ab}^{OM11e} = \rho \int_{\Omega^{e}} N_{a} N_{b} d\,\Omega \end{array} $$
((55))
$$\begin{array}{*{20}l} &{} M_{0ab}^{OM22e} = \rho \int_{\Omega^{e}} N_{a} N_{b} d\,\Omega \end{array} $$
((56))
$$\begin{array}{*{20}l} &{} C_{0ab}^{OM11e} = A_{M} \rho\int_{\Omega^{e}} N_{a} N_{b} d\,\Omega + A_{K} (2\mu+\lambda) \int_{\Omega^{e}} N_{a,x} N_{b,x} d\,\Omega + A_{K} \mu \int_{\Omega^{e}} N_{a,y} N_{b,y} d\,\Omega \end{array} $$
((57))
$$\begin{array}{*{20}l} &{} C_{0ab}^{OM22e} = A_{M} \rho\int_{\Omega^{e}} N_{a} N_{b} d\,\Omega + A_{K} (2\mu+\lambda) \int_{\Omega^{e}} N_{a,y} N_{b,y} d\,\Omega + A_{K} \mu \int_{\Omega^{e}} N_{a,x} N_{b,x} d\,\Omega \end{array} $$
((58))
$$\begin{array}{*{20}l} &{} C_{0ab}^{OM12e} = A_{K} \mu \int_{\Omega^{e}} N_{a,y} N_{b,x} d\,\Omega + A_{K} \lambda \int_{\Omega^{e}} N_{a,x} N_{b,y} d\,\Omega \end{array} $$
((59))
$$\begin{array}{*{20}l} &{} C_{0ab}^{OM21e} = A_{K} \mu \int_{\Omega^{e}} N_{a,x} N_{b,y} d\,\Omega + A_{K}\lambda \int_{\Omega^{e}} N_{a,y} N_{b,x} d\,\Omega \end{array} $$
((60))
$$\begin{array}{*{20}l} &{} K_{0ab}^{OM11e} = (2\mu+\lambda) \int_{\Omega^{e}} N_{a,x} N_{b,x} d\,\Omega + \mu \int_{\Omega^{e}} N_{a,y} N_{b,y} d\,\Omega \end{array} $$
((61))
$$\begin{array}{*{20}l} &{} K_{0ab}^{OM22e} = (2\mu+\lambda) \int_{\Omega^{e}} N_{a,y} N_{b,y} d\,\Omega + \mu \int_{\Omega^{e}} N_{a,x} N_{b,x} d\,\Omega \end{array} $$
((62))
$$\begin{array}{*{20}l} &{} K_{0ab}^{OM12e} = \mu \int_{\Omega^{e}} N_{a,y} N_{b,x} d\,\Omega + \lambda \int_{\Omega^{e}} N_{a,x} N_{b,y} d\,\Omega \end{array} $$
((63))
$$\begin{array}{*{20}l} &{} K_{0ab}^{OM21e} = \mu \int_{\Omega^{e}} N_{a,x} N_{b,y} d\,\Omega + \lambda \int_{\Omega^{e}} N_{a,y} N_{b,x} d\,\Omega \end{array} $$
((64))

The arrays of this form, with the first subscript taken as 0, correspond to the first and second terms in (30). The expressions corresponding to the first and second terms in (31) are similar, except that the integration takes place in an element in Ω I rather than in Ω and that the first subscript of M and K is taken as m instead of 0 (m=1,…,P).

The contribution from the \({B^{m}_{I}}\) term (denoted by the superscript BI) to the damping matrices is (m=1,…,P) is:

$$\begin{array}{*{20}l} & C_{mab}^{BI11e} = \frac{2\mu+\lambda}{c_{x}} \tilde a_{m-1} \int_{\Gamma_{I}} N_{a}N_{b} d\,\Gamma \end{array} $$
((65))
$$\begin{array}{*{20}l} & A_{mab}^{BI11e} =-\frac{2\mu+\lambda}{c_{x}} a_{m-1} \int_{\Gamma_{I}} N_{a}N_{b} d\,\Gamma \end{array} $$
((66))
$$\begin{array}{*{20}l} & C_{mab}^{BI22e} = \frac{\mu}{c_{y}} \tilde a_{m-1} \int_{\Gamma_{I}} N_{a}N_{b} d\,\Gamma \end{array} $$
((67))
$$\begin{array}{*{20}l} & A_{mab}^{BI22e} =-\frac{\mu}{c_{y}} a_{m-1} \int_{\Gamma_{I}} N_{a}N_{b} d\,\Gamma, \end{array} $$
((68))

The contribution from the \({B^{m}_{I}}\) term to the stiffness matrixes is:

$$\begin{array}{*{20}l} & K_{mab}^{BI11e} = \frac{2\mu+\lambda}{c_{x}} \tilde b_{m-1} \int_{\Gamma_{I}} N_{a}N_{b} d\,\Gamma \end{array} $$
((69))
$$\begin{array}{*{20}l} & B_{mab}^{BI11e} =-\frac{2\mu+\lambda}{c_{x}} b_{m-1} \int_{\Gamma_{I}} N_{a}N_{b,x} d\,\Gamma - (2\mu+\lambda) \int_{\Gamma_{I}} N_{a}N_{b,x}d\,\Gamma \end{array} $$
((70))
$$\begin{array}{*{20}l} & K_{mab}^{BI12e} = \lambda \int_{\Gamma_{I}} N_{a}N_{b,y}d\,\Gamma \end{array} $$
((71))
$$\begin{array}{*{20}l} & K_{mab}^{BI21e} = \mu \int_{\Gamma_{I}} N_{a}N_{b,y}d\,\Gamma \end{array} $$
((72))
$$\begin{array}{*{20}l} & K_{mab}^{BI22e} = \frac{\mu}{c_{y}} \tilde b_{m-1} \int_{e} N_{a}N_{b} d\,\Gamma \end{array} $$
((73))
$$\begin{array}{*{20}l} & B_{mab}^{BI22e} =-\frac{\mu}{c_{y}} b_{m-1} \int_{e} N_{a}N_{b} d\,\Gamma - \mu \int_{\Gamma_{I}} N_{a}N_{b,x}d\,\Gamma, \end{array} $$
((74))

Similarly, the contribution from the \({B^{m}_{E}}\) term (denoted by the superscript BE) to the damping matrices is (m=0,…,P−1):

$$\begin{array}{*{20}l} & G_{mab}^{BE11e} =-\frac{2\mu+\lambda}{c_{x}} \tilde a_{m} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma \end{array} $$
((75))
$$\begin{array}{*{20}l} & C_{mab}^{BE11e} = \frac{2\mu+\lambda}{c_{x}} a_{m} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma \end{array} $$
((76))
$$\begin{array}{*{20}l} & G_{mab}^{BE22e} =-\frac{\mu}{c_{y}} \tilde a_{m} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma \end{array} $$
((77))
$$\begin{array}{*{20}l} & C_{mab}^{BE22e} = \frac{\mu}{c_{y}} a_{m} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma, \end{array} $$
((78))

and with the LK termination:

$$\begin{array}{*{20}l} & C_{Pab}^{BE11e} = \rho c_{L} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma \end{array} $$
((79))
$$\begin{array}{*{20}l} & C_{Pab}^{BE22e} = \rho c_{T} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma, \end{array} $$
((80))

The contribution from the \({B^{m}_{E}}\) term to the stiffness matrix is:

$$\begin{array}{*{20}l} & K_{mab}^{BE11e} = \frac{2\mu+\lambda}{c_{x}} b_{m} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma \end{array} $$
((81))
$$\begin{array}{*{20}l} & H_{mab}^{BE11e} =-\frac{2\mu+\lambda}{c_{x}} \tilde b_{m} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma +\left(2\mu+\lambda\right) \int_{\Gamma_{E}} N_{a}N_{b,x}d\,\Gamma \end{array} $$
((82))
$$\begin{array}{*{20}l} & K_{mab}^{BE12e} =- \lambda \int_{\Gamma_{E}} N_{a}N_{b,y}d\,\Gamma \end{array} $$
((83))
$$\begin{array}{*{20}l} & K_{mab}^{BE21e} =- \mu \int_{\Gamma_{E}} N_{a}N_{b,y}d\,\Gamma \end{array} $$
((84))
$$\begin{array}{*{20}l} & H_{mab}^{BE22e} =-\frac{\mu}{c_{y}} \tilde b_{m} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma + \mu \int_{\Gamma_{E}} N_{a}N_{b,x}d\,\Gamma \end{array} $$
((85))
$$\begin{array}{*{20}l} & K_{mab}^{BE22e} = \frac{\mu}{c_{y}} b_{m} \int_{\Gamma_{E}} N_{a}N_{b} d\,\Gamma. \end{array} $$
((86))

The different contributions (denoted OM, BI and BE) to the matrices M, C, k, A, B, G and H are summed, and the resulting matrices are assembled by the standards FE assembly method, to form the system (47)–(54).

Computational aspects

Parameters of the absorbing condition

The parameters a m , b m , \(\tilde a_{m}\) and \(\tilde b_{m}\) in the recursive relations (15) and (16) are defined in a manner equivalent to their definition in the case of an absorbing boundary condition for elasticity defined on a single boundary (see for example [21]). Based on the analysis presented there we have:

$$ {a}_j= \cos {\theta}_j,\kern1em {b}_j=\frac{\overset{2}{ \sin }{\theta}_j}{T \cos {\theta}_j},\kern1em {\tilde{a}}_j= \cos {\overset{\sim }{\theta}}_j,\kern1em {\overset{\sim }{b}}_j=\frac{\overset{2}{ \sin }{\overset{\sim }{\theta}}_j}{T \cos {\overset{\sim }{\theta}}_j}\kern1em . $$
((87))

where T is the time period of interest and the “angles of incidence” θ j and \(\tilde {\theta }_{j}\) control the accuracy and stability of the particular double absorbing boundary. In previous studies we have used optimized parameter values (see, e.g., [14,21]), which enabled the control of both propagating and evanescent modes, and yielded excellent accuracy. In the numerical experiments described below we restrict ourselves to the case \(\theta _{j}=\tilde {\theta }_{j}=0\), which is the simplest possible set of parameters.

In the present scheme, using optimized parameters, which include small values of θ j , gave rise to some instability difficulties.

Discretization in time

The system (47)–(52) was solved using a fully implicit Newmark method, where all the degrees of freedom were obtained at once in each time step: u i in Ω and all the auxiliary variables in Ω L . The standard second-order implicit scheme with β=0.25 and γ=0.5 turned out to be unstable for this problem. Therefore, the procedure carried out in the acoustic case [13] was repeated here: we modify the Newmark parameters to be β L =0.36, γ L =0.7 inside the layer and along an additional row of elements in Ω I which is adjacent to Γ I (in this row of elements we calculate the derivative u i,x and use it in calculations performed in the layer). Throughout the rest of the domain Ω I , we take β L =0.25, γ L =0.5. The damped Newmark scheme is of first-order accuracy; however, this fact is not of much concern, since the layer may be regarded as a purely numerical construction, while in the region of interest the scheme has second-order accuracy in time. In theory, it is possible that the reduced accuracy could affect the reflection coefficient of the DAB, although no such effect has been observed in [13].

Stability

The stability of each of the examples described herein was investigated by two methods. The first, direct method consisted of running the scheme for very long periods of time (typically 500,000 time steps) and observing the behavior of the displacement at a point located inside the solution domain Ω. When the solution grew larger than a predefined value (typically 1), the solution was noted as unstable. This method gives an immediate indication of instability; however, it may judge a method to be stable when in fact divergence might occur at very large times (larger than 500,000 time steps).

We therefore used another method to determine stability, by solving an associated eigenvalue problem. In this regard, we have two options: to investigate the stability either of the semi-discrete problem (i.e., after FE discretization and before time discretization), or of the fully-discrete problem. The fully discrete problem is, of course, the problem that we actually solve and therefore its stability is crucial, but it is also of interest to investigate the stability of the semi-discrete problem.

We start by considering the semi-discrete problem. We rewrite the system (47)–(52) as a general second order system:

$$ \hat{\boldsymbol{M}} \ddot{\boldsymbol{u}} + \hat{\boldsymbol{C}} \dot{\boldsymbol{u}} + \hat{\boldsymbol{K}} \boldsymbol{u} =\hat {\boldsymbol{f}} $$
((88))

where \(\hat {\boldsymbol {M}}\) is the global mass matrix consisting of \(\boldsymbol {M}^{11}_{m}\) and \(\boldsymbol {M}^{22}_{m}, \hat {\boldsymbol {C}}\) is the global damping matrix consisting of \(\boldsymbol {C}^{11}_{m}, \boldsymbol {C}^{22}_{m}, \boldsymbol {A}^{11}_{m}, \boldsymbol {A}^{22}_{m}, \boldsymbol {G}^{11}_{m}\) and \(\boldsymbol {G}^{22}_{m}\), and \(\hat {\boldsymbol {K}}\) is the global stiffness matrix consisting of \(\boldsymbol {K}^{11}_{m}\), \(\boldsymbol {K}^{12}_{m}\), \(\boldsymbol {K}^{21}_{m}\), \(\boldsymbol {K}^{22}_{m}\), \(\boldsymbol {B}^{11}_{m}\), \(\boldsymbol {B}^{22}_{m}\), \(\boldsymbol {H}^{11}_{m}\) and \(\boldsymbol {H}^{22}_{m}\). The stability of this semi-discrete system is determined from the solution of the corresponding quadratic eigenvalue problem

$$ \left(\Lambda^{2}\hat{\boldsymbol{M}} + \Lambda \hat{\boldsymbol{C}} + \hat{\boldsymbol{K}}\right) \boldsymbol{u} = \mathbf{0} $$
((89))

where stability is established if all the eigenvalues Λ satisfy Re(Λ)≤0.

Now we consider the fully-discrete scheme. The fully discrete form of (88) is given in [31]:

$$\begin{array}{*{20}l} & \boldsymbol{M}^{*}= \hat{\boldsymbol{M}} + \gamma \Delta t \hat{\boldsymbol{C}} + \beta \Delta t^{2} \hat{\boldsymbol{K}} \end{array} $$
((90))
$$\begin{array}{*{20}l} & \boldsymbol{C}^{*}= -2\hat{\boldsymbol{M}}+\left(1-2\gamma \right) \Delta t \hat{\boldsymbol{C}} + \left(1/2-2\beta+\gamma \right) \Delta t^{2} \hat{\boldsymbol{K}} \end{array} $$
((91))
$$\begin{array}{*{20}l} & \boldsymbol{K}^{*}= \hat{\boldsymbol{M}} - \left(1-\gamma \right)\Delta t \hat{\boldsymbol{C}} + \left(1/2+ \beta-\gamma \right) \Delta t^{2} \hat{\boldsymbol{K}} \end{array} $$
((92))
$$\begin{array}{*{20}l} & \boldsymbol{M}^{*} \boldsymbol{u}_{n+1} + \boldsymbol{C}^{*} \boldsymbol{u}_{n} + \boldsymbol{K}^{*} \boldsymbol{u}_{n-1} + \Delta t^{2} \boldsymbol{f}^{*}_{n} = \mathbf{0} \end{array} $$
((93))

where u n is the value of u at time step n, and \(\boldsymbol {f}^{*}_{n}\) is some expression for the external force vector at time step n. In this case we solve the eigenvalue problem

$$ \left(\Lambda^{*2} \boldsymbol{M}^{*} + \Lambda^{*} \boldsymbol{C}^{*} + \boldsymbol{K}^{*}\right) \boldsymbol{u} = \mathbf{0} $$
((94))

and require for stability that all the eigenvalues Λ satisfy |Λ |≤1.

Our analysis based on (89) and (94), and our numerical experiments, yielded the following results for the DAB scheme described above. The semi-discrete formulation is long-time unstable, namely there are eigenvalues of (89) with a (small) positive real part. However, the fully-discrete problem is stable; this is obtained by both running the scheme for long times and through the eigenvalue analysis of (94). These facts give us a view of the stability properties of DAB at the various levels. At the continuous level, DAB for elastodynamics is believed to be stable. In [14] a well-posedness proof was provided for the DAB scheme, albeit for the acoustic problem, and in [20] stability was proved for an elastodynamics formulation using the same high-order ABC that the DAB in the present paper is based on. Since, as noted above, the semi-discrete problem is found to be unstable, we conclude that the FE formulation destabilizes the DAB scheme. Luckily, the dissipative time discretization that we employ regains stability. Of course, it would have been better to find a FE formulation that would maintain the stability of the continuous DAB formulation.

Achieving robust stability is still an issue for the present scheme. Unfortunately, the scheme loses its long-time stability when either of the following changes is applied: (a) replacing of the periodic boundary conditions on the north and south boundaries by traction-free conditions; (b) using non-zero (but not too large) coefficients for the low-order terms in order to capture efficiently evanescent waves; (c) using small computational parameter values \(a_{m}^{*}\) and/or \(b_{m}^{*}\), such as the optimal or quasi-optimal parameter values proposed in [14]; (d) taking much softer material parameters, i.e., much slower medium wave speeds. All these instabilities occur in a homogeneous elastic medium as well as in a layered one, but do not occur in the acoustic case. Research is underway to overcome these difficulties.

Accuracy

In order to measure the accuracy of the DAB method in the elastic setting, the solution in the truncated waveguide was compared to a reference solution u h,ref, which was constructed in such a manner that the waves generated in Ω I do not encounter a truncation boundary at x E . This was done in one of two ways: (a) We solve the same problem in a waveguide in which the truncation boundary is located much farther from the source of the wave, so that waves reflected from the far truncation boundary do not reach back to x=x E during the simulation. Thus, the difference between the solution u h computed in Ω and the reference solution u h,ref is approximately the DAB error. (b) If extending the waveguide significantly is computationally prohibitive (as in the case of the numerical example of the thin layer given below), one can extend the waveguide only slightly, and apply a high order version of the DAB method at the end of the extended domain. In such a case there will be reflections from the truncated boundary in the extended waveguide that will penetrate to the left of x E , but such reflections will be small, and the obtained difference between u h computed in Ω and the reference solution u h,ref will be an upper bound of the approximate DAB error.

We thus define the error measure

$$ E = \frac{\left\|\boldsymbol{u}^{h}-\boldsymbol{u}^{h,ref}\right\|_{\Omega_{I}\times[0,T]}} {\left\|\boldsymbol{u}^{h,ref}\right\|_{\Omega_{I}\times[0,T]}}\, $$
((95))

where \(\|\cdot \|_{\mathcal {M}}\) is the l 2 norm calculated on the manifold , and T is the simulation time. Note that the errors are measured outside of Ω L . In one case we shall also consider the evolution of the relative error in time, to which end we define

$$ (t) = \frac{\left\|\boldsymbol{u}^{h}-\boldsymbol{u}^{h,ref}\right\|_{\Omega_{I}}(t)}{\sqrt{A\left(\Omega_{I}\right)}}\, $$
((96))

where A(Ω I ) is the area of the domain Ω I .

In all of the numerical examples presented in Section “3”, the reference solution was calculated using method (b) outlined above. In the case of a homogeneous medium we also calculated the reference solution using method (a), thus satisfying ourselves that method (b) is indeed a legitimate method for error estimation.

Results and discussion

Initial pulse in a homogeneous medium

The discrete DAB scheme was tested on a numerical example of a waveguide with b=3, x I =10. The mesh of the solved problem was composed of 60×200 square elements, with n L elements across the width of Ω L , where the side length of each square element was h=0.05. Since the DAB layer extended n L elements beyond x E , the total length of the waveguide varied according to x E =x I +h·n L . We take λ=1, μ=1 (which corresponds to Poisson’s ratio ν=0.25), and ρ=1. The time step is taken as Δ t=0.005.

The C 1-continuous initial conditions used in the all the test runs (unless specified otherwise) of the first example problem were

$$\begin{array}{*{20}l} & u_{x0}(x,y)= \begin{cases} \left[(x-8.5)^{2}-1\right]^{2}\sin\left(\frac{2\pi y}{b}\right) & \text{for}\quad 7.5\leqslant x\leqslant 9.5 \\ 0 & \text{otherwise} \end{cases} \end{array} $$
((97))
$$\begin{array}{*{20}l} & u_{y0}(x,y)=0 \end{array} $$
((98))
$$\begin{array}{*{20}l} & v_{x0}(x,y)=0 \end{array} $$
((99))
$$\begin{array}{*{20}l} & v_{y0}(x,y)=0. \end{array} $$
((100))

The DAB scheme was tested for accuracy by running the scheme for different values of P and n L , and comparing the results to a waveguide 2.5 times longer than x E . The results without damping (ξ=0) are shown in Figures 3 and 4 for two discretization densities: the one mentioned above (Figure 4) and one which is twice coarser in space and time (Figure 3).

Figure 3
figure 3

Initial pulse in a homogeneous medium, no Rayleigh damping, discretization with h=0.1 and Δ t=0.01: error E as a function of the DAB order P, for different values of the DAB thickness.

Figure 4
figure 4

Initial pulse in a homogeneous medium, no Rayleigh damping, discretization with h=0.05 and Δ t=0.005: error E as a function of the DAB order P, for different values of the DAB thickness.

As can be seen from the accuracy plots, the high-order DAB acts as an accurate absorbing boundary, reducing the L 2 error in time and space as P increases, up to the discretization level. With P=0, which amounts to the use of the LK condition on Γ E , the error is large (around 20%), but with P=6 the error is reduced by an order of magnitude. The error also decreases as n L is increased.

The cause for the dependence of the error on n L is the presence of evanescent modes in the model, which are uncontrolled due to our simple choice of the parameters θ j =~ θ j =0.

Let us compare the error curves corresponding to n L =2,4,6 in Figure 3 (the coarse-mesh results) to the error curves corresponding to n L =4,8,12 in Figure 4 (the fine-mesh results), respectively. The three pairs of values of n L correspond to the same DAB thickness (in terms of physical distance), since the elements of the coarse mesh are twice as large as those of the fine mesh (the h ratio is 2). The error curves behave similarly in terms of reduction with P, except that the errors corresponding to the fine mesh are much smaller. For example, for large P and for the thickest DAB, the coarse-mesh model yields an error of about 7·10−3 whereas the fine-mesh model yields an error of about 2·10−3. The reduction factor of 3.5 is a bit smaller than the theoretical reduction factor of 4 of the discretization error (expected from the second-order accuracy in space and time). What the theory does not take into account is the discretization error associated with the DAB equations themselves, which apparently contributes slightly to the total error.

One may wonder why the error “saturates” at a certain level for large values of P, rather than approaching zero as P; after all, the reference solution uses the same mesh density and time-step size as the solution in Ω. The explanation is as follows. One should distinguish between three different exact ABCs associated with P: (a) the exact ABC at the continuous level, denoted B u=0; (b) the ABC obtained by discretizing the continuous exact ABC, denoted B h u h=0; and (c) the exact ABC at the discrete level, denoted G h u h=0. It should be noted that B h and G h are not the same. In other words, taking the exact ABC at the continuous level and discretizing it is not the same as taking the exact ABC at the discrete level. The reference solution is equivalent to a solution of the truncated problem with the ABC G h u h=0. On the other hand, when we use our numerical scheme and take P, our solution satisfies B h u h=0. The difference between the two solutions is of the order of the discretization error.

This was confirmed by a direct calculation in [16].

Therefore, as we increase P, we cannot go much below the discretization error. This explains why the error “saturates” in the error graph of Figure 4 and all the subsequent figures.

Next we investigate the effect of damping on the accuracy of the method. Damping coefficients A K and A M that are characteristic of soil need to be selected. We select values for these coefficients from the literature, specifically in the range of values suggested for soils in [32]. We compare the magnitude of the solution vector at a point in the middle of Γ I with and without mass-proportional and stiffness-proportional damping for n L =4 and P=10 as shown in Figure 5. It is clear that the added damping reduces somewhat the magnitude of the solution, and smooths out its oscillations; the larger the damping the smaller the oscillations. However, the small stiffness-proportional coefficient A K from the range of values specified in [32] has only a small effect on the solution for relatively short times.

Figure 5
figure 5

Initial pulse in a homogeneous medium, with Rayleigh damping: |u| at a single point for various values of the damping parameters.

In Figure 6 we show the accuracy of the absorbing boundary for the damped case for n L =4 for different values of P. The damping reduces somewhat the error of the DAB scheme, although for the suggested values of the damping the reduction is small. For large values of P the stiffness-proportional damping increases the error somewhat when compared to the error of the scheme without stiffness-proportional damping and with the same mass-proportional damping.

Figure 6
figure 6

Initial pulse in a homogeneous medium, with Rayleigh damping: error E as a function of the DAB order P, for different values of the damping parameters, with n L =4.

In the sequel we examine other problems, in which we shall only consider the undamped case, i.e., with A M =0, A K =0.

Initial pulse in a two-layer medium

In this example we use an overall geometry, initial and boundary conditions that are similar to the ones used in the first example. In the present case, however, the medium is divided into two horizontal layers with different material properties. The interface between the layers is the line y=b/2=1.5. In the upper layer (y>1.5), which we denote medium 1, we set λ 1=1 and μ 1=1. In the lower layer (y<1.5), which we denote medium 2, we set λ 2=1, while the value of μ 2 will vary in the following experiments. We take ρ=1 everywhere. The procedure was carried out both for the solved waveguide and for the reference waveguide.

Snapshots comparing the solution at various times to the solution in the reference waveguide for μ 2=0.75, P=10 and n L =4 are shown in Figure 7. Each snapshot consists of two subplots: the upper one depicts the reference solution \(u_{i}^{h,ref}\) over the extended domain, and the lower one is the computed solution \({u_{i}^{h}}\) in the truncated domain. The x component of the solutions is shown in the left column and the y component is shown in the right column. The DAB layer is indicated by showing the vertical boundary Γ I .

Figure 7
figure 7

Pulse propagation in a two-layer medium, snapshots of solution: (a)(b) u x and u y at t=0; (c)(d) u x and u y at t=2; (e)(f) u x and u y at t=4; (g)(h) u x and u y at t=10.

The snapshots show good agreement between the reference and the actual solutions. The obtained accuracy in this case is illustrated better by the error plot given in Figure 8. It can be seen that the behavior of the error in this case does not differ much from the behavior in the homogeneous case.

Figure 8
figure 8

Two-layer medium, μ 2=0.75: error E as a function of the DAB order P, for different values of the DAB thickness.

In Figures 9 and 10 we show the error as a function of P for different values of the μ 2 parameter for n L =4 and n L =12, respectively. It can be seen that the presence of two layers with different material properties in the domain does not change the DAB error drastically when we follow the guidelines outlined in the section “Layered media” for the selection of the appropriate constants in the definition of the boundary condition. For large value of P, the error is larger in the case where the material properties of the two layers differ significantly, namely the case μ 2/μ 1=3, than in the other cases. We also note that by enlarging the DAB layer from n L =4 to n L =12, the accuracy is improved, since some of the evanescent waves are better represented.

Figure 9
figure 9

Two layer medium: error E as a function of the DAB order P, for different values of the material properties, with n L =4.

Figure 10
figure 10

Same as Figure 9, with n L =12.

Persistent couple

We consider a problem involving the same waveguide and the same mesh as in the previous examples b=3, x I =10, with homogeneous material properties λ=1, μ=1, but where instead of initial conditions driving the solution, we impose the following two persistent point forces F 1 and F 2 inside Ω I :

$$\begin{array}{*{20}l} & \boldsymbol{F}_{1}(x=9.5, y=1, t) = (1,0) \end{array} $$
((101))
$$\begin{array}{*{20}l} & \boldsymbol{F}_{2}(x=9.5, y=2, t) = (-1,0), \end{array} $$
((102))

thus a couple is imposed inside the waveguide, near Γ I . The plot of the space-time error as a function of the DAB order P is shown in Figure 11 for different values of n L .

Figure 11
figure 11

Persistent couple problem: error E as a function of the DAB order P, for different values of the DAB thickness.

It can be seen from the plot that the present problem is in a certain sense more difficult than the problem with the initial driving pulse, as the error lines drop sharply only after P=4 in this case, and the error level to which the graphs reach for large P is somewhat higher than in the case of the initial pulse. Nevertheless, a decrease of almost two orders of magnitude in the error is obtained for all values of n L .

Three layers

In the last example we consider a setup of a uniform medium through which passes a thin horizontal layer made from a different material. We consider a domain with a width of 9 and a length of 2.5 plus a DAB layer. The mesh is composed of square elements with the side length of 0.05. We apply a DAB layer with a thickness of 4 elements. The thin horizontal layer is located at 4<y<5. As before, we apply periodic boundary conditions on Γ S and Γ N .

We define \(r = \sqrt {(x-0.5)^{2}+(y-3.5)^{2}}\), a=0.3, and introduce the following source function:

$$\begin{array}{*{20}l} & u_{x}(x,y,t=0) = \left\{\begin{array}{ll} 2\left(\frac{r}{a}\right)^{3} - 3\left(\frac{r}{a}\right)^{2} + 1 & \text{ for}\quad 0\leqslant r\leqslant a \\ 0 & \text{ otherwise} \end{array}\right. \end{array} $$
((103))
$$\begin{array}{*{20}l} & u_{y}(x,y,t=0) = 0. \end{array} $$
((104))

In the thin layer, we set λ=μ=1, whereas in the surrounding media we set λ=μ=2. For this problem, Figure 12 shows the L 2 spatial error as a function of time, e(t), for three values of P. It is apparent that the LK ABC (which corresponds to P=0) yields a large error in this case. The DAB error converges very quickly with P to an error significantly smaller.

Figure 12
figure 12

Three layer problem: error e(t) for three values of DAB order P.

In Figure 13 we compare three solutions. We fix the Poisson’s ratio at ν=0.25 by taking λ=μ in all the phases. In all cases the background medium has the properties λ=μ=2. The left subplots correspond to a thin layer with λ=μ=1, the middle subplots correspond to a thin layer with λ=μ=2, i.e., a homogeneous medium, and the right subplot correspond to a thin layer with λ=μ=4. We take ρ=1 everywhere. This creates a hierarchy for the wave velocities c L and c T : in the layer with λ=μ=1 (left subplots) these velocities are smaller by a factor of \(\sqrt {2}\) than those in the surrounding regions, while in the layer with λ=μ=4 (right subplots) they are larger by a factor of \(\sqrt {2}\) than the velocities in the surrounding regions. Since the wave source is located below the thin layer, the latter serves as a wave decelerator or accelerator, as the case may be. This effect is clearly apparent in the figure; see, e.g., Figure 13(e).

Figure 13
figure 13

A three-layer medium, snapshots of solution. Solutions in three media are compared, where Poisson’s ratio is held fixed at ν=0.25 in all the phases. In all cases the background medium has the properties λ=μ=2. The left subplots correspond to a thin layer with λ=μ=1, the middle subplots correspond to a thin layer with λ=μ=2, i.e., a homogeneous medium, and the right subplot correspond to a thin layer with λ=μ=4. Shown are snapshots of: (a) u x at t=0; (b) u y at t=0; (c) u x at t=1; (d) u y at t=1; (e) u x at t=2; (f) u y at t=2; (g) u x at t=3; (h) u y at t=3; (i) u x at t=4; (j) u y at t=4; (k) u x at t=5; (l) u y at t=5.

The solution is non-symmetric despite the periodic boundary conditions, since the wave source is not centrally located but lies slightly below the horizontal layer. The solution in the heterogeneous media is characterized by much more reverberation. No clear spurious reflection of waves is observed. It should be remarked that each plot has its own scaled color map; therefore wave intensities cannot be compared between subplots in this figure.

The errors corresponding to these solutions (not shown here), with respect to a reference solution obtained with a longer domain, have the same characteristics as those shown in the previous examples.

Conclusions

In this paper the Double Absorbing Boundary (DAB) method for solving unbounded domain problems was applied to elastodynamics problems, in conjunction with the FE method. This method shares some of the properties of using a high-order ABC on a single boundary and of PML, but has some relative advantages with respect to both.

Several important extensions will be considered in our future work. Three of them are:

  • Some stability issues still need to be resolved. In particular, when the periodic boundary conditions used in this paper were replaced by some physical conditions, like traction-free conditions, numerical stability was lost. It seems that this difficulty, which is currently under investigation, is not associated with the original DAB formulation, at the continuous level, but with the semi-discrete formulation.

  • The behavior of DAB in the presence of corners joining two straight artificial boundaries will be investigated. The use of DAB on a boundary with corners is expected to be as straight forward as the case with PML (roughly speaking, a “cross product” of the x and y formulations), and thus to be free of the difficulties associated with high order ABCs on a single boundary in the presence of corners.

  • Extension to an anisotropic medium also seems possible, although it would require a more involved adaptation. It would be especially interesting to see how DAB behaves in those cases of anisotropy where standard PMLs are unstable.

References

  1. Hagstrom T (1999) Radiation Boundary Conditions for the Numerical Simulation of Waves. Acta Numerica 8: 47–106.

    Article  MathSciNet  Google Scholar 

  2. Givoli D (2004) High-Order Local Non-Reflecting Boundary Conditions: A Review. Wave Motion 39: 319–326.

    Article  MATH  MathSciNet  Google Scholar 

  3. Givoli D (2008) Computational Absorbing Boundaries. In: Marburg S Nolte B (eds)Computational Acoustics of Noise Propagation in Fluids, Chapter 5, 145–166.. Springer, Berlin.

    Google Scholar 

  4. Bermudez A, Hervella-Nieto L, Prieto A, Rodriguez R (2010) Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems. Archives Comput Meth Engng 17: 77–107.

    Article  MATH  Google Scholar 

  5. Bérenger JP (1994) A Perfectly Matched Layer for the Absorption of Electromagnetic Waves. J. Comput. Phys. 114: 185–200.

    Article  MATH  MathSciNet  Google Scholar 

  6. Collino F (1993) High Order Absorbing Boundary Conditions for Wave Propagation Models. Straight Line Boundary and Corner Cases. In: Kleinman R et al. (eds)Proc. 2nd Int. Conf. on Mathematical & Numerical Aspects of Wave Propagation, 161–171.. SIAM, Delaware.

    Google Scholar 

  7. Engquist B, Majda A (1979) Radiation Boundary Conditions for Acoustic and Elastic Wave Calculations. Comm Pure Appl Math 32: 313–357.

    Article  MATH  MathSciNet  Google Scholar 

  8. Bayliss A, Turkel E (1980) Radiation Boundary Conditions for Wave-Like Equations. Comm Pure Appl Math 33: 707–725.

    Article  MATH  MathSciNet  Google Scholar 

  9. Rabinovich D, Givoli D, Bécache E (2010) Comparison of High-order Absorbing Boundary Conditions and Perfectly Matched Layers in the Frequency Domain. Int J Num Meth Biomed Engng (Formerly Commun Numer Meth Engng) 26: 1351–1369.

    Article  MATH  Google Scholar 

  10. Asvadurov S, Druskin V, Guddati M, Knizherman L (2003) On optimal finite difference approximation of PML. SIAM J Numer Anal 41: 287–305.

    Article  MATH  MathSciNet  Google Scholar 

  11. Hagstrom T, Warburton T (2004) A New Auxiliary Variable Formulation of High-Order Local Radiation Boundary Conditions: Corner Compatibility Conditions and Extensions to First Order Systems. Wave Motion 39: 327–338.

    Article  MATH  MathSciNet  Google Scholar 

  12. Rabinovich D, Givoli D, Bielak J, Hagstrom T (2011) A Finite Element Scheme with a High Order Absorbing Boundary Condition for Elastodynamics. Comput. Meth. Appl. Mech. Engng 200: 2048–2066.

    Article  MATH  MathSciNet  Google Scholar 

  13. Hagstrom T, Givoli D, Rabinovich D, Bielak J (2014) The Double Absorbing Boundary Method. J. Comput. Phys 259: 220–241.

    Article  MathSciNet  Google Scholar 

  14. Baffet D, Hagstrom T, Givoli D (2014) Double Absorbing Boundary Formulations for Acoustics and Elastodynamics. SIAM J Sci Comput 36: A1277–A1312.

    Article  MATH  MathSciNet  Google Scholar 

  15. Hagstrom T, Warburton T (2009) Complete Radiation Boundary Conditions: Minimizing the Long Time Error Growth of Local Methods. SIAM J. Numer. Anal 47: 3678–3704.

    Article  MATH  MathSciNet  Google Scholar 

  16. Hagstrom T, Mar-Or A, Givoli D (2008) High-Order Local Absorbing Conditions for the Wave Equation: Extensions and Improvements. J Comput Phys 227: 3322–3357.

    Article  MATH  MathSciNet  Google Scholar 

  17. Bécache E, Givoli D, Hagstrom T (2010) High Order Absorbing Boundary Conditions for Anisotropic and Convective Wave Equations. J. Comput. Phys 229: 1099–1129.

    Article  MATH  MathSciNet  Google Scholar 

  18. Hagstrom T, Bécache E, Givoli D, Stein K (2012) Complete Radiation Boundary Conditions for Convective Waves. Commun. Comput. Phys 11: 610–628.

    MathSciNet  Google Scholar 

  19. Mar-Or A, Givoli D (2009) High Order Global-Regional Model Interaction: Extension of Carpenter’s Scheme. Int. J. Numer. Meth. Engng 77: 50–74.

    Article  MATH  Google Scholar 

  20. Baffet D, Bielak J, Givoli D, Hagstrom T, Rabinovich D (2012) Long-Time Stable High-Order Absorbing Boundary Conditions for Elastodynamics. Comput. Meth. Appl. Mech. Engng241–244: 20–37.

    Article  MathSciNet  Google Scholar 

  21. Rabinovich D, Givoli D, Hagstrom T, Bielak J (2013) Stress-Velocity Complete Radiation Boundary Conditions. J. Comput. Acoust., 21: 1350003–1-38.

    Article  MathSciNet  Google Scholar 

  22. Lysmer J, Kuhlemeyer RL (1969) Finite Dynamic Model for Infinite Media. J Eng Mech Div ASCE 95: 859–877.

    Google Scholar 

  23. Bamberger A, Chalindar B, Joly P, Roberts JE, Teron JL (1988) Absorbing Boundary Conditions for Rayleigh Waves. SAIM J Sci Stat Comput 9: 1016–1049.

    Article  MATH  MathSciNet  Google Scholar 

  24. Bao HS, Bielak J, Ghattas O, Kallivokas LF, O’Hallaron DR, Shewchuk JR, Xu JF (1998) Large-Scale Simulation of Elastic Wave Propagation in Heterogeneous Media on Parallel Computers. Comput. Meth. Appl. Mech. Engng 152: 85–102.

    Article  MATH  Google Scholar 

  25. Bielak J, Ghattas O, Kim EJ (2005) Parallel Octree-Based Finite Element Method for Large-Scale Earthquake Ground Motion Simulation. Comput Model Engng Sci 10: 99–112.

    MATH  MathSciNet  Google Scholar 

  26. Day SM, Graves R, Bielak J, Dreger D, Larsen S, Olsen KB, Pitarka A, Ramirez-Guzman L (2008) Model for Basin Effects on Long-Period Response Spectra in Southern California. Earthquake Spectra 24: 257–277.

    Article  Google Scholar 

  27. Duru K (2014) A Perfectly Matched Layer for the Time-Dependent Wave Equation in Heterogeneous and Layered Media. J. Comput. Phys 257: 757–781.

    Article  MathSciNet  Google Scholar 

  28. Hanasoge SM, Komatitsch D, Gizon L (2010) An Absorbing Boundary Formulation for the Stratified, Linearized, Ideal MHD Equations Based on an Unsplit, Convolutional Perfectly Matched Layer. Astronomy & Astrophys 522: A87/1–A87/8.

    Article  Google Scholar 

  29. Collino F, Tsogka C (2001) Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media. Geophys 66: 294–307.

    Article  Google Scholar 

  30. Hagstrom T (1240) High-Order Radiation Boundary Conditions for Stratified Media and Curvilinear Coordinates. J. Comput. Acoust., 20: 002–1240018.

  31. Hughes TJR (1987) The Finite Element Method. Prentice Hall, Englewood Cliffs, N.J.

    MATH  Google Scholar 

  32. Ju S-H, Ni S-H (2007) Determining Rayleigh Damping Parameters of Soils for Finite Element Analysis. Int J Numer. Anal Meth Geomech 31: 1239–1255.

    Article  MATH  Google Scholar 

Download references

Acknowledgments

This work was supported by the US-Israel Binational Science Foundation (BSF), grant number 890020 (Technion number 2011303). The work of DG was also supported by the fund provided through the Lawrence and Marie Feldman Chair in Engineering; that of TH by NSF grant DMS-1418871 and ARO grant W911NF-09-1-0344; and that of JB by the U.S. National Science Foundation Award No. NSF OCI-0749227. Any conclusions or recommendations in the paper are those of the authors and do not necessarily represent the views of the BSF, NSF or ARO.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Dan Givoli.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

DG and TH developed the theory of the DAB method for elastodynamics, JB provided input on practical and computational aspects as well as on the structure and content of the manuscript, DR performed the coding, carried out the numerical investigations, and wrote most of the manuscript. All authors read and approved the final manuscript.

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rabinovich, D., Givoli, D., Bielak, J. et al. The double absorbing boundary method for elastodynamics in homogeneous and layered media. Adv. Model. and Simul. in Eng. Sci. 2, 3 (2015). https://doi.org/10.1186/s40323-015-0026-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40323-015-0026-8

Keywords