A time-domain finite element boundary integral approach for elastic wave scattering

The response of complex scatterers, such as rough or branched cracks, to incident elastic waves is required in many areas of industrial importance such as those in non-destructive evaluation and related fields; we develop an approach to generate accurate and rapid simulations. To achieve this we develop, in the time domain, an implementation to efficiently couple the finite element (FE) method within a small local region, and the boundary integral (BI) globally. The FE explicit scheme is run in a local box to compute the surface displacement of the scatterer, by giving forcing signals to excitation nodes, which can lie on the scatterer itself. The required input forces on the excitation nodes are obtained with a reformulated FE equation, according to the incident displacement field. The surface displacements computed by the local FE are then projected, through time-domain BI formulae, to calculate the scattering signals with different modes. This new method yields huge improvements in the efficiency of FE simulations for scattering from complex scatterers. We present results using different shapes and boundary conditions, all simulated using this approach in both 2D and 3D, and then compare with full FE models and theoretical solutions to demonstrate the efficiency and accuracy of this numerical approach.


Introduction
There is considerable scientific interest in both 2D and 3D elastic wave interaction with complex scatterers involving several disciplines. For instance, in seismology, measurements and simulations are used to investigate the reflected seismic wavefield from random rough surfaces for time-lapse monitoring techniques [19,33] and for inversion purposes [21]. In Non-destructive evaluation (NDE) and ultrasonics, the scattering behavior from rough defects can increase the measurement uncertainty and possibly alter the inspection results [25,45] with consequent implications on safety criteria. Motivated by these applications, amongst others, we develop an efficient, and general, tool to model elastic wave scattering from complex scatterers.
Naturally, there are pre-existing approaches and to compute the elastic wave scattering in the above applications many different modelling procedures have been developed. Analytical formulae are, of course, attractive and capable of handling a variety of scattering problems but only within certain restricted ranges of geometry or frequency. Classical approaches such as separation of variables (SEP) [43] provide exact mathematical formulae to calculate scattered waves but only from scatterers with regular and simple geometries, such as a spherical voids, and cannot handle complex branched cracks or irregular shapes. In terms of frequency one can employ methods such as the Kirchhoff approximation (KA) which is based on the assumption of an infinite tangential plane can predict scattering signals from defects with irregular shapes, but these approaches can generate unacceptable errors when the roughness is large or with grazing incidence/scattering angles [36,37].
Given the restrictions on analytical results there is natural interest in numerical approaches such as finite difference (FD) [7], finite element (FE) [4,8,30] and boundary element (BE) methods [13,14,17,41] that can be implemented in situations where the analytical solutions or the approximationbased methods are not reliable, but with a penalty in terms of computational cost; highly irregular or rough surfaces require accurate and highly resolved meshing and consequent increases in the numbers of unknowns. Comparison between the various methods is not straightforward, Frehner et al. [8] shows that FE is considered to be more accurate and robust than FD when modelling elastic wave scattering from defects with complex shapes. Compared with meshing the whole domain using FE, BE only needs to mesh the boundary of the scatterer and hence gives a smaller matrix system. However, the nature of FE results in a sparse matrix to be solved while the matrix from BE is not sparse. To further complicate comparison, and the use of specific methods, there are demands from industry for the use of standard FE packages that satisfy agreed industy standards. As a result, commercial FE software such as Abaqus (Dassault Systemes Simulia Corp., Providence, RI) and Ansys (Ansys Inc., PA) have been widely applied in industry and are standard numerical tools. There are also other approaches that take advantage of software advances, for instance, researchers at Imperial College have developed a GPU based time domain FE solver Pogo [11] to significantly accelerate FE simulations.
Despite these advances, and particularly in 3D and in unbounded domains, computational costs can become prohibitive for FE. A powerful approach is to use a so-called hybrid method [9,12,24,26,31,38,42] to reduce the computation effort. A hybrid method [26] in general involves a fast calculation of the global wave fields for a simplified material, and a more involved computation of a small local volume enclosing the irregular defects, complex material, or topography thereby concentrating the computational effort where it is most required. Amongst all the various different hybrid methods available [3,5,18,42,44], the FE-BI method has become of recent interest for large unbounded domains [15,16,24,35,38,42]; here one first numerically computes the local scattering in a small FE domain with the aid of an absorbing region [2,28,29] to prevent undesirable reflections due to the finite domain, and then calculate the scattered waves via a boundary integral in a sequential order.
The FE-BI method can be implemented either in the frequency domain or in the time domain depending upon the application. In the frequency domain this has been extensively studied and implemented for acoustic wave scattering problems (see, for example, the book by Ihlenburg [12]), with somewhat less work on elastic wave scattering in solids. Solid media support both compressional (P) and shear (S) modes, and also surface/interfacial waves such as Rayleigh and Stoneley waves. An elastodynamic finite element local scattering model has been developed to compute the scattering matrix, particularly on the ultrasonic array applications in the far field [38,42], and similar concepts have been applied  to simulate the elastic wave scattering of defects confined in waveguides [22,24,35]. However, much less work has been done on the coupled FE-BI for elastic waves in the time domain: A time-domain solution would provide simulated waveforms convenient for data processing in many applications, such as ultrasonic imaging [45] and seismic full wave form inversion [40]. In addition, the time domain FE scheme does not need to invert any global matrix because the global mass and damping matrices can both be diagonalized [11].
In terms of hybrid approaches [9] implemented it to compute the transient scattering signal using Auld's reciprocity principle and required a bespoke code because the FE calculation is based on a fictitious domain method with a new family of mixed finite elements [2]. Shen and Giurgiutiu [35] develop a hybrid FE method in the time domain with a specific application for guided wave structural health monitoring and more generally a hybrid method [31] has been developed with implementation for the time domain FE explicit scheme to solve more generic scattering problems; this approach provides an analytical hybrid solution to link transducer and defect responses computed from a numerical method. The hybrid method as described in, say, Rajagopal et al. [31] is shown in Fig. 1  potential is applied to calculate the incident wave propagating to the defect box. In the second step, the interaction between the incident wave and the scatterer is computed again by the time-domain FE inside the defect box. The excitation of the FE model is realized by applying forces to an excitation line (2D) several elements away from the defect. In the third step, the scattering field is collected by a monitoring box around the scatterer, and the hybrid scheme is used to calculate the receiving signals at the transducer.
Here we describe a robust and efficient time domain approach based on the hybrid concept, but we make specific modifications that incorporate the boundary integral approach showing that this provides considerable advantages over the previous hybrid method such as the work in [31]. Specifically, we use a reformulated FE equation to calculate the required forcing signals of the excitation nodes. The excitation and the monitoring nodes for the local FE box can both be on the surface of the scatterer that can have an arbitrary shape; hence the size of the FE model is reduced, and simultaneously the computational effort for the boundary integral is also minimized. In addition, the boundary integral is represented as a superposition of retarded time traces contributed from boundary values; the algorithm is completely implemented in the time domain and avoids any loop using Fourier transforms from moving between time domain and frequency formulations as in [31]. Finally, the new approach is flexible and capable of modelling more complex scenarios, including scattering from rough surfaces in a half space, and from scatterers with different boundary conditions. We show through several numerical examples that the new method achieves very good accuracy, and also high efficiency vis-a-vis the full FE model or the analytical solution.
This paper is organized as follows: Sect. 2 introduces the methodology, which includes the theoretical formulae and numerical implementations. Section 3 presents numerical examples in 2D and 3D, and concluding remarks are made in Sect. 4.

Methodology
Following the basic steps of the hybrid concept, as shown in Fig. 1, the modelling procedure we use is shown schematically in Fig. 2. Firstly, the transmitter is modelled to obtain the incident wave displacement field, for example, by the Rayleigh integral [10]. A time-domain FE reformulation is then used to calculate the required forcing signals applied on the excitation nodes, which can be on the scatterer surface. Only the elements attached at the excitation nodes are required in the second step; by using the forces obtained from the previous step as an input, a standard 2nd-order FE explicit scheme is implemented to compute the scattering displacement inside a small box. An absorbing region with a thickness of around one wavelength is added to eliminate unwanted reflections from the boundary. In the last step, the displacement signals on the scatterer surface are recorded, and substituted into the time-domain boundary integral to calculate the scattering signals at the receiver.

Finite element formulation
The 2nd-order elastodynamic FE equation for one element is [23]: where u e is the displacement vector, M e is the mass matrix, K e is the stiffness matrix and C e is the damping matrix. b e is the nodal force caused by the boundary traction and f e is the nodal force caused by the applied body force. (3) In NDE, the scatterer is normally a crack with a stress-free boundary condition, implying that the assembly of the nodal forces from the boundary tractions is zero. By adapting Eq.
(1) and assembling the required matrices, the following two equations can be obtained: In Eqs. (2) and (3), M att and K att are the local mass and stiffness matrices assembled only from elements attached at the excitation nodes. M all , K all and C all are corresponding matrices for the whole FE box including the absorbing region, and f ex (t) are forces at excitation/boundary nodes. M att and K att are subsets of M all and K all and the partitions are shown in Fig. 3. By substituting the displacement field of the incident wave u in (t) at the nodes of the attached elements into Eq. (2), one can calculate the forces at the excitation nodes f ex (t) using the central difference in the time domain.
where u in (t −δt), u in (t) and u in (t +δt) are the incident wave displacement vectors at the previous, current and next time steps, respectively. We choose to apply forces to excite the FE box rather than simply imposing the incident displacement. The reason is that we need to record the superposition of the incident and the scattering field on the boundary of the scatterer for the boundary integral. Applying the boundary displacement would result in only recording the incident field.
If the boundary condition is not stress-free, the boundary traction is not zero and hence the force f ex (t) needs to be subtracted by an additional term b ex (t): (5) in which B is the strain-displacement matrix, E is the matrix containing the elastic constants [23], n is the normal vector at the boundary, s is the shape function and Γ is the boundary surface for integration. The shape function varies per element type and it can be found in a standard textbook and references [23,32]. For a linear triangular element the shape function is , where and η are the local coordinates for a master element.
In the second step, f ex (t) obtained from Eq. (2) are taken into Eq. (3) as an input to calculate the displacement field u(t + δt) [11] in the local FE box: This equation is a standard 2nd-order explicit scheme to solve elastic wave problems as implemented in available software packages, such as Pogo and Abaqus. The acceptable thickness of the absorbing layers is approximately 1λ max using the stiffness reduction matrix (SRM) technique [28]. Here λ max denotes the maximum wavelength, and it normally equals to the compressional wavelength λ p . Assuming that the maximum dimension of the scatterer is L s , then in general the dimension of the FE box is just L s + 2λ max . Such a small computational region makes the implementation of FE in 3D feasible. In the frequency domain, it is possible to impose a non-reflecting boundary technique [39] and hence no absorbing region is required. However, currently there is no such solution in the time domain for the 2nd-order FE explicit scheme. In summary, the force f ex (t) is first calculated according to the known incident wave field using Eq.
(2), and these forces are then used in Eq.

Boundary integral
Once the surface displacements are obtained from the local FE simulation, a time-domain boundary integral can be applied to compute the scattered waves according to the Huygens' principle. The general boundary integral formulae in the time domain is [1]: where u i (r, τ ) and σ i j (r, τ ) represent the displacement and stresses at a point r at a time τ on the surface of the defect. R is the vector indicating the position of the observation point, and n j is the unit normal vector pointing outside the surface. u G i;k (|R − r|, t − τ ) and σ G i j;k (|R − r|, t − τ ) are the time dependent elastodynamic Green's function in the unbounded domain. The vector notations related with the boundary integral are depicted in Fig. 4.
It is desirable in engineering to use simplified expressions for the numerical implementation of the time-domain boundary integral. Here we start alternatively from the frequency domain and then transform the expression back to the time domain. In the frequency domain, the general boundary integral is: where the expressions of the Green's function σ G i j;k (|R − r|, ω) and u G i;k (|R − r|, ω) can be found in [1]. To analytically manipulate the stress Green's function [25], the far field condition that k α |R − r| 1 (α = p, s) is assumed, and terms involving 1/|R − r| 2 are neglected. Hence Eq. (8) in a 3D isotropic medium can be expressed as: where T α (t − D/c α ) and V α (t − D/c α ) represent travelling waves which are called surface retarded tractions and velocities.
In reality, the displacement U is the master variable which can be computed directly from the FE as shown in Eq. (6). The velocity V (t − D/c α ) is therefore calculated as the first order derivative of U(t − D/c α ). The traction T (t − D/c α ) on the other hand needs to be estimated by the following steps: (1) Locate the elements which share the same boundary node.
(2) Recover the stresses at the boundary elements using Σ = EBu(t). (3) Take the average value of the recovered stresses of these elements to estimate the stress at the boundary node, and then calculate the corresponding traction.
Step (3) is critical when a linear triangular/tetrahedral element is implemented to mesh irregular geometries. This is because the stress value is constant across such an element, and the averaging procedure can improve the accuracy of the recovered stress at the node.
The boundary of the defect can be discretized into small triangular facets in 3D as shown in Fig. 4 to perform the integration numerically: where D m = R − r m , and r m is the centre of the mth facet.
The boundary values are delayed in the time domain for integration, and they are linearly interpolated from the output of the FE computation.
Since |R| |Δr|, the vector between the observing point and the point at the crack surface can be approximated as: Note that in Eq. (13) we expand D with the first-order approximation relative to the centre of each facet rather than the usual way with respect to the origin of the vector r)(e.g. D ≈ R −R· r). Hence the boundary integral is more accurate when the distance R is not large enough for the requirement of the far field. Next a general form F(t − D/c α ) is used to represent either T (t − D/c α ) or V (t − D/c α ). By replacing D with the approximation from Eq. (13), F(t − D/c α ) can be expanded using a Taylor series: The integration of F(t − D/c α ) across one facet is therefore: where I m refers to the area of the mth facet. In Eq. (15) the first term assumes that the boundary values are constant across this facet, which can be approximated as being at the centre point. In the numerical implementation the scattered wave at the centre is calculated as an average of those at the facet nodes, so this first term is equivalent to applying a trapezoidal rule to the integration. Higher order Taylor expansion terms which refer to the derivatives of the retarded time traces account for the time variation across the facet. However, it will be shown through several numerical examples that integration based on simply keeping the first term is sufficiently accurate with the mesh density that has been tested.
In 2D the boundary integration formula has a slightly different form from that in 3D because of a different Green's function, and for 2D a scaling factor 2iπ k α needs to be included in the frequency domain formulation: Transforming the above equation to the time domain yields a convolution between the scaling factor and the boundary formula. tapered plane wave in the following simulations. Note that the proposed method can also model the scattering with an arbitrary incident wave field produced from a real transducer.

2D examples
The first example in 2D is the plane wave scattering from a side-drilled hole (SDH) with a diameter of 4 mm (≈2.6λ p ) as shown in Fig. 5a. The corresponding full FE region has a dimension of 51 × 51 mm 2 (≈33λ p × 33λ p ). The SRM technique [28] is adopted here to reduce the thickness of the absorbing region to about 1.5 mm (≈1λ p ). Linear triangular elements (equivalent to CPE3 in Abaqus) are applied to mesh the domain automatically. The element size is λ p /60, sufficiently small for the convergence requirement of both scattered P and S waves [6]. An excitation line is placed 20 mm (≈13λ p ) above the SDH to produce a plane P wave propagating along the negative y-direction, and the scattering signals are recorded by a monitoring circle 20 mm away from the SDH. The full FE model is solved using the Abaqus explicit solver.
In contrast, the size of the FE box (8×8 mm 2 ≈ 5λ p ×5λ p ) used in the FE-BI method is much smaller by comparing Fig.  5a, b, and the excitation nodes are located at the SDH surface in the FE box denoted by red dots. The forcing applied at the defect surface is calculated using Eq. (4). After running the FE explicit scheme in the box, the boundary displacements are used to calculate the scattering signals via the boundary integral Eq. (11), at the monitoring circle. The first term in Eq.
(11) is neglected due to the stress-free boundary condition. The quantity of the boundary velocity in Eq. (11) must be the total field (i.e., incident field plus scattering field). Figure 6a, c show good agreements between the full FE model and the FE box for the scattered P and S waves at a near grazing angle (θ s = 80 • ). The scattering amplitude, which here is defined as the ratio between the peaks between the scattering and the incident signals, with θ s ranging from 0 • to 360 • is shown in Fig. 6b, d, respectively. The mean absolute error (MAE) from all scattering angles is used to measure the deviation and it can be expressed as: where M is the number of the scattering angles, A(θ n s ) and A re f (θ n s ) are the scattering amplitude calculated using the FE-BI method and the full FE method. For the SDH case, the MAE of the scattering amplitude are calculated as 1.5 and 3.2% for P and S waves, respectively.
The element size used in the full FE model is λ p /60, which is sufficiently small as it easily satisfies established convergence criteria [6]. On the contrary, in the FE-BI box the element size is λ p /30, much larger than that used in the full FE simulation. It implies a somewhat more relaxed convergence requirement using the FE-BI method than the full FE model. The reason is that the local FE simulation mitigates the error caused by the mesh dispersion, when modelling the wave propagating from the source to the defect and vice-versa.
Giving the full FE model the same element size of λ p /30 would cause noticeable errors for the scattered S waves due to the mesh dispersion. The number of nodes in the local FE box is around 25,000, which is significantly smaller than that for the full FE model (4.1 million), a reduction of approximately 170 times. In addition, the total time for running the full FE model is around 32 min, while for the small FE box it only takes 20 s. A huge improvement of the computational efficiency is therefore achieved. The second example is the scattering of a plane P wave from a rough surface as shown in Fig. 7. To approximate the plane wave scattering from an infinitely long surface a Gaussian tapered plane wave in [34] is adopted here. A spatial Gaussian window is added to the plane wave so that the amplitude of the incident wave impinging on the rough surface gradually reduces to zero at the ends, in order to avoid the edge effects when performing the boundary integral along the surface with a finite length. The total length of the surface is 24 mm (≈15λ p ), based on which the half beam width is approximately 4 mm (≈2.6λ p ). One realization of the Gaussian surface with RMS = λ p /4 and correlation length λ 0 = λ p /2 is generated using the spectral method [37]. The dimension of the full FE model shown in Fig. 7a is 51 × 25.5 mm 2 (≈33λ p × 16λ p ). An excitation line with the same length of the surface is placed 4 mm (≈2.6λ p ) above the rough surface. The incident wavefront distorts when the ideal Gaussian plane wave propagates towards the surface, and hence in the full FE model we put the excitation line as close as possible to the rough surface. The scattered waves are recorded by a semi-circle of the nodes 20 mm away from the centre of the surface. Some redundant waves are noticed propagating to the positive y-direction. These waves are generated by the excitation line due to the vibration of the source nodes in a free space.
For comparison, Fig. 7b shows the corresponding FE box which is 29 × 4 mm 2 (≈18.7λ p × 2.6λ p ). Specifically, the excitation is realized by the nodes on the rough surface with applied forces calculated from Eq. (4). Computed surface displacements from the local FE box are substituted into the boundary integral Eq. (11) to obtain the scattering signals at the same monitoring points in the full FE model. For the rough surface, the scattered waveforms become much more complicated especially for the mode converted S waves, due to the increased diffuse field as shown in Fig. 8c. However, good agreements of the scattered P and S waveforms calculated from the full FE model and the FE-BI model are still found. In addition, the FE-BI method also accurately predicts the scattering amplitude across all the scattering angles, and the values of the MAE are 3.8 and 5.1% for P and S waves, respectively. The number of elements for the full FE model and the local FE box in Fig. 7a, b are 2.1 million and 47,000, and the computational effort is greatly reduced by the new method.

3D examples
In 3D the plane wave scattering from a spherical void and an inclusion with the same shape are simulated. The com-puted results from the proposed method are then compared with those calculated from the theoretical formulae [43]. It is currently very difficult to simulate a full 3D FE model due to the limit of computational power. The 3D FE box has dimensions of 6.2 × 6.2 × 7.2 mm 3 (≈4λ p × 4.7λ p × 4λ p )as shown in Fig. 9, with the absorbing region being 1.5 mm (≈1λ p ) thick. The spherical void with a diameter of 1.2 mm (≈0.77λ p ) is created by subtracting a 3D solid sphere from the FE box. Linear tetrahedral elements Assuming a plane wave propagating in the negative zdirection, the forcing signals excited on the defect surface can be calculated using Eq. (4). After executing the FE explicit scheme in the local box, the calculated surface displacements are substituted into Eq. (11) to calculate the scattering signals from 0 • to 360 • . The distance between the observing point and the centre of the defect is approximately 100λ p . Figure  9a shows the mesh around the void, and the local scattering wavefield is shown in Fig. 9b. Figure 9c shows a good agreement for the scattered P wave signals when θ s = 30 • , using the FE-BI method and the theoretical solution [43]. From Fig. 9d, the MAE of the scattering amplitude between the two approaches is 2.1%. A small deviation of 4.2% can be seen when θ s is around 180 • , corresponding to the transmission direction.
The same geometry and bulk material (Aluminium) are used for the simulation of the scattering from a spherical inclusion, filled with Alumina (Young's modulus, 390 GPa; density, 3950 kg/m 3 ; and Poisson ratio, 0.22). The impedance ratio for the compressional wave between the bulk medium and the inclusion is hence 0.40. In practice, the spherical inclusion and the cubic box without the sphere are meshed separately, and then joined to form the whole mesh shown in Fig. 10a. Note that the two separate volumetric meshes should produce the same boundary mesh at the interface to make the volumetric mesh compatible.
We do not apply forces on the surface of the inclusion directly in the FE box here because it is found to produce unwanted waves which will pollute the scattering field. For example, exciting the lower half spherical surface would give waves travelling in both positive and negative z-directions, which is caused by the vibration of source nodes in the free space. The waves travelling along the positive z-direction is unwanted and will contaminate the recorded waves on the inclusion surface. Hence here an excitation plane 1 mm above the inclusion is applied instead. The forces are given to the plane and generate the required incident field travelling downward along the negative z-direction. Figure 10b shows the scattering field from the inclusion, and the surface nodes to record boundary displacements. Elements associated with these nodes are also selected to post-process the tractions, and an average is performed to estimate the tractions at the surface of the inclusion. After taking the boundary velocity and traction into Eq. (11), the scattering signals at the observing points can be calculated. A good match between the theory and the simulation can be found in Fig. 10c, d. A small error of the amplitude at the transmission direction (θ s = 0 • ) might be caused by the fact that the linear tetrahedral element is not sufficiently accurate to estimate the boundary stress. A common way to improve the accuracy is to increase the mesh density, especially around the surface of the inclusion. Several alternative ways may be useful to recover the stresses more accurately based on the current mesh density [20]. For example, the nodal point forces (NPF) method as developed in [27] can be potentially implemented here.

Conclusions
We have presented an accurate and high-performance numerical method to solve transient elastic wave scattering prob-lems, in particular for 2D or 3D complex scatterers. It couples the FE and the boundary integral and is completely implemented in the time domain. The time-domain FE equation is reformulated to calculate the required forcing signals applied on the excitation nodes according to arbitrary incident waves. These excitation nodes can be the same as the boundary nodes of the scatterer. Then a small local FE box is executed using the standard explicit scheme with the computed forces as an excitation. The boundary displacement of the scatterer is obtained and integrated in the time domain to calculate the scattering signals for different modes.
The method proposed in this article significantly improves the efficiency of time-domain FE simulations, enabling a direct comparison of the waveforms at different scattering angles obtained numerically and experimentally. The incident wave can be in principle a wide band signal without loss of efficiency. In addition, it can potentially be applied to many research areas involving the use of the scattering signals in solids, such as full waveform inversion, ultrasonic tomography, TOFD for sizing an internal defect, and a variety of elastic wave imaging algorithms.