SAFE-PML approach for modal study of waveguides with arbitrary cross sections immersed in inviscid fluid

Journal of Sound and Vibration 406 (2017) 181–196

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

SAFE-PML approach for modal study of waveguides with arbitrary cross sections immersed in inviscid fluid Peng Zuo, Zheng Fan n School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore

a r t i c l e i n f o

abstract

Article history: Received 9 February 2017 Received in revised form 30 May 2017 Accepted 1 June 2017 Handling Editor: L.G. Tham

Ultrasonic guided wave is an important non-destructive tool for large area inspections of immersed structures as well as fluid characterizations. In this paper, a numerical tool is developed for the modal study of immersed waveguides with arbitrary cross sections, by coupling the Semi-Analytical Finite Element (SAFE) method with Perfectly Matched Layer (PML). The model is first validated on waveguides with regular cross sections with analytical solutions. It is then applied to immersed waveguides with rectangular cross sections and L-shaped cross sections, showing the potential of guided waves for NDT applications and fluid characterizations. & 2017 Elsevier Ltd All rights reserved.

Keywords: Immersed waveguides SAFE-PML Inviscid fluid

1. Introduction Understanding properties of ultrasonic guided waves in immersed waveguides is important for various industries. On one hand, ultrasonic guided waves can be used as an efficient tool for Non-Destructive Testing (NDT) over a large area in immersed waveguides (such as immersed plates or underwater pipelines), as they can propagate for a long distance from a single generation position. On the other hand, the interaction between guided waves and the surrounding fluid make it possible to sense the presence and nature of the fluid [1–3]. In both cases, understanding modal properties, especially the attenuation of guided waves is essential. The modal properties in regular waveguides (such as plates, rods and pipes) immersed in fluid can be extracted through analytical methods, i.e. Transfer Matrix Method [4], Global Matrix Method [5] and dispersion equations [6–8]. For waveguides with irregular cross sections, Semi-Analytical Finite Element (SAFE) method has been widely used [9–14]. With such method, only the cross section needs to be meshed by finite element, and the wave is assumed to propagate harmonically along the propagation direction. More recently, the SAFE method has been extended to immersed/embedded waveguides [1,15–17]. Hladky-Hennion et al. [15] proposed radiation elements to model the infinite fluid as a finite domain. Fan et al. [1] used an absorbing region with artificial damping properties to absorb leaking waves thus simulating an infinite extent of the fluid medium. In this method, by coupling the SAFE method and the absorbing region modeling techniques (SAFE-AR), leaky waves can be absorbed and the attenuation can be extracted from the complex wavenumber when the absorbing region is long enough. An alternative approach to model leaky waves has been proposed by Mazzotti et al. [16] through coupling the SAFE method with a regularized 2.5D Boundary Element Method (BEM). In this method, the infinite fluid is described by BEM and n

Corresponding author. E-mail address: [email protected] (Z. Fan).

http://dx.doi.org/10.1016/j.jsv.2017.06.001 0022-460X/& 2017 Elsevier Ltd All rights reserved.

182

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

manipulated as a single, wavenumber and frequency dependent, finite element of infinite extension. More recently, Perfectly Matched Layer (PML) were introduced [17,18] in the SAFE method for the modelling of embedded waveguides, whose impedance perfectly match with the neighbor medium and it can absorb almost all propagation waves of non-tangential angles-of-incidence and any desired frequencies with relatively low computational cost. The concept of PML was first created in the context of electromagnetic waves in 1994 [19] and then it was viewed as the result of complex coordinate stretching [20,21]. Since then, the same ideas were also applied to other wave equations, such as Helmholtz equations [22], Euler equations [23] and even elastodynamic wave equations [24–27]. Coupling SAFE method with PML approach (SAFEPML) is a further extension of the PML method in elastodynamic wave equations. By combining of the SAFE and the PML method, the leaky waves can be damped, while the attenuation in the axial direction can be extracted and represented by the imaginary part of the wavenumber. However, most of the studies on SAFE-PML method focus on embedded waveguides, which only contains solid medium. This paper develops coupled SAFE-PML [28] method and applies it to waveguides immersed in an inviscid fluid. It starts with mathematical derivations of the coupled SAFE-PML equations, followed by the validation of two examples which can be solved by the matrix method program, including a steel plate with water-loaded on both sides and a cylindrical steel bar immersed in water. The validated method is then used to study modal properties on an immersed rectangular bar and an L-shaped bar immersed in water, showing the potential of guided waves for NDT applications and fluid characterizations.

2. Numerical model 2.1. Mathematical framework of the SAFE-PML method for immersed waveguides In this section, the SAFE-PML equations for waveguides immersed in an inviscid fluid have been developed. The PML method can be viewed as a result of complex coordinate stretching, where the ordinary solutions of wave equations in real coordinate space are mapped to a complex coordinate space in which the wave field decays exponentially. Through coordinate stretching, the wave is attenuated within a finite length, so that the infinite medium can be truncated into a finite domain [29]. The SAFE method starts with the three-dimensional approach and reduces the problem to two dimension by assuming harmonic wave propagating along the axial direction. Therefore, the SAFE method only needs to solve the equations in the cross section of the waveguide and hence, the coupled SAFE-PML method just needs to do the complex coordinate stretching in the cross section (as shown in Fig. 1). The model is composed of two parts: infinite inviscid fluid and a solid waveguide with arbitrary cross section immersed in the fluid. The cross section is parallel to the (x1, x2) plane and the wave propagates in the x3 direction. The stretched coordinates in the waveguide are defined as

Fig. 1. Schematics of (a) a solid waveguide with arbitrary cross section immersed in an infinite inviscid fluid and (b) two dimensional model with PML for calculation of modal properties of the waveguide, where d x and d x denote the start position of the PML in the x1 and x2 direction, respectively; h x and h x2 2 1 1 represent the length of the PML in the x1 and x2 direction, respectively; Ωs and Ωf indicate the solid domain and the fluid domain, respectively; ∂Ωinner is interface between the solid waveguide and the fluid and ∂Ωouter denotes the outside boundary of the PML.

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196 x1 ⎧ γ1( ξ )dξ ⎪ x˜1( x1) = 0 ⎪ x 2 ⎨ , ⎪ x˜2( x2) = 0 γ2( ξ )dξ ⎪ ⎩ x˜3( x3) = x3

183

∫

∫

(1)

where x˜1, x˜2, and x˜3 denote the stretched coordinates. γ1( x1), γ2( x2) are non-zero, continuous, complex-valued coordinate stretching functions, also called PML functions, which satisfy

γ1( x1) = 1 for x1 ≤ d x1 and Im{ γ1} > 0 for x1 > d x1,

(2)

γ2( x2) = 1 for x2 ≤ d x2 and Im{ γ2} > 0 for x2 > d x2.

(3)

Therefore, this change of variable yields for any function f˜ :

⎧ ∂f˜ 1 ∂f ⎪ = γ1 ∂x1 ⎪ ∂x˜1 ⎪ ⎪ ∂f˜ 1 ∂f ⎨ , = γ2 ∂x2 ⎪ ∂x˜2 ⎪ ∂f ⎪ ∂f˜ = ⎪ ∂x3 ⎩ ∂x˜3

(4)

where f˜ ( x˜1, x˜2, x˜3, t ) = f ( x1, x2, x3, t ). 2.1.1. Wave equations in the fluid With the assumption of harmonic wave propagating in the x˜3 direction, the acoustic pressure in the fluid can be written as ˜ p( x˜1, x˜2 , x˜3, t ) = P ( x˜1, x˜2)e I( kx3 − ωt ) ,

(5)

in which I and k represent the imaginary unit and wavenumber, respectively; ω = 2πf denotes the angular frequency and f being the frequency. Considering a compressive fluid with compressibility coefficient K, the acoustic pressure p is related to the displacement u f through the constitutive equations written as

˜ ·u f , p = − K∇

(6)

˜ denotes the gradient operator with respect to the stretched coordinates. The differential equations of motion in the where ∇ fluid domain are

˜ p, ρ f ω2u f = ∇

(7)

f

where ρ is the density of the fluid. Combining Eqs. (6) and (7), the differential equations of motion in the fluid can be written as

˜ ·( K ∇ ˜ p) + ρ f ω2p = 0. ∇

(8)

Using some intermediary transformations, Eq. (8) can be rewritten as

⎡ ⎛ ⎞2 2 ⎛ 1 ⎞⎛ 1 ⎞ ∂P ⎤ 1 ∂ P ⎜⎜ ⎟⎟⎜⎜ ⎟⎟′ ⎥ − Kk 2P + ρ f ω2P = 0, ⎟⎟ K + 2 γ x ∂ ⎝ γi ⎠⎝ γi ⎠ ∂xi ⎥⎦ ⎠ ⎝ i i i=1 ⎣ 2

∑ ⎢⎢ K ⎜⎜

(9)

in which the apostrophe means the derivative with respect to the independent variable. 2.1.2. Wave equations in the solid waveguide In the solid waveguide, the real coordinate space is coincident with the complex coordinate space (x˜1 = x1, x˜2 = x2 and x˜3 = x3), as the coordinate stretching functions are equal to one (γ1 = 1 and γ2 = 1). Thus equation derivation for the solid waveguide is placed in the real coordinate space for simplicity. The wave propagation in the solid waveguide is also assumed to be harmonic in the x3 direction, thus the displacement in the solid waveguide can be written as

uis ( x1, x2 , x3, t ) = Uis( x1, x2)e I( kx3 − ωt ) , s Ui

(10)

where is the displacement in the cross section of the waveguide. The differential equations of motion in the solid waveguide can be presented as

184

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

⎡ ∂ 2ujs ⎤ ⎥ + ρ s ω2uis = 0 for i = 1, 2, 3, ⎢ Cimjl ⎢ ∂x m∂xl ⎥⎦ j , m, l = 1 ⎣ 3

∑

(11)

s

where ρ is the density of the solid; Cimjl are the components of the elastic constant tensor, which can be related to the component of the stiffness matrix Cαβ (1 ≤ α , β ≤ 6) by the conventional notation: α = i = m if i = m, α = 9 − ( i + m) otherwise; β = j = l if j = l , β = 9 − ( j + l ) otherwise. Using some intermediary transformations and Einstein's convention of summation for economy of presentation, Eq. (11) can be written as

Cimjl

∂ 2U js ∂x m∂xl

(

)

+ I Ci3jm + Cimj3 k

∂U js ∂x m

− Ci3j3k 2U js + ρ s ω2δijU js = 0,

(12)

with summation over the indices j = 1, 2, 3 and m , l = 1, 2. Eq. (12) is a quadratic eigenvalue problem, where the eigensolutions are the wavenumber k and mode shape of the waveguide at the chosen angular frequency ω. In order to cast it into a general linear form, a new variable should be introduced [18,30–32]

V js = kU js.

(13)

Thus, the governing equations of the solid waveguide can be cast into

Cimjl

∂ 2U js ∂x m∂xl

(

∂V s

) ∂x j

+ I Ci3jm + Cimj3

− kCi3j3V js + ρ s ω2δijU js = 0,

(14)

m

ρ s ω2δijV js = kρ s ω2δijU js,

(15)

with summation over j = 1, 2, 3 and m , l = 1, 2, and δij is the Kronecker delta. 2.1.3. Fluid loading interface conditions For inviscid fluid, continuity of the normal displacement as well as stress is applied as the interface boundary condition. In the fluid domain, the following relations can be obtained

˜ p) = ρ f ω2K n˜ f ·us , n˜ f ·( K ∇

(16)

f

s

where n˜ is the outward unit vector of the fluid domain on the interface in the complex stretched coordinate; u is the displacement of the interface calculated in the solid domain. On the solid-fluid interface, the PML function should be equal to 1, which means that the real coordinate is coincident with the complex coordinate. Therefore n˜ f = n f is applied, where n f is the outward unit vector in the real coordinate. With the consideration of n3f = 0, the interface condition can be cast into

⎛ f ∂P ⎞ f ∂P f K ⎜ n1 + n2 ⎟ = ρ f ω2K n1f U1s + n2 U2s , ∂x2 ⎠ ⎝ ∂x1

(

)

(17)

In the solid waveguide domain, the boundary conditions at the interface can be written as

T = − pn s

(18)

where ns is the outward unit vector of the solid domain, related to the outward unit vector of the fluid domain by ns = − n f ; T is the stress vector at the interface which can be presented by

Ti = Cimjl

∂U js ∂xl

nms + ICimj3V jsnms = − Pnis ,

(19)

with summation over the indices j = 1, 2, 3 and m , l = 1, 2. 2.2. FEM formalism and implementation of the SAFE-PML formulas The corresponding finite element eigenvalue formalism with eigenvalue

∇·( c∇u + α u − r ) − au − β·∇u + daλ u − eaλ2u = f ,

λ is [33] (20)

in which all matrix coefficients can take real or complex values, and u represents the set of variables to be determined. The generalized Neumann boundary condition and Dirichlet boundary condition are expressed in the same FEM code as

n·( c∇u + α u − r ) + qu = g , u = h.

(21) (22)

For the fluid domain, Eq. (9) is compared with Eq. (20), the acoustic pressure needs to be chosen as the finite element

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

185

variable, u = P , and the matrix coefficients become

⎤ ⎡ ⎛ ⎞2 ⎥ ⎢ ⎜ 1⎟ K 0 ⎥ ⎢ ⎜⎝ γ ⎟⎠ 1 ⎥; a = − ρ f ω2; e = K; c = ⎢ a ⎢ ⎛ 1 ⎞2 ⎥ ⎥ ⎢ 0 ⎟ ⎜ K⎜ ⎟ ⎢⎣ ⎝ γ2 ⎠ ⎥⎦ ⎡ ⎛ ⎞⎛ ⎞ ⎤ ⎢ K ⎜⎜ 1 ⎟⎟⎜⎜ 1 ⎟⎟′ ⎥ ⎢ ⎝ γ1 ⎠⎝ γ1 ⎠ ⎥ ⎡ ⎤ ⎥; α = r = ⎢ 0⎥; da = f = 0. β = ⎢ ⎣ 0⎦ ⎞ ⎛ ⎞ ⎛ ⎢ 1 1 ⎥ ⎢ K ⎜⎜ ⎟⎟⎜⎜ ⎟⎟′⎥ ⎣ ⎝ γ2 ⎠⎝ γ2 ⎠ ⎦

(23)

⎡⎣ U s U s U s V s V s V s ⎤⎦T 1 2 3 1 2 3

For the solid domain, the new variables u = should be chosen as the finite element variable and the matrix coefficients can been obtained from Ref. [30]. Regarding to the fluid loading conditions at the interface between the solid waveguide and the fluid, Neumann conditions are used at the interface with

q = 0 and g = ρ f ω2K ( n1U1s + n2U2s ),

(24)

in the fluid domain and T q = 0 and g = ⎡⎣ Pn1 Pn2 0 0 0 0⎤⎦ ,

(25)

in the solid domain, where n1 and n2 are the component of the unit vector normal to the interface which are predefined in the package and can be called directly (here, n = n f is assumed. If n = − n f , the sign of the interface conditions will change). On the outer boundaries, the Neumann boundary conditions with q ¼ 0 and g ¼ 0 correspond to sound hard boundaries, which means the normal derivative of the acoustic pressure is zero; and sound soft boundaries, where the acoustic pressure vanishes, can also be implemented by use of the Dirichlet boundary condition u = 0 on the corresponding boundaries. 2.3. PML function, PML length and mode filtering In SAFE-PML method, the absorption efficiency of leaky waves in the PML strongly depends on the choice of the PML function and the PML length. Numerous PML functions have been tried in the literature. For example, a constant function has generally been chosen as a PML function [34,35]. However it may suffer from the artificial discontinuity of the PML function at the interface between the physical domain and the PML. A continuous function, e.g. γ = 1 + IΦ/ω , where Φ is a parabolic function inside the PML and ω is the angular frequency of the propagating wave, is usually chosen as the PML function for the problems in the time domain [36–38]. More recently, an exponential PML function is investigated and used to study the wave propagation in buried pipes [39]. In our cases, a continuous parabolic function for both the real and imaginary parts of the PML function is therefore used which has demonstrated good performance in problems in the frequency domain [18,28]:

⎧1 if x1 ≤ d x1 ⎪ ⎪ 2 ⎛ x −d ⎞ γ1( x1) = ⎨ ; x1 ⎪ 1 + γ^ ⎜⎜ 1 ⎟⎟ if x1 > d x 1 1 ⎪ h x1 ⎠ ⎝ ⎩

(26)

⎧1 if x2 ≤ d x2 ⎪ ⎪ 2 ⎛ x −d ⎞ γ2( x2) = ⎨ , x2 ⎪ 1 + γ^ ⎜⎜ 2 ⎟⎟ if x2 > d x 2 2 ⎪ h x2 ⎠ ⎝ ⎩

(27)

where γ^1 = a1 + Ib1( b1 > 0) and γ^2 = a2 + Ib2 ( b2 > 0) quantify the PML absorption and will be given explicitly in each case study. As leaky waves grow exponentially in the transverse directions, placing the PML close to the waveguide can reduce the effect of the exponential growth of the leaky modes. In order to estimate the length of the PML, a simplified 2D plane wave propagation model can be used to approximately predict the length of the PML with the formula [28]

h≥

6.9 , kleakb

(28)

186

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

where h is the length of the PML; b is the imaginary part of the PML function and kleak represents the wavenumber of the longitudinal wave in the fluid. When the PML function is given, b can be determined. The length of the PML can be obtained by calculating the minimum wavenumber of the longitudinal wave in the frequency range of interest. In the calculation, three kinds of modes can be found in the immersed waveguide system. They are trapped mode, leaky mode and radiation mode [17]. As the name indicates, trapped modes propagate in the waveguide with energy concentrating in the waveguide. Leaky modes propagate along the waveguide with some energy leaking into the surrounding fluid. Radiation modes resonate mainly in the fluid domain, and they are of less interest in practical applications. In this paper, a filtering condition is applied in post-processing to identify and remove the radiation modes:

P¯Ωs

> η,

P¯Ωf

(29)

where η is a user-defined parameter, identifying the criterion of the model; P¯ represents the average energy flux along the wave propagation direction x3 and is defined as

P¯ =

∫S Px3 ds ∫S ds

, (30)

in which S denotes the area of the cross section. Px3 represents the absolute value of the energy flux. In the solid waveguide, the energy flux is defined as

⎡ ⎛ Iω ⎞ ⎤ s s s + u2s *σ32 + u3s *σ33 Px3 = − Re⎢ ⎜ ⎟ u1s *σ31 ⎥, ⎣⎝ 2 ⎠ ⎦

(

)

(31)

where the asterisk stands for the complex conjugate, and the component of the stress tensor can be calculated by

⎛ ∂u s ⎞ s = C55⎜ Iku1s + 3 ⎟; σ31 ∂x1 ⎠ ⎝ ⎛ ∂u s ⎞ s σ32 = C44⎜ 3 + Iku2s ⎟; ⎝ ∂x2 ⎠ s σ33 = C13

∂u1s ∂u s + C23 2 + C33Iku3s . ∂x1 ∂x2

(32)

In the fluid, the energy flux can be calculate by the following formula

⎡ ⎛ Iω ⎞ f * ⎤ Px3 = − Re⎢ ⎜ ⎟u3 σ33f ⎥, ⎦ ⎣⎝ 2 ⎠ where

f σ33

= − p and the displacement

u3f =

Ikp ρ f ω2

(33) f u3

can be calculated by

.

(34)

All of the quantities used in the post-processing can be easily extracted from the eigensolutions of the FEM solver.

3. Validation of the SAFE-PML method Two validation cases of immersed waveguides with regular cross sections are chosen to calculate using the SAFE-PML method, and the solutions from analytical methods are used for comparison. 3.1. A steel plate with water-loaded on both sides The first example is a steel plate with water-loaded on both sides. This is a representative geometry that has been studied by many other researchers [40–43]. It has a simple geometry, and therefore can be easily solved by the matrix method [44]. 3.1.1. Model description To model the propagation of wave along an infinitely wide plate, a narrow strip of the structure is used with periodic boundary conditions [30] representing continuity of displacements and stresses (or pressure in the water) between the two edges, as shown in Fig. 2. The narrow strip of the plate is 0.2 mm in width and 1 mm in thickness. The thickness of the water is 1 mm, which is attached to the plate on both sides. According to Eq. (28), the thickness of

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

187

Fig. 2. Schematic of the SAFE-PML model used for a steel plate with water-loaded on both sides.

the PML is dependent on the PML function and the minimum wavenumber of the longitudinal wave in the water at the frequencies of interest. In the simulation, the PML function is γ^2 = 3 + 12I and 29 frequencies will be solved from 100 kHz to 1500 kHz with increments of 50 kHz. The minimum wavenumber of the longitudinal wave can be obtained at f = 100 kHz : kleak = min( kL ) = 418.88 m−1, where kL is the wavenumber of the longitudinal wave in the water. Thus, based on Eq. (28), a rough estimation of the length of the PML can be obtained: h x2 ≥ 1.37 mm . In the simulation, the length of the PML is chosen to be 1.5 mm. The material properties are given in Table 1. Water loading interface conditions are implemented at the interface between the plate and the water according to Eqs. (24) and (25). In addition, Neumann boundary conditions with q ¼ 0 and g ¼ 0 are also imposed on the outer surface of the PML. The whole geometry is meshed by 868 triangular elements of the first order. These elements are automatically generated by the software and the number of degrees of freedom is 3660. In the post-processing, η = 2 is used to identify and filter the radiation modes.

3.1.2. Results Fig. 3 presents dispersion curves for the phase velocity and the attenuation, which are computed from the real and imaginary part of the wavenumber, respectively. Excellent agreements can be observed between the results from the SAFEPML method and analytical solutions. Three modes including A0 mode, SH0 mode and S0 mode can be found in the frequencies and some obvious features can be obtained. Since the A0 mode has a significant out of plane displacement leading to a strong coupling between the solid and the f luid, the attenuation of the A0 mode is relatively larger than that of the S0 mode which means a large part of energy of the A0 wave leaks into the water. In contrast, the SH0 mode is non-leaky over the range of the frequencies. Additionally, the A0 mode is not continuous when it crosses the longitudinal velocity of the water, causing a significant increase in attenuation at that point. Table 1 Mechanical properties of materials used in the models. Material

Steel Water

(kg/m3)

Longitudinal velocity (m/s )

Shear velocity (m/s )

7932 1000

5960.5 1500

3260 …

Density

188

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

Fig. 3. (a) Phase velocity and (b) attenuation dispersion curves for a steel plate with water-loaded on both sides. The dots represent the solutions from the SAFE-PML method, and the solid lines are the solutions obtained from the matrix method program [44]. The horizontal dashed lines in (a) indicate the velocity of the water.

3.2. A cylindrical steel bar immersed in water The second validation case includes a cylindrical steel bar immersed in water. This problem can also be solved by the matrix method program and has been used as a validation case for the SAFE model in the literature [1]. 3.2.1. Model description The schematic of the second validation case is shown in Fig. 4. The steel bar is 2 mm in diameter and immersed in water which is modeled by a square with length of 4 mm. The PML functions in this model are γ^ = 3 + 12I and γ^ = 3 + 12I, and the 1

2

model is solved for 39 frequencies from 100 kHz to 2000 kHz with increments of 50 kHz. Based on Eq. (28), the minimum length of the PML is h x1 = h x2 ≥ 1.37 mm . Thus a layer of 1.5 mm is used as the PML attached to the water. The material properties are listed in Table 1. Water loading interface conditions to describe the interaction between the cylindrical bar and the surrounding water are implemented based on Eqs. (24) and (25). Neumann boundary conditions with q ¼ 0 and g ¼ 0 are imposed on the outside boundaries of the PML. The whole geometry is meshed by 2206 triangular elements of the first order and the number of degrees of freedom is 6326. In the post-processing, η = 3 is used to identify and filter the radiation modes.

Fig. 4. Schematic of the SAFE-PML model for a cylindrical steel bar immersed in water.

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

189

Fig. 5. (a) Phase velocity and (b) attenuation dispersion curves for a cylindrical steel bar immersed in water. The dots represent the solutions from the SAFE-PML method, and the solid lines are the solutions obtained from the matrix method program [44]. The horizontal dash lines in (a) indicate the longitudinal velocity of the water. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

3.2.2. Results Fig. 5 presents the dispersion curves for leaky guided waves in the cylindrical bar. Excellent agreements can be found between the SAFE-PML method and the analytical solutions in the dispersion curves for both phase velocity and attenuation. It can be seen that the F(1,1) mode goes across the longitudinal velocity of the water around 250 kHz. When the phase velocity of the F(1,1) mode is slower than that of the water, no energy leaks into the water; while when it goes across the velocity of the water, some energy leaks into the water, and the F(1,1) mode changes from non-leaky to leaky [5]. Fig. 6 presents the mode shape of the guided modes at 1000 kHz and the corresponding displacements on the x2 axis in the cross section, which can give an intuitive interpretation of the absorption efficiency of the PML and attenuation feature of the modes shown in Fig. 5(b). The arrows in the mode shapes represent displacements in the cross section and the color indicates the energy flux in the x3 direction. According to the displacements, it is obvious that leaky waves can be totally damped in a very short distance when they go into the PML, showing the PML works properly, being an efficient tool to attenuate leaky waves. It also can be seen that the F(1,1) mode have antisymmetric mode shape with energy concentrates in the two opposite corners where the radial displacements dominate and strongly couple into the water, leading to a relatively high attenuation. As no displacement coupling occurs between the cylindrical bar and the water in the T(0,1) mode (shown in Fig. 6(c) and (d)), there is no displacement in the surrounding water in this mode, thus the wave propagates without attenuation. The symmetrical longitudinal mode L(0,1) is dominant in the axial direction, together with most of the energy concentrating in the center of the bar, thus this mode is weakly coupled with the water, showing low attenuation at this frequency. By comparison to the SAFE-AR model in Ref. [1], it is obviously that, only 1.5 mm PML can effectively damp the leaky waves over the frequency region, which is about one fourth of the length of the absorbing region in the SAFE-AR model, therefore leading to a significant reduction in the computational cost as well as the number of the radiation modes. From the two validation cases, it can be concluded that the SAFE-PML method is an efficient method to study the properties of leaky guided waves in immersed waveguides.

4. SAFE-PML method in waveguides with irregular cross sections In this section, the SAFE-PML method is applied to a rectangular steel bar and an L-shaped steel bar, where no analytical solutions can be found and small material damping components are considered in the latter case. 4.1. A rectangular steel bar immersed in water Rectangular bars are typical structures with a wide application in f luid characterization [1,45–47], therefore the study of modal properties of the rectangular bar is practically important. However, only particular guided modes in such waveguide were studied in previous work, while no effort has been attempted to understand modal properties of all guided modes due to lack of efficient models. 4.1.1. Modal description The geometry of the model is shown in Fig. 7. A rectangular steel bar of 1.1 mm × 2.2 mm immersed in water is studied. The surrounding water is modeled as a square shape with length of 4 mm. The PML functions in the simulation are chosen

190

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

Fig. 6. Mode shapes (on the left) and displacements (on the right) on the x2 axis from the ordinate origin to the endpoint of the PML for the cylindrical bar immersed in water at 1000 kHz: (a) mode shapes and (b) displacements of the F(1,1) mode with k = 2338.35 + 24.20I m−1; (c) mode shapes and (d) displacements of the T(0,1) mode with k = 1927.36 m−1; (e) mode shapes and (f) displacements of the L(0,1) mode with k = 1255.66 + 8.21I m−1. Ur and Uθ represent the radial displacement and the circumferential displacement in the cross section, respectively; U3 is the displacement in the axial direction. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

as γ^1 = 3 + 12I and γ^2 = 3 + 12I. The SAFE-PML model is solved for 39 frequencies from 50 to 1000 kHz with increments of 25 kHz. The minimum length of the PML can be calculated based on Eq. (28): h x1 = h x2 ≥ 2.75 mm . In this study, 3 mm is used as the length of the PML. All of the material properties are given in Table 1.

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

191

Fig. 7. Schematic of the SAFE-PML model for a rectangular steel bar immersed in water.

Similarly as in previous examples, fluid loading conditions are imposed at the interface between the rectangular bar and the water. Neumann boundary conditions with q ¼ 0 and g ¼ 0 are used on the outside boundaries of the PML. The whole geometry is meshed with 2064 triangular elements and the number of degrees of freedom is 4872. The mode filter parameter η = 3 is used in the post-processing to pick up the guided mode propagating in the rectangular bar. 4.1.2. Results Fig. 8 presents the dispersion curves for leaky guided waves in the rectangular bar. The modes have been named as in Ref. [48] according to the convention for modes in a rectangular bar in vacuum. It can be seen that the first order longitudinal mode L1 is little dispersive and almost non-leaky at a relatively low frequency range (0 − 700 kHz ). These properties make this mode a potential candidate for NDT application in immersed rectangular bars. The two features have also been observed in the first order torsional mode T1, the difference being that the T1 mode experiences almost no attenuation in a relatively narrow frequency range (0 − 200 kHz). In this range, the non-circular cross section can cause the fluid to be trapped at the corners of the cross section and thus will affect the propagation of the torsional wave. Therefore, this mode can sense the presence and nature of the surrounding fluid and can be used as a sensitive tool for fluid characterizations [1].

Fig. 8. (a) Phase velocity and (b) attenuation dispersion curves for a rectangular steel bar immersed in water. The horizontal dash lines in (a) indicate the longitudinal velocity of the water.

192

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

Fig. 9. Mode shapes for the rectangular bar immersed in water at 500 kHz for (a) B1x mode with k = 1646.23 + 85.56I m−1, (b) B1y mode with k = 1328.72 + 16.03I m−1, (c) T1 mode with k = 1292.22 + 20.46I m−1 and (d) L1 mode with k = 606.40 + 1.05I m−1. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

The two bending modes (B1x and B1y ) are similar to the F(1,1) mode in the cylindrical bar. They have non-leaky feature when the phase velocity of the mode is slower than that of the water. However, as their phase velocities go above the velocity of the water, the attenuation increases rapidly due to energy leaking into the water. The corresponding mode shape of the guided modes at 500 kHz are presented in Fig. 9. The arrows represent displacements in the cross section, and the color indicates the energy flux in the x3 direction. The bending modes (B1x and B1y ) dominate in the x1 and x2 direction and are strongly coupled with the water in the dominant vibration direction, thus leading to a lot of energy radiating into the water. It also can be seen that the torsional mode (T1) cause displacements in the water and has a relatively high attenuation at the frequency, which is in contrast to the case in the cylindrical bar where the torsional mode is non-attenuative. In the longitudinal mode (L1), as the dominant displacement is in the x3 direction, the vibration in the x1 and x2 direction becomes weak, so that it is slightly attenuated at the frequency.

4.2. An L-shaped steel bar immersed in water The second case is to study the modal properties of an L-shaped steel bar immersed in water. Small material damping in the steel is considered in the model, including longitudinal bulk wave attenuation βL = 0.003 Np/wavelength and shear bulk wave attenuation βS = 0.008 Np/wavelength , which has been reported by using SAFE-BEM in the literature [49]. In this case, the elastic constant tensor becomes a complex one combining elastic properties and viscous properties of the material, which can be easily calculated from the assumed rheological behavior of the material [49].

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

193

Fig. 10. Schematic of the SAFE-PML model for an L-shaped steel bar immersed in water.

4.2.1. Modal description Fig. 10 presents the schematic of the SAFE-PML model for the L-shaped bar. The L-shaped bar with dimension 30 × 20 × 4 mm is immersed in water which is modeled by a square with length of 40 mm. The PML function in this model are γ^ = 3 + 12I and γ^ = 3 + 12I, and the model is solved for 19 frequencies from 10 kHz to 100 kHz with increments of 1

2

5 kHz. The minimum length of the PML can be calculated based on Eq. (28): h x1 = h x2 ≥ 13.73 mm . In this case, a layer of 15 mm is used as the PML attached to the water. The material elastic properties are listed in Table 1. Fluid loading conditions are imposed at the interface between the L-shaped bar and the water. Neumann boundary conditions with q ¼ 0 and g ¼ 0 are used on the outside boundaries of the PML. The whole geometry is meshed by 3210 triangular elements with the number of degrees of freedom being 12522. η = 3 is used in the post-processing to pick up guided modes propagating in the L-shaped bar. 4.2.2. Results Fig. 11 shows the dispersion curves for leaky guided waves in the L-shaped bar. Some low order modes have been named as in Ref. [49]. It is obviously that the m1 mode (longitudinal mode) is non-dispersive and almost non-leaky at low frequency range (0–40 kHz), and is therefore particularly useful for NDT applications in this structure. The m2 mode (first pseudo-flexural mode) experience weakly dispersive over the whole frequency range, and it changes from non-leaky to

Fig. 11. (a) Phase velocity and (b) attenuation dispersion curves for an L-shaped steel bar immersed in water. The horizontal dash lines in (a) indicate the longitudinal velocity of the water.

194

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

Fig. 12. Mode shapes for the L-shaped bar immersed in water at 30 kHz for (a) m4 mode with k = 187.53 + 0.11I m−1, (b) m3 mode with k = 174.20 + 0.12I m−1, (c) m2 mode with k = 125.98 + 2.10I m−1, (d) h1 mode with k = 83.93 + 1.14I m−1, (e) h2 mode with k = 74.07 + 3.03I m−1 and (f) m1 mode with k = 36.25 + 0.07I m−1. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

leaky with attenuation increasing monotonically. For the m3 (second pseudo-flexural mode) and m4 (pseudo-torsional mode) modes, they are non-leaky over a relatively large frequency range as their phase velocities are slower than that of the water. However, as the phase velocity is very close to the water velocity, corresponding to a transition zone, they would not be suitable for practical applications [49]. Comparing to the SAFE-BEM model, it can be found that the phase velocity dispersion curves agree very well with that predicted from SAFE-BEM method in Ref. [49]. For the attenuation dispersion curves, most of the modes also agree well with the prediction from SAFE-BEM method, although some discrepancies are found in the m3 and m4 mode at high frequency

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

195

region (90–100 kHz). The possible reason is that in such frequency range, the phase velocity of these modes are close to the longitudinal velocity of the water, causing inaccuracy in the SAFE-PML method. The same problem is also encountered in Ref. [18] and enhancing the absorption strength of the PML may improve the results. Fig. 12 presents the mode shape of the guided modes at 30 kHz. The arrows represent displacements in the cross section, and the color indicates the energy flux in the x3 direction. It can be observed that the m3 and m4 modes have similar characteristics on the mode shape and energy distribution with energy flux concentrating on the end of the L-shaped bar. As the h1 and h2 modes are high order modes, both of them show complex mode shape and energy flux distribution. In addition, the m1 mode is a longitudinal mode dominating in the x3 direction, thus most of the energy concentrates in the L-shaped bar and the attenuation should be small due to the weak coupling between the bar and the water.

5. Conclusions In this work, the SAFE-PML method has been developed for waveguides immersed in an infinite inviscid fluid. The governing equations are derived by combining the SAFE method with the PML method, taking account of the solid-fluid interaction conditions. As it is implemented in a commercial FE package, the SAFE-PML method provides a fast and efficient way to investigate the modal properties of the immersed waveguides where energy can leak into the surrounding fluid. The method has been first validated on two immersed waveguides with regular cross sections, including a steel plate with water-loaded on both sides and a cylindrical steel bar immersed in water, showing good agreements with analytical solutions. The method is then applied to an immersed rectangular bar and an L-shaped bar immersed in water, illustrating the potential of the guided waves in NDT and fluid characterizations. It should be noted that, the main drawback of the SAFE-PML method is to select optimal PML parameters, including the PML function and the length of the PML. It is clear that an optimal PML function can reduce the length of the PML leading to a relatively low computational cost as well as the number of the radiation modes. Even though the PML function used in this paper shows good performance, and the length of the PML can be predicted by a 2D plane wave propagation model, the selection of the optimal PML parameters is still an open problem, which requires extra effort in future work.

Acknowledgements This work was supported by the Singapore Maritime Institute under SMI Simulation and Modelling R&D Programme.

References [1] Z. Fan, M.J.S. Lowe, M. Castaings, C. Bacon, Torsional waves propagation along a waveguide of arbitrary cross section immersed in a perfect fluid, J. Acoust. Soc. Am. 124 (4) (2008) 2002–2010. [2] T.K. Vogt, M.J.S. Lowe, P. Cawley, Measurement of the material properties of viscous liquids using ultrasonic guided waves, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 51 (6) (2004) 737–747. [3] F. Cegla, P. Cawley, M.J.S. Lowe, Material property measurement using the quasi-Scholte mode – a waveguide sensor, J. Acoust. Soc. Am. 117 (3) (2005) 1098–1107. [4] A.H. Nayfeh, P.B. Nagy, General study of axisymmetric waves in layered anisotropic fibers and their composites, J. Acoust. Soc. Am. 99 (2) (1996) 931–941. [5] B.N. Pavlakovic, Leaky Guided Ultrasonic Waves in NDT (Ph.D. thesis), Imperial College, London, 1998. [6] R.N. Thurston, Elastic waves in rods and clad rods, J. Acoust. Soc. Am. 64 (1) (1978) 1–37. [7] J.A. Simmons, E. Drescher-Krasicka, H. Wadley, Leaky axisymmetric modes in infinite clad rods. I, J. Acoust. Soc. Am. 92 (2) (1992) 1061–1090. [8] M. Viens, Y. Tsukahara, C.K. Jen, J.D.N. Cheeke, Leaky torsional acoustic modes in infinite clad rods, J. Acoust. Soc. Am. 95 (2) (1994) 701–707. [9] T. Hayashi, W.J. Song, J.L. Rose, Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example, Ultrasonics 41 (3) (2003) 175–183. [10] L. Gavrić, Finite element computation of dispersion properties of thin-walled waveguides, J. Sound Vib. 173 (1) (1994) 113–124. [11] X. Yu, M. Ratassepp, P. Rajagopal, Z. Fan, Anisotropic effects on ultrasonic guided waves propagation in composite bends, Ultrasonics 72 (2016) 95–105. [12] P. Manogharan, P. Rajagopal, K. Balasubramaniam, Longitudinal guided waves confined in radius filler regions of composite joints, J. Acoust. Soc. Am. 140 (1) (2016) 334–343. [13] P. Rajagopal, R.K. Pattanayak, Ultrasonic guided waves in elliptical annular cylinders, J. Acoust. Soc. Am. 138 (3) (2015) EL336–EL341. [14] A. Ramdhas, R.K. Pattanayak, K. Balasubramaniam, P. Rajagopal, Antisymmetric feature-guided ultrasonic waves in thin plates with small radius transverse bends from low-frequency symmetric axial excitation, J. Acoust. Soc. Am. 134 (3) (2013) 1886–1898. [15] A.C. Hladky-Hennion, P. Langlet, M. De Billy, Finite element analysis of the propagation of acoustic waves along waveguides immersed in water, J. Sound Vib. 200 (4) (1997) 519–530. [16] M. Mazzotti, I. Bartoli, A. Marzani, E. Viola, A coupled SAFE-2.5D BEM approach for the dispersion analysis of damped leaky guided waves in embedded waveguides of arbitrary cross-section, Ultrasonics 53 (7) (2013) 1227–1241. [17] F. Treyssède, K.L. Nguyen, A.S. Bonnet-BenDhia, C. Hazard, Finite element computation of trapped and leaky elastic waves in open stratified waveguides, Wave Motion 51 (7) (2014) 1093–1107. [18] K.L. Nguyen, F. Treyssède, C. Hazard, Numerical modeling of three-dimensional open elastic waveguides combining semi-analytical finite element and perfectly matched layer methods, J. Sound Vib. 344 (2015) 158–178. [19] J.P. Bérenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114 (2) (1994) 185–200. [20] C.M. Rappaport, Perfectly matched absorbing boundary conditions based on anisotropic lossy mapping of space, IEEE Microw. Guid. Wave Lett. 5 (3) (1995) 90–92. [21] F.L. Teixeira, W.C. Chew, General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media, IEEE Microw. Guid. Wave Lett. 8 (6) (1998) 223–225.

196

[22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49]

P. Zuo, Z. Fan / Journal of Sound and Vibration 406 (2017) 181–196

I. Harari, M. Slavutin, E. Turkel, Analytical and numerical studies of a finite element PML for the Helmholtz equation, J. Comput. Acoust. 8 (01) (2000) 121–137. F.Q. Hu, On absorbing boundary conditions for linearized Euler equations by a perfectly matched layer, J. Comput. Phys. 129 (1) (1996) 201–219. W.C. Chew, Q.H. Liu, Perfectly matched layers for elastodynamics: a new absorbing boundary condition, J. Comput. Acoust. 4 (04) (1996) 341–359. F.D. Hastings, J.B. Schneider, S.L. Broschat, Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation, J. Acoust. Soc. Am. 100 (5) (1996) 3061–3069. U. Basu, A.K. Chopra, Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation, Comput. Methods Appl. Mech. Eng. 192 (11) (2003) 1337–1375. S. François, M. Schevenels, G. Lombaert, G. Degrande, A two-and-a-half-dimensional displacement-based PML for elastodynamic wave propagation, Int. J. Numer. Methods Eng. 90 (7) (2012) 819–837. P. Zuo, X.D. Yu, Z. Fan, Numerical modeling of embedded solid waveguides using SAFE-PML approach using a commercially available finite element package, NDTE Int. 90 (2017) 11–23. F.L. Teixeira, W.C. Chew, Complex space approach to perfectly matched layers: a review and some new developments, Int. J. Numer. Model. Electron. Netw. Device Fields 13 (2000) 441–455. M.V. Predoi, M. Castaings, B. Hosten, C. Bacon, Wave propagation along transversely periodic structures, J. Acoust. Soc. Am. 121 (4) (2007) 1935–1944. M. Castaings, M. Lowe, Finite element model for waves guided along solid systems of arbitrary section coupled to infinite solid media, J. Acoust. Soc. Am. 123 (2) (2008) 696–708. T. Hayashi, D. Inoue, Calculation of leaky lamb waves with a semi-analytical finite element method, Ultrasonics 54 (6) (2014) 1460–1469. COMSOL, User’s Guide and Introduction, COMSOL MULTIPHYSICS. URL 〈http://www.comsol.com/http://www.comsol.com/〉. (Accessed 20 October 2016). Y.O. Agha, F. Zolla, A. Nicollet, S. Guenneau, On the use of PML for the computation of leaky modes, COMPEL 27 (2008) 95–109. A. Pelat, S. Felix, V. Pagneux, A coupled modal-finite element method for the wave propagation modeling in irregular open waveguides, J. Acoust. Soc. Am. 129 (3) (2011) 1240–1249. I. Harari, U. Albocher, Studies of FE/PML for exterior problems of time-harmonic elastic waves, Comput. Methods Appl. Mech. Eng. 195 (29) (2006) 3854–3879. F. Collino, C. Tsogka, Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics 66 (1) (2001) 294–307. A. Bermúdez, L. Hervella-Nieto, A. Prieto, R. Rodrıguez, An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems, J. Comput. Phys. 223 (2) (2007) 469–488. W. Duan, R. Kirby, P. Mudge, T.H. Gan, A one dimensional numerical approach for computing the eigenmodes of elastic waves in buried pipelines, J. Sound Vib. 384 (2016) 177–193. A.H. Nayfeh, D. Chimenti, Propagation of guided waves in fluid-coupled plates of fiber-reinforced composite, J. Acoust. Soc. Am. 83 (5) (1988) 1736–1743. F. Ahmad, N. Kiyani, F. Yousaf, M. Shams, Guided waves in a fluid-loaded transversely isotropic plate, Math. Probl. Eng. 8 (2) (2002) 151–159. G. Belloncle, H. Franklin, F. Luppé, J. Conoir, Normal modes of a poroelastic plate and their relation to the reflection and transmission coefficients, Ultrasonics 41 (3) (2003) 207–216. G. Belloncle, F. Luppé, H. Franklin, J.M. Conoir, Influence of the slow wave on the relation between the angular resonances and the leaky guided modes properties for a poroelastic plate embedded in water, Ultrasonics 42 (1) (2004) 511–514. B.N. Pavlakovic, M.J.S. Lowe, DISPERSE: A System for Generating Dispersion Curves, User’s Manual, 2003. H.H. Bau, Torsional wave sensor – a theory, J. Appl. Mech. Trans. ASME 53 (4) (1986) 846–848. C.C. Smit, E.D. Smith, The analysis and results of a continuous wave ultrasonic densitometer, J. Acoust. Soc. Am. 104 (3) (1998) 1413–1417. C.L. Shepard, B.J. Burghard, M.A. Friesel, B.P. Hildebrand, X. Moua, A.A. Diaz, C.W. Enderlin, Measurements of density and viscosity of one-and twophase fluids with torsional waveguides, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46 (3) (1999) 536–548. W.B. Fraser, Stress wave propagation in rectangular bars, Int. J. Solids Struct. 5 (4) (1969) 379–397. M. Mazzotti, A. Marzani, I. Bartoli, Dispersion analysis of leaky guided waves in fluid-loaded waveguides of generic shape, Ultrasonics 54 (1) (2014) 408–418.