Towards a combined perfectly matching layer and infinite element formulation for unbounded elastic wave problems

This paper presents a framework for implementing a novel perfectly matching layer and infinite element (PML+IE) combination boundary condition for unbounded elastic wave problems in the time domain. To achieve this, traditional hexahedral finite elements are used to model wave propagation in the inner domain and IE test functions are implemented in the exterior domain. Two alternative implementations of the PML formulation are studied: the case with constant stretching in all three dimensions and the case with spatially dependent stretching along a single direction. The absorbing ability of the PML+IE formulation is demonstrated by the favourable comparison with the reflection coefficient for a plane wave incident on the boundary achieved using a finite-element-only approach where stress free boundary conditions are implemented at the domain edge. Values for the PML stretching function parameters are selected based on the minimisation of the reflected wave amplitude and it is shown that the same reduction in reflection amplitude can be achieved using the PML+IE approach with approximately half of the number of elements required in the finite-element-only approach.


Introduction
In numerical modelling and simulation, where resources are limited by memory and processing power [1], it is useful to study a small spatial domain containing the region of interest and simply assume that the outer domain extends to infinity [2]; for example, in aerospace problems, the infinite domain could be atmospheric air while the region of interest relates only to the localised flow around an airplane wing [3]. In cases such as this, it is necessary to employ a boundary condition on the exterior of the computational domain to prevent outgoing In the general, three-dimensional formulation, the domain of interest F is modelled by standard FEs (shaded here) and bordered in all directions by elements modelled using PML+IE combinations; this domain is denoted I . The schematic on the right represents a cross-section in any plane of the three-dimensional case.
waves reflecting back into the region of interest. A number of methods exist to deal with these reflections within numerical simulations such as absorbing boundary conditions (ABCs) [4], perfectly matching layers (PMLs) [5] and infinite elements (IEs) [6]. A review of the use of such boundary conditions within the finite element method (FEM) for time-harmonic acoustics is carried out in [7][8][9][10][11].
The PML technique (first presented by Berenger [5]) is based on the use of an absorbing layer, with the matching medium designed to absorb without reflection and prevent any wave travelling back into the computational domain. While Berenger dealt with electromagnetic waves, Chew and Liu [12] proved that there exists a fictitious elastic PML half-space in solids which completely absorbs elastic waves in spite of the coupling between compressional and shear wave modes. In the same year, Lyons et al. [13] demonstrated the accuracy and potential of the PML in a finite element (FE) formulation. The work of Liu and Tao [14] and Qi and Geers [15] extended the PML to simulate acoustic wave propagation in absorptive media. The PML was also shown to be effective in numerically solving the Helmholtz equation in [16,17] and was developed further for time harmonic elastodynamics in [18,19]. Methods for modelling the elastic wave equation with PMLs in both the time and frequency domain are examined in [20][21][22][23][24][25][26]. Examination of PMLs in the specific case of a forced, open, elastic waveguide has been carried out in [25], where PMLs were implemented in the radial direction (the direction orthogonal to the wave propagation direction) along the length of the waveguide to model a buried bar. The modal analysis that this allowed investigated the roll of this radial PML on the wavenumber spectrum, the attenuation as a function of frequency, and the form of the mode shapes.
An alternative approach to truncating the FE mesh is to employ IEs at the boundary. These IEs essentially extend the element domain to infinity, and are based on the shape functions used in the interior elements which are multiplied by an appropriate decay function to achieve the desired behaviour at infinity. The IE method was discussed in detail by Bettess [6] while Astley provides a review of IE formulations with various element types and assesses their accuracy in [27]. IEs have been applied to the acoustic wave equation in the frequency domain [28][29][30][31][32] and in the time domain [33][34][35], as well as to the elastic wave equation [36,37]. Transient IEs have also been explored [38,39].
Building on this existing body of work, the authors developed a combined PML+IE method for unbounded acoustic problems in the frequency domain and assessed its performance for a spherical resonator [40]. This approach was extended in [41] to consider the scalar wave equation in the time domain. In the present paper, a method for studying elastodynamic waves in a three-dimensional, homogeneous volume with the PML+IE formulation implemented at the boundary (see Figure 1) is presented. This work is believed to be the first FE implementation of a combined PML+IE formulation for the vector elastic wave equation in the time domain. A numerical study of a three-dimensional, homogeneous volume placed in a vacuum with PML+IE boundary conditions applied at the right-most face of the domain (the face x 1 = L 1 shown in Figure 2) is conducted, allowing an empirical comparison of the method to a FE-only approach where stress free boundary conditions are employed at this boundary. In contrast to the elastodynamic wave problem considered in [25], which is developed specifically for waveguides where Dirichlet boundary conditions are implemented on the outer edge of the PML domain and the system is solved using an implicit scheme, in this paper the focus is on examining the time-domain reflection coefficient from the end of a three-dimensional volume in the direction of wave propagation, owing to the PML+IE formulation. Here, an IE formulation is implemented in conjunction with the PML at the open end of the three-dimensional volume and standard stress-free boundary conditions are The interior domain F is of fixed length L 1 and is meshed using standard FEs, depicted here by the shaded region. The exterior domain I is implemented at the right end of the three-dimensional volume (the open, semi-infinite part of the domain) and uses the PML+IE combination. The face X at x 1 = X 1 is a notional face that will later be allowed to tend to infinity. applied along its length. In addition, this paper opts for an explicit FE formulation, choosing to compromise some numerical accuracy in return for increased computational efficiency. This will facilitate the use of the method in three-dimensional problems with PML+IE boundary conditions at all domain boundaries.
We begin with the elastodynamic wave equation and a Fourier transform in time is taken in order to introduce appropriate coordinate stretching transformations. The variational formulation is used to introduce IEs under the assumption that the material is locally isotropic. The problem is then discretised with traditional hexahedral FEs modelling the inner domain and IE test functions employed in the exterior domain. Two scenarios are considered: the case with constant PML stretching in all three directions and the case with spatially dependent stretching in one direction. To arrive at an explicit time-domain formulation, diagonalisation of the velocity coefficient matrix is implemented. A reflection coefficient is defined and used to compare the new PML+IE formulation with an extended FE-domain-only implementation (where stress-free boundary conditions are implemented in place of the PML+IE formulation). Values for the PML stretching function parameters are then selected based on the minimisation of the reflected wave amplitude and it is shown that to achieve the same reduction in reflected wave amplitude as observed in the PML+IE case, twice the number of elements (and, thus, twice the runtime) is required in the FE-only implementation.

Geometry and governing equations
Consider the problem of an open-ended, homogeneous, three-dimensional volume with the geometry shown in Figure 2, where X is a notional face at x 1 = X 1 which will later be allowed to tend to infinity, F is the inner domain, modelled by conventional FEM techniques, and I is the outer domain, modelled by PML and IE combinations. The elements which lie on the left-hand boundary of the domain are excited by some time-limited input function (depicted by the arrows in Figure 2). The governing elastodynamic equations are [42] where ρ is density, v is velocity, σ ij denotes the stress tensor and with C ijkl the stiffness tensor. We now define a transformation where Here α j , β j are functions which can be used to fine tune the stretching function and ω is the angular frequency. By writing (1) and (2) in terms of the stretched coordinatesx j (assuming v i (x, 0) = 0 and σ ij (x, 0) = 0), and taking a Fourier transform in time we can write Stress-free boundary conditions are implemented on all faces except X where the Sommerfeld radiation condition is employed, where = O 1/X 2 1 and k = ω/c is the wavenumber (here c can be chosen as either the wave speed of a compressional or shear wave). Thus, by multiplying both sides of (5) by a test function w, integrating over the whole domain = F ∪ I and applying the divergence theorem, we can eventually write Note, to examine the case of a clamped domain (that is, where zero velocity boundary conditions are applied along the length of the domain), we would instead multiply (6) by our test function w, integrate and apply the divergence theorem. This would provide a different start point for the forthcoming analysis and so is beyond the scope of this paper, but would make for an interesting extension in future studies.
To discretise equation (8), . Note that in these expressions, subscript j refers to the node in the FE discretisation and φ j denotes the basis functions given by The test function w is also replaced by a series of test functions θ i of compact support given by The functional form of such basis functions will be given in the following sections. Assuming an isotropic stiffness tensor, equation (8) gives rise to the equations where λ and μ denote the Lamé constants.
In order to progress, some decisions must be made about the stretching function. Two scenarios will be considered: the first has constant stretching in all three directions, that is, s 1 = s 2 = s 3 ≡ s; the second scenario has stretching in only one direction, retaining some spatial dependency, where s 2 = s 3 ≡ 1 and s 1 is a function of x 1 . Note, however, that in both cases s is still a function of frequency.

Constant stretching
To simplify the integrals which appear in (11)-(13), we choose first to implement a constant stretching approach, where s j is independent of x j , that is, s 1 = s 2 = s 3 ≡ s. We define the basis functions as functions of a local coordinate system (ξ , η, ζ ) centred on the node j where L 1 is the fixed length of the interior domain. Next, a parameterisation in terms of these local coordinates of the FEs and IEs is introduced. In I , this parameterisation is only used in the x 2 and x 3 directions. The mappings are shown in Figures 3 and 4 and the test functions are chosen as where the choice of w i is based on the Astley-Leis IE [27]. Now with the mappings shown in Figure 3, we assume that instead where x i for i = 1, 2, 3 denote the element sidelengths. Hence, the N i basis functions have value 1 at node i and value 0 at all other nodes, and form a continuous function whose support is in the elements of which node i is a vertex.  It is also worth noting that, in this work, the assumption has been made that the stress components at a local node i , given byγ ni for n = 1, . . . , 6, are equal to the stress components at local node j , given byγ nj for n = 1, . . . , 6, when nodes i and j belong to the same element. Therefore,γ ni =γ nj ≡ψ n for i , j ∈ e , where e is either a FE or IE. Consequently, the stress componentγ ni at a local node i , is replaced by the stress component for the element under consideration,ψ n . Now, for computational speed, we want to derive an explicit scheme to solve the discretised elastodynamic equations. Currently, the discretised system of (11)-(13) lead to an implicit set of algebraic equations in the unknowns. Deriving an explicit scheme would require the inversion of a very large coefficient matrix which could only be conducted numerically, would be computationally expensive and demand increased memory requirements. An alternative approach is to approximate this matrix by a diagonal one whose inversion is trivial and, hence, facilitates the derivation of a computationally efficient explicit integration scheme. This approximation is made by summing the entries within each row of the velocity coefficient matrix and replacing the diagonal entry with this sum, setting all other entries to zero. This is predicated on the assumption that there are no sharp changes in the velocities and so adjacent nodes have very similar values and hence one can approximate the value at one node by the value at its neighbour. Specifically, in this constant stretching scenario, this approximation is subject to the conditions α μ/ρc 2 and α (c p /c s ) 2 /4, where c s is the shear wave speed and c p is the longitudinal wave speed. Physically, from the first condition we observe that, if c is set equal to the shear wave speed c s , then α 1. This is in agreement with the second condition when c s ≈ c p /2, which is a typical ratio between the longitudinal and shear wave speeds in many engineering materials. Through the diagonalisation of the velocity coefficient matrix, a frequency domain approximation for each FE and IE can eventually be reached [43]. The left-hand side of (11)-(13) for each element then become (after taking inverse Fourier transforms) and where and where · denotes the ceiling function. Similarly, for each IE, the right-hand side of (11)-(13) become and where α, β are the stretching function parameters from equation (4) (which have no spatial dependence for this constant stretching case). To obtain an explicit form, Euler's method is used and, for example, from (20) we can write and from (24), we obtain To calculate the velocities at a global level, the FE and IE approximations can be recombined n is the number of FEs that share global node n as a vertex, b (I) n is the number of IEs that share global node n as a vertex, F (F) is a vector whose entries are the combination of stress terms used in (20)- (22) when global node n is local node i in a FE, and F (I) is a vector whose entries are the combination of stress terms used in (24)- (26) when global node n is local node i in an IE. Thus, the evolution of the velocity can be written in the form where In a similar manner, the FE stress equations which are derived from (6) yield the explicit forms For each IE we have and

Spatially dependent stretching
As attention is restricted in this paper to the reflection coefficient from the x 2 -x 3 plane (see Figure 2), then the PML will only stretch along the x 1 -axis and so we can set s 2 = s 3 = 1. The basis functions, test functions and mappings are still defined as in the constant stretching case and the expressions for the FE given by (20)-(22) remain unchanged. A lengthy derivation given in [43] eventually gives rise to the IE time-domain approximations and Here where we choose and whereᾱ is some constant that can be used to fine tune the PML (see [43] for further details and justification). Note that m and n must be chosen carefully so that I * 1 and I * 2 converge to non-zero constants as X 1 → ∞. As before, there exists a suitable parameter regime which justifies the diagonalisation of the velocity coefficient matrix (which is necessary to arrive at these explicit forms of A). It transpires that I * 2 is subject to two conditions: I * 2 ρL 1 c/3μ and I * 2 4L 1 μ/3c(λ + 2μ). By non-dimensionalising with respect to the spatial variable so thatx 1 = x 1 /L 1 , equation (45) can be rewritten Here we see that the factor L 1 /c is common to both sides of the second inequality, and so the integral in I * 2 is constrained by the material properties of the volume (similar to the constant stretching case) and, hence, must be much smaller than 4c 2 s /3c 2 p . On further inspection, it is shown in [43] that the integral term in equation (45) can be approximated by a function of the free parameterᾱ. By choosingᾱ > 1, we can ensure that both inequalities hold in the materials of interest (for example, steels and other common engineering metals).
Once again, by implementing Euler's method, an explicit form for each IE can be obtained, for example, equation (41) becomes

Now recombining the elements in order to calculate the velocities at a global level gives
where b (F) n is the number of FEs that share global node n as a vertex, b (I) n is the number of IEs that share global node n as a vertex, F (F) is a vector whose entries are the combination of stress terms used in (20)- (22) when global node n is local node i in a FE, and F (I) is a vector whose entries are the combination of stress terms used in (41)- (43) when global node n is local node i in an IE. Thus, we have Letting then the velocity can be written in the form The stress equations can be calculated on an element by element basis, therefore, for each FE, the stress components are given by (29)- (34). For each IE, we can write and where

Results
Having derived PML+IE formulations with constant stretching and spatially dependent stretching, each case was implemented in an explicit FE code using Fortran (see Appendix A for the implementation). A comparison is made between the PML+IE formulation and the FE-only implementation (where the same elastodynamic problem is modelled using hexahedral FEs with stress-free boundary conditions implemented at x 1 = L 1 ) using a reflection coefficient p refl , defined as the ratio of the returning reflected wave to the incident wave measured at a node in the centre of the volume. It is given by where, for a node at position x * = (a x 1 , b x 2 , b x 3 ), |p(t)| , where t i is the number of timesteps taken for the incident wave to first reach this node and t r is the number of timesteps taken for the reflected wave to return to this node. The p max values are calculated over a window in time of size 2δt n to ensure that the arrival time of the maximum amplitude of the wavefront is accurately captured. We then use this reflection coefficient to find values of the stretching function parameters that minimise the amplitude of the reflected wave. In this work, a steel block (c = 5690 m/s, ρ = 7700 kg/m 3 , λ = 97.3 GPa, μ = 76 GPa) that has n 1 = 51 nodes in length and n 2 = n 3 = 101 nodes in width and height, is considered, where the nodes are spaced equally ( x 1 = x 2 = x 3 = 1 × 10 −5 m) and the element side length is chosen to be less than the wavelength/10 (which is shown to be sufficient for satisfactory convergence when choosing first-order linear elements [44])

Constant stretching
In the case of constant stretching in the PML function, there are two parameters to optimise, namely α and β from (4). To minimise the reflection at the domain boundary, we study the effects on the reflected amplitude at a fixed point in the spatial domain as these stretching function parameters are varied. The effect of α is assessed in Figure 5. As α appears in the IE weighting W (I) in (28), it is assumed that it should be of the order x 1 /L 1 so that W (I) is similar to W (F) . However, to justify diagonalisation of the velocity coefficient matrix we know that the condition α (λ + 2μ)/4μ must hold and so small values close to unity are tested. Figure 5 shows that the reflection coefficient increases as α increases and so the best choice for the stretching function parameter α is given by a value close to 1 (we choose α = 1.001). The effect of β is assessed in Figure 6. As β appears in the velocity update in (27), it is assumed that the coefficient of p (t) n should be between 0 and 1. It is clear in Figure  6 that smaller values of β produce the lowest reflection coefficients and so we now choose β = 10.
To compare the reflected wave amplitude in the FE-only case and the PML+IE case, the first velocity component p along a horizontal line in the centre of the test volume at a fixed point in time (chosen as the point immediately after reflection from the end of the open end of the volume occurs) is plotted in Figure 7. It can be seen that the reflected wave has a greater amplitude for the case where stress-free boundary conditions are implemented, demonstrating that the PML+IE formulation (with α = 1.001 and β = 10) is successful in reducing the reflection from the boundary. For the same set of parameters, the effect of the domain length in the x 1 direction on the reflection coefficient was studied and it was shown that to achieve a similar reduction in reflected amplitude in the case where stress-free boundary conditions are applied at x 1 = L 1 , the domain had to be more than doubled. Note that the time taken to run the simulations scales approximately linearly with the  length of the domain in the x 1 direction, as is expected. Importantly, there is little to no difference in runtime between the FE-only formulation and the PML+IE formulation and so the PML+IE formulation can produce a reflection coefficient equal to that of the FE-only formulation by using less than half the number of nodes (memory) and requires only half the runtime.

Spatially dependent stretching
In order to evaluate the integrals present in the formulation of the time-domain IE approximation in the case of non-constant stretching, values m = 1 and n = 1/4 are used in (48) and (49) [43]. Therefore, there is only one degree of freedom to explore, namely the stretching function parameterᾱ. The effect ofᾱ on the reflection coefficient has been assessed by exploring values close to the singularity atᾱ = 1 without breaching the conditions I * 2 ρL 1 c/3μ and I * 2 4L 1 μ/3c(λ + μ) (see [43]). It transpires that values very close to one in fact produce larger reflection coefficients and soᾱ = 2 is chosen for the remaining parametric studies. Figure 8 shows the first velocity component p recorded along a horizontal line in the centre of the test volume at a fixed point in time in the case whereᾱ = 2. Similar to the results presented in the constant stretching case, it can be seen that the reflected wave has a greater amplitude for the FE-only case (around 1) than for the PML+IE case (close to 0.67), demonstrating that the PML+IE formulation is successful in reducing the reflection from the boundary by around 33%. Similar to the constant stretching case, withᾱ = 2, the PML+IE formulation can produce a reflection coefficient equal to that of the FE-only formulation by using less than half the number of nodes and taking around the half the time to run the simulation [43].

Conclusion
An IE approach has been successfully combined with a PML to produce a new boundary condition for unbounded wave problems in the time domain and a novel formulation for an explicit FE approach using this new PML+IE combination has been presented. Two approaches were considered: the first where the PML stretching function has constant coefficients, and the second where a spatial dependency is retained. Although it would be useful to make comparisons between the PML+IE formulation and IE-only or PML-only formulations, it is not possible to simply switch off one of these components in this formulation. Instead, we compare the reflected wave amplitude with that observed in the case that standard stress-free boundary conditions are applied at the open end of the spatial domain. It was found that in both the constant stretching and spatially dependent stretching cases, the new combined PML+IE is successful in reducing the reflection coefficient by approximately 33%. It was observed that to achieve the same reduction in the reflected wave amplitude using the FE-only approach (with stress-free boundary conditions implemented at the open end) would require double the memory and twice the CPU time. Note here that both formulations produced similar results for the example studied. However, it may be in other scenarios that one approach outperforms the other and this would be an interesting avenue for future work.
Of course, in an ideal world, the new PML+IE formulation presented here would not only reduce the amplitude of the wave reflected from the boundary, but eradicate it completely. In this work, our inability to remove the reflected waveform entirely can in part be attributed to the approximations made in the diagonalisation of the velocity coefficient matrix, which is necessary to compute an efficient explicit FE scheme. However, now that the authors have developed a framework which combines PML and IE boundary conditions, it would be interesting future work to examine the results when an implicit scheme is implemented and quantitatively compare the performance of the method with that of other ABCs. set stretching parameterᾱ if I * 2 >> 4L 1 μ 3c(λ+2μ) or I * 2 >> ρL 1 c 3μ then STOP end if end if Initialise velocities, nodal weightings and stresses at all nodes Define the mass of a finite element as in equation (28) if constant stretching then Define the mass of an infinite element as in equation (28)  where N1 is the number of nodes in the x 1 direction Assign weights to all the nodes of the finite elements using the local coordinate index as shown in Figure 3 end for Assign weights to all the nodes of the infinite elements using the local coordinate index as shown in Figure 4 end for end for for timesteps from 1 to 1400 do Increase the time by δt for I = 1, . . . , NNOD do where NNOD is the total number of nodes in the spatial domain if constant stretching then Update velocities at node I using equation (27)) else if spatially dependent stretching then Update velocities at node I using equation (51)) end if for I = 1, . . . , NNOD do Set stresses to zero end for end for Apply an initial sinusoidal wave in the x 1 direction only for the first 35 timesteps for K = 1, . . . , N3 − 1 do for J = 1, . . . , N2 − 1 do for I = 2, . . . , N1 do Set elastic constants to be equal throughout the nodes this allows scope for heterogeneous materials end for Convert global node numbering to a local node index for the finite elements for I = 1, . . . , N1 − 1 do Calculate stresses on finite elements using equations (29)-(34) end for Convert global node numbering to a local node index for the infinite elements if constant stretching then Calculate stresses on infinite elements using equations (35)-(39) else if spatially dependent stretching then Calculate stresses on infinite elements using equations (52)-(56) end if end for end for Find a strand in the middle of the spatial domain and find the maximum velocity in the x 1 direction and the position at which this occurs on the strand Find the reflection coefficient at a specified node in the middle of the spatial domain and write to a file end for Find the runtime and write this to a file END