Comptes Rendus Mécanique

Three-dimensional hybrid asynchronous perfectly matched layer


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 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 displacementbased 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 stressdisplacement 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 C. R. Mécanique -2020, 348, no 12, 1003-1030 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.

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 xi : R → C.
The complex coordinates are defined by In the above equation, ω denotes the circular frequency and f p i 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 e i 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 p i on the attenuation in the PML as well as on the wave reflection at the interface between a nondissipative 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 frequencyindependent 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].

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 Equations ( 2)-( 4) constitute the strong form of the 1D PML medium in the frequency domain.Parameters σ and ε are the complex values denoting the stress and the strain, respectively, ρ is the density, and E is Young's modulus for P waves.For S waves, by replacing Young's modulus E by the Sijia Li et al. 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 λ(x), and the elastic constitutive relationship remains intact.
We adopt the scaling and attenuation functions f e = 0, f p = γ, respectively, where γ is a constant positive real value.It gives the expression of the constant complex-valued stretching function The deformation equation and the elastic constitutive relationship are used to replace the stress term in the wave propagation equation as follows: Introducing λ into the equation of motion, The dependence of the complex coordinates on the factor iω 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: 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(ω 0 t − kx)), the expression of the wave number k can be obtained: where v denotes the velocity of P waves or S waves.The expression of the rightward propagating wave in the 1D PML medium is 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: It can be seen that the wave frequency ω 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

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 .
Based on the continuity of displacements and equilibrium of stresses at the interface, we can write From the continuity of displacements at interface (15), we have 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 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 Finally, taking into account (17), we derive the remarkable property of the PML for all frequencies: 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 γ.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 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 realvalued 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]: where γ 0 is a user-tunable scalar parameter.In fact, the larger the value of γ 0 , the larger the discretization error.Namely, more reflected spurious waves are produced at the interface with larger γ 0 values as is shown in the following numerical investigations.The logarithmic decrement of the PML domain δ is obtained by integrating (11) along the thickness of the PML: We define the attenuation coefficient R attenution from the logarithmic decrement: For instance, if we want to reach a logarithmic decrement target of δ = 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 attenution 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 γ 0 can be obtained:

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.

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: where C i j kl are the components of the elastic constitutive tensor.Then, we introduce the following notation for the PML region.Ω PML is the region of the PML bounded by prescribed traction force on Γ N PML 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: 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: 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: 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 U = t 0 u dτ.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 straindeformation relationship identical to Basu's formulation, which involves the strain rate and the time integral of the strain tensor given by All the matrices involved in ( 28) and ( 30) depend on scaling functions f e i (x i ) and attenuation functions f p i (x i ).Their expressions are presented in Appendix A.

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 υ be the test function belonging to an appropriate space.The weak formulation is obtained by integrating over the PML domain: Here, taking into account the scaling and damping functions, the expressions of the modified strain tensors are as follows: The internal force is expressed as follows:

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 where 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: The inertial system matrices M, C, K, and K are assembled from their respective element-level matrices.The element-level matrices are obtained by adopting a quadrature formula in every hexahedral element: Taking into account (35), the internal force term P e int can be written as where the matrices Bee , Bep , and Bpp 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 σ represents the six-component vector of stresses and Σ and Σ 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: Using the assumptions given in (44), the third equation of the system in (30) leads to the expression of the strain εn+1 at the end of the time step: The matrices Fε , F Q , B ε , and B 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 P e n+1 can be written in terms of the element velocity and displacement vectors ( Ue n+1 and U e n+1 ) as well as a term denoted as P(ε e n , E e n , Σ e n , Σe n ) depending only on known quantities at the beginning of the time step.The element-wise internal force P e n+1 is written as where the matrix B is defined as a function of the previous matrices Bee , Bep , and Bpp and the time step as follows: The known part of the internal force at the beginning of the time step is given by 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 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 : 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.

Weak form for subdomain coupling
Let Ω be a bounded domain belonging to R 3 with a regular boundary.The time interval of interest is J = [0, T ].The domain Ω is divided into two parts Ω 1 and Ω 2 as shown in Figure 3: The parameter Γ I denotes the interface between the two subdomains, subdomain Ω 1 representing the non-dissipative medium (the domain of interest) and subdomain Ω 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 Ω 1 is characterized by its density ρ 1 ; Young's modulus E 1 ; Poisson's coefficient ν 1 ; the body force b 1 ; the Dirichlet prescribed displacement on Γ D 1 , u D 1 ; and the traction force at the Neumann condition on The subdomain Ω 2 is characterized by its density ρ 2 ; Young's modulus E 2 ; Poisson's coefficient ν 2 ; the body force b 2 ; the Dirichlet prescribed displacement on Γ D 2 , u D 2 ; and the traction force at the Neumann condition on Before writing the weak form of the coupled problem in Ω divided into two parts Ω 1 and Ω 2 , we classically introduce test functions v 1 and v 2 , which belong to appropriate spaces: , where d is the space dimension equal to 3. The solutions belong to the following spaces: 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 Γ 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.Next, using a dual Schur formulation, we can write the principle of virtual power for transient dynamics.We find the solution u 1 (t ) ∈ W 1 , u 2 (t ) ∈ W 2 , and λ(t ) ∈ Q for which the following weak form is satisfied: 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 where L 1 and L 2 are the constraint matrices of Boolean type in the case of matching meshes at the interface Γ I 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 Ω 1 is integrated independently in time by a second-orderaccurate Newmark explicit time integration scheme.Moreover, the subdomain Ω 2 is treated by an extended third-order-accurate Newmark implicit time integration scheme [21] or a secondorder explicit time integration method following the central difference scheme.

Multi-time-step implicit PML
As illustrated in Figure 3, an explicit time integrator with a fine time step ∆t 1 imposed by the Courant-Friedrichs-Lewy (CFL) condition [31] is adopted for the subdomain Ω 1 .An implicit time integrator with a large time step ∆t 2 is used for the subdomain Ω 2 because the implicit scheme is unconditionally stable with ∆t 2 = m∆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 ∆t 2 = [t 0 ; t m ], while the equilibrium of subdomain 1 is attained at the end of every fine time step ∆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: The first equation is the discrete in space equation of motion of the subdomain Ω 1 defined at the end of the fine time step ∆t 1 = [t j −1 ; t j ].The second equation is the discrete in space equation of motion of the subdomain Ω 2 , corresponding to the PML medium, defined at the end of the large time step ∆t 2 = [t 0 ; t m ].The third equation is the discrete velocity continuity at the interface.The subdomain Ω 1 is integrated in time by a Newmark explicit scheme (β 1 = 0 and γ 1 = 1/2) with a lumped mass matrix M 1 .We define U j−1,p 1 as the predictor displacement and Uj−1,p 1 as the predictor acceleration, which are classically introduced in approximate Newmark formulas: 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: Regarding the subdomain Ω 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, β 2 and γ 2 are the usual Newmark parameters related to the classical constant average acceleration scheme.
C. R. Mécanique -2020, 348, no 12, 1003-1030 These are equal to 1/4 and 1/2, respectively; α 2 is an additional parameter required for the thirdorder extended Newmark scheme, which is equal to 1/12.Velocities, displacements, and time integrals of displacement are expressed as functions of predictors as follows: Here, the predictors U0,p 2 , U 0,p 2 , and U 0,p 2 are defined as follows: 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: where the effective stiffness matrix M2 in the PML subdomain is defined as

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 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 U m 2 , the mid-step velocity Um−1/2 2 , and the integral in time of the displacement U m 2 are known using the classical central difference scheme given by Um−1/2 In (69), the damping term, given by (C 2 + C2 ) Um−1/2

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 nondiagonal 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 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].

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 Ufree,j 1 and Ufree,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 λ.The kinematic quantities of the subdomain Ω 2 at t j Ufree,j 2 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: Here, the interface operator and the right-hand-side vector are defined as follows: The parameter M2 is the effective stiffness matrix used in the implicit computation.The explicit computation follows along the same lines by replacing M2 with the lumped mass matrix M 2 .The interface operator 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 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.

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 nondissipative elastic medium to a PML medium.The soil subdomain is assumed to be linear elastic with the following material characteristics: ρ 1 = 1700 kg/m 3 , E 1 = 10 MPa, and ν 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 γ 0 on the accuracy, the model is composed of a soil subdomain of 300 m and a PML subdomain with R attenuation being equal to 0.01 and m PML being equal to 2. Different PML lengths from 10 m to 300 m are investigated, leading to different γ 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:   Non-harmonic waves are investigated by considering a Ricker incident wave defined as follows: 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  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 E F < λ min /20, where λ min is the minimal wavelength.
On the left of Figure 7, the value of the attenuation coefficient γ 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 γ 0 .It is also clearly observed that with larger length and smaller γ 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 attenution 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  end of the model.The reflections from the end can be easily controlled by R attenution .In terms of the reflections at the interface, an appropriate value of the length should be applied to obtain satisfactory results.

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 ∆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 (∆t 2 = 50∆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    [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.

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 semiinfinite 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.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 attenution value of 0.01 and m 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 ∆t 1 = ∆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 ∆t 1 = 0.025 s.The PML subdomain is integrated with a large time step ∆t 2 = m∆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 (∆t 2 = m∆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-orderaccurate 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 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   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 complexcoordinate-stretched equations but also large time steps can be adopted in the PML subdomain for reducing the computational time.

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 3 and ν 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 with the same material properties as those of the two soil layers.The PML is designed by using the following parameters: R attenution 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 3 , ν 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   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 ∆t 1 = 0.0025 s and a time step ratio m = 3.In other words, the time step ∆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 ∆t 1 = 0.025 s and the same time step ratio m = 3, this demonstrates that PML accuracy  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.

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 C. R. Mécanique -2020, 348, no 12, 1003-1030 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.The scalar values from the right-hand side of the equation of motion in (28) are as follows:  The components of the above derivative matrix are given for an index i = 1, 2, 3 without the summation convention: In addition, Bep and Bpp are defined similarly by replacing Ñ ee I i with Ñ ep I i and Ñ pp I i , respectively.The Voigt notation is adopted for the stress and strain tensors, thereby giving the following vectors: The constitutive relationship for an isotropic elastic medium is where D is the material constitutive matrix expressed as follows: Furthermore, additional matrices have to be defined for the strain-deformation relationship given in (45).We express the Fε and B ε matrices depending on shape function derivatives as well as scaling and attenuation functions: with the matrices Moreover, the component in the B ε matrix is given by N l I i = F l i i N I ,i for i = 1, 2, 3. Finally, F εQ and B Q are defined similarly by replacing F ϵ with F Q , where F Q = F p F l .

Figure 1 .
Figure 1.Wave propagation from elastic medium to PML medium.

Figure 2 .
Figure 2. Evolution of the damping functions in PML subdomain.
defining the decomposition of the boundary conditions into Dirichlet and Neumann conditions.In addition, g N denotes the C. R. Mécanique -2020, 348, no 12, 1003-1030

2
)Here, u 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 C. R. Mécanique -2020, 348, no 12, 1003-1030

Figure 4 .
Figure 4.The algorithm for multi-subdomain coupling at the initialization stage and over a large time step.

Figure 5 .
Figure 5. Numerical model of a semi-infinite elastic 3D bar subjected to horizontal displacement.

Figure 6 .
Figure 6.Waveform and Fourier transform of the Ricker wavelet.

Figure 7 .
Figure 7. γ 0 and maximum numerical reflection coefficient as a function of PML length.

Figure 8 .
Figure 8.Time history of wave propagation at the observation point C with a PML length of 200 m.

Figure 9 .
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).

Figure 10 .
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.

Figure 11 .Table 3 .
Figure 11.Displacements recorded at the observation point using different time step ratios.

Figure 12 .
Figure 12.Time histories of kinetic and internal energies computed using different time step ratios.

Table 4 .Table 5 .
Energy errors using different time step ratios Kinetic energy (%) Internal energy (%) Normalized CPU time for different time step ratios m = 1 m = 3 m = 5

Figure 13 .
Figure 13.Rigid foundation on a layered soil.The quarter model of a PML-truncated semiinfinite heterogeneous medium subjected to a uniform force.

Figure 14 .Table 6 .
Figure 14.Snapshots of displacement magnitudes at different times.Table 6. Displacement errors using different time step ratios and time steps

Figure 15 .
Figure 15.Displacements recorded at the observation point using different time step ratios and time steps.

Table 7 .
Normalized CPU time for different time step ratios and time steps m = 1, ∆t 1 = 0.025 m = 3, ∆t 1 = 0.025 m = 3, ∆t 1 = 0 j = [1 + f e i (x i )][1 + f e j (x j )], f ep i j = [1 + f e i (x i )] f p j (x j ) + [1 + f e j (x j )] 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

Table 2 .
Reflections and CPU times using different time step ratios for E/E co-simulations 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 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 C. R. Mécanique -2020, 348, no 12, 1003-1030