Time domain numerical modelling of guided wave excitation in fluid-filled pipes

Acoustic waves have been widely used to inspect pipeline defects, including leakage, blockage and corrosion etc. The time-domain excitation of guided waves in the fluid/solid coupled pipeline system has rarely been studied theoretically. This requires the incorporation of a source term in the coupled system. This article introduces a finite element based numerical model to study the excitation of guided waves in the coupled system with a sound source either in the fluid or on the pipe wall. A compatible pair of time-domain Perfectly Matched Layers (PMLs) have been proposed to absorb pipe end reflections in the elastic wall and in the fluid respectively. These fluid/ solid PML numerical formulations are coupled and the implementation of the coupling is introduced. The numerical model is validated against an enlarged model without PML, and excellent agreement has been achieved. The numerical model shows that the dominant excited wave mode is a fluid type wave mode when the sound source is a fluid source or a radial elastic line source. However, the dominant excited wave mode is a structure type wave mode when the sound source is an axial elastic line source.


Introduction
Acoustic waves are widely used to inspect pipeline defects in nondestructive testing and condition monitoring applications. This can be based on a pulse-echo principle [1], an acoustic correlation approach [2] or a pulse reflectometry technique [3]. In pulse-echo and pulse reflectometry applications, the incident sound wave is normally excited by a transient sound source and time domain analysis is important to facilitate direct comparison between scattered signals from theoretical predictions and experimental measurements [1,3]. The time-domain excitation of guided waves in fluid-filled pipes is rarely investigated theoretically in the literature. In this article, a finite element based numerical technique is proposed to study the excitation of guided waves in fluid-filled pipes. Two consistent time-domain Perfectly Matched Layers (PMLs) have been proposed to close the modelling domain in the elastic pipe wall and in the fluid respectively. The elastic PML incorporates a mixed displacement-strain un-split formulation, whereas the fluid PML incorporates a mixed velocity potential-velocity un-split formulation. The two PMLs are structured in a way that facilitates the implementation of fluid/solid coupling in the absorbing regions. The numerical technique is validated by comparing PML solutions with solutions obtained from an enlarged model without PML. On this basis, the exited wave modes from a sound source either in the fluid or on the pipe wall are calculated and compared.
Analytical approaches have been proposed to study wave propagation in fluid-filled pipes. Fuller and Fahy investigated dispersion and energy distribution properties of guided waves in a fluid-filled pipe [4]. A thin-shell theory has been used which allows the characteristic equation to be formed based on the expansion of the determinant of a third order matrix. Sinha et al. studied axisymmetric wave propagation along fluid-loaded pipes [5,6]. The wave equations are in the framework of the linear theory of elasticity, and the displacements are expressed in terms of the Bessel functions with unknown coefficients. The classical perfect-slip boundary conditions have been used at the solid-fluid interface, and modal solutions are obtained when the determinant of a six order matrix vanishes. Lowe [7] developed a general method for modal solutions of multi-layered elastic waveguides using the global matrix technique. This delivers Lamb wave solutions in elastic isotropic waveguides. In addition, the global matrix technique is extended to include cylindrical layers [8]. These analytical approaches [4][5][6][7][8] require an iterative numerical root-search procedure to be used to find the roots of the characteristic equations, which limits the computational efficiency of these approaches.
A few approximate analytical solutions have also been proposed which avoids the numerical root-searching procedure. Muggleton et al. proposed explicit expressions for wavenumber solutions of two fundamental axisymmetric wave modes in a fluid-filled buried pipe [9,10]. The slip bonding condition has been used on the pipe/soil interface. Gao et al. improved the wavenumber solution for the fundamental fluid dominant mode in the fluid filled buried pipe [11]. The perfect bonding condition between the pipe and the soil is adopted, which is considered to be a more appropriate boundary condition between two elastic domains. These approximate wavenumber solutions [9][10][11] are valid below the ring frequency of the pipe because of thin shell approximations.
The semi-analytical finite element (SAFE) method is a highly efficient numerical approach to calculate modal solutions of fluid-filled pipes. The SAFE method adopts a variable separation technique, and an analytical plane wave solution is used in the axial direction of the fluid-filled pipe, so that it is only required to mesh the cross-section of the fluid-filled pipe. Kalkowski et al. proposed a variational formulation to examine fluid-filled pipes that can be immersed in fluid or buried in soil [12]. Duan and Kirby introduced a weighted residual formulation based on the Galerkin's principle to study the coupled problem [13]. The conventional slip boundary between the fluid and the pipe has been used to combine equations in the form of a standard eigenvalue problem [12,13]. The frequency domain PML technique is used to close the computational domain for buried or immersed pipe, and the exponential stretching function proposed by Duan et al. [14] is seen to work particularly well for these modal solutions.
The SAFE technique is only valid for a waveguide that is uniform in the axial direction. This method cannot be directly used to study excitation or scattering problems. Williams et al. investigated scattering by a single flange in a fluid-filled pipe [15]. The waveguide is divided into three uniform sections, and modal expansions are used in each section of the waveguide. Point collocation and mode matching techniques are used to link modal solutions on the interfaces between these regions. The incident sound wave is assumed to be a single mode with known coefficient, and the reflection and transmission coefficients are generated.
The above SAFE, point collocation and mode matching techniques are used in the frequency domain. The pipeline system is assumed to be infinitely long in these models, which implies that the length of the pipe is much longer than the wavelength in practice. To study the excitation of guided waves in infinitely long fluid-filled pipes in the time domain, the modelling domain has to be closed at both ends of the pipeline system. A time domain PML method can be used to absorb outgoing waves near the ends of the pipe. Although the frequency domain (time harmonic) PML technique has been used in SAFE applications [12][13][14], this cannot be readily extended to the time domain. In frequency domain models the systematic matrices are formulated in the complex space. Thus the complex PML stretching functions and coordinates can be directly incorporated into these equations without further modification.
However, the systematic matrices are formulated in the real space in standard time marching models. As a result, the complex PML stretching functions and coordinates cannot be directly included and executed in real space models. A common technique to address this is to derive the PML governing equations in the frequency domain, and convert the frequency domain equations to the time domain through an inverse Fourier transform. This conversion can deliver a set of real space equations at the cost of introducing auxiliary fields or convolution operations in the PML region. Komatitsch and Tromp proposed a split-field PML formulation where the gradient operator is split in terms of components perpendicular and parallel to the interface between the normal and PML region, and the displacement is further split into four parts [16]. Komatitsch and Martin later used the so named convolutional or complex frequency-shifted technique to improve the efficiency of the PML at grazing incidence [17]. The equations are presented in the first-order velocity-stress mixed formulation, which cannot be readily used in numerical schemes which are based on the second-order wave equation involving only displacement. Furthermore, the complex frequency-shifted PML technique suffers from degrading absorption of propagating waves at low frequencies.
Kucukcoban and Kallivokas proposed a displacement-stress mixed two dimensional formulation that is compatible with standard displacement only numerical codes for the regular unstretched region, and can also be executed with the Newmark time marching scheme [18,19]. The wave equation is presented in the displacement-stress mixed form, and a separate kinematic equation is proposed which links the stretched strain with displacement. The two governing equations can be solved with the support of auxiliary variables. This formulation retains the advantage of unsplit fields and requires no convolutional or complicated time integration schemes. Fathi et al. extended the two-dimensional formulation to three dimensions, and discussed several time integration schemes [20]. Pled and Desceliers reviewed the latest developments on the PML method for modelling elastic wave propagation in unbounded domains [21]. An alternative displacement-strain mixed formulation is proposed which is close to the displacement-stress mixed formulation [18][19][20], and an extended Newmark time marching scheme is introduced. It should be noted that the coordinate stretching function used in Ref. [21] is a simplified version of the standard stretching function, with the real part fixed to be one. This could reduce the efficiency of the PML technique.
The time-domain PML has also been used in the fluid domain. Gao proposed an unsplit complex frequency-shifted technique to study twodimensional wave propagation in the fluid [22]. Three non-physical auxiliary variables are introduced, and the final time domain equations are first-order and second order equations which can be executed without convolution operations. Assi and Cobbold studied transient wave propagation in a two-dimensional unbounded fluid-solid medium [23]. The fluid domain equations are presented in the pressure-velocity mixed formulation, and non-physical auxiliary variables are introduced in the solid domain. However, the numerical implementation of the governing equations is not introduced and the equations are executed in a commercial FEM package COMSOL.
In this article, a time domain fluid/solid coupled finite element model is proposed to study the excitation of guided waves in fluid-filled pipes. In the solid domain, the computational region is closed by an unsplit displacement-strain mixed PML formulation. This is similar to the work of Pled and Desceliers [21]. However, the formulation proposed here is based on the standard two-parameter stretching function, which is more efficient than the single-parameter stretching function based formulation especially when evanescent waves are present. In the fluid domain, a new unsplit velocity potential-velocity mixed formulation is proposed, which is compatible to the elastic domain PML. The fluid/solid coupling on the boundary interfaces are discussed in both the regular and PML regions. An extended Newmark time marching algorithm is used to execute the time domain model. The article is structured as follows: in section 2, the theoretical development of the weak form equations, the coupling and the time marching scheme are detailed. In section 3 the numerical model is validated, and predictions are given for excited sound waves from a sound source either in the fluid or around the circumference of the pipe. Conclusions are then drawn in section 4.

Theoretical development
The structure of the fluid-solid coupled system is shown in Fig. 1. The pipe is assumed to be infinitely long. PMLs are placed both in the fluid and in the elastic pipe wall to close the infinite domain. The fluid and solid PMLs both start from the same axial position and terminate at the same axial position. The sound source can be either in the fluid or on the  The governing equations are presented based on a general threedimensional model, however, the system studied in this article can be simplified as a two-dimensional axisymmetric model, and the simplification is detailed in this section. The traction free boundary condition has been applied on the exterior surface of the pipe.

Elastic domain
The equation of motion for elastic wave propagation in the pipe wall can be written as [24]: or in the matrix form Here, σ qj (q,j = 1,2 or 3) is a symmetric stress tensor of rank two, ρ is density, t is time, u q is the displacement in the q direction, and f q is a body force. The repeated indices indicate summation throughout this article. The cylindrical coordinate system is used here which facilitates the simplification of the model for axisymmetric problems. For computational reasons, it is practical to rewrite stress in the vector form, and the divergence of stress is given in the matrix form in Appendix A1.
To derive the weak form of Eq. (1), a Galerkin based method is used by taking the inner product of Eq. (1) with an arbitrary weighting function v q . The resulting product is then integrated over the elastic computational domain Ω s (excluding PMLs). Integration by parts and the divergence theorem gives: Here, ε qj is the strain tensor, and σ qj = σ jq has been used in the development of Eq. (2). The double dot over a variable denotes the second order partial derivative with respect to time. Γ s is the piecewise smooth boundary of the elastic domain Ω s , and n s = [n s1 n s2 n s3 ] is the unit outward normal vector to Γ s . The stress-strain constitutive relation is given as: The double dot product of the stiffness and strain tensor is defined by the summation over repeated subscripts k and l in Eq. (3). The matrix formulation of Eq. (3) is given in Appendix A2.
In the PML region the coordinate x j is replaced by a stretched coordinate x j , defined as where ξ j (x j ) is a non-zero, continuous and complex-valued coordinate stretching function defining the scaling and damping profiles in the x j direction. This definition requires that dx j = ξ j dx j , and The stretching function is usually chosen to be frequency dependent in the time domain: where, h j and d j are real valued scaling and damping functions, ω is the radian frequency and i = ̅̅̅̅̅̅ ̅ − 1 √ . The particular form of iω in Eq. (5) facilitates the easy implementation of the inverse Fourier transform from frequency domain equations to time domain equations. Note that the real part of the stretching function, h j , can be simplified and fixed to be one. However, using a value larger than one in the PML improves the attenuation of evanescent waves [18].
The elastic equation of motion in the transformed coordinate system can be written in the frequency domain as: The over-hat indicates the Fourier transform of the variable, and the external body force has been set to zero since it is in the PML. It is convenient to start the derivation from the frequency domain, because ξ j is a complex variable, and it will be shown here that a number of procedures can be used to ensure that when transformed from frequency domain to time domain the final systematic equation involves only real variables which can be executed in standard time marching programmes.
Multiply Eq. (6) with ζ, where ζ = ξ 1 ξ 2 ξ 3 . It is clear that ξ j is a function of x j only, and thus ζ ξ j is independent of x j . The governing equation can be rewritten as: Here, D R 3 , D ∈ 3 and D ⊓ 3 are given in Appendix A3-A5 respectively. The stress-strain constitutive relation in the stretched coordinate system is given as: The strain-displacement kinematic relation in the stretched coordinate system is given as [18,21]: where û is the displacement vector and the gradient of û in the cylindrical coordinate system is given in Appendix A6. The stretching tensor J is given in Appendix A7. It is possible to substitute Eqs. (8) and (9) into Eq. (7), and obtain a single equation with only displacements as unknowns. However, the resulting time-domain equation will be onerous and no clear procedure can be used to solve such an equation. In this article, the strain tensor ε is introduced as an auxiliary variable, and the final systematic equation will be in a displacement-strain mixed form.
Eq. (9) can be pre and post-multiplied by iωJ T and J, respectively. Rearranging and grouping similar terms, the following kinematic equation can be obtained [18,19]: where J = H + 1 iω D has been used in the development of Eq. (10). H and D are given in appendix A8-A9.
The weak form of Eq. (7) is obtained similarly by using the Galerkin's principle. The inner product with an arbitrary weighting function v q is taken. The resulting product is integrated over the elastic PML domain Ω sp . Finally, integration by parts and the divergence theorem gives: Here, Γ sp is the piecewise smooth boundary of the elastic PML domain Ω sp .
The kinematic Eq. (10) is weighted by an arbitrary weighting function w. The resulting product is integrated over the elastic PML domain Ω sp , and the weak form equation can be written as: Now Eq. (8) is substituted into Eq. (11), and the inverse Fourier transform can be applied to Eq. (11&12) to get time domain equations.
To facilitate the inverse Fourier transform to Eq. (11&12), we introduce additional auxiliary fields, such that: , Here, an over dot denotes differentiation with respect to time, and F − 1 denotes an inverse Fourier transform.
By applying the inverse Fourier transform to Eq. (11&12) and substituting Eq. (13) into the transformed equations, we have and ∫ Here, the equations can be simplified by assuming that no coordinate stretching is introduced in the circumferential direction of the pipe, so that h 2 = 1, and d 2 = 0. The sub-indexes 1, 2 and 3 corresponds to r, θ and z respectively. Furthermore, it is assumed that the sound source is placed either at the centre of the pipe or around the circumference of the pipe. The general three-dimensional Eq. (14&15) can then be simplified as an axisymmetric model, with only radial and axial displacements considered in this article.
The finite element method proceeds by discretising the displacement and strain using standard shape functions in the vector form, so that and where N u and N e are global trial (or shape) functions, u q and ε k are nodal values of displacement and strain expressed in the column vector form. In addition, N is a row vector of length m s , and m s is the number of nodes in the elastic region. Φ is a row vector of length m sp , and m sp is the number of nodes in the elastic PML region. N u is a 2 × 2m s matrix, and N e is a 4 × 4m sp matrix. Isoparametric elements are used throughout this article, which requires that the weighting functions are equal to the shape functions so that v q = [N u ] T , and w = [N e ] T . Eqs. 2 and (14) &15) can be combined in the matrix form as: Here The submatrices in these equations are given in Appendix A10-A29. Note that the coordinate stretching starts outwards from the interface between the elastic domain Ω s and the PML domain Ω sp , and thus D R 3 is an identity matrix, D ∈ 3 and D ㄇ 3 are zero matrices on the Ω s / Ω sp interface. This implies the surface integrations in Eq. (22) can be cancelled out on the interface between the two regions, because the unit outward normal vectors point towards each other.
Furthermore, a traction free boundary condition has been applied on the exterior surface of the pipe. The residual stresses are also assumed to be zero on the far end of the PML surface, because of damping in the PML region. Surface integrations on these surfaces can thus be avoided. However, the solid pipe is coupled to the fluid on the coupling interfaces Γ fsi and Γ fsp (shown in Fig. 1), and thus surface integrations must be carried out on the fluid/solid interfaces. Eq. (22) can be rewritten as: The surface integrations in Eq. (23) will be considered later in section 2.3.
All the matrices and vectors in Eq. (18) are real, and Eq. (18) can be executed in an extended Newmark time marching procedure, which will be detailed in section 2.4.

Fluid domain
The equation of motion for wave propagation in a stationary, uniform fluid can be written as [25]: where ϕ is the velocity potential, c is the speed of sound, and f 0 is a body force in the fluid.
Eq. (24) is multiplied by a weighting function v and integrated over the fluid domain Ω f (excluding PMLs). After applying Green's theorem, the weak form equation is given as: ∫ Here, Γ f is the piecewise smooth boundary of the fluid domain Ω f , and n f = [n f1 n f2 n f3 ] is the unit outward normal vector to Γ f . Now we move to the PML region, and the same coordinate stretching function defined in Eqs. (4) and (5) are used here. The wave equation in the transformed coordinate system can be written in the frequency domain as: The over-hat indicates the Fourier transform of the variable, and the external body force has been set to zero since it is in the PML region.
The following auxiliary vector field is now introduced: Here, ψ j is physically interpreted as the fluid particle velocity in the PML region.
The weak form of Eq. (28) is obtained by taking the inner product with an arbitrary weighting function v. The resulting product is integrated over the fluid PML domain Ω fp . Finally, integration by parts and the divergence theorem gives: Here, Γ fp is the piecewise smooth boundary of the fluid PML domain Ω fp .
Both sides of Eq. (27) are now multiplied by ξ j . The weak form of the constitutive equation is obtained by weighting Eq. (27) with an arbitrary function w and integrating over the fluid PML domain Ω fp : To facilitate the inverse Fourier transform to Eq. (29&30), we introduce additional auxiliary fields, such that: Clearly, and By applying the inverse Fourier transform to Eq. (29&30) and substituting Eq. (31) into the transformed equations, we have ∫ and ∫ The equations are now simplified by assuming that no coordinate stretching is introduced in the circumferential direction, so that h 2 = 1, and d 2 = 0. Furthermore, it is assumed that the fluid sound source is placed at the centre of the pipe. The general three-dimensional Eq. (32&33) can then be simplified as an axisymmetric model, with only radial and axial particle velocities considered in this article.
The finite element method proceeds by discretising the velocity potential and velocity using standard shape functions in the vector form, so that and where N ϕ and N ψ are global trial (or shape) functions, ϕ q and ψ k are nodal values of velocity potential and velocity expressed in the column vector form. In addition, N ϕ is a row vector of length m f , and m f is the number of nodes (or degrees of freedom) in the fluid region. Φ is a row vector of length m fp , and m fp is the number of nodes in the fluid PML region. Isoparametric elements are used throughout this article, which requires that the weighting functions are equal to the shape functions so Eqs. 25 and (32) &33) can be combined in the matrix form as: The submatrices in these equations are given in Appendix B1-B15. Similar to solid domain surface integrations given in Eq. (22&23), it is only necessary to carry out the fluid domain surface integrations on the solid/fluid interfaces Γ fsi and Γ fsp in Eq. (40), and the details are given in the next section 2.3. Furthermore, it can be seen in Appendix B15 that D ⊓ 2 simplifies to a zero matrix since d 2 = 0. H f is thus a zero matrix, and it is kept in Eq. (36) for ease of extension to a threedimensional model only.

Fluid/solid coupling
The coupling between the fluid and solid domains must be considered in surface integrations in Eq. (23&40). The surface integrations require boundary conditions to be specified on the fluid/solid interfaces.
The integrations on the coupling interface Γ fsi will be considered first. It is assumed that the normal displacements are continuous, i.e., Here, u f is the fluid particle displacement. The velocity potential is defined by: where v f is the fluid particle velocity.
Substitute Eq. (41&42) into Eq. (40). The surface integration on Γ fsi in Eq. (40) is given as: ∫ The normal solid traction force and fluid pressure are also continuous on the Γ fsi interface: The boundary condition Eq. (44) is also valid on the Γ fsp interface, and Eq. (45) can be written in the frequency domain as: It is clear that the force term F 0 in the fluid domain involves solid particle velocity ∂u1 ∂t and displacement u 1 , whereas the force term F s in the solid domain involves fluid potential derivative ∂ϕ ∂t and potential ϕ. It is thus a coupled problem with integrations shown on the coupling interfaces Γ fsi and Γ fsp . Note that these ordinary and first order derivative terms can be readily obtained in a Newmark time marching algorithm which will be shown in the next section.
To improve the condition of the matrices in the coupled system, the fluid potential is normalised by: Here, f c is the centre frequency of the sound source, and the bar over a variable denotes normalisation.
The final systematic Eq. (18) in the solid domain can be normalised as: where The final systematic Eq. (36) in the fluid domain can be normalised as: where F 0 is given in Eq. (51). Eq. (55&57) can be implemented using an incremental version Newmark time marching algorithm detailed in the next section. The coupled equations are calculated based on an iterative 'staggered scheme' in each time step [26]. In this article, the fluid domain Eq. (57) is solved first with known initial conditions. The calculated potential ϕ and potential derivative ∂ϕ ∂t are then substituted into Eq. (56), and the solid domain Eq. (55) can then be solved. This procedure is iterated again in each time step, and in the second iteration the updated fluid/solid domain solution will be used on the coupling interfaces for the calculation of F 0 and F s . The numerical solution will be converged after only two iterations in each time step [26].

Extended newmark time marching
Eq. (55) and Eq. (57) are identical in form, and thus can be solved using the same time marching algorithm. The wave equations are discretised in the time domain with a time interval Δt, defined as Δt = T/N. Here T is the total modelling duration, and N is a positive integer.
The general time domain discretised equation can be written as: Here, a n is denoted as discretised temporal approximation of U in Eq. (55), or P in Eq. (57) at time t n , where t n = nΔt for each n = 0, 1, 2,⋯N − 1. I n and J n are corresponding discretised auxiliary fields. The initial conditions are assumed to be static conditions, i.e., a 0 = 0, ȧ 0 = 0, ä 0 = 0, I 0 = 0, J 0 = 0 and f 0 = 0.
The auxiliary fields I n and J n are expressed as [23]: The Newmark algorithm is unconditionally stable when β 3 ≥ 1/2 [27].
Considering Eq. (60&61), the incremental quantities can be rewritten as: Eq. (63) is the first step of the time marching algorithm. All the incremental quantities in Eq. (62) can be updated by substituting the solution of Eq. (63) into Eq. (62). Now quantities at the time step t n+1 are evaluated by a n+1 = a n + Δa n , ȧ n+1 =ȧ n + Δȧ n and so on [23]. Eq. (63) thus permits the evaluation of quantities from the current time step t n to the next time step t n+1 , and the time marching procedure is now complete.

Numerical results and discussions
The time domain model detailed in section 2 can be executed with careful specification of an incident sound wave. It is preferable to choose a frequency range within which the incident wave mode is relatively non-dispersive so that defect location can be accurately estimated in non-destructive testing and condition monitoring applications [1][2][3]. A modal analysis shall be carried out to study dispersion properties of guided waves in fluid-filled pipes prior to the numerical excitation study, and the SAFE model of Duan and Kirby [13] is used here. The pipe is chosen to be a 3 inch steel pipe with an outer radius of 44.65 mm and an inner radius of 39 mm respectively [28]. The properties of the steel pipe are c L = 5960 m/s, c T = 3260 m/s, and ρ = 8030 kg/ m 3 . Here, c L and c T are the longitudinal and shear bulk wave velocity of the pipe respectively. The pipe is filled with water with properties given by c = 1480 m/s, and ρ 0 = 1000 kg/m 3 .

Guided wave modal properties
Phase velocity, energy velocity and energy distribution properties of guided waves in the water-filled pipe are shown in Fig. 2. The energy distribution (ratio) is defined as the sum of the sound power in the fluid divided by the sum of the sound power in the pipe. The labelling of wave modes follows the principle introduced by Duan and Kirby [13]. Here, FS(0,1) is a fluid type mode at low frequencies below 15 kHz, and energy transfers from the fluid to the pipe wall at high frequencies. SF(0,1) is a structural type mode at low frequencies below 15 kHz, and energy transfers from the pipe wall to the fluid at high frequencies. T (0,1) is torsional wave mode in the pipe wall that is uncoupled to the fluid, and thus all the energy of this mode is in the pipe wall. SF(0,2), SF(0,3) and SF(0,4) are cut-on wave modes with energy mainly in the pipe wall above but close to cut-on frequencies. The frequency dependent energy transferring properties of fluid/pipe coupled wave modes are also discussed in Ref. [13] for immersed and buried pipes. It is clear that some wave modes are particularly dispersive around frequencies such as 18 kHz, 26 kHz etc, and these frequencies shall be avoided in non-destructive testing and condition monitoring applications.

Convergence and validation
We will now examine the time domain model presented in section 2. The numerical model was programmed in MATLAB® and executed on a desktop with a 3.5 GHz AMD CPU and an 16 GB RAM. The incident sound wave is chosen to be a narrow frequency bandwidth 10 cycle Hanning windowed pulse given as: Here, δ is the Dirac delta function, and r 0 and z 0 are the radial and axial location of the sound source respectively. f c is the centre frequency of the incident pulse. t d is the sound source time delay and N c is the number of cycles. In this article, the time delay is fixed to be zero and the number of cycles is 10, i.e., t d = 0 and N c = 10. The force term defined in Eq. (64) is a scalar function if the sound source is in the fluid. However, the force term is a vector function if the sound source is on the pipe wall. In the latter case, the force is generally applied in the radial or axial direction. This can be done by substituting Eq.  (23). Note that the second-order Ricker wavelet (or Mexican hat wavelet) is often used as the incident function in seismic applications [16][17][18][19]. However, the frequency bandwidth of the Ricker wavelet is wider than the 10 cycle Hanning windowed pulse used here. To avoid wave dispersion, the Ricker wavelet is rarely used as the incident wave in pipeline inspection applications.
The complex exponential stretching function proposed by Duan et al. [14] is used here. The radial coordinate stretching function is used to study pipes buried in elastic media or immersed in fluid. Although the numerical model proposed in section 2 can be used to model buried or immersed pipes, the additional modal, convergence studies and the interpretation of associated trapped or leaky wave modes require further investigation and won't be covered here for the sake of space. No coordinate stretching is applied in the radial direction, so that h 1 = 1, and d 1 = 0. In the axial direction of the pipe, the following exponential stretching functions are used: and where z = |z − z e |/h. z e is the entrance axial position of the PML and the thickness of the PML is h. α and β are real valued constants. In this article, the scaling parameter is fixed to be three, i.e., α = 3, and the damping parameter is defined as: where the ceil operator indicates rounding the element in the square bracket to the nearest integer greater than or equal to that element. R 0 is the reflection coefficient and fixed to be R 0 = 1*10 − 10 in this article. A commonly used method to validate a PML technique is to compare the PML truncated solution with an enlarged domain solution with fixed boundaries [18][19][20]. This method is used here to validate the accuracy of the PML model with a sound source located at the centre of the pipe in the fluid, i.e., r 0 = 0 and z 0 = 0. The total length of the pipe is 12 m in the PML model, including a PML thickness of h = 0.05 m at each end of the pipe. In the enlarged domain model, the total length of the pipe is 40 m. The centre frequency of the Hanning windowed pulse is 10 kHz. At this frequency, SF(0,1) has an energy velocity of 5013.6 m/s and FS(0,1) has an energy velocity of 1262.1 m/s. Eight-noded quadratic elements are used to mesh the fluid and pipe with an element density of at least 9 nodes per wavelength for the shortest bulk wave up to 40 kHz. In the PML model, 60,000 quadratic elements are generated in the fluid region, and 6000 are generated in the pipe region respectively. In the enlarged domain model, 160,000 quadratic elements are generated in the fluid  region, and 16,000 elements are generated in the pipe region respectively. The marching time step Δt is 0.002 ms. The total modelling duration N⋅Δt is 8.19 ms in both models. Fig. 3(a-c) show the comparison of time domain velocity potential solutions from both models at three receiver locations. The velocity potential presented here is normalised according to Eq. (54), and the normalised velocity potential is used throughout Section 3. The three receivers are all located on the central axis of the pipe, i.e., r = 0. Excellent agreement between the PML model and the enlarged model can be observed. The reflected sound waves from both pipe ends in the PML model have been totally absorbed. This shows the efficiency and accuracy of the fluid/solid coupled PML solution proposed in this article. Two modes have been excited with the sound source in the fluid. The fluid type mode FS(0,1) is the dominant mode, however, the structure type mode SF(0,1) can also be observed with an amplitude that is only around 1.2% of the amplitude of FS(0,1). In Fig. 3(a), the two wave modes are partially superimposed as it takes time for them to separate. The velocity of each wave mode can be simply verified by calculating the ratio of receiver distance to the travelling time of the corresponding pulse between two receivers. However, the accuracy of this method is influenced by the dispersion characteristics within the frequency bandwidth of the sound source.
The influence of PML thickness on the absorption of pipe end reflections is then examined and shown in Fig. 4. The receiver location is chosen to be z = 5 m, which is close to the PML entrance. The sound source and receiver are both located on the central axis of the pipe, i.e., r = 0. The out-going travelling waves at this receiver location have been shown in Fig. 3(c). Fig. 4 shows that only out-going waves can be observed when PML thickness h = 0.05 m or h = 0.1 m. However, both out-going and partially reflected waves from pipe end can be observed when PML thickness h = 0.01 m or h = 0.03 m. The radius of the pipe and the shortest wavelength of the propagating waves have been used as references for setting up the PML thickness. For all the four case studies shown in Fig. 4, the pipe end reflection from SF(0,1) has been totally absorbed. SF(0,1) has a faster wave speed than FS(0,1), however, the amplitude of SF(0,1) is relatively small and it is easier for PML to damp this wave mode. The pipe end reflection from FS(0,1) can be clearly seen in Fig. 4 when the PML thickness h = 0.01 m. In this case, the amplitude of the reflected FS(0,1) is less than 8% of the amplitude of the out-going FS(0,1), indicating majority energy of the reflected FS(0,1) has been absorbed by the PML. The PML thickness is then increased to h = 0.03 m, and the amplitude of the reflected FS(0,1) is reduced to be less than 0.5% of the amplitude of the out-going FS(0,1). However, partial shown in Ref. [14]. The PML thickness is fixed to be h = 0.05 m in the following studies. A more accurate method to calculate frequency-dependent wave velocities and identify wave modes is to carry out a two-dimensional Fourier transform (2 d FFT) to time-domain signals [29]. 128 receivers are placed on the central axis of the pipe with an even interval of 0.03 m between two adjacent receivers. The location of the first receiver is r = 0, and z = 1 m. The location of the last receiver is r = 0, and z = 4.81 m. The time domain velocity potentials at all the 128 receiver points are calculated in the PML model, and a two-dimensional Fourier transform is then applied to the assembled data stored in the form of a matrix. This transfers the data from the space-time domain to the wavenumber-frequency domain, shown in Fig. 5. To identify wave modes, the wavenumber solution obtained from the SAFE model [13] is plotted on Fig. 5 for comparison. It can be seen that the 2 d FFT solution matches the SAFE solution very well. This further validates the PML model proposed in this article. Note that the wavenumber resolution of the Fourier transformed signal shown in Fig. 5 can be further improved or reduced by adding or removing sampling points in the space domain. Fig. 6(a-f) show the fluid velocity potential/pipe radial displacement distribution at different times. Fig. 6(a-e) are solutions obtained from the PML model, and Fig. 6(f) is the solution obtained from the enlarged model. Excellent agreement has been observed between sound fields in the unstretched region in the PML and enlarged models, however, not all solutions from the enlarged model have been presented, for the sake of space. Fig. 6(e-f) show zoomed-in solutions from both models at the same time t = 5.6 ms, and can thus be compared directly. This time is chosen because majority part of the FS(0,1) wave has entered the PML, while the remaining part is still outside of the PML. This helps to examine whether any sound wave has been reflected by the PML. It can be seen in Fig. 6(e) that the PML (with a thickness of only 0.05 m) has absorbed all energy of the FS(0,1) wave that has entered the PML. The sound field in the unstretched region is identical to that shown in Fig. 6 (f) from the enlarged model. Fig. 6(a-f) also show there is a phase difference between the fluid velocity potential field and the pipe radial velocity field, which is seen as the slight misalignment between pulses in the axial direction. This is consistent with Eqs. 41 and 42, as the gradient of the potential gives radial velocity which is continuous across the fluid/pipe interface. In Fig. 6(a), both SF(0,1) and FS(0,1) can be observed. The structural type SF(0,1) mode has a wavelength of 0.52 m at 10 kHz, and a 10 cycle Hanning windowed pulse spans over a width of 5.2 m. The first and last 2 cycles of SF(0,1) are not clearly visible in Fig. 6(a), due to that the amplitude of this mode is very small compared to the FS(0,1) mode. At t = 2 ms and afterwards, SF(0,1) has entered the PML region and is absorbed, and thus this mode cannot be seen in Fig. 6(b-f). The fluid type mode FS(0,1) has a wavelength of 0.134 m at 10 kHz, and a 10 cycle Hanning windowed pulse spans over a width of 1.3 m. The pulse spanning width is actually increasing as this mode propagates along the pipe because of the dispersive nature of this wave mode.

Time domain wave propagation in fluid-filled pipes
The time domain PML model has been validated in Figs. 3-6 for a sound source located in the fluid. Now the PML model is used to study guided waves excited by an axisymmetric line source located around the outer circumference of the pipe at r = 0.04465 m, and z = 0. The incident sound is the 10 cycle Hanning windowed pulse with a centre frequency of 10 kHz. The structural excitation force is applied in the radial direction. Fig. 7. (a & b) show the time domain signals received in the fluid and on the pipe wall respectively. It can be seen that the dominant excited wave mode is the FS(0,1) mode. The velocity potential amplitude of SF(0,1) is around 4.1% of the amplitude of FS(0,1). Fig. 7 (b) show the displacements of these two wave modes. It is clear that FS (0,1) has a dominant radial displacement, while SF(0,1) has a dominant axial displacement. This implies that the FS(0,1) mode can be better detected by measuring the fluid potential field or the pipe radial displacement. The SF(0,1) mode can be better detected by measuring the pipe axial displacement in the low frequency range studied here. Fig. 8 (a & b) show the fluid velocity potential/pipe displacement distribution at t = 1 ms with a radial line source located around the outer circumference of the pipe at r = 0.04465 m, and z = 0. At this time, both FS(0,1) and SF(0,1) can be observed in the PML model. It is clear that the axial displacement of SF(0,1) is stronger than its radial displacement, which is consistent with Fig. 7 (b). The axial displacement of FS(0,1) is slightly distorted in the range between 0 and 1 m, due to the superimposition between FS(0,1) and the trailing part of SF(0,1). The axial displacement of FS(0,1) is also decreasing from the inner radius of the pipe to the outer radius.
To study the relationship between the excited wave modes and the direction of the sound source force term in Eq. (23), the excitation force is now applied in the axial direction of the pipe. Fig. 9. (a & b) show the time domain signals received in the fluid and on the pipe wall respectively, with an axial line source located around the outer circumference of the pipe at r = 0.04465 m, and z = 0. Fig. 9(a) shows that the velocity potential amplitude of the FS(0,1) mode is around 9% of the amplitude of the SF(0,1) mode. Note that SF(0,1) is a structural type mode, with most of the energy in the pipe wall. Comparing Fig. 9(a) with Fig. 7(a), the maximum fluid potential excited by the axial line source is much smaller than that excited by the radial line source. Fig. 9(b) shows that only the displacements of SF(0,1) can be observed. The displacements of FS(0,1) are too small compared to SF(0,1), that this mode is almost invisible. Fig. 9. (a & b) show that the dominant excited wave mode is SF (0,1) when the excitation force is applied in the axial direction of the pipe. Fig. 10 (a & b) show the fluid velocity potential/pipe displacement distribution at t = 1 ms with an axial line source located around the outer circumference of the pipe at r = 0.04465 m, and z = 0. It can be seen that the only visible mode is SF(0,1), and FS(0,1) is invisible because of the relative small amplitude of this wave mode. The SF(0,1) mode has a dominant displacement in the axial direction of the pipe. This mode is an extensional type wave mode, with most of its energy in the pipe wall in the low frequency range, which can also be seen in Fig. 2 (c). The radial displacement of this wave mode is relatively small, and consequently the fluid potential amplitude is small following the continuity of radial displacement across the fluid/pipe wall interface.

Conclusions
In this article, a numerical time-domain PML truncated model is proposed to study the excitation of guided waves in fluid filled pipes. In the elastic domain, a displacement-strain mixed formulation is proposed in the PML, whereas in the fluid domain, a velocity potential-velocity mixed formulation is proposed in the PML. The strain in the solid and particle velocity in the fluid are auxiliary variables introduced to guarantee that complex variables can be eliminated and convolution  operations can be avoided, following an inverse Fourier transform. The final systematic equations in the elastic and fluid domains are presented in an identical form and can be executed with the same extended Newmark time marching algorithm. Furthermore, the PML formulations in the elastic and fluid domains are compatible with each other, and the fluid/solid coupling in the PML can be implemented efficiently. This is done by substituting the solid particle velocity and displacement into the fluid domain on the coupling interface, and substituting the fluid potential and its derivative into the solid domain on the coupling interface through an iterative 'staggered scheme'.
The accuracy and efficiency of the PML model is examined and validated. The thickness of the PML is chosen to be 0.05 m in the fluid and solid domains. Time domain solutions are obtained. The twodimensional Fourier transform of time domain signals delivers frequency-wavenumber solution which can be compared to SAFE based modal solution, and excellent agreement has been observed. The PML model is then compared to an enlarged numerical model without PML, and identical solutions have been achieved in the unstretched region. The PML model is used to study the excitation of guided waves with a sound source in the fluid or on the pipe wall. In the low frequency range studied here, the dominant excited wave mode is the fluid type mode FS (0,1) when the sound source is in the fluid or when a radial line source is applied on the outer circumference of the pipe wall. However, the dominant excited wave mode is the structural type SF(0,1) mode when the sound source is an axial line source applied around the circumference of the pipe.

Author statement
W Duan is responsible for development of equations, finite element programming in Matlab, validation of the model, data analysis and drafting of the article.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.