Three-dimensional hybrid asynchronous perfectly matched layer for wave propagation in heterogeneous semi-infinite media

Sijia Li, Michael Brun, Irini Djeran-Maigre, Sergey Kuznetsov

► To cite this version:

Sijia Li, Michael Brun, Irini Djeran-Maigre, Sergey Kuznetsov. Three-dimensional hybrid asynchronous perfectly matched layer for wave propagation in heterogeneous semi-infinite media. Comptes Rendus. Mécanique, 2020, 348 (12), pp.1003-1030. 10.5802/crmeca.59. hal-04102425

HAL Id: hal-04102425

https://hal.science/hal-04102425

Submitted on 1 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Sijia Li, Michael Brun, Irini Djeran-Maigre and Sergey Kuznetsov

Three-dimensional hybrid asynchronous perfectly matched layer for wave propagation in heterogeneous semi-infinite media

Volume 348, issue 12 (2020), p. 1003-1030

Published online: 19 January 2021

https://doi.org/10.5802/crmeca.59

<figure><img src="image_2.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_3.png" /></figure>

Les Comptes Rendus. Mécanique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org
e-ISSN : 1873-7234
Short paper / *Note*

# Three-dimensional hybrid asynchronous
perfectly matched layer for wave propagation
in heterogeneous semi-infinite media

Sijia Li*, a, b, Michael Brun c, Irini Djeran-Maigre a and Sergey Kuznetsov d

a Univ Lyon, INSA-Lyon, GEOMAS, 69621 Villeurbanne, France

b Structural Mechanics and Coupled Systems Laboratory, Conservatoire National des
Arts et Métiers, 292 rue Saint-Martin, 75141 Paris Cedex 03, France

c LEM3 CNRS UMR 7239, Université de Lorraine, 7 rue Félix Savart, 57070 METZ,
France

d Institute for Problems in Mechanics, 119526 Moscow, Russia

*E-mails*: sijia.li1@insa-lyon.fr (S. Li), Michael.Brun@insa-lyon.fr (M. Brun), irini.djeran-maigre@insa-lyon.fr (I. Djeran-Maigre), kuzn-sergey@yandex.ru (S. Kuznetsov)

**Abstract.** This paper presents an efficient hybrid asynchronous three-dimensional (3D) perfectly matched layer (PML) for modeling unbounded domains. The proposed unsplit explicit or implicit 3D PML formulation is implemented in the framework of a heterogeneous asynchronous time integrator. It is fully versatile in terms of time integrators and time step sizes according to partitions while conserving classical finite element formulations in the elastic domain without complex-valued stretched coordinates. Examples of a semi-infinite bar, Lamb's test, and a soil-structure interaction problem with PML-truncated semi-infinite heterogeneous media are investigated to illustrate the efficiency of the proposed PML in terms of accuracy and CPU time.

**Keywords.** Three-dimensional perfectly matched layers (PMLs), Unbounded domains, Wave propagation, Transient analysis, Subdomain coupling, Heterogeneous asynchronous time integrator (HATI).

Manuscript received 16th July 2020, revised 3rd November 2020, accepted 17th November 2020.

## 1. Introduction

One of the critical issues regarding the numerical simulation of wave propagation problems in unbounded domains using the finite element method is finding a suitable technique to simulate infinite media. The simplest way is to consider a very large extended numerical mesh. However, this approach involves a high computational time, in particular when long-time simulation is of interest. Hence, non-reflecting boundary conditions are required at the boundary of

<footnote>* Corresponding author.</footnote>
the truncated domain for mimicking infinite or semi-infinite media. Several kinds of artificial boundaries have been developed for numerical methods to avoid spurious waves being reflected at the boundary. They include infinite elements (Bettess [1], Houmat [2]), absorbing boundary conditions (Enquist et al. [3]), absorbing layer methods (Kosloff and Kosloff [4], Semblat et al. [5], Rajagopal et al. [6], Zafati et al. [7]), and perfectly matched layers (PMLs).

The PML proposed by Bérenger [8] for absorbing electromagnetic waves, and shortly after interpreted by Chew [9] in terms of complex coordinate stretching, has been increasingly used for dealing with infinite media in the context of finite difference, finite element, and spectral element methods. The PML media, constructed by applying complex-valued coordinate stretching to the elastic wave equation, provides the same attenuation for all frequencies and non-reflecting features in the continuous setting for all angles of incidence at the interface. This makes it more efficient than other absorbing layers. The first PML adapted to elastodynamic equations was formulated according to a velocity-stress form in the framework of the finite difference method by using a split procedure for the components of velocities with respect to the interface (Chew and Liu [10], Collino and Tsogka [11]). The physical variable was partitioned into two components that were orthogonal and parallel to the boundary, leading to an increase in the number of unknown parameters. Then, an unsplit formulation, known as C-PML, requiring the computation of convolution integrals, was developed by Wang et al. [12]. Matzen [13] extended the C-PML approach to the finite element method. For time-domain elastodynamics, Basu and Chopra [14, 15] also proposed an unsplit formulation in two-dimensional (2D) problems, which was displacement-based for straightforward finite element implementation. Here, convolution integrals were not computed, but additional quantities such as integrals of stress and strain in time were required. Later, Basu [16] extended this method to the three-dimensional (3D) case in explicit dynamics, which was later implemented in finite element software such as LS-DYNA and DIANA [17, 18]. Recently, Brun et al. [19] have implemented Basu and Chopra's formulation for 2D problems using a dual subdomain coupling approach. The issue of computing integrals of stress and strain in time was circumvented by Kucukcoban and Kallivokas [20] by introducing a combined stress-displacement formulation at the expense of increase in system size. This approach was finally extended to 3D problems by Fathi et al. [21]. The present paper presents an efficient hybrid (different time integrators) asynchronous (different time steps) 3D PML for modeling unbounded domains through a standard displacement-based finite element method, which is well suitable for finite element implementation. The proposed unsplit 3D explicit/implicit PML formulation is implemented in the framework of a heterogeneous asynchronous time integrator (HATI) [22, 23], which employs the dual approach with Lagrange multipliers for subdomain coupling. This enables the PML to be treated independently using an explicit or implicit scheme with large time steps while conserving classical finite element formulations in the elastic domain to optimize computational efficiency.

In this paper, we first consider the strong form of 1D wave propagation in a PML medium and present the classical design equation of the PML. This enables us to choose appropriate parameters for the PML. The frequency-independent absorbing capabilities of the PML turn out to be very similar to the case of absorbing layers based on Kosloff damping [4], but with the advantage of being reflectionless at the interface between the non-dissipative interior domain and the PML. Next, the strong and weak forms of 3D PML proposed by Basu [16] are rewritten. For this purpose, first, the complex stretching function is changed so as to avoid introducing the characteristic length. Second, a convenient expression of the internal force is developed for integration using implicit or explicit time integration schemes. In Section 4, the weak formulation of the coupled problem, including the interior and PML subdomains, is presented according to a dual coupling approach requiring the introduction of Lagrange multipliers. The versatility of the HATI framework enables the interior domain to be handled by the classical finite element
formulation. Moreover, the PML is dealt with by complex-valued stretched coordinates and an appropriate time integration scheme with its own time step to be chosen independent of the time-stepping procedure adopted in the interior subdomain. Finally, various examples, including a 3D semi-infinite bar, Lamb's test, and soil-structure interaction (SSI) with different soil layers, are investigated to illustrate the efficiency of the HATI formulation in terms of accuracy and CPU time.

## 2. Design of a perfectly matched layer

The PML model is formulated by introducing the complex-valued stretching functions into the classical elastodynamic equations in a frequency domain. The main idea is to replace the real coordinate $x_i$, which denotes the x, y, and z coordinates for the index equal to 1, 2, and 3, by the complex coordinates $\tilde{x}_i : \mathbb{R} \to \mathbb{C}$.

The complex coordinates are defined by

$$ \frac{\partial \tilde{x}_i}{\partial x_i} = \lambda_i(x_i) = 1 + f_i^e(x_i) + \frac{f_i^p(x_i)}{i\omega}. \quad (1) $$

In the above equation, $\omega$ denotes the circular frequency and $f_i^p$ is the attenuation function, which is positive real-valued as a function of $x_i$. The attenuation function serves to attenuate the propagating waves in the $x_i$ direction, whereas the scaling function $f_i^e$ attenuates the evanescent waves by stretching the coordinate variable $x_i$ [24]. It has to be noted that this expression is different from that adopted by Basu and Chopra [14, 15] and Basu [16] so as to avoid introducing a characteristic length of the problem under consideration.

In the following, we focus on the propagating waves in the 1D case to design the PML attenuation performance. To this end, we study the effect of the damping function $f_i^p$ on the attenuation in the PML as well as on the wave reflection at the interface between a non-dissipative elastic medium and the dissipative PML. Indeed, the design of the PML aims at damping out all the incident waves from the domain of interest while minimizing the spurious waves reflected at the boundary of the truncated domain. For this purpose, the strong form of wave propagation in 1D PML media is investigated. It will be shown that the attenuation formulation proposed by Kosloff and Kosloff [4] shares the same absorbing and frequency-independent capabilities of the PML as underlined by Carcione and Kosloff [25]. In addition, the non-reflecting characteristic at the interface between a non-dissipative elastic medium and a PML medium is shown by considering the continuous problem of wave propagation in contrast to the Kosloff medium, which in theory is not reflectionless at the interface [26].

### 2.1. 1D Wave propagation in PML medium

The governing equations in elastodynamics are modified by complex coordinate stretching. The displacement in one-dimensional (1D) PML medium is governed by the modified equations

$$ \frac{1}{\lambda(x)} \frac{\partial \sigma}{\partial x} = -\rho \omega^2 u, \quad (2) $$

$$ \varepsilon = \frac{1}{\lambda(x)} \frac{\partial u}{\partial x}, \quad (3) $$

$$ \sigma = E\varepsilon. \quad (4) $$

Equations (2)–(4) constitute the strong form of the 1D PML medium in the frequency domain. Parameters $\sigma$ and $\varepsilon$ are the complex values denoting the stress and the strain, respectively, $\rho$ is the density, and $E$ is Young's modulus for P waves. For S waves, by replacing Young's modulus $E$ by the
shear modulus G, the same equations hold. It can be seen that the equation of motion and the deformation equation have been modified by the introduction of a complex-valued stretching function $\lambda(x)$, and the elastic constitutive relationship remains intact.

We adopt the scaling and attenuation functions $f^e = 0$, $f^p = \gamma$, respectively, where $\gamma$ is a constant positive real value. It gives the expression of the constant complex-valued stretching function

$$ \lambda = 1 - \frac{i\gamma}{\omega}. \qquad (5) $$

The deformation equation and the elastic constitutive relationship are used to replace the stress term in the wave propagation equation as follows:

$$ E \frac{1}{\lambda^2} \frac{\partial^2 u}{\partial x^2} = -\rho \omega^2 u. \qquad (6) $$

Introducing $\lambda$ into the equation of motion,

$$ E \frac{\partial^2 u}{\partial x^2} = \rho(\omega i)^2 u + 2\rho(\omega i) u\gamma + \rho u\gamma^2. \qquad (7) $$

The dependence of the complex coordinates on the factor $i\omega$ allows an easy application of the inverse Fourier transform when expressing the PML in the time domain. We apply the inverse Fourier transform to obtain the wave propagation in the time domain of the PML medium as follows:

$$ E \frac{\partial^2 u}{\partial x^2} = \rho \ddot{u} + 2\rho \dot{u} \gamma + \rho u \gamma^2. \qquad (8) $$

Here, it can be remarked that the above equation of motion is the same as the Kosloff damping formulation originally proposed by Kosloff and Kosloff [4] before the seminal paper of Bérenger [6] on PML for electromagnetic waves. As a result, it will be shown in the following that the same attenuation capabilities will be derived for the Kosloff medium as those related to the PML medium.

By introducing the harmonic solution $u(x, t) = u_0 \exp(i(\omega_0 t - kx))$, the expression of the wave number $k$ can be obtained:

$$ k = \frac{\omega_0}{\nu} \left( 1 - i \frac{\gamma}{\omega_0} \right), \qquad (9) $$

where $\nu$ denotes the velocity of P waves or S waves. The expression of the rightward propagating wave in the 1D PML medium is

$$ u(x, t) = u_0 \exp(i(\omega_0 t - k_0 x)) \exp\left(-\frac{\gamma x}{\nu}\right), \quad k_0 = \frac{\omega_0}{\nu}. \qquad (10) $$

Using the previous expression in (10), the absorbing capability of PML is given in the form of logarithmic decrement as a function of thickness and constant attenuation scalar of the PML:

$$ \delta = \ln \left( \frac{|u(x)|}{|u(x + \Delta x)|} \right) = \frac{\gamma \Delta x}{\nu}. \qquad (11) $$

It can be seen that the wave frequency $\omega_0$ has no effect on the absorbing ability of the PML with regard to the logarithmic decrement, which means that all waves with all frequencies can be attenuated in the same way. Thus, PML turns out to be independent of frequency. The velocity of P waves is higher than the velocity of S waves for the same medium. In other words, based on the above relationships, to achieve the same logarithmic decrement, the layer thickness for damping out S waves should be smaller than that related to P waves. Therefore, the velocity of P waves $v_p$ is adopted for the design of the absorbing layer.

It should be remarked again that the PML medium has the same wave propagation form in the time domain as in the Kosloff medium because both media have the same equation of motion as given in (8). However, the behavior of the PML medium at the interface with the non-dissipative medium turns out to be better than that of the Kosloff medium. In comparison with the Kosloff
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Wave propagation from elastic medium to PML medium.</figcaption></figure>

medium, the change in the deformation equation renders the PML reflectionless at the interface between the elastic medium and the PML medium as demonstrated in the following.

## 2.2. Wave propagation from elastic medium to PML medium

The wave propagation problem from an elastic medium to a PML medium is considered below in the case of 1D harmonic waves. As shown in Figure 1, three components have to be taken into account: the incident wave $u_1$, the transmitted wave $u_2$, and the reflected wave $u_R$.

$$u_1(x, t) = A \exp \left[ i\omega_0 \left( t - \frac{x}{v_1} \right) \right], \quad (12)$$

$$u_2(x, t) = T \exp \left[ i\omega_0 \left( t - \frac{x}{v_2} \right) \right] \exp \left[ -\frac{\gamma x}{v_2} \right], \quad (13)$$

$$u_R(x, t) = R \exp \left[ i\omega_0 \left( t + \frac{x}{v_1} \right) \right]. \quad (14)$$

Based on the continuity of displacements and equilibrium of stresses at the interface, we can write

$$u_2(x = 0, t) = u_1 + u_R(x = 0, t), \quad (15)$$

$$\frac{E_2}{\lambda} \partial_x u_2(x = 0, t) = E_1 (\partial_x u_1 + \partial_x u_R)(x = 0, t). \quad (16)$$

From the continuity of displacements at interface (15), we have

$$T = A + R. \quad (17)$$

By substituting (12)–(14) in (16) and using the definition of the complex stretching function given in (6), the continuity of stresses can be expressed as

$$\frac{E_2}{1 + \frac{\gamma}{i\omega_0}} \left[ \left( -\frac{i\omega_0}{v_2} \right) T + \left( -\frac{\gamma}{v_2} \right) T \right] = E_1 \left( \frac{i\omega_0}{v_1} R - \frac{i\omega_0}{v_1} A \right). \quad (18)$$

Assuming the same Young's modulus ($E_1 = E_2$) and wave velocities ($v_1 = v_2$) in both media, the continuity of stresses at the interface can be simplified as

$$T = A - R. \quad (19)$$

Finally, taking into account (17), we derive the remarkable property of the PML for all frequencies:

$$T = A, \quad R = 0. \quad (20)$$

This means that the incident wave is equal to the transmitted wave; no reflected wave will be produced at the interface. In other words, in theory, the PML is completely reflectionless, and this is true for all non-zero constant attenuation parameters $\gamma$. As reported before, the difference between the PML and Kosloff media lies in the interface behavior, which is reflectionless in the case of the PML in contrast to the Kosloff medium owing to the introduction of the complex-valued stretching function into the deformation equation. Nonetheless, this property
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Evolution of the damping functions in PML subdomain.</figcaption></figure>

is only valid in the theoretical derivation. Indeed, although there is no reflection at the interface analytically, spatial discretization introduces spurious reflections at the interface. Therefore, optimal PML parameters need to be applied for minimizing these numerical reflections. The real-valued positive functions should be monotonically increasing and vanish at the interface so that the contrast is minimized in the discrete setting between the physical domain and the unphysical PML. Classically, the damping function $f^p$ is written as a polynomial of degree $m_{PML}$ as follows [14–16, 20, 21, 27]:

$$ \gamma(x) = \gamma_0 \left(\frac{x-x_0}{L}\right)^{m_{PML}}, \qquad (21) $$

where $\gamma_0$ is a user-tunable scalar parameter. In fact, the larger the value of $\gamma_0$, the larger the discretization error. Namely, more reflected spurious waves are produced at the interface with larger $\gamma_0$ values as is shown in the following numerical investigations.

The logarithmic decrement of the PML domain $\delta$ is obtained by integrating (11) along the thickness of the PML:

$$ \delta = \int_{x_0}^{x_0+L} \frac{\gamma_0}{v_{\rho}} \left(\frac{x-x_0}{L}\right)^{m_{PML}} dx = \frac{\gamma_0 L}{(m_{PML}+1) v_{\rho}}. \qquad (22) $$

We define the attenuation coefficient $R_{attenuation}$ from the logarithmic decrement:

$$ R_{attenuation} = \left( \frac{|u(x+L)|}{|u(x)|} \right)^2 = e^{-2\delta}. \qquad (23) $$

For instance, if we want to reach a logarithmic decrement target of $\delta = \ln(10)$, this means that 90% of the amplitude of the incident wave is absorbed from the interface to the end of the PML. Next, the attenuation also occurs for the reflection process from the end of the PML toward the interface. Hence, the incident wave is attenuated by 99%, and the attenuation coefficient $R_{attenuation}$ is theoretically equal to 1% before the space and time discretization.
Finally, we can propose the general formula to design the PML based on the 1D harmonic wave problem presented in the PML medium. After choosing the $R_{attenuation}$ value, the total thickness $L$ and the power $m_{PML}$ of the damping function $\gamma_0$ can be obtained:

$$ \gamma_0 = \frac{(m_{PML} + 1)}{2L} \times v_{\rho} \times \ln\left(\frac{1}{R_{attenuation}}\right). \qquad (24) $$

### 3. Three-dimensional PML

In this section, the discrete formulation of the PML for 3D elastodynamics is presented, leading to an efficient method for calculating the internal force in the PML domain. The main steps of the PML development are resumed. The details about matrices related to the derivatives of the shape functions of the hexahedral element, combined with attenuation and scaling functions of the PML, can be found in Appendix A.

#### 3.1. Strong form of the three-dimensional PML

As shown in Section 2.1, the PML formulation is obtained by modifying the governing equations defined in the frequency domain. The frequency-domain equations for PML are obtained by applying the complex-valued stretching functions related to three directions:

$$ \sum_j \frac{1}{\lambda_j(x_j)} \frac{\partial \sigma_{ij}}{\partial x_j} = -\omega^2 \rho u_i, \qquad (25) $$

$$ \sigma_{ij} = \sum_{k,l} C_{ijkl} \varepsilon_{kl}, \qquad (26) $$

$$ \varepsilon_{ij} = \frac{1}{2} \left[ \frac{1}{\lambda_j(x_j)} \frac{\partial u_i}{\partial x_j} + \frac{1}{\lambda_i(x_i)} \frac{\partial u_j}{\partial x_i} \right], \qquad (27) $$

where $C_{ijkl}$ are the components of the elastic constitutive tensor.

Then, we introduce the following notation for the PML region. $\Omega_{PML}$ is the region of the PML bounded by $\Gamma_{PML} = \Gamma_{PML}^D + \Gamma_{PML}^N$, where $\Gamma_{PML}^D \cap \Gamma_{PML}^N = \emptyset$, defining the decomposition of the boundary conditions into Dirichlet and Neumann conditions. In addition, $g_N$ denotes the prescribed traction force on $\Gamma_{PML}^N$ and $J = [0, T]$ is the time interval of interest.

Thanks to the introduction of the stretching functions expressed in (1), the inverse Fourier transform can be easily applied to the previous frequency-domain equations, leading to the following equations in the time domain:

$$ \underline{\underline{\text{div}}}(\underline{\underline{\sigma}} \tilde{F}^{ee} + \underline{\underline{\Sigma}} \tilde{F}^{ep} + \underline{\underline{\Sigma}} \tilde{F}^{pp}) = \rho f_M \underline{\underline{u}} + \rho f_C \underline{\underline{u}} + \rho f_K \underline{\underline{u}} + \rho f_H \underline{\underline{u}}, \qquad (28) $$

$$ \underline{\underline{\sigma}} = \underline{\underline{C}} : \underline{\underline{\varepsilon}}, \qquad (29) $$

$$ \begin{aligned} \underline{\underline{F}}^{eT} \underline{\underline{\varepsilon}} \underline{\underline{F}}^{e} + (\underline{\underline{F}}^{pT} \underline{\underline{\varepsilon}} \underline{\underline{F}}^{e} + \underline{\underline{F}}^{eT} \underline{\underline{\varepsilon}} \underline{\underline{F}}^{p}) + \underline{\underline{F}}^{pT} \underline{\underline{E}} \underline{\underline{F}}^{p} \\ = \frac{1}{2} [\underline{\underline{F}}^{eT}(\underline{\underline{\nabla u}}) + (\underline{\underline{\nabla u}})^T \underline{\underline{F}}^{e}] + \frac{1}{2} [\underline{\underline{F}}^{pT}(\underline{\underline{\nabla u}}) + (\underline{\underline{\nabla u}})^T \underline{\underline{F}}^{p}]. \end{aligned} \qquad (30) $$

The first equation of the above system is the equation of motion in the PML, which is complemented by clamped Dirichlet conditions and zero traction forces at the Neumann conditions:

$$ \begin{cases} \underline{\underline{u}} = 0 & \text{on } \Gamma_{PML}^D \\ (\underline{\underline{\sigma}} \tilde{F}^{ee} + \underline{\underline{\Sigma}} \tilde{F}^{ep} + \underline{\underline{\Sigma}} \tilde{F}^{pp}) \underline{\underline{n}} = 0 & \text{on } \Gamma_{PML}^N. \end{cases} \qquad (31) $$

In addition to the stress tensor, the time-domain PML involves the time integral of the stress tensor and the time integral of the time-integral stress tensor, which are defined as follows:

$$ \underline{\underline{\Sigma}} = \int_{0}^{t} \underline{\underline{\sigma}} dt, \quad \underline{\underline{\Sigma}} = \int_{0}^{t} \underline{\underline{\Sigma}} dt. \qquad (32) $$
It is also noted that the equation of motion is now a third-order differential equation with four fields: the classical displacement; velocity and acceleration fields, which are complemented by the time integral of the displacement expressed as $\underline{U} = \int_{0}^{t} \underline{u} d\tau$. The same equation of motion has been obtained by Basu [16], where one multiplicative factor arises from a slightly different choice of the stretching function as previously discussed. The second equation represents the classical constitutive relationship for an elastic medium. The third equation is the PML strain-deformation relationship identical to Basu's formulation, which involves the strain rate and the time integral of the strain tensor given by

$$ \underline{E} = \int_{0}^{t} \underline{\epsilon} dt. \qquad (33) $$

All the matrices involved in (28) and (30) depend on scaling functions $f_i^e(x_i)$ and attenuation functions $f_i^p(x_i)$. Their expressions are presented in Appendix A.

## 3.2. Weak form of the three-dimensional PML

The space discretization is displacement-based, following a standard finite element formulation. The space and time discretization is summarized in the following before presenting the time coupling of a hybrid multi-time-step PML with the physical domain. Let $\underline{v}$ be the test function belonging to an appropriate space. The weak formulation is obtained by integrating over the PML domain:

$$ \int_{\Omega} \rho f_M \underline{v} \cdot \dot{\underline{u}} d\Omega + \int_{\Omega} \rho f_C \underline{v} \cdot \dot{\underline{u}} d\Omega + \int_{\Omega} \rho f_K \underline{v} \cdot \underline{u} d\Omega + \int_{\Omega} \rho f_H \underline{v} \cdot \underline{U} d\Omega + \int_{\Omega} \underline{\tilde{\epsilon}}^{ee}: \underline{\sigma} d\Omega \\ + \int_{\Omega} \underline{\tilde{\epsilon}}^{ep}: \underline{\Sigma} d\Omega + \int_{\Omega} \underline{\tilde{\epsilon}}^{pp}: \underline{\tilde{\Sigma}} d\Omega = \int_{\Gamma_N} \underline{v} \cdot (\underline{\sigma} \underline{\tilde{F}}^{ee} + \underline{\Sigma} \underline{\tilde{F}}^{ep} + \underline{\tilde{\Sigma}} \underline{\tilde{F}}^{pp}) \cdot \underline{n} d\Gamma. \qquad (34) $$

Here, taking into account the scaling and damping functions, the expressions of the modified strain tensors are as follows:

$$ \left\{ \begin{aligned} \underline{\tilde{\epsilon}}^{ee} &= \frac{1}{2} \left[ (\nabla \underline{v}) \underline{\tilde{F}}^{ee} + \underline{\tilde{F}}^{eeT} (\nabla \underline{v})^T \right] \\ \underline{\tilde{\epsilon}}^{ep} &= \frac{1}{2} \left[ (\nabla \underline{v}) \underline{\tilde{F}}^{ep} + \underline{\tilde{F}}^{epT} (\nabla \underline{v})^T \right] \\ \underline{\tilde{\epsilon}}^{pp} &= \frac{1}{2} \left[ (\nabla \underline{v}) \underline{\tilde{F}}^{pp} + \underline{\tilde{F}}^{ppT} (\nabla \underline{v})^T \right]. \end{aligned} \right. \qquad (35) $$

The internal force is expressed as follows:

$$ p_{\text{int}}^{e} = \int_{\Omega} \underline{\tilde{\epsilon}}^{ee}: \underline{\sigma} \, d\Omega + \int_{\Omega} \underline{\tilde{\epsilon}}^{ep}: \underline{\Sigma} \, d\Omega + \int_{\Omega} \underline{\tilde{\epsilon}}^{pp}: \underline{\tilde{\Sigma}} \, d\Omega. \qquad (36) $$

## 3.3. Finite element discretization

In the following, we consider the space discretization for a classical eight-node hexahedral element with linear shape functions. The approximation of the displacement is given as follows: $u_e(x, y, z) = N(x, y, z)U_e$ of size 24 × 1, where $U_e$ gathers the nodal displacements of the eight nodes; the matrix $N(x, y, z)$ of size 3 × 24 contains the nodal shape functions $N(x, y, z) = [N_1 \mathbf{I} \ N_2 \mathbf{I} \ \dots \ N_8 \mathbf{I}]$, where $\mathbf{I}$ is the 3 × 3 identity matrix.

Introducing the finite element discretization into the weak form of the equation of motion in (34), the semi-discrete equation of motion can be derived:

$$ \mathbf{M}\ddot{\mathbf{U}} + \mathbf{C}\dot{\mathbf{U}} + \mathbf{K}\mathbf{U} + \bar{\mathbf{K}}\bar{\mathbf{U}} + \mathbf{P}_{\text{int}} = \mathbf{P}_{\text{ext}}. \qquad (37) $$

The inertial system matrices $\mathbf{M}$, $\mathbf{C}$, $\mathbf{K}$, and $\bar{\mathbf{K}}$ are assembled from their respective element-level matrices. The element-level matrices are obtained by adopting a quadrature formula in every
hexahedral element:

$$ \mathbf{M}^{\mathrm{e}} = \int_{\Omega^{\mathrm{e}}} \rho f_{\mathrm{M}} \mathbf{N}^{\mathrm{T}} \mathbf{N} \mathrm{d}\Omega, \qquad (38) $$

$$ \mathbf{C}^{\mathrm{e}} = \int_{\Omega^{\mathrm{e}}} \rho f_{\mathrm{C}} \mathbf{N}^{\mathrm{T}} \mathbf{N} \mathrm{d}\Omega, \qquad (39) $$

$$ \mathbf{K}^{\mathrm{e}} = \int_{\Omega^{\mathrm{e}}} \rho f_{\mathrm{K}} \mathbf{N}^{\mathrm{T}} \mathbf{N} \mathrm{d}\Omega, \qquad (40) $$

$$ \bar{\mathbf{K}}^{\mathrm{e}} = \int_{\Omega^{\mathrm{e}}} \rho f_{\mathrm{H}} \mathbf{N}^{\mathrm{T}} \mathbf{N} \mathrm{d}\Omega. \qquad (41) $$

Taking into account (35), the internal force term $\mathbf{P}_{\mathrm{int}}^{\mathrm{e}}$ can be written as

$$ \mathbf{P}_{\mathrm{int}}^{\mathrm{e}} = \int_{\Omega^{\mathrm{e}}} \bar{\mathbf{B}}^{\mathrm{ee}T} \hat{\boldsymbol{\sigma}} \mathrm{d}\Omega + \int_{\Omega^{\mathrm{e}}} \bar{\mathbf{B}}^{\mathrm{ep}T} \hat{\boldsymbol{\Sigma}} \mathrm{d}\Omega + \int_{\Omega^{\mathrm{e}}} \bar{\mathbf{B}}^{\mathrm{pp}T} \hat{\boldsymbol{\Sigma}} \mathrm{d}\Omega, \qquad (42) $$

where the matrices $\bar{\mathbf{B}}^{\mathrm{ee}}$, $\bar{\mathbf{B}}^{\mathrm{ep}}$, and $\bar{\mathbf{B}}^{\mathrm{pp}}$ depend on the derivatives of the shape functions and the scaling and attenuation functions of the PML. Their expressions are presented in Appendix A. In (42), we use the Voigt notation, where $\hat{\boldsymbol{\sigma}}$ represents the six-component vector of stresses and $\hat{\boldsymbol{\Sigma}}$ and $\hat{\boldsymbol{\Sigma}}$ are the successive time integrals of stresses.

For the time-stepping procedure over the time step $[t_n; t_{n+1}]$, additional relationships are assumed:

$$ \hat{\mathbf{E}}_{n+1} = \hat{\mathbf{E}}_n + \hat{\boldsymbol{\varepsilon}}_n \Delta t, \quad \hat{\boldsymbol{\Sigma}}_{n+1} = \hat{\boldsymbol{\Sigma}}_n + \hat{\boldsymbol{\sigma}}_n \Delta t, \quad \hat{\boldsymbol{\Sigma}}_{n+1} = \hat{\boldsymbol{\Sigma}}_n + \hat{\boldsymbol{\Sigma}} \Delta t, \qquad (43) $$

$$ \hat{\boldsymbol{\varepsilon}}_{n+1} = \frac{\hat{\boldsymbol{\varepsilon}}_{n+1} - \hat{\boldsymbol{\varepsilon}}_n}{\Delta t}. \qquad (44) $$

Using the assumptions given in (44), the third equation of the system in (30) leads to the expression of the strain $\hat{\boldsymbol{\varepsilon}}_{n+1}$ at the end of the time step:

$$ \hat{\boldsymbol{\varepsilon}}_{n+1} = \frac{1}{\Delta t} \left[ \frac{1}{\Delta t} \hat{\mathbf{F}}^{\mathrm{e}} \hat{\boldsymbol{\varepsilon}}_n - \hat{\mathbf{F}}^{\mathrm{Q}} \hat{\mathbf{E}}_n + \mathbf{B}^{\mathrm{e}} \dot{\mathbf{U}}_{n+1} + \mathbf{B}^{\mathrm{Q}} \mathbf{U}_{n+1} \right]. \qquad (45) $$

The matrices $\hat{\mathbf{F}}^{\mathrm{e}}$, $\hat{\mathbf{F}}^{\mathrm{Q}}$, $\mathbf{B}^{\mathrm{e}}$, and $\mathbf{B}^{\mathrm{Q}}$ depending on the derivatives of shape functions as well as scaling and attenuation functions are defined in Appendix A. It has to be noted that the above strain–deformation relationship in the 3D case is the same as that in the 2D case [15].

Here, we propose a convenient expression of the internal force, which is different from the computation of the strain terms in the paper by Basu [16] for the 3D PML. The internal force is decomposed into two parts. The first part of the internal force contains only known quantities at time $t_n$, whereas the second part contains the unknown quantities at time $t_{n+1}$. Thus, the element-wise internal force vector $\mathbf{P}_{n+1}^{\mathrm{e}}$ can be written in terms of the element velocity and displacement vectors ($\dot{\mathbf{U}}_{n+1}^{\mathrm{e}}$ and $\mathbf{U}_{n+1}^{\mathrm{e}}$) as well as a term denoted as $\mathbf{P}(\boldsymbol{\varepsilon}_n^{\mathrm{e}}, \mathbf{E}_n^{\mathrm{e}}, \boldsymbol{\Sigma}_n^{\mathrm{e}}, \tilde{\boldsymbol{\Sigma}}_n^{\mathrm{e}})$ depending only on known quantities at the beginning of the time step. The element-wise internal force $\mathbf{P}_{n+1}^{\mathrm{e}}$ is written as

$$ \mathbf{P}_{n+1}^{\mathrm{e}} = \left( \int_{\Omega_e} \bar{\mathbf{B}}^{\mathrm{T}} \frac{1}{\Delta t} \mathrm{D}\mathbf{B}^{\mathrm{e}} \mathrm{d}\Omega \right) \dot{\mathbf{U}}_{n+1} + \left( \int_{\Omega_e} \bar{\mathbf{B}}^{\mathrm{T}} \frac{1}{\Delta t} \mathrm{D}\mathbf{B}^{\mathrm{Q}} \mathrm{d}\Omega \right) \mathbf{U}_{n+1} + \mathbf{P}(\boldsymbol{\varepsilon}_n^{\mathrm{e}}, \mathbf{E}_n^{\mathrm{e}}, \boldsymbol{\Sigma}_n^{\mathrm{e}}, \tilde{\boldsymbol{\Sigma}}_n^{\mathrm{e}}), \qquad (46) $$

where the matrix $\bar{\mathbf{B}}$ is defined as a function of the previous matrices $\bar{\mathbf{B}}^{\mathrm{ee}}$, $\bar{\mathbf{B}}^{\mathrm{ep}}$, and $\bar{\mathbf{B}}^{\mathrm{pp}}$ and the time step as follows:

$$ \bar{\mathbf{B}} = \bar{\mathbf{B}}^{\mathrm{ee}} + \Delta t \bar{\mathbf{B}}^{\mathrm{ep}} + \Delta t^2 \bar{\mathbf{B}}^{\mathrm{pp}}. \qquad (47) $$

The known part of the internal force at the beginning of the time step is given by

$$ \mathbf{P}(\boldsymbol{\varepsilon}_n^{\mathrm{e}}, \mathbf{E}_n^{\mathrm{e}}, \boldsymbol{\Sigma}_n^{\mathrm{e}}, \tilde{\boldsymbol{\Sigma}}_n^{\mathrm{e}}) = \left( \int_{\Omega_e} \bar{\mathbf{B}}^{\mathrm{T}} \frac{1}{\Delta t^2} \mathrm{D}\boldsymbol{\varepsilon}_n^{\mathrm{e}} \mathrm{d}\Omega \right) - \left( \int_{\Omega_e} \bar{\mathbf{B}}^{\mathrm{T}} \frac{1}{\Delta t} \mathrm{D}\mathbf{F}^{\mathrm{Q}} \mathbf{E}_n^{\mathrm{e}} \mathrm{d}\Omega \right) \\ + \left( \int_{\Omega_e} \bar{\mathbf{B}}^{\mathrm{ep}T} \boldsymbol{\Sigma}_n^{\mathrm{e}} \mathrm{d}\Omega \right) + \left( \int_{\Omega_e} \bar{\mathbf{B}}^{\mathrm{pp}T} \tilde{\boldsymbol{\Sigma}}_n^{\mathrm{e}} \mathrm{d}\Omega \right) + \left( \int_{\Omega_e} \bar{\mathbf{B}}^{\mathrm{pp}T} \Delta t \boldsymbol{\Sigma}_n^{\mathrm{e}} \mathrm{d}\Omega \right). \qquad (48) $$
To express the part of the internal force that has to be computed at the end of the time step given by the first two terms in (46), the element-level matrices are defined as

$$ \tilde{\mathbf{C}}^{\mathbf{e}} = \int_{\Omega_{\mathbf{e}}} \tilde{\mathbf{B}}^{\mathbf{T}} \frac{1}{\Delta t} \mathbf{DB}^{\mathbf{e}} \mathrm{d}\Omega, \quad (49) $$

$$ \tilde{\mathbf{K}}^{\mathbf{e}} = \int_{\Omega_{\mathbf{e}}} \tilde{\mathbf{B}}^{\mathbf{T}} \frac{1}{\Delta t} \mathbf{DB}^{\mathbf{Q}} \mathrm{d}\Omega. \quad (50) $$

These two matrices can be viewed as one viscous matrix operating on velocities and one additional stiffness matrix operating on displacements. The proposed method for calculating the internal force can be employed in explicit and implicit time integration. Finally, after assembling the element matrices given in (38)–(41) and (49)–(50), the space and time discrete equation of motion is obtained at the end time $t_{n+1}$:

$$ \mathbf{M}\ddot{\mathbf{U}}_{n+1} + (\mathbf{C} + \tilde{\mathbf{C}})\dot{\mathbf{U}}_{n+1} + (\mathbf{K} + \tilde{\mathbf{K}})\mathbf{U}_{n+1} + \bar{\mathbf{K}}\bar{\mathbf{U}}_{n+1} + \mathbf{P}(\varepsilon_{\mathbf{n}}, \mathbf{E}_{\mathbf{n}}, \mathbf{\Sigma}_{\mathbf{n}}, \tilde{\mathbf{\Sigma}}_{\mathbf{n}}) = \mathbf{F}_{\text{ext}}. \quad (51) $$

The above equation is third-order in time, requiring a specific time integration procedure. Here, a dual approach is preferred for coupling the elastic domain and PML. Indeed, the elastic domain and the PML are integrated in time by using the powerful and flexible framework of HATI to enable choosing in each partition the appropriate time integrator with its own time step while conserving the classical finite element formulation in other subdomains.

# 4. Hybrid asynchronous coupling between the physical domain and PML

## 4.1. Weak form for subdomain coupling

Let $\Omega$ be a bounded domain belonging to $\mathbb{R}^3$ with a regular boundary. The time interval of interest is $J = [0, T]$. The domain $\Omega$ is divided into two parts $\Omega_1$ and $\Omega_2$ as shown in Figure 3: $\Omega_1 \cap \Omega_2 = \emptyset$ and $\partial\Omega_1 \cap \partial\Omega_2 = \Gamma_I$. The parameter $\Gamma_I$ denotes the interface between the two subdomains, subdomain $\Omega_1$ representing the non-dissipative medium (the domain of interest) and subdomain $\Omega_2$ denoting the PML medium. Subdomain 1 is related to a linear elastic behavior and subdomain 2 is related to the PML region that was presented previously. The subdomain $\Omega_1$ is characterized by its density $\rho_1$; Young's modulus $E_1$; Poisson's coefficient $\nu_1$; the body force $\underline{\mathbf{b}}_1$; the Dirichlet prescribed displacement on $\Gamma_1^D$, $\underline{u}_1^D$; and the traction force at the Neumann condition on $\Gamma_1^N$, $g_1^N$. The subdomain $\Omega_2$ is characterized by its density $\rho_2$; Young's modulus $E_2$; Poisson's coefficient $\nu_2$; the body force $\underline{\mathbf{b}}_2$; the Dirichlet prescribed displacement on $\Gamma_2^D$, $\underline{u}_2^D$; and the traction force at the Neumann condition on $\Gamma_2^N$, $g_2^N$.

Before writing the weak form of the coupled problem in $\Omega$ divided into two parts $\Omega_1$ and $\Omega_2$, we classically introduce test functions $\underline{\mathbf{v}}_1$ and $\underline{\mathbf{v}}_2$, which belong to appropriate spaces: $\underline{\mathbf{v}}_1 \in W_1^* = \{\underline{\mathbf{v}}_1 \in (H^1(\Omega_1))^d \text{ and } \underline{\mathbf{v}}_1 = 0 \text{ on } \Gamma_1^D\}$ and $\underline{\mathbf{v}}_2 \in W_2^* = \{\underline{\mathbf{v}}_2 \in (H^1(\Omega_2))^d \text{ and } \underline{\mathbf{v}}_2 = 0 \text{ on } \Gamma_2^D\}$, where $d$ is the space dimension equal to 3. The solutions belong to the following spaces: $\underline{u}_1(t) \in W_1$, $W_1 = \{\underline{u}_1 \in (H^1(\Omega_1))^d \text{ and } \underline{u}_1 = u_1^D \text{ on } \Gamma_1^D\}$, $\underline{u}_2(t) \in W_2$, and $W_2 = \{\underline{u}_2 \in (H^1(\Omega_2))^d \text{ and } \underline{u}_2 = u_2^D \text{ on } \Gamma_2^D\}$. According to the dual Schur approach, the introduction of the Lagrange multipliers allows us to glue the velocities of the two subdomains at the interface $\Gamma_I$ [28, 29]. The Lagrange multipliers belong to the adapted dual trace space related to the interface between the two subdomains denoted by $Q$. All the space variables considered are assumed to be sufficiently smooth and regular.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Domain Ω divided into two subdomains Ω₁ and Ω₂.</figcaption></figure>

Next, using a dual Schur formulation, we can write the principle of virtual power for transient dynamics. We find the solution $\underline{u}_1(t) \in W_1$, $\underline{u}_2(t) \in W_2$, and $\underline{\lambda}(t) \in Q$ for which the following weak form is satisfied: $\forall \underline{v}_1 \in W_1^*, \forall \underline{v}_2 \in W_2^*$ and $\forall \underline{\mu} \in Q$.

$$ 
\begin{aligned}
& \int_{\Omega_1} \rho_1 \underline{v}_1 \cdot \underline{\ddot{u}}_1 d\Omega + \int_{\Omega_1} \underline{\varepsilon}(\underline{v}_1) : \underline{\underline{\sigma}}_1 d\Omega + \int_{\Omega_2} \rho f_M \underline{v}_2 \cdot \underline{\ddot{u}}_2 d\Omega + \int_{\Omega_2} \rho f_C \underline{v}_2 \cdot \underline{\dot{u}}_2 d\Omega \\
& + \int_{\Omega_2} \rho f_K \underline{v}_2 \cdot \underline{\underline{u}}_2 d\Omega + \int_{\Omega_2} \rho f_H \underline{v}_2 \cdot \underline{\underline{U}}_2 d\Omega \\
& + \int_{\Omega_2} \underline{\varepsilon}^{ee} : \underline{\underline{\sigma}}_2 d\Omega + \int_{\Omega_2} \underline{\varepsilon}^{ep} : \underline{\underline{\Sigma}}_2 d\Omega + \int_{\Omega_2} \underline{\varepsilon}^{pp} : \underline{\underline{\Sigma}}_2 d\Omega \\
& + \int_{\Gamma_I} \underline{v}_1 \cdot \underline{\lambda} d\Gamma + \int_{\Gamma_I} \underline{v}_2 \cdot \underline{\lambda} d\Gamma + \int_{\Gamma_I} \underline{\mu} \cdot (\underline{\dot{u}}_1 - \underline{\dot{u}}_2) d\Gamma \\
& = \int_{\Omega_1} \underline{v}_1 \cdot \underline{b}_1 d\Omega + \int_{\Gamma_1^N} \underline{v}_1 \cdot \underline{g}_1^N d\Gamma + \int_{\Gamma_2^N} \underline{v}_2 \cdot (\underline{\sigma} \underline{\tilde{F}}^{ee} + \underline{\underline{\Sigma}} \underline{\tilde{F}}^{ep} + \underline{\underline{\Sigma}} \underline{\tilde{F}}^{pp}) n d\Gamma.
\end{aligned}
\quad (52) $$

Then, we follow along the classical lines of finite element discretization. At the interface between the subdomains, the continuity of velocities is imposed by the condition

$$ \mathbf{L}_1 \dot{\mathbf{U}}_1 + \mathbf{L}_2 \dot{\mathbf{U}}_2 = 0, \quad (53) $$

where $\mathbf{L}_1$ and $\mathbf{L}_2$ are the constraint matrices of Boolean type in the case of matching meshes at the interface $\Gamma_1$ as presented in previous works [19,26]. From the weak form of the global problem in (52), semi-discrete equations can be derived in space corresponding to the two equations of motion related to the two subdomains, which are completed by a kinematic condition. In the following, the hybrid integration of this set of equations is carried out to propose a 3D hybrid asynchronous PML. For time discretization, the GC method proposed by Gravouil and Combescure is employed [28, 30], which belongs to the HATI framework. Applying continuity of velocities at the interface, it was demonstrated that the coupling GC method is stable for all time integrators (implicit and explicit) belonging to the Newmark family [29], where their own time steps depend on subdomains. Different time integrators with their own time steps
can be adopted depending on the considered subdomain, rendering the proposed framework very useful for coupling complex PML formulations, while conserving classical finite element formulations and time integrators in other subdomains.

In the following, the subdomain $\Omega_1$ is integrated independently in time by a second-order-accurate Newmark explicit time integration scheme. Moreover, the subdomain $\Omega_2$ is treated by an extended third-order-accurate Newmark implicit time integration scheme [21] or a second-order explicit time integration method following the central difference scheme.

## 4.2. Multi-time-step implicit PML

As illustrated in Figure 3, an explicit time integrator with a fine time step $\Delta t_1$ imposed by the Courant-Friedrichs-Lewy (CFL) condition [31] is adopted for the subdomain $\Omega_1$. An implicit time integrator with a large time step $\Delta t_2$ is used for the subdomain $\Omega_2$ because the implicit scheme is unconditionally stable with $\Delta t_2 = m\Delta t_1$, where m is the time step ratio between the two subdomains. In this way, hybrid (different schemes) asynchronous (different time steps depending on subdomains) absorbing layers can be obtained. The equilibrium of subdomain 2 is attained at time $t_m$ at the end of the large time step $\Delta t_2 = [t_0; t_m]$, while the equilibrium of subdomain 1 is attained at the end of every fine time step $\Delta t_1 = [t_{j-1}; t_j]$. The gluing of the velocities at the interface is defined at the fine time scale.

Finally, the weak form given in (52) together with the velocity continuity equation in (53) and the expression of the interface terms as a function of Lagrange multipliers can be expressed in the following discrete form in space and time:

$$ M_1 \ddot{U}_1^j + K_1 U_1^j = F_1^{\text{ext},j} - L_1^T \lambda^j, \quad (54) $$

$$ M_2 \ddot{U}_2^m + (C_2 + \tilde{C}_2) \dot{U}_2^m + (K_2 + \tilde{K}_2) U_2^m + \bar{K}_2 \bar{U}_2^m + P_2(\epsilon_0, E_0, \Sigma_0, \bar{\Sigma}_0) = F_2^{\text{ext},m} - L_2^T \hat{\lambda}^m, \quad (55) $$

$$ L_1 \ddot{U}_1^j + L_2 \dot{U}_2^j = 0. \quad (56) $$

The first equation is the discrete in space equation of motion of the subdomain $\Omega_1$ defined at the end of the fine time step $\Delta t_1 = [t_{j-1}; t_j]$. The second equation is the discrete in space equation of motion of the subdomain $\Omega_2$, corresponding to the PML medium, defined at the end of the large time step $\Delta t_2 = [t_0; t_m]$. The third equation is the discrete velocity continuity at the interface. The subdomain $\Omega_1$ is integrated in time by a Newmark explicit scheme ($\beta_1 = 0$ and $\gamma_1 = 1/2$) with a lumped mass matrix $M_1$. We define $\ddot{U}_1^{j-1,\mathbf{P}}$ as the predictor displacement and $\dot{U}_1^{j-1,\mathbf{P}}$ as the predictor acceleration, which are classically introduced in approximate Newmark formulas:

$$ \ddot{U}_1^{j-1,\mathbf{P}} = U_1^{j-1} + \Delta t \dot{U}_1^{j-1} + (\frac{1}{2} - \beta_1) \Delta t^2 \ddot{U}_1^{j-1}, \quad (57) $$

$$ \dot{U}_1^{j-1,\mathbf{P}} = \ddot{U}_1^{j-1} + \Delta t(1 - \gamma_1)\ddot{U}_1^{j-1}. \quad (58) $$

The classical approximate Newmark formulas in terms of displacements and velocities at the end of the time step $t_j$ are expressed in acceleration form as follows:

$$ \ddot{U}_1^j = U_1^{j-1,\mathbf{P}} + \beta_1 \Delta t^2 \ddot{U}_1^j, \quad (59) $$

$$ \dot{U}_1^j = \ddot{U}_1^{j-1,\mathbf{P}} + \gamma_1 \Delta t \dot{U}_1^j. \quad (60) $$

Regarding the subdomain $\Omega_2$, we use an implicit third-order extended Newmark scheme as proposed by Fathi et al. [21]. For the implicit third-order extended Newmark scheme, $\beta_2$ and $\gamma_2$ are the usual Newmark parameters related to the classical constant average acceleration scheme.
These are equal to 1/4 and 1/2, respectively; $\alpha_2$ is an additional parameter required for the third-order extended Newmark scheme, which is equal to 1/12. Velocities, displacements, and time integrals of displacement are expressed as functions of predictors as follows:

$$ \bar{\mathbf{U}}_{2}^{\mathrm{m}}=\bar{\mathbf{U}}_{2}^{0, \mathrm{p}}+\alpha_{2} \Delta t^{3} \dot{\mathbf{U}}_{2}^{\mathrm{m}}, \qquad (61) $$

$$ \mathbf{U}_{2}^{\mathrm{m}}=\mathbf{U}_{2}^{0, \mathrm{p}}+\beta_{2} \Delta t^{2} \ddot{\mathbf{U}}_{2}^{\mathrm{m}}, \qquad (62) $$

$$ \dot{\mathbf{U}}_{2}^{\mathrm{m}}=\dot{\mathbf{U}}_{2}^{0, \mathrm{p}}+\gamma_{2} \Delta t \ddot{\mathbf{U}}_{2}^{\mathrm{m}}. \qquad (63) $$

Here, the predictors $\dot{\mathbf{U}}_{2}^{0, \mathrm{p}}$, $\mathbf{U}_{2}^{0, \mathrm{p}}$, and $\bar{\mathbf{U}}_{2}^{0, \mathrm{p}}$ are defined as follows:

$$ \bar{\mathbf{U}}_{2}^{0, \mathrm{p}}=\bar{\mathbf{U}}_{2}^{0}+\Delta t \mathbf{U}_{2}^{0}+\frac{\Delta t^{2}}{2} \dot{\mathbf{U}}_{2}^{0}+\left(\frac{1}{6}-\alpha_{2}\right) \Delta t^{3} \ddot{\mathbf{U}}_{2}^{0}, \qquad (64) $$

$$ \mathbf{U}_{2}^{0, \mathrm{p}}=\mathbf{U}_{2}^{0}+\Delta t \dot{\mathbf{U}}_{2}^{0}+\left(\frac{1}{2}-\beta_{2}\right) \Delta t^{2} \ddot{\mathbf{U}}_{2}^{0}, \qquad (65) $$

$$ \dot{\mathbf{U}}_{2}^{0, \mathrm{p}}=\dot{\mathbf{U}}_{2}^{0}+(1-\gamma_{2}) \Delta t \ddot{\mathbf{U}}_{2}^{0}. \qquad (66) $$

Introducing the above approximate Newmark formulas into (54)–(55) leads to the equations of motion whose unknown parameters to be solved are the acceleration terms written as follows:

$$ \left\{ \begin{array}{l} \mathbf{M}_{1} \ddot{\mathbf{U}}_{1}^{j}=\mathbf{F}_{1}^{\mathrm{ext}, j}-\mathbf{K}_{1} \mathbf{U}_{1}^{j-1, p}-\mathbf{L}_{1}^{\mathrm{T}} \lambda^{j}, \\ \tilde{\mathbf{M}}_{2} \ddot{\mathbf{U}}_{2}^{m}=-\mathbf{P}_{2}(\varepsilon_{0}, \mathbf{E}_{0}, \Sigma_{0}, \tilde{\Sigma}_{0})-(\mathbf{C}_{2}+\tilde{\mathbf{C}}_{2}) \dot{\mathbf{U}}_{2}^{0, p}-\left(\mathbf{K}_{2}+\tilde{\mathbf{K}}_{2}\right) \mathbf{U}_{2}^{0, p}-\bar{\mathbf{K}}_{2} \bar{\mathbf{U}}_{2}^{0, p}-\mathbf{L}_{2}^{\mathrm{T}} \lambda^{m}, \\ \mathbf{L}_{1} \dot{\mathbf{U}}_{1}^{j}+\mathbf{L}_{2} \dot{\mathbf{U}}_{2}^{j}=0, \end{array} \right. \qquad (67) $$

where the effective stiffness matrix $\tilde{\mathbf{M}}_2$ in the PML subdomain is defined as

$$ \tilde{\mathbf{M}}_{2}=\mathbf{M}_{2}+\gamma_{2} \Delta t(\mathbf{C}_{2}+\tilde{\mathbf{C}}_{2})+\beta_{2} \Delta t^{2}(\mathbf{K}_{2}+\tilde{\mathbf{K}}_{2})+\alpha_{2} \Delta t^{3} \bar{\mathbf{K}}_{2}. \qquad (68) $$

### 4.3. Multi-time-step explicit PML

The explicit version of the 3D PML is also proposed using the expression of the internal force in (46). The discrete equations of motion become

$$ \left\{ \begin{array}{l} \mathbf{M}_{1} \ddot{\mathbf{U}}_{1}^{j}=\mathbf{F}_{1}^{\mathrm{ext}, j}-\mathbf{K}_{1} \mathbf{U}_{1}^{j-1, p}-\mathbf{L}_{1}^{\mathrm{T}} \lambda^{j}, \\ \mathbf{M}_{2} \ddot{\mathbf{U}}_{2}^{m}=-\mathbf{P}_{2}(\varepsilon_{0}, \mathbf{E}_{0}, \Sigma_{0}, \tilde{\Sigma}_{0})-(\mathbf{C}_{2}+\tilde{\mathbf{C}}_{2}) \dot{\mathbf{U}}_{2}^{m-1 / 2}-\left(\mathbf{K}_{2}+\tilde{\mathbf{K}}_{2}\right) \mathbf{U}_{2}^{m}-\bar{\mathbf{K}}_{2} \bar{\mathbf{U}}_{2}^{m}-\mathbf{L}_{2}^{\mathrm{T}} \lambda^{m}, \\ \mathbf{L}_{1} \dot{\mathbf{U}}_{1}^{j}+\mathbf{L}_{2} \dot{\mathbf{U}}_{2}^{j}=0. \end{array} \right. \qquad (69) $$

Here, it has to be noted that the known part of the internal force in the PML remains unchanged with respect to the previous implicit formulation. In contrast to the previous implicit case, the displacement $\mathbf{U}_2^m$, the mid-step velocity $\dot{\mathbf{U}}_2^{m-1/2}$, and the integral in time of the displacement $\bar{\mathbf{U}}_2^m$ are known using the classical central difference scheme given by

$$ \dot{\mathbf{U}}_{2}^{m-1 / 2}=\dot{\mathbf{U}}_{2}^{m-3 / 2}+\Delta t \ddot{\mathbf{U}}_{2}^{m-1}, \qquad (70) $$

$$ \mathbf{U}_{2}^{m}=\mathbf{U}_{2}^{m-1}+\Delta t \dot{\mathbf{U}}_{2}^{m-1 / 2}, \qquad (71) $$

$$ \bar{\mathbf{U}}_{2}^{m}=\bar{\mathbf{U}}_{2}^{m-1}+\Delta t \mathbf{U}_{2}^{m}. \qquad (72) $$

In (69), the damping term, given by $(\mathbf{C}_2 + \tilde{\mathbf{C}}_2)\dot{\mathbf{U}}_2^{m-1/2}$, is classically taken into account at the mid-step velocity. This time lag in the velocity is introduced to avoid solving the system by a non-diagonal matrix when the central difference scheme is adopted [32]. In fact, it can be observed that (69) only involves the mass matrix, which can be lumped so as to speed up time stepping. It is important to note that only the mass matrix $\mathbf{M}_2$ is lumped in our proposed explicit version of the 3D PML and the other PML matrices are kept consistent. This approach is called the “intermediate formulation” in Basu’s work [16].
### 4.4. Interface problem

In the following, we deal with both implicit and explicit versions of the 3D PML. Following along the lines of the coupling GC method [28, 30], the kinematic quantities are divided into two parts: free and linked quantities. The free quantities, denoted by $\dot{\mathbf{U}}_1^{\text{free},j}$ and $\dot{\mathbf{U}}_2^{\text{free},j}$, are calculated by taking into account internal and external forces without considering interface forces. The linked quantities are obtained from the interface loads given by the Lagrange multiplier vector $\boldsymbol{\lambda}$. The kinematic quantities of the subdomain $\Omega_2$ at $t_j$ $\dot{\mathbf{U}}_2^{\text{free},j}$ are interpolated between the free quantities at the beginning and at the end of the large time step. It can be shown [28, 30] that the velocity continuity at the interface, given in (56), leads to a reduced-size interface problem whose unknown parameters are the Lagrange multipliers at the fine time scale:

$$ H\boldsymbol{\lambda}^j = \mathbf{b}^j. \qquad (73) $$

Here, the interface operator and the right-hand-side vector are defined as follows:

$$ \left\{ \begin{array}{l} \mathbf{H} = \gamma_1 \Delta t \mathbf{L}_1 \mathbf{M}_1^{-1} \mathbf{L}_1^T + \gamma_2 \Delta t \mathbf{L}_2 \tilde{\mathbf{M}}_2^{-1} \mathbf{L}_2^T, \\ \mathbf{b}^j = \mathbf{L}_1 \dot{\mathbf{U}}_1^{\text{free},j} + \mathbf{L}_2 \dot{\mathbf{U}}_2^{\text{free},j}. \end{array} \right. \qquad (74) $$

The parameter $\tilde{\mathbf{M}}_2$ is the effective stiffness matrix used in the implicit computation. The explicit computation follows along the same lines by replacing $\tilde{\mathbf{M}}_2$ with the lumped mass matrix $\mathbf{M}_2$. The interface operator $\mathbf{H}$ is called the Steklov-Poincaré operator, which can be viewed as the condensed effective stiffness matrix of degrees of freedom belonging to the interface between the two subdomains. The right-hand-side vector $\mathbf{b}^j$ only depends on the free velocities computed in both subdomains without considering the interface forces. It can be seen as a predictor value projected to the degrees of freedom belonging to the interface. The detailed computational stages are presented in Figure 4.

## 5. Numerical examples

### 5.1. Numerical test of a semi-infinite 3D elastic bar

The numerical model of a semi-infinite elastic bar subjected to horizontal displacement at the free end is set up as shown in Figure 5. It simulates the propagation of P waves from a non-dissipative elastic medium to a PML medium. The soil subdomain is assumed to be linear elastic with the following material characteristics: $\rho_1 = 1700 \text{ kg/m}^3$, $E_1 = 10 \text{ MPa}$, and $v_1 = 0.24$. The velocity of P waves $C_p$ is 83 m/s. The same material characteristics are applied in the PML subdomain.

To investigate the effect of $\gamma_0$ on the accuracy, the model is composed of a soil subdomain of 300 m and a PML subdomain with $R_{\text{attenuation}}$ being equal to 0.01 and $m_{\text{PML}}$ being equal to 2. Different PML lengths from 10 m to 300 m are investigated, leading to different $\gamma_0$ values on the basis of the general design formula given in (24). The observation point C is located 20 m from the left end of the model. The simulation is conducted by using a homogeneous time step in both subdomains. The results are compared to the reference results using an extended mesh.

To distinguish the difference between the PML results and the reference results, the error in the PML solution is computed with respect to the reference results from the extended mesh as follows:

$$ \text{Error (\%)} = \frac{\max_n |u_p(t_n) - u_{\text{ref}}(t_n)|}{\max_n |u_{\text{ref}}(t_n)|}. \qquad (75) $$

Here, $u_{\text{ref}}$ and $u_p$ are the displacements of the extended mesh model and the PML model, respectively. In fact, one part of the error arises from the reflections at the interface between the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The algorithm for multi-subdomain coupling at the initialization stage and over a large time step.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Numerical model of a semi-infinite elastic 3D bar subjected to horizontal displacement.</figcaption></figure>

soil and the PML subdomain because of spatial discretization. The other part arises from the end of the PML model depending on $R_{attenuation}$. Therefore, the error in the PML solution is the maximum reflection between the two parts with respect to the amplitude of the incident wave. It is called the maximum numerical reflection coefficient.

Non-harmonic waves are investigated by considering a Ricker incident wave defined as follows:

$$ \text{Ric}(t, t_p, t_s) = A \left( 2\pi^2 \frac{(t-t_s)^2}{t_p^2} - 1 \right) \exp \left( -\pi^2 \frac{(t-t_s)^2}{t_p^2} \right). \quad (76) $$

In Figure 6, the Ricker wave plotted in the time and frequency domains has three parameters: the fundamental period $t_p$, the time shift $t_s$, and the amplitude $A$. The chosen values are as follows: $t_p = 3$ s, $t_s = 3$ s, and $A = 1$. Thus, the finite element size of the eight-node hexahedral
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Waveform and Fourier transform of the Ricker wavelet.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** $\gamma_0$ and maximum numerical reflection coefficient as a function of PML length.</figcaption></figure>

elements in the longitudinal direction composing the mesh (displayed in Figure 5) has to be designed to accurately reproduce the propagation of the input Ricker wave. To achieve sufficient accuracy, the finite element size for linear finite elements should respect the following relationship: $L_{EF} < \lambda_{min}/20$, where $\lambda_{min}$ is the minimal wavelength.

On the left of Figure 7, the value of the attenuation coefficient $\gamma_0$, given in the design equation of the PML in (24), is plotted as a function of the PML length for $R_{attenuation}$ equal to 0.01 and $m_{PML}$ equal to 2. On the right of Figure 7, the maximum numerical reflection is plotted as a function of the PML length. It is obvious that the greater the length of the PML subdomain, the smaller the necessary value of $\gamma_0$. It is also clearly observed that with larger length and smaller $\gamma_0$, the maximum numerical reflection coefficient decreases and better accuracy can be achieved. In fact, the larger the length, the more the elements existing in the PML to describe the attenuation of the waves. Therefore, less spurious reflection is produced at the interface between the soil and the PML subdomain.

The time history of wave propagation at the observation point C with a PML length of 200 m is shown in Figure 8. The first reflection from the interface is 0.99%. The second reflection from the end of the PML subdomain defined by an $R_{attenuation}$ value of 0.01 is 0.78%. In short, the principle behind the design of the PML subdomain is to control the reflections from the interface and the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Time history of wave propagation at the observation point C with a PML length of 200 m.</figcaption></figure>

end of the model. The reflections from the end can be easily controlled by $R_{attenuation}$. In terms of the reflections at the interface, an appropriate value of the length should be applied to obtain satisfactory results.

## 5.2. Numerical test of a semi-infinite 3D elastic bar with a refined mesh in the middle

In the previous semi-infinite elastic bar, a detailed mesh is introduced locally into the middle of the first elastic soil subdomain. The same observation point C is considered. The aim is to investigate, through a very simple test, a representative case of a more complex SSI problem. In this problem, a refined mesh is introduced into the soil mesh so as to model site effects and vibrations of the structure in a localized near-field area. To simulate the refined mesh, the size of the elements in the middle region of the soil subdomain is reduced by a factor of 100. As a consequence, the time step satisfying the CFL condition, which is imposed by the mechanical properties and the finite element size in the elastic soil subdomain, is reduced to $\Delta t_1 = 0.0005$ s. This value is 50 times smaller than that in the previous test. Moreover, the finite element size at the PML interface is kept unchanged. In the fully explicit calculation, the time step in the PML subdomain is substantially reduced and leads to additional computational time. As a result, it is more suitable to use hybrid asynchronous time integration, which is beneficial to optimize the computational time with a fine time step in the soil subdomain and an independent large time step in the other subdomain. The reference results are obtained by fully explicit computation with the extended mesh using the fine time step imposed by the detailed mesh.

Regarding HATI computations, we investigate time step ratios $m$ from 1 to 100 for the implicit version of the PML and from 1 to 50 for the explicit version. The last value corresponds to the CFL condition in the PML ($\Delta t_2 = 50\Delta t_1 = 0.025$ s). In Figure 9(a), the results obtained for the implicit PML are compared to the reference results in terms of displacement at the observation point. It is shown that a high time step ratio can be employed ($m = 50$) while maintaining the reflected spurious wave at a low level. The advantage of the implicit PML is that it enables selecting a time step higher than the CFL condition with a time step ratio equal to 100. In Figure 9(b), the case of the explicit PML is also compared to the reference results. This figure highlights very good behavior of the multi-time-step explicit PML even for a large time step ($m = 50$). The accuracy and wall-clock times of both PML versions are listed in Tables 1 and 2.

In regard to CPU time for different time step ratios, as the time step ratio increases, the CPU time decreases significantly. In other words, the proposed 3D hybrid asynchronous PML turns out to be efficient in terms of accuracy and CPU time thanks to the versatility of the employed
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Time history of wave propagation at the observation point C: (a) multi-time-step implicit PML (E/I); (b) multi-time-step explicit PML (E/E).</figcaption></figure>

<table><caption>**Table 1.** Reflections and CPU times using different time step ratios for E/I co-simulations</caption><thead><tr><th></th><th>Fully explicit</th><th>E/I (m=5)</th><th>E/I (m=10)</th><th>E/I (m=30)</th><th>E/I (m=50)</th><th>E/I (m=100)</th></tr></thead><tbody><tr><th>Reflection (%)</th><td>0.98</td><td>0.98</td><td>0.95</td><td>1.01</td><td>1.39</td><td>2.37</td></tr><tr><th>CPU time (s)</th><td>170.47</td><td>45.60</td><td>28.34</td><td>15.33</td><td>11.97</td><td>9.63</td></tr></tbody></table>

<table><caption>**Table 2.** Reflections and CPU times using different time step ratios for E/E co-simulations</caption><thead><tr><th></th><th>Fully explicit</th><th>E/E (m = 5)</th><th>E/E (m = 10)</th><th>E/E (m = 30)</th><th>E/E (m = 50)</th></tr></thead><tbody><tr><th>Reflection (%)</th><td>0.98</td><td>0.94</td><td>0.90</td><td>1.56</td><td>2.33</td></tr><tr><th>CPU time (s)</th><td>170.47</td><td>40.76</td><td>24.50</td><td>13.03</td><td>10.64</td></tr></tbody></table>

HATI framework. Finally, it should be noted that both PML versions are stable. It is recalled that the explicit version of the PML in Basu's work [16] is unstable when the mass matrix is lumped and the other PML matrices are kept consistent (“intermediate formulation”). Here, the proposed explicit 3D PML, also implemented by the “intermediate formulation”, is stable and endowed with multi-time-step capabilities as shown in HATI computations.

### 5.3. 3D Lamb's test

To evaluate the effectiveness of the hybrid asynchronous PML, Lamb's test has been simulated. In Lamb's test, the concentrated load applied to the surface of an infinite half-space medium generates three types of waves propagating through the soil: P, S, and Rayleigh waves [33]. Consequently, Lamb's test can be considered as a good test for assessing the performance of the PML. Non-harmonic waves are investigated by considering a Ricker incident wave as detailed in Section 5.1. The chosen parameters are $t_p = 3$ s, $t_s = 3$ s, and $A = 2$ MN.

As illustrated in Figure 10, the numerical model is a quarter model of a PML-truncated semi-infinite homogeneous medium subjected to a concentrated force. It is composed of bounded soil (subdomain 1) of size 50 m and a PML (subdomain 2) of thickness 50 m. The same material characteristics are adopted as in the previous numerical model of the semi-infinite elastic bar.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** 3D Lamb's test modeled using PML. The quarter model of a PML-truncated semi-infinite homogeneous medium subjected to a concentrated force.</figcaption></figure>

The size of the eight-node hexahedral element (5 m × 5 m × 5 m) has been taken into account so as to control the inherent wave dispersion. The reference results are computed from an extended mesh. The PML design employs the following parameters: $R_{\text{attenuation}}$ value of 0.01 and $m_{\text{PML}}$ value of 2. A recording point is located on the surface of the subdomain soil at 20 m from each symmetric side.

In the case of a homogeneous time step, the time step satisfying the CFL condition imposed by the mechanical properties in the soil subdomain and the finite element size is applied in both subdomains. Namely, $m = 1$ and $\Delta t_1 = \Delta t_2 = 0.025$ s. We can observe that the displacement obtained by the PML agrees with the reference results. The reflected spurious wave is 1.15% in the X and Z directions and 2.82% in the Y direction as listed in Table 3. In this paper, the bounded soil subdomain and the PML are limited to a size of 50 m, corresponding only to 1/5 of the P wavelength. By using a larger model, better accuracy can be achieved. Using the GC method, the classical second-order Newmark explicit time integration scheme is conserved in the soil subdomain without introducing complex-coordinate-stretched equations in the interior domain. Moreover, thanks to the versatility and stability of the HATI adopted in this paper, it is possible to use a larger time step in the PML domain as is done in the following.

In the case of heterogeneous time steps, the subdomain soil is integrated with a fine time step $\Delta t_1 = 0.025$ s. The PML subdomain is integrated with a large time step $\Delta t_2 = m\Delta t_1$ to reduce the computational time. The time histories of displacements in the three directions at the observation point with different time step ratios $m (\Delta t_2 = m\Delta t_1)$ equal to 1, 3, and 5 are shown in Figure 11. The errors in comparison to the reference results are presented in Table 3. It can be noted that the different curves are quite close and reflections increase as the time step ratio increases. In the X and Z directions, in comparison to the reference results, the amplitude of the spurious wave varies from 1.15% to 5.31% with respect to the amplitude of the incident wave. In the Y direction, the maximum reflection increases from 2.82% to 7.70%. The observed decrease in accuracy as the time step ratio increases can be explained by the following points. First, it is considered that the loss of accuracy can be mainly due to the additional first-order-accurate assumptions made in (43)–(44) for performing the time integration of the complex strain–deformation relation in (30). For a bigger time step, more numerical errors are introduced
<figure><img src="image_14.png" /><figcaption>Figure 11. Displacements recorded at the observation point using different time step ratios.</figcaption></figure>

<table><caption>Table 3. Displacement errors using different time step ratios</caption><thead><tr><th></th><th>Displacement X (%)</th><th>Displacement Y (%)</th><th>Displacement Z (%)</th></tr></thead><tbody><tr><td>m = 1</td><td>1.15</td><td>2.82</td><td>1.15</td></tr><tr><td>m = 3</td><td>3.18</td><td>5.24</td><td>3.18</td></tr><tr><td>m = 5</td><td>5.31</td><td>7.70</td><td>5.31</td></tr></tbody></table>

due to the approximation. Second, the GC coupling algorithm is known to be dissipative when heterogeneous time steps are used between the subdomains, generating spurious waves at the interface. It has been demonstrated that for the GC method, when adopting the same time step, the second order of accuracy is achieved. When different time steps are adopted, the first order of accuracy is achieved due to a slight spurious dissipation at the interface [28, 30].

The kinetic and internal energies of the soil subdomain are computed for different time step ratios *m* as shown in Figure 12. The errors in comparison to the reference results are computed using (75) and are presented in Table 4. It can be observed that the errors are small for different time step ratios and that the errors increase with time step ratio. The CPU times for different time step ratios *m* are listed in Table 5 in a normalized form divided by the CPU time of the homogeneous time step. It shows that with increase in time step ratio, the CPU times decrease
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Time histories of kinetic and internal energies computed using different time step ratios.</figcaption></figure>

<table><caption>**Table 4.** Energy errors using different time step ratios</caption><thead><tr><th></th><th>Kinetic energy (%)</th><th>Internal energy (%)</th></tr></thead><tbody><tr><td><i>m</i> = 1</td><td>0.343</td><td>0.126</td></tr><tr><td><i>m</i> = 3</td><td>0.554</td><td>0.213</td></tr><tr><td><i>m</i> = 5</td><td>0.772</td><td>0.300</td></tr></tbody></table>

<table><caption>**Table 5.** Normalized CPU time for different time step ratios</caption><thead><tr><th></th><th><i>m</i> = 1</th><th><i>m</i> = 3</th><th><i>m</i> = 5</th></tr></thead><tbody><tr><td>Normalized CPU time</td><td>1</td><td>0.334</td><td>0.189</td></tr></tbody></table>

significantly, highlighting the significance of hybrid asynchronous time integration. This implies that by using explicit/implicit co-computation, not only the classical Newmark explicit time integration scheme can be conserved in the soil subdomain without introducing complex-coordinate-stretched equations but also large time steps can be adopted in the PML subdomain for reducing the computational time.

### 5.4. 3D rigid foundation on a layered heterogeneous elastic half-space

The classical SSI problem of a rigid foundation on a heterogeneous half-space is investigated as displayed in Figure 13. The load is defined by a Ricker wave with the same parameters as in Lamb's test. Three different subdomains are considered: the soil medium (subdomain 1), the PML medium (subdomain 2), and rigid foundation (subdomain 3).

The soil subdomain is assumed to be linear elastic and is composed of two layers. The thickness of each layer is 25 m with surface dimensions 50 m × 50 m. A recording point is located on the surface of the subdomain soil at 20 m from each symmetric side to assess the efficiency of the PML layers for modeling an infinite heterogeneous half-space medium. The common material parameters of the soil layers are $ρ_1 = 1700$ kg/m³ and $v_1 = 0.24$. The second layer is characterized by a Young's modulus value that is twice greater than that of the first layer (10 MPa). Similarly, to match the soil subdomain, the interface between the layers has to be taken into account in the PML subdomain around the soil with a thickness of 50 m, leading to two PMLs
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Rigid foundation on a layered soil. The quarter model of a PML-truncated semi-infinite heterogeneous medium subjected to a uniform force.</figcaption></figure>

with the same material properties as those of the two soil layers. The PML is designed by using the following parameters: $R_{attenuation}$ value of 0.01 and $m_{PML}$ value of 2.

The rigid foundation on the soil is characterized by a thickness of 5 m, surface dimensions of 10 m × 10 m, $ρ_3 = 1700$ kg/m³, $v_3 = 0.24$, and $E_3 = 1000$ MPa. The last value is 100 times greater than the Young's modulus value in the soil subdomain. Consequently, the time step satisfying the CFL condition imposed by the mechanical properties of the rigid foundation is 0.0025 s, which is 10 times smaller than the time step required in the soil subdomain. If the same explicit time integration scheme is adopted for the soil subdomain and for the rigid foundation subdomain, the time step in the soil subdomain is reduced, which leads to additional computational time. As a result, it is of great importance to couple the soil subdomain and the rigid foundation by a coupling algorithm and adopt an implicit time integration scheme for the rigid foundation subdomain. Finally, by using the subdomain coupling strategy, three different subdomains are coupled within the multi-time-step explicit/implicit co-simulation. The soil medium is integrated with a time step $Δt_1$. The rigid foundation is integrated using a classical second-order Newmark implicit scheme. The PML is integrated using an extended third-order Newmark implicit scheme. The time step is $Δt_3 = Δt_2 = mΔt_1$, where *m* denotes the time step ratio.

First, we consider the case with the time step ratio *m* = 1 and the time step $Δt_1$ = 0.025 s satisfying the CFL condition imposed by the mechanical properties of the soil subdomain. The displacements recorded at point C are compared with the reference results obtained from an extended mesh. Here, the coupling is only between different time integrators because the same time step size is adopted in all the subdomains. From Figure 15 and Table 6, it can be observed that good agreement is achieved in comparison to the reference results. The reflected spurious wave is 1.11% in the X and Z directions and 2.97% in the Y direction. The snapshots of displacement magnitudes at different times are displayed in Figure 14. The first snapshot at time 3.875 s shows the propagation of the maximum peak of Ricker incident waves. The second snapshot at time 4.35 s shows that the maximum peak of the Ricker incident waves begins to be absorbed in the PML region, which is followed by an additional smaller peak produced in the soil subdomain. The third snapshot at time 5.05 s shows that the maximum peak has been absorbed by the PML and that the absorption of the smaller peak begins. The last snapshot at 5.9 s reveals that the absorption of waves in the PML region is almost completed. No obvious reflections can
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Snapshots of displacement magnitudes at different times.</figcaption></figure>

<table><caption>**Table 6.** Displacement errors using different time step ratios and time steps</caption><thead><tr><th></th><th>Displacement X (%)</th><th>Displacement Y (%)</th><th>Displacement Z (%)</th></tr></thead><tbody><tr><td>m = 1, Δt<sub>1</sub> = 0.025</td><td>1.11</td><td>2.97</td><td>1.11</td></tr><tr><td>m = 3, Δt<sub>1</sub> = 0.025</td><td>6.55</td><td>5.26</td><td>6.55</td></tr><tr><td>m = 3, Δt<sub>1</sub> = 0.0025</td><td>1.51</td><td>1.87</td><td>1.51</td></tr></tbody></table>

be observed at the PML interface or from the boundaries, thereby indicating very satisfactory performance of the PML.

In the case of heterogeneous time steps, two examples are demonstrated. The first is carried out by using a time step ratio $m = 3$ and a time step $\Delta t_1 = 0.025$ s. As presented in Table 6, the accuracy decreases with increasing time step ratio, which is similar to Lamb's test. The other computation is carried out with a finer time step size $\Delta t_1 = 0.0025$ s and a time step ratio $m = 3$. In other words, the time step $\Delta t_1$ is taken as the CFL critical time step in the rigid foundation, corresponding to the time step of fully explicit computation. As shown in Figure 15, the PML results are in good agreement with the reference results, achieving an error of 1.51% in the X and Z directions and 1.87% in the Y direction. In comparison to the first case with the time step size $\Delta t_1 = 0.025$ s and the same time step ratio $m = 3$, this demonstrates that PML accuracy
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Displacements recorded at the observation point using different time step ratios and time steps.</figcaption></figure>

<table><caption>**Table 7.** Normalized CPU time for different time step ratios and time steps</caption><thead><tr><th></th><th>m = 1, Δt<sub>1</sub> = 0.025</th><th>m = 3, Δt<sub>1</sub> = 0.025</th><th>m = 3, Δt<sub>1</sub> = 0.0025</th></tr></thead><tbody><tr><td>Normalized CPU time</td><td>1</td><td>0.287</td><td>3.53</td></tr></tbody></table>

depends on the size of the time step. In short, the size of the time step has an important effect on the accuracy of the PML. The smaller the time step in the PML, the better the accuracy.

In regard to the computational time normalized by the CPU time for the case with the time step size $Δt_1 = 0.025$ s and the time step ratio $m = 1$, which is given in Table 7, an important reduction in computational time can be obtained by using the time step size $Δt_1 = 0.025$ s and the time step ratio $m = 3$. For a finer time step size $Δt_1 = 0.0025$ s corresponding to the time step of fully explicit computation, more time steps should be calculated in the numerical simulation, thereby resulting in a longer computational time.

## 6. Conclusion

A novel 3D PML, which is suitable for finite element implementation, has been proposed in this paper for transient elastodynamics. The displacement-based PML, making use of the unsplit
formulation and the convenient expression of the internal force in the PML domain, is integrated into the heterogeneous (different time integrators) asynchronous (different time steps) time integrator framework. This enables the PML to be dealt with independently by an explicit or an implicit scheme with large time steps while conserving classical finite element formulations in the elastic domain for optimizing computational efficiency.

A bar modeled by the 3D finite element method is demonstrated as the first example, which is followed by more advanced examples such as Lamb's test and SSI problems with PML-truncated semi-infinite homogeneous and heterogeneous media. In addition to the straightforward implementation features in finite element and spectral element software, the proposed 3D hybrid asynchronous PML turns out to be accurate, stable, and efficient in terms of CPU time thanks to the versatility of the employed HATI framework. It is also shown that the stiffest subdomain in the SSI problem, such as a rigid foundation on a multi-layered soil domain, can be dealt with by implicit time integration without impacting the choice of time integrators and time steps in multi-layered soil and PML media. Future work is in progress for applying the developed 3D hybrid asynchronous PML to complex SSI problems such as wave barriers in heterogeneous media.

## Appendix A.

In the following, we summarize the matrices expressed as a function of scaling and attenuation functions $f^e$ and $f^p$ related to the PML.

$$ \underline{F}^e = \begin{bmatrix} 1 + f_1^e(x_1) & 0 & 0 \\ 0 & 1 + f_2^e(x_2) & 0 \\ 0 & 0 & 1 + f_3^e(x_3) \end{bmatrix}, \quad \underline{F}^p = \begin{bmatrix} f_1^p(x_1) & 0 & 0 \\ 0 & f_2^p(x_2) & 0 \\ 0 & 0 & f_3^p(x_3) \end{bmatrix}, \\ \underline{\tilde{F}}^{ee} = \begin{bmatrix} f_{23}^{ee} & 0 & 0 \\ 0 & f_{13}^{ee} & 0 \\ 0 & 0 & f_{12}^{ee} \end{bmatrix}, \quad \underline{\tilde{F}}^{ep} = \begin{bmatrix} f_{23}^{ep} & 0 & 0 \\ 0 & f_{13}^{ep} & 0 \\ 0 & 0 & f_{12}^{ep} \end{bmatrix}, \quad \underline{\tilde{F}}^{pp} = \begin{bmatrix} f_{23}^{pp} & 0 & 0 \\ 0 & f_{13}^{pp} & 0 \\ 0 & 0 & f_{12}^{pp} \end{bmatrix}, \quad (77) $$

where

$$ \begin{cases} f_{ij}^{ee} = [1 + f_i^e(x_i)][1 + f_j^e(x_j)], \\ f_{ij}^{ep} = [1 + f_i^e(x_i)]f_j^p(x_j) + [1 + f_j^e(x_j)]f_i^p(x_i), \\ f_{ij}^{pp} = f_i^p(x_i)f_j^p(x_j). \end{cases} \quad (78) $$

The scalar values from the right-hand side of the equation of motion in (28) are as follows:

$$ \left\{ \begin{aligned} f_M &= [1 + f_1^e(x_1)][1 + f_2^e(x_2)][1 + f_3^e(x_3)], \\ f_C &= [1 + f_1^e(x_1)][1 + f_2^e(x_2)]f_3^p(x_3) + [1 + f_1^e(x_1)][1 + f_3^e(x_3)]f_2^p(x_2) \\ &\quad + [1 + f_2^e(x_2)][1 + f_3^e(x_3)]f_1^p(x_1), \\ f_K &= f_1^p(x_1)f_2^p(x_2)[1 + f_3^e(x_3)] + f_2^p(x_2)f_3^p(x_3)[1 + f_1^e(x_1)] \\ &\quad + f_1^p(x_1)f_3^p(x_3)[1 + f_2^e(x_2)], \\ f_H &= f_1^p(x_1)f_2^p(x_2)f_3^p(x_3). \end{aligned} \right. \quad (79) $$

Then, we introduce the element-wise finite element discretization for the weak form of the internal force terms in (42). The matrices containing shape function derivatives of an eight-node
hexahedral element combined with the previous scaling and attenuation functions are expressed as follows:

$$ \tilde{\mathbf{B}}_{\mathrm{I}}^{\mathrm{ee}} = \begin{bmatrix} \tilde{N}_{11}^{ee} & 0 & 0 \\ 0 & \tilde{N}_{12}^{ee} & 0 \\ 0 & 0 & \tilde{N}_{13}^{ee} \\ \tilde{N}_{12}^{ee} & \tilde{N}_{11}^{ee} & 0 \\ \tilde{N}_{13}^{ee} & 0 & \tilde{N}_{11}^{ee} \\ 0 & \tilde{N}_{13}^{ee} & \tilde{N}_{12}^{ee} \end{bmatrix}, \qquad (80) $$

$$ \tilde{\mathbf{B}}^{\mathrm{ee}} = [\tilde{\mathbf{B}}_{1}^{\mathrm{ee}} \quad \tilde{\mathbf{B}}_{2}^{\mathrm{ee}} \dots \quad \tilde{\mathbf{B}}_{8}^{\mathrm{ee}}]. \qquad (81) $$

The components of the above derivative matrix are given for an index $i = 1,2,3$ without the summation convention:

$$ \tilde{N}_{1i}^{ee} = \tilde{F}_{ii}^{ee} N_{i,i}, \quad \tilde{N}_{1i}^{ep} = \tilde{F}_{ii}^{ep} N_{i,i}, \quad \tilde{N}_{1i}^{pp} = \tilde{F}_{ii}^{pp} N_{i,i}. \qquad (82) $$

In addition, $\tilde{\mathbf{B}}^{\mathrm{ep}}$ and $\tilde{\mathbf{B}}^{\mathrm{pp}}$ are defined similarly by replacing $\tilde{N}_{1i}^{ee}$ with $\tilde{N}_{1i}^{ep}$ and $\tilde{N}_{1i}^{pp}$, respectively.

The Voigt notation is adopted for the stress and strain tensors, thereby giving the following vectors:

$$ \hat{\boldsymbol{\sigma}} = \begin{cases} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \end{cases}, \quad \hat{\boldsymbol{\varepsilon}} = \begin{cases} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \end{cases}. \qquad (83) $$

The constitutive relationship for an isotropic elastic medium is

$$ \hat{\boldsymbol{\sigma}} = \mathbf{D} \hat{\boldsymbol{\varepsilon}}, \qquad (84) $$

where $\mathbf{D}$ is the material constitutive matrix expressed as follows:

$$ \mathbf{D} = \begin{bmatrix} k + 4\mu/3 & k - 2\mu/3 & k - 2\mu/3 \\ k - 2\mu/3 & k + 4\mu/3 & k - 2\mu/3 \\ k - 2\mu/3 & k - 2\mu/3 & k + 4\mu/3 \\ \mu & 0 & 0 \\ 0 & \mu & 0 \\ 0 & 0 & \mu \end{bmatrix}. \qquad (85) $$

Furthermore, additional matrices have to be defined for the strain-deformation relationship given in (45). We express the $\hat{\mathbf{F}}^{\epsilon}$ and $\mathbf{B}^{\epsilon}$ matrices depending on shape function derivatives as well as scaling and attenuation functions:

$$ \hat{\mathbf{F}}^{\epsilon} = \begin{bmatrix} (F_{11}^{\epsilon})^2 & 0 & 0 & 0 & 0 & 0 \\ 0 & (F_{22}^{\epsilon})^2 & 0 & 0 & 0 & 0 \\ 0 & 0 & (F_{33}^{\epsilon})^2 & 0 & 0 & 0 \\ 0 & 0 & 0 & F_{11}^{\epsilon} F_{22}^{\epsilon} & 0 & 0 \\ 0 & 0 & 0 & 0 & F_{11}^{\epsilon} F_{33}^{\epsilon} & 0 \\ 0 & 0 & 0 & 0 & 0 & F_{22}^{\epsilon} F_{33}^{\epsilon} \end{bmatrix}, \qquad (86) $$
$$ \mathbf{B}^{\epsilon} = \begin{bmatrix} F_{11}^{\epsilon} N_{I1}^{l} & 0 & 0 \\ 0 & F_{22}^{\epsilon} N_{I2}^{l} & 0 \\ 0 & 0 & F_{33}^{\epsilon} N_{I3}^{l} \\ F_{11}^{\epsilon} N_{I2}^{l} & F_{22}^{\epsilon} N_{I1}^{l} & 0 \\ F_{22}^{\epsilon} N_{I3}^{l} & 0 & F_{33}^{\epsilon} N_{I1}^{l} \\ 0 & F_{22}^{\epsilon} N_{I3}^{l} & F_{33}^{\epsilon} N_{I2}^{l} \end{bmatrix}, \qquad (87) $$

with the matrices

$$ \underline{F}^l = \left[ \underline{F}^p + \frac{\underline{F}^e}{\Delta t} \right]^{-1}, \quad \underline{F}^e = \underline{F}^e \underline{F}^l. \qquad (88) $$

Moreover, the component in the $\mathbf{B}^\epsilon$ matrix is given by $N_{Ii}^l = F_{ii}^l N_{I,i}$ for $i = 1, 2, 3$. Finally, $\hat{\mathbf{F}}^{\epsilon Q}$ and $\mathbf{B}^Q$ are defined similarly by replacing $\underline{F}^e$ with $\underline{F}^Q$, where $\underline{F}^Q = \underline{F}^p \underline{F}^l$.

## Supplementary data

Supporting information for this article is available on the journal's website under https://doi.org/10.5802/crmeca.59 or from the author.

## References

[1] P. Bettess, "Infinite elements", *Int. J. Numer. Meth. Eng.* **11** (1977), p. 53-64.

[2] A. Houmat, "Mapped infinite p-element for two-dimensional problems of unbounded domains", *Comput. Geotech.* **35** (2008), p. 608-615.

[3] B. Enquist B, A. Majda, "Absorbing boundary conditions for the numerical simulation of waves", *Math. Comput.* **31** (1977), p. 629-651.

[4] D. Kosloff, R. Kosloff, "Absorbing boundaries for wave propagation problems", *J. Comput. Phys.* **63** (1986), p. 363-376.

[5] J. F. Semblat, L. Lenti, A. Gandomzadeh, "A simple multi-directional absorbing layer method to simulate elastic wave propagation in unbounded domains", *Int. J. Numer. Meth. Eng.* **85** (2011), p. 1543-1563.

[6] P. Rajagopal, M. Drozdz, E. A. Skelton, "On the use of the absorbing layers to simulate the propagation of elastic waves in unbounded isotropic media using commercially available finite element packages", *NDT and E Int.* **51** (2012), p. 30-40.

[7] E. Zafati, M. Brun, I. Djeran-Maigre, F. Prunier, "Design of an efficient multi-directional explicit/implicit Rayleigh absorbing layer for seismic wave propagation in unbounded domain using a strong form formulation", *Int. J. Numer. Meth. Eng.* **106** (2015), p. 83-112.

[8] J. P. Bérenger, "A perfectly matched layer for the absorption of electromagnetic waves", *J. Comput. Phys.* **114** (1994), p. 185-200.

[9] W. C. Chew, W. H. Weedon, "A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates", *Microw. Opt. Technol. Lett.* **7** (1994), p. 599-604.

[10] W. C. Chew, Q. H. Liu, "Perfectly matched layers for elastodynamics: a new absorbing boundary condition", *J. Comput. Acoust.* **4** (1996), p. 341-359.

[11] F. Collino, C. Tsogka, "Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media", *Geophysics* **66** (2001), p. 294-307.

[12] T. Wang, X. Tang, "Finite-difference modeling of elastic wave propagation: a nonsplitting perfectly matched approach", *Geophysics* **68** (2003), p. 1749-1755.

[13] R. Matzen, "An efficient finite element time-domain formulation for the elastic second-order wave equation: A non split complex frequency shifted convolutional PML", *Int. J. Numer. Meth. Eng.* **88** (2011), p. 951-973.

[14] U. Basu, A. Chopra, "Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation", *Comput. Meth. Appl. Mech. Eng.* **192** (2003), p. 1337-1375.

[15] U. Basu, A. Chopra, "Perfectly matched layers for transient elastodynamics of unbounded domains", *Int. J. Numer. Meth. Eng.* **59** (2004), p. 1039-1074.

[16] U. Basu, "Explicit finite element perfectly matched layer for transient three-dimensional elastic waves", *Int. J. Numer. Meth. Eng.* **77** (2009), p. 151-176.

[17] LS-DYNA, LS-DYNA Keyword User's Manual, Livermore Software Technology Corporation, 2019.
[18] J. Mamie, W. P. Kikstra, DIANA User's Manual, DIANA FEA BV, 2016, Version 10.1.

[19] M. Brun, E. Zafati, I. Djeran-Maigre, F. Prunier, "Hybrid asynchronous perfectly matched layer for seismic wave propagation in unbounded domains", *Finite Elem. Anal. Des.* **122** (2016), p. 1-15.

[20] S. Kucukcoban, L. F. Kallivokas, "Mixed perfectly-matched-layers for direct transient analysis in 2D elastic heterogeneous media", *Comput. Meth. Appl. Mech. Eng.* **200** (2011), p. 57-76.

[21] A. Fathi, B. Poursartip, L. F. Kallivokas, "Time-domain hybrid formulations for wave simulations in three-dimensional PML-truncated heterogeneous media", *Int. J. Numer. Meth. Eng.* **101** (2015), p. 165-198.

[22] M. Brun, A. Gravouil, A. Combescure, "Two FETI-based heterogeneous time step coupling methods for Newmark and $\alpha$-schemes derived from the energy method", *Comput. Meth. Appl. Mech. Eng.* **283** (2015), p. 130-176.

[23] A. Gravouil, A. Combescure, M. Brun, "Heterogeneous asynchronous time integrators for computational structural dynamics", *Int. J. Numer. Meth. Eng.* **102** (2015), p. 202-232.

[24] Q. H. Liu, "Perfectly matched layers for elastic waves in cylindrical and spherical coordinates", *J. Acoust. Soc. Am.* **105** (1999), p. 2075-2084.

[25] J. M. Carcione, D. Kosloff, "Representation of matched-layer kernels with viscoelastic mechanical models", *Int. J. Numer. Anal. Model.* **10** (2013), p. 221-232.

[26] S. Li, M. Brun, I. Djeran-Maigre, S. Kuznetsov, "Hybrid asynchronous absorbing layers based on Kosloff damping for seismic wave propagation in unbounded domains", *Comput. Geotech.* **109** (2019), p. 69-81.

[27] I. Harari, U. Albocher, "Studies of FE/PML for exterior problems of time-harmonic elastic waves", *Comput. Meth. Appl. Mech. Eng.* **195** (2006), p. 3854-3879.

[28] A. Gravouil, A. Combescure, "A multi-time-step explicit-implicit method for non-linear structural dynamics", *Int. J. Numer. Meth. Eng.* **50** (2001), p. 199-225.

[29] N. M. Newmark, "A method of computation for structural dynamics", *J. Eng. Mech. Div. (ASCE)* **85** (1959), p. 67-94.

[30] A. Combescure, A. Gravouil, "A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis", *Comput. Meth. Appl. Mech. Eng.* **191** (2002), p. 1129-1157.

[31] T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*, Prentice-Hall, Englewood Cliffs, NJ, 1987.

[32] T. Belytschko, W. K. Liu, B. Moran, *Nonlinear Finite Elements for Continua and Structures*, Wiley, New York, 2000.

[33] H. Lamb, *Proceedings of the 38th Royal Society of London*, vol. 72, Royal Society of London, London, 1903, 128-130 pages.