 Research Article
 Open Access
 Published:
The double absorbing boundary method for elastodynamics in homogeneous and layered media
Advanced Modeling and Simulation in Engineering Sciences volume 2, Article number: 3 (2015)
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 highorder 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 highorder ABC on two parallel artificial boundaries, which are a small distance apart, and thus form a thin nonreflecting 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 timedependent elastic waves in homogeneous and layered media. The equations are written in secondorder 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 absorbingboundary type schemes have been those based on the use of a highorder 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 [14].
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 adhoc damping and stretching profiles.
The first highorder 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]. Highorder 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 highorder 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 highorder ABC scheme, the order of the ABC is simply an input parameter. The high derivatives that initially appear when designing a highorder ABC are eliminated by introducing auxiliary variables ϕ _{ j } on the boundary.
Recently the two approaches — highorder 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 highorder ABC proposed by Hagstrom and Warburton [11].
As mentioned above, each of the two classes of techniques has relative advantages. One major disadvantage of highorder 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 highorder 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 wellresolved. By contrast, in the case of highorder 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 highorder 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 highorder ABC, but enjoys some of the advantages that each of them lacks. In the new method, called the Double Absorbing Boundary (DAB) method, a highorder 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 highorder 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 linearlyincreasing cost. The numerical properties of DAB, like accuracy, stability and sensitivity to discretization, are similar to those of a highorder 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 wellposedness proof was provided for the DAB scheme for the acoustics problem written in secondorder form. The energy method was employed to obtain uniformintime estimates of the norm of the solution and the auxiliary functions, thus establishing the wellposedness and asymptotic stability of the DAB formulation. In addition, in [14] the DAB was applied to problems in 2D isotropic elastodynamics, written in firstorder conservation form. The problem was discretized using the LaxWendroff finite difference scheme.
Although DAB is a general approach, and in principle can be used with any highorder ABC applied on the double boundary, in [13,14] the ABCs of the form proposed by Hagstrom and Warburton (HW) 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 uniformintime error estimate for both propagating and decaying waves. The HW 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 [1618], and where exterior sources are present (nesting) [19].
Recently, an HW type ABC was applied to problems in elastodynamics [20,21]. This is the first known highorder ABC for elastodynamics that is longtime stable. The key to obtaining stability turns out to be the use of the LysmerKuhlemeyer (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., [2426]. In the stable highorder ABC for elastodynamics, all the recursive relations except the last one are scalar in nature (as in the HW 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 continuouslystratified 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 unboundeddomain problem of elastodynamics
We consider a twodimensional semiinfinite 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.
In the waveguide we consider the equation of elastodynamics with mass and stiffnessproportional damping, i.e.,
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 bodyforce 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 massproportional and stiffnessproportional damping coefficients, which together define the Rayleigh damping that exists in the system.
We define the wave speeds c _{ L } and c _{ T } as
where λ and μ are the Lamé constants.
Boundary conditions are specified on the three boundaries: on Γ _{ N } and Γ _{ S } we impose periodic boundary conditions:
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
Initial conditions are also given, i.e.,
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 semiinfinite domain by introducing the artificial boundary Γ _{ E }, located at x=x _{ E } and spanning 0≤y≤b. Slightly to the west of Γ _{ E } we set an interface denoted Γ _{ I }, located at x=x _{ I }, with 0≤y≤b. 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 subdomains: 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 semiinfinite problem in that domain. Thus, Ω _{ L } will act as an absorbing or nonreflecting 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:
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:
On Γ _{ S } and Γ _{ N } (or more precisely the parts of these boundaries that are in the layer, i.e., for x _{ I }≤x≤x _{ E } and 0≤y≤b) we apply the same (periodic) boundary conditions as for u _{ i }, i.e.,
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 loworder terms:
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 highorder 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 HW 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 LysmerKuhlemeyer (LK) termination condition on Γ _{ E }:
Figure 2 illustrates the “ladder” structure of the DAB scheme, namely the flow of information on the two boundaries bounding the DAB 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:
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 nth layer, and denoting c=c _{ x }=c _{ y }, we obtain for each layer n:
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
all the recursive relations (21) for the different layers can be replaced by the following recursive relation written out for a reference layer:
Now, when we consider the termination conditions for the different layers, we take into account that
Therefore, on the interface between the layers we obtain from the second termination condition (18) that:
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:
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:
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
where we have
We rewrite the recursive relations (15) and (16) in the following form:
On the boundary Γ _{ I }:
On the boundary Γ _{ E }:
and (17) and (18) serve as the termination conditions:
in which case, by substituting (34)–(37) in the terms \({B^{m}_{E}}\), \({B^{m}_{I}}\) in (32) and (33), these terms become:
and (38) and (39) are substituted directly into (33), with m=P:
Here and elsewhere, the xderivatives on the boundaries Γ _{ I },Γ _{ E } are calculated in a onesided 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.
Semidiscrete 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 finitedimensional approximations
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 globallevel 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:
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 elementlevel 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 elementlevel 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:
For m=1,…,P−1
and also:
with the initial conditions
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):
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:
The contribution from the \({B^{m}_{I}}\) term to the stiffness matrixes is:
Similarly, the contribution from the \({B^{m}_{E}}\) term (denoted by the superscript BE) to the damping matrices is (m=0,…,P−1):
and with the LK termination:
The contribution from the \({B^{m}_{E}}\) term to the stiffness matrix is:
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:
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 secondorder 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 firstorder 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 secondorder 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 semidiscrete problem (i.e., after FE discretization and before time discretization), or of the fullydiscrete 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 semidiscrete problem.
We start by considering the semidiscrete problem. We rewrite the system (47)–(52) as a general second order system:
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 semidiscrete system is determined from the solution of the corresponding quadratic eigenvalue problem
where stability is established if all the eigenvalues Λ satisfy Re(Λ)≤0.
Now we consider the fullydiscrete scheme. The fully discrete form of (88) is given in [31]:
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
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 semidiscrete formulation is longtime unstable, namely there are eigenvalues of (89) with a (small) positive real part. However, the fullydiscrete 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 wellposedness 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 highorder ABC that the DAB in the present paper is based on. Since, as noted above, the semidiscrete 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 longtime stability when either of the following changes is applied: (a) replacing of the periodic boundary conditions on the north and south boundaries by tractionfree conditions; (b) using nonzero (but not too large) coefficients for the loworder 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 quasioptimal 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
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
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
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).
As can be seen from the accuracy plots, the highorder 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 coarsemesh results) to the error curves corresponding to n _{ L }=4,8,12 in Figure 4 (the finemesh 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 coarsemesh model yields an error of about 7·10^{−3} whereas the finemesh 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 secondorder 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 timestep 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 massproportional and stiffnessproportional 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 stiffnessproportional coefficient A _{ K } from the range of values specified in [32] has only a small effect on the solution for relatively short times.
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 stiffnessproportional damping increases the error somewhat when compared to the error of the scheme without stiffnessproportional damping and with the same massproportional damping.
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 twolayer 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 }.
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.
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.
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 }:
thus a couple is imposed inside the waveguide, near Γ _{ I }. The plot of the spacetime error as a function of the DAB order P is shown in Figure 11 for different values of n _{ L }.
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 {(x0.5)^{2}+(y3.5)^{2}}\), a=0.3, and introduce the following source function:
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.
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).
The solution is nonsymmetric 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 highorder 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 tractionfree 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 semidiscrete 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.
 2
Givoli D (2004) HighOrder Local NonReflecting Boundary Conditions: A Review. Wave Motion 39: 319–326.
 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.
 4
Bermudez A, HervellaNieto L, Prieto A, Rodriguez R (2010) Perfectly Matched Layers for TimeHarmonic Second Order Elliptic Problems. Archives Comput Meth Engng 17: 77–107.
 5
Bérenger JP (1994) A Perfectly Matched Layer for the Absorption of Electromagnetic Waves. J. Comput. Phys. 114: 185–200.
 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.
 7
Engquist B, Majda A (1979) Radiation Boundary Conditions for Acoustic and Elastic Wave Calculations. Comm Pure Appl Math 32: 313–357.
 8
Bayliss A, Turkel E (1980) Radiation Boundary Conditions for WaveLike Equations. Comm Pure Appl Math 33: 707–725.
 9
Rabinovich D, Givoli D, Bécache E (2010) Comparison of Highorder 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.
 10
Asvadurov S, Druskin V, Guddati M, Knizherman L (2003) On optimal finite difference approximation of PML. SIAM J Numer Anal 41: 287–305.
 11
Hagstrom T, Warburton T (2004) A New Auxiliary Variable Formulation of HighOrder Local Radiation Boundary Conditions: Corner Compatibility Conditions and Extensions to First Order Systems. Wave Motion 39: 327–338.
 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.
 13
Hagstrom T, Givoli D, Rabinovich D, Bielak J (2014) The Double Absorbing Boundary Method. J. Comput. Phys 259: 220–241.
 14
Baffet D, Hagstrom T, Givoli D (2014) Double Absorbing Boundary Formulations for Acoustics and Elastodynamics. SIAM J Sci Comput 36: A1277–A1312.
 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.
 16
Hagstrom T, MarOr A, Givoli D (2008) HighOrder Local Absorbing Conditions for the Wave Equation: Extensions and Improvements. J Comput Phys 227: 3322–3357.
 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.
 18
Hagstrom T, Bécache E, Givoli D, Stein K (2012) Complete Radiation Boundary Conditions for Convective Waves. Commun. Comput. Phys 11: 610–628.
 19
MarOr A, Givoli D (2009) High Order GlobalRegional Model Interaction: Extension of Carpenter’s Scheme. Int. J. Numer. Meth. Engng 77: 50–74.
 20
Baffet D, Bielak J, Givoli D, Hagstrom T, Rabinovich D (2012) LongTime Stable HighOrder Absorbing Boundary Conditions for Elastodynamics. Comput. Meth. Appl. Mech. Engng241–244: 20–37.
 21
Rabinovich D, Givoli D, Hagstrom T, Bielak J (2013) StressVelocity Complete Radiation Boundary Conditions. J. Comput. Acoust., 21: 1350003–138.
 22
Lysmer J, Kuhlemeyer RL (1969) Finite Dynamic Model for Infinite Media. J Eng Mech Div ASCE 95: 859–877.
 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.
 24
Bao HS, Bielak J, Ghattas O, Kallivokas LF, O’Hallaron DR, Shewchuk JR, Xu JF (1998) LargeScale Simulation of Elastic Wave Propagation in Heterogeneous Media on Parallel Computers. Comput. Meth. Appl. Mech. Engng 152: 85–102.
 25
Bielak J, Ghattas O, Kim EJ (2005) Parallel OctreeBased Finite Element Method for LargeScale Earthquake Ground Motion Simulation. Comput Model Engng Sci 10: 99–112.
 26
Day SM, Graves R, Bielak J, Dreger D, Larsen S, Olsen KB, Pitarka A, RamirezGuzman L (2008) Model for Basin Effects on LongPeriod Response Spectra in Southern California. Earthquake Spectra 24: 257–277.
 27
Duru K (2014) A Perfectly Matched Layer for the TimeDependent Wave Equation in Heterogeneous and Layered Media. J. Comput. Phys 257: 757–781.
 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.
 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.
 30
Hagstrom T (1240) HighOrder 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.
 32
Ju SH, Ni SH (2007) Determining Rayleigh Damping Parameters of Soils for Finite Element Analysis. Int J Numer. Anal Meth Geomech 31: 1239–1255.
Acknowledgments
This work was supported by the USIsrael 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 DMS1418871 and ARO grant W911NF0910344; and that of JB by the U.S. National Science Foundation Award No. NSF OCI0749227. 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
Affiliations
Corresponding author
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/.
About this article
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/s4032301500268
Received:
Accepted:
Published:
Keywords
 Double absorbing boundary
 Absorbing boundary condition
 Highorder
 Auxiliary variables
 Elastic waves
 Elastodynamics
 Layered media
 Finite elements