A Spectral Method Algorithm for Modeling the Dispersion of Non-Axisymmetric Modes in Fluid-Filled Elastic Tubes

Abstract

In order to design a low-noise water-filled pipeline system, it is necessary to obtain knowledge of the dispersion characteristics of axial propagation modes in different water-filled elastic tubes. In this work, an algorithm is developed based on the spectral method, which has previously been used to solve the dispersion of axisymmetric modes in cylindrical structures but has not yet been applied to non-axisymmetric modes. The algorithm can obtain the dispersion characteristics, modal displacement, and stress distribution of axial propagation modes in a fluid-filled elastic multi-layer tube. The algorithm behaves well both at low and ultrasonic frequencies, and it is suitable for any tube dimensions, wall thickness and layers. The results of a water-filled PMMA tube obtained using the spectral method algorithm were verified using a COMSOL simulation, while the dispersion curves of the same tube from the literature were found to be missing some low-order modes. In addition, the dispersion curves of a water-filled three-layer tube are given. The spectral method algorithm has the advantages of fast calculation speed, less computational resources consumed, accurate results, and no modal omission.

1. Introduction

Fluid-filled pipeline systems are widely used in engineering. Pumps, valves, and other pipeline system components often vibrate during operation. Vibrational energy propagates in the fluid and pipe walls in the form of waves, producing a negative effect on the system performance and the surrounding environment. Axially propagating waves contain axisymmetric and non-axisymmetric modes related to frequencies. The accurate and efficient calculation of dispersion curves of theses modes and the corresponding displacement and stress distribution is a key step in studying the dispersion characteristics of propagation modes in a fluid-filled elastic tube, which lays a theoretical foundation for the subsequent design of low-noise pipeline systems and non-destructive pipeline testing tools.

A great deal of research has been carried out on the dispersion characteristics of axially propagating modes in hollow cylindrical shells and fluid-filled tubes. Jacobi [1] took the ideal liquid cylinder as a research object and obtained the dispersion curves of low-frequency axisymmetric modes under various non-dissipative boundary conditions: rigid walls, pressure-release walls, an infinite liquid boundary, liquid walls, and thin solid walls. In order to study propagation modes over a wider frequency range, Lin and Morgan [2] used the approximate equations governing the motion of an elastic tube wall to derive the characteristic equation of the modes. Gazis [3,4] described the motion of cylindrical shells using the linear elastic theory proposed by Pochhammer [5] and Chree [6] and obtained the eigenmode characteristic equation of an infinitely long isotropic hollow cylinder. Del Grosso [7] calculated the dispersion relation of axially propagating modes in liquid-filled tubes with arbitrary wall thicknesses based on the exact longitudinal and transverse wave equations. Kumar [8] obtained the dispersion equation of axially propagating modes in a liquid-filled cylindrical shell using the exact three-dimensional linear elastic equation and the fluid motion equation and discussed the characteristics of the bending mode ($n=1$). Using the theory of Del Grosso [7], Lafleur and Shields [9] calculated the dispersion relation of axisymmetric propagating modes in a water-filled elastic tube. Using linear elasticity and classical perfect slip boundary conditions at the solid–fluid boundary, Sinha et al. [10,11] analyzed the behavior of axisymmetric longitudinal waves in cylindrical shells with different internal and external fluids and also carried out an experimental verification. Easwaran and Munjal [12] used the same method to study the lowest-order acoustic modes in a liquid-filled impedance tube. Greenspon and Singer [13] replaced the elastic constants in the three-dimensional elastic equation with complex viscoelastic material parameters and studied the axisymmetric modes and non-axisymmetric modes in a liquid-filled viscoelastic tube. Berliner and Solecki [14,15] also derived the dispersion equation of wave propagation in a transversely isotropic elastic tube from the three-dimensional elastic equation and discussed the dispersion characteristics of the non-axisymmetric mode $n=1$. Baik [16] applied Del Grosso’s [7] theory to viscous liquids to predict the decay of axisymmetric modes in liquid-filled tubes.

The common aspect of the above research is that the modal dispersion characteristics were obtained by the so-called “root-finding method”. As for a multi-layer wave guide, the traditional ways of constructing system matrices are the transfer matrix [17,18] technique and the global matrix [19] technique. When the thickness of the multi-layer model is large or the solving frequencies are high, the transfer matrix technique is prone to numerical instability [20], and although the global matrix technique does not have this problem, its matrix size becomes larger as the number of model layers increases, and the solving process of the characteristic equation becomes extremely slow.

The spectral method has been widely used in fluid dynamics since the 1970s [21,22], but it was not until 2004 that it was first introduced to acoustic engineering by Adamou [23] for solving the dispersion equation of modes in elastic waveguide structures. Research has shown that the spectral method is easier to implement in a numerical program than the traditional root-finding method, provides faster calculation speeds, and is more suitable for dealing with curving, damping, and inhomogeneous and anisotropic model problems. Using the spectral method, Karpfinger et al. [24,25] developed a set of algorithms for solving the dispersion characteristics of axisymmetric modes in multi-layer cylindrical structures.

To date, there has been no research on the dispersion characteristics of non-axisymmetric modes in fluid-filled elastic tubes using the spectral method. In this paper, a new algorithm based on the spectral method is proposed, which obtains the dispersion curves of the axial propagation modes in fluid-filled elastic tubes. It can also calculate the displacement distribution and stress distribution of the axial propagation modes in the radial direction. The algorithm proposed here works well in spite of tube dimensions, wall thickness, layers and frequency ranges. The results of the algorithm are verified using a COMSOL simulation. We found that there existed some modes omission of the dispersion curves presented in the literature [26].

2. Theory

2.1. Wave Equations

In this section, using linear elastic theory, the governing wave propagation equations in the elastic tube wall and the fluid are derived, respectively. The corresponding boundary conditions at different interfaces are given.

Figure 1 shows an infinite, homogeneous, and isotropic fluid-filled elastic tube model. The inner radius of the tube is $r_1$, the outer radius is $r_2$, and the cylindrical coordinate system $\left(r,\theta ,z\right)$ is established with the axis of the cylindrical tube as the z-axis. Waves propagating along the axial direction in a fluid-filled elastic tube can be divided into longitudinal waves and torsional waves, both of which have axisymmetric modes and non-axisymmetric modes [27]. In some references, longitudinal and torsional waves are defined as waves with only axisymmetric modes, and waves with non-axisymmetric modes are collectively called flexural waves. In this paper, the positive integer $n$ is used to represent the circumferential order of the mode, and the axisymmetric mode satisfies $n=0$.

For the elastic tube model in Figure 1, assuming that the outer surface satisfies the stress-free boundary condition, the wave equation governing the propagation of elastic waves [27] is

$$\mu {\nabla }^2\overrightarrow{U}+(\lambda +\mu )\nabla \nabla \cdot \overrightarrow{U}=\rho \frac{{\partial }^2\overrightarrow{U}}{\partial t^2}$$

(1)

where $\overrightarrow{U}$ represents the displacement field, which is a function of the three cylindrical coordinates and time. The bulk wave velocities in the elastic material are determined by the density $\rho$ and the Lamé constants $\lambda$ and $\mu$. If $c_1$ represents the dilatational longitudinal bulk wave velocity and $c_2$ is the shear bulk wave velocity, then

$$c_1=\sqrt{\frac{\lambda +2\mu }{\rho }}$$

(2)

$$c_2=\sqrt{\frac{\mu }{\rho }}$$

(3)

The displacement field $\overrightarrow{U}$ can be described using the dilatational scalar potential $\Phi$ and the equivoluminal vector potential $\overrightarrow{\Psi }$:

$$\overrightarrow{U}=\nabla \Phi +\nabla \times \overrightarrow{\Psi }$$

(4)

Substituting Equation (4) into Equation (1) gives

$${\nabla }^2\Phi =\frac{1}{c_1^2}\frac{{\partial }^2\Phi }{\partial t^2}$$

(5)

$${\nabla }^2\overrightarrow{\Psi }=\frac{1}{c_2^2}\frac{{\partial }^2\overrightarrow{\Psi }}{\partial t^2}$$

(6)

Using elasticity theory, the potentials $\Phi$ and $\overrightarrow{\Psi }$ have the following form in cylindrical coordinates:

$${\nabla }^2\Phi =\left(\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\theta }^2}+\frac{{\partial }^2}{\partial z^2}\right)\Phi$$

(7)

$$\begin{array}{ll}{\nabla }^2\overrightarrow{\Psi }& ={\nabla }^2\left({\phi }_r{\overrightarrow{e}}_r+{\phi }_{\theta }{\overrightarrow{e}}_{\theta }+{\phi }_z{\overrightarrow{e}}_z\right)\\ & =\left({\nabla }^2{\phi }_r-\frac{1}{r^2}{\phi }_r-\frac{2}{r^2}\frac{\partial {\phi }_{\theta }}{\partial \theta }\right){\overrightarrow{e}}_r+\left({\nabla }^2{\phi }_{\theta }-\frac{1}{r^2}{\phi }_{\theta }+\frac{2}{r^2}\frac{\partial {\phi }_r}{\partial \theta }\right){\overrightarrow{e}}_{\theta }+{\nabla }^2{\phi }_z{\overrightarrow{e}}_z\end{array}$$

(8)

Gazis [3] presented expressions for the potentials satisfying Equations (5) and (6):

$$\begin{array}{l}\Phi =\varphi \cos n\theta e^{i\left(k_{z}z-\omega t\right)}\\ {\phi }_z={\phi }_1\sin n\theta e^{i\left(k_{z}z-\omega t\right)}\\ {\phi }_r={\phi }_2\sin n\theta e^{i\left(k_{z}z-\omega t\right)}\\ {\phi }_{\theta }={\phi }_3\cos n\theta e^{i\left(k_{z}z-\omega t\right)}\end{array}$$

(9)

where the integer $n$ is known as the circumferential order of a wave mode and $k_z$ is the axial wavenumber. $\varphi$, ${\phi }_1$, ${\phi }_2$, and ${\phi }_3$ are unknown functions related only to $r$.

Let

$$\begin{array}{l}{h}_1={\phi }_2+{\phi }_3\\ h_2={\phi }_2-{\phi }_3\\ h_3={\phi }_1\end{array}$$

(10)

Substituting Equation (9) into Equations (5) and (6), and dropping $e^{i\left(k_{z}z-\omega t\right)}$, after a series of mathematical transformations, the equation of a wave propagating in an infinitely long, homogeneous, and isotropic elastic tube wall can be obtained as

$$\begin{array}{l}\left(\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{n^2}{r^2}+\frac{{\omega }^2}{c_1^2}\right)\varphi =k_z^2\varphi \\ \left(\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{{\left(n-1\right)}^2}{r^2}+\frac{{\omega }^2}{c_2^2}\right)h_1=k_z^2{h}_1\\ \left(\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{{\left(n+1\right)}^2}{r^2}+\frac{{\omega }^2}{c_2^2}\right)h_2=k_z^2{h}_2\\ \left(\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{n^2}{r^2}+\frac{{\omega }^2}{c_2^2}\right)h_3=k_z^2{h}_3\end{array}$$

(11)

The wave equation in the fluid is

$${\nabla }^2{\Phi }_{fl}=\frac{1}{c_{fl}^2}\frac{{\partial }^2{\Phi }_{fl}}{\partial t^2}$$

(12)

where the subscript $fl$ represents the fluid and $c_{fl}$ is the sound velocity of the fluid in the tube.

The displacement potential function ${\Phi }_{fl}$ is

$${\Phi }_{fl}={\varphi }_{fl}\cos n\theta e^{i\left(k_{z}z-\omega t\right)}$$

(13)

Substituting Equation (13) into Equation (12), the equation of an acoustic wave propagating in the fluid is obtained as

$$\left(\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{n^2}{r^2}+\frac{{\omega }^2}{c_{fl}^2}\right){\varphi }_{fl}=k_z^2{\varphi }_{fl}$$

(14)

2.2. Boundary Conditions

In order to obtain the dispersion equation of axial propagation modes in a fluid-filled elastic tube, it is necessary to set reasonable boundary conditions.

If the elastic tube wall is composed of multi-layer materials, at the solid 1–solid 2 interface, the continuous condition to be satisfied is that the displacement components $u_r$, $u_{\theta }$, $u_z$ and the surface stress components ${\tau }_{rr}$, ${\tau }_{r\theta }$, ${\tau }_{rz}$ in the solid 1 are equal to the corresponding components in the solid 2, as shown in Equation (15). Subscripts $s1$ and $s2$ represent the solid 1 and the solid 2, respectively, and $r=r_i$ represents the interface.

$$\begin{array}{l}{\left[u_{rs1},u_{\theta s1},u_{zs1}\right]}_{r=r_i}={\left[u_{rs2},u_{\theta s2},u_{zs2}\right]}_{r=r_i}\\ {\left[{\tau }_{rrs1},{\tau }_{r\theta s1},{\tau }_{rzs1}\right]}_{r=r_i}={\left[{\tau }_{rrs2},{\tau }_{r\theta s2},{\tau }_{rzs2}\right]}_{r=r_i}\end{array}$$

(15)

For the solid–fluid interface, if the viscosity of the fluid cannot be ignored, the continuous condition at the interface is that the displacement components and the surface stress components in the solid are equal to the corresponding components in the viscous fluid, as shown in Equation (16). Subscripts $s$ and $f$ represent the solid and the viscous fluid, respectively.

$$\begin{array}{l}{\left[u_{rs},u_{\theta s},u_{zs}\right]}_{r=r_i}={\left[u_{rf},u_{\theta f},u_{zf}\right]}_{r=r_i}\\ {\left[{\tau }_{rrs},{\tau }_{r\theta s},{\tau }_{rzs}\right]}_{r=r_i}={\left[{\tau }_{rrf},{\tau }_{r\theta f},{\tau }_{rzf}\right]}_{r=r_i}\end{array}$$

(16)

If the viscosity of the fluid is negligible, the boundary condition at the interface is a perfect slip boundary condition. That is, the radial component of the displacement of the solid and the fluid are equal, the circumferential component and the axial component of the displacement are discontinuous, and the surface stress components are still equal. The shear stress of the non-viscous fluid does not exist, that is, ${\tau }_{r\theta f}={\tau }_{rzf}=0$, as shown in Equation (17):

$$\begin{array}{l}{\left[u_{rs}\right]}_{r=r_i}={\left[u_{rf}\right]}_{r=r_i}\\ {\left[{\tau }_{rrs},{\tau }_{r\theta s},{\tau }_{rzs}\right]}_{r=r_i}={\left[{\tau }_{rrf},0,0\right]}_{r=r_i}\end{array}$$

(17)

The outer tube wall satisfies the stress-free boundary condition, as shown in Equation (18). The subscript $s$ represents the tube wall and $r=r_2$ represents the outer tube wall.

$${\left[{\tau }_{rrs},{\tau }_{r\theta s},{\tau }_{rzs}\right]}_{r=r_2}={\left[0,0,0\right]}_{r=r_2}$$

(18)

In addition, for an infinitely long tube, the inner and outer wall of the tube should meet:

$$\nabla \cdot \overrightarrow{\Psi }=0$$

(19)

The above Equations (15)–(19) constitute the boundary conditions that a fluid-filled elastic tube needs to be satisfied under different conditions. Expressions for these physical quantities in the above boundary conditions are given below.

2.2.1. Physical Quantities in the Tube Wall

From elastic theory, the displacement vector of the elastic tube wall is $\overrightarrow{U}=u_r{\overrightarrow{e}}_r+u_{\theta }{\overrightarrow{e}}_{\theta }+u_z{\overrightarrow{e}}_z$ [27], where

$$u_r=\frac{\partial }{\partial r}\Phi +\frac{1}{r}\frac{\partial }{\partial \theta }{\phi }_z-\frac{\partial }{\partial z}{\phi }_{\theta }$$

(20)

$$u_{\theta }=\frac{1}{r}\frac{\partial }{\partial \theta }\Phi +\frac{\partial }{\partial z}{\phi }_r-\frac{\partial }{\partial r}{\phi }_z$$

(21)

$$u_z=\frac{\partial }{\partial z}\Phi +\frac{1}{r}\frac{\partial }{\partial r}\left(r{\phi }_{\theta }\right)-\frac{1}{r}\frac{\partial }{\partial \theta }{\phi }_r$$

(22)

Substituting Equation (9) into Equations (20)–(22) and using $h_1,h_2,h_3$ to represent ${\phi }_1,{\phi }_2,{\phi }_3$ accruing to Equation (10), yields the displacement field:

$$u_r=\left(\frac{\partial }{\partial r}\varphi -\frac{1}{2}{\widehat{h}}_1+\frac{1}{2}{\widehat{h}}_2+\frac{n}{r}{h}_3\right)\cos n\theta \cdot e^{i\left(k_{z}z-\omega t\right)}$$

(23)

$$u_{\theta }=\left(-\frac{n}{r}\varphi +\frac{1}{2}{\widehat{h}}_1+\frac{1}{2}{\widehat{h}}_2-\frac{\partial }{\partial r}{h}_3\right)\sin n\theta \cdot e^{i\left(k_{z}z-\omega t\right)}$$

(24)

$$\begin{array}{ll}{\widehat{u}}_z& =i{k}_z{u}_z\\ & =\left[-k_z^2\varphi +\frac{1}{2}\left(\frac{\partial }{\partial r}+\frac{1-n}{r}\right){\widehat{h}}_1-\frac{1}{2}\left(\frac{\partial }{\partial r}+\frac{1+n}{r}\right){\widehat{h}}_2\right]\cos n\theta \cdot e^{i\left(k_{z}z-\omega t\right)}\end{array}$$

(25)

where

$$\begin{array}{l}{\widehat{h}}_1=i{k}_z{h}_1\\ {\widehat{h}}_2=i{k}_z{h}_2\end{array}$$

(26)

The stress tensor components in the elastic tube wall can be obtained using the following equations [27]:

$${\tau }_{rr}=\lambda \left(\frac{\partial u_r}{\partial r}+\frac{u_r}{r}+\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{\partial u_z}{\partial z}\right)+2\mu \frac{\partial u_r}{\partial r}$$

(27)

$${\tau }_{r\theta }=\mu \left(\frac{1}{r}\frac{\partial u_r}{\partial \theta }+\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}\right)$$

(28)

$${\tau }_{rz}=\mu \left(\frac{\partial u_r}{\partial z}+\frac{\partial u_z}{\partial r}\right)$$

(29)

Substituting Equations (23)–(25) into equations for ${\tau }_{rr}$, ${\tau }_{r\theta }$, and ${\tau }_{rz}$, one obtains the stress field:

$$\begin{array}{ll}{\tau }_{rr}& =\{[\lambda (\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{n^2}{r^2}-k_z^2)+2\mu \frac{{\partial }^2}{\partial r^2}]\varphi -\mu \frac{\partial }{\partial r}{\widehat{h}}_1\\ & +\mu \frac{\partial }{\partial r}{\widehat{h}}_2+2\mu (-\frac{n}{r^2}+\frac{n}{r}\frac{\partial }{\partial r})h_3\}\cos n\theta \cdot e^{i\left(k_{z}z-\omega t\right)}\end{array}$$

(30)

$$\begin{array}{ll}{\tau }_{r\theta }& =[2\mu (\frac{n}{r^2}-\frac{n}{r}\frac{\partial }{\partial r})\varphi +\frac{1}{2}\mu (\frac{\partial }{\partial r}+\frac{n-1}{r}){\widehat{h}}_1+\frac{1}{2}\mu (\frac{\partial }{\partial r}-\frac{n+1}{r}){\widehat{h}}_2\\ & +\mu (-\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{n^2}{r^2})h_3]\sin n\theta \cdot e^{i\left(k_{z}z-\omega t\right)}\end{array}$$

(31)

$$\begin{array}{ll}{\widehat{\tau }}_{rz}& =i{k}_z{\tau }_{rz}\\ & =[-2\mu k_z^2\frac{\partial }{\partial r}\varphi +\frac{\mu }{2}(\frac{{\partial }^2}{\partial r^2}+\frac{1-n}{r}\frac{\partial }{\partial r}+\frac{n-1}{r^2}+k_z^2){\widehat{h}}_1\\ & +\frac{\mu }{2}(-\frac{{\partial }^2}{\partial r^2}-\frac{n+1}{r}\frac{\partial }{\partial r}+\frac{n+1}{r^2}-k_z^2){\widehat{h}}_2-\mu k_z^2\frac{n}{r}{h}_3]\cos n\theta \cdot e^{i\left(k_{z}z-\omega t\right)}\end{array}$$

(32)

In addition, in the cylindrical coordinate system,

$$\nabla \cdot \overrightarrow{\Psi }=\frac{1}{r}\frac{\partial \left(r{\phi }_r\right)}{\partial r}+\frac{1}{r}\frac{\partial {\phi }_{\theta }}{\partial \theta }+\frac{\partial {\phi }_z}{\partial z}$$

(33)

Using ${\widehat{h}}_1$, ${\widehat{h}}_2$, ${\widehat{h}}_3$ to represent $i{k}_z\nabla \cdot \overrightarrow{\Psi }$, one obtains

$$i{k}_z\nabla \cdot \overrightarrow{\Psi }=\left[\frac{1}{2}\left(\frac{\partial }{\partial r}+\frac{1-n}{r}\right){\widehat{h}}_1+\frac{1}{2}\left(\frac{\partial }{\partial r}+\frac{1+n}{r}\right){\widehat{h}}_2-k_z^2{h}_3\right]\sin n\theta \cdot e^{i\left(k_{z}z-\omega t\right)}$$

(34)

2.2.2. Physical Quantities in the Fluid

The displacement components of the non-viscous fluid are expressed as

$$u_{rf}=\frac{\partial {\Phi }_{fl}}{\partial r}=\frac{\partial }{\partial r}{\varphi }_{fl}\cos n\theta e^{i\left(k_{z}z-\omega t\right)}$$

(35)

$$u_{\theta f}=\frac{1}{r}\frac{\partial {\Phi }_{fl}}{\partial \theta }=-\frac{n}{r}{\varphi }_{fl}\sin n\theta e^{i\left(k_{z}z-\omega t\right)}$$

(36)

$${\widehat{u}}_{zf}=i{k}_z{u}_{zf}=i{k}_z\frac{\partial {\Phi }_{fl}}{\partial z}=-k_z^2{\varphi }_{fl}\cos n\theta e^{i\left(k_{z}z-\omega t\right)}$$

(37)

The stress tensor components of the inviscid fluid are

$${\tau }_{rrf}=-\lambda \frac{{\omega }^2}{c_{fl}^2}{\varphi }_{fl}\cos n\theta \cdot e^{i\left(k_{z}z-\omega t\right)}$$

(38)

$${\tau }_{r\theta f}={\tau }_{rzf}=0$$

(39)

The wave equations in Equation (11), Equation (14), and the boundary conditions in Equations (15)–(19) constitute an eigenvalue problem. The traditional root-finding method is to select the appropriate general solution form of the wave equations, and to substitute the general solution into the appropriate boundary conditions to obtain a set of uniform linear geometric equations. The determinant of the coefficient matrix is set equal to zero, and the so-called dispersion equation is obtained. Finally, the root of the dispersion equation is found using numerical methods, and the dispersion relationship $\omega \left(k_z\right)$ of modes is obtained.

The general solution of the wave equations in the cylindrical coordinate system contains various forms of Bessel functions, and it is often difficult to separate and determine the root of different modes. If factors such as fluid viscosity and solid loss need to be considered, the dispersion equation becomes more difficult to solve. The spectral method introduced in the next section can help effectively avoid these problems.

3. The Spectral Method Algorithm for a Fluid-Filled Elastic Tube

3.1. The Spectral Method Theory

The advantage of the spectral method is that it does not require calculation of the root of the Bessel functions, but instead uses the differential matrix to solve the eigenvalue problem. The core idea is to approximate unknown functions to be solved in the differential equation through global interpolation using high-order orthogonal polynomials. As for the eigenvalue problem of a fluid-filled elastic tube, the differentiation occurs on a bounded interval, $0\le r\le r_2$; the appropriate choice is to use Chebyshev polynomials to construct the differential matrix according to Boyd [28]. The choice of polynomials for different unknown function $f\left(x\right)$ is listed in Appendix A,

Table A1. The choice of basis functions for different unknown functions.

If	Basis Set Is
$f\left(x\right)$ is periodic	Fourier series
$f\left(x\right)$  $\mathrm{is} \mathrm{periodic} \& \mathrm{symmetric} \mathrm{about} x=0$	Fourier cosine
$f\left(x\right)$  $\mathrm{is} \mathrm{periodic} \& \mathrm{antisymmetric} \mathrm{about} x=0$	Fourier sine
$x\in \left[a,b\right]$$\& f\left(x\right)$ is non-periodic	Chebyshev polys, Legendre polys.

.

The differential matrix $D_x^{\left(n\right)}$ is defined as follows according to Adamou [23]:

$$\left(\begin{array}{c}{f}_1^{\left(n\right)}\\ f_2^{\left(n\right)}\\ ⋮\\ f_N^{\left(n\right)}\end{array}\right)\approx \underset{D_x^{\left(n\right)}}{\underbrace{\left(\begin{array}{cccc}{D}_{11}^{\left(n\right)}& D_{12}^{\left(n\right)}& \cdots & D_{1N}^{\left(n\right)}\\ D_{21}^{\left(n\right)}& \ddots & & ⋮\\ ⋮& & \ddots & ⋮\\ D_{N1}^{\left(n\right)}& \cdots & \cdots & D_{NN}^{\left(n\right)}\end{array}\right)}}\left(\begin{array}{c}{f}_1\\ f_2\\ ⋮\\ f_N\end{array}\right)$$

(40)

where $f_i$ is the value of the function $f(x)$ evaluated at the interpolation point $x_i$, $i=1,2,...N$, $f_i{}^{(n)}$ is the $n$th derivative of $f(x_i)$, and the $N\times N$ matrix $D_x^{\left(n\right)}$ is the differential matrix.

The MATLAB program chebdif [29] generates interpolation points on the interval $-1\le x\le 1$:

$$x_i=\cos \left(\frac{\left(i-1\right)\pi }{N-1}\right)$$

(41)

Here, interpolation points for $r$ are

$$r_i=\frac{(r_2-r_1)x_i+r_1+r_2}{2}$$

(42)

where $r_2$ denotes the outer radius, $r_1$ denotes the inner radius.

$D_x^{\left(n\right)}$ generated by chebdif [29] approximates $d^m/d{x}^m$, but here requires the approximation of $d^m/d{r}^m$, which is given by

$$D_r^{\left(n\right)}=\left(\frac{2}{r_2-r_1}\right)D_x^{\left(n\right)}$$

(43)

Wave equations and boundary conditions only deal with $r$ derivatives, so the spectral method algorithm in the next section uses $D^{\left(n\right)}$ to represent $D_r^{\left(n\right)}$ for simplicity.

Adamou [23] gave a heuristic argument about the theoretical accuracy of the spectral method. For $f\left(x\right)$ on the Chebshev interval $-1\le x\le 1$, it requires at least two interpolation points per wavelength to accurately resolve $f\left(x\right)$. If $N$ is large enough, the interpoint spacing is asymptotically $\frac{\pi }{N}$ according to Equation (41). Assuming the wavelength of $f\left(x\right)$ is $\lambda$, and it is roughly constant across the domain, then the condition for accurate resolution is

$$\frac{\pi }{N}\le \frac{\lambda }{2}$$

(44)

A more formal demonstration based on the convergence rates of Chebyshev series can be found in the study of Gottlieb and Orszag [30].

3.2. The Spectral Method Algorithm

Wave equations and the corresponding boundary conditions in Section 2 are expressed in terms of the differential matrix as follows.

The wave Equation (11) in the elastic tube wall can be written as

$$\underset{\mathrm{L}}{\underbrace{\left(\begin{array}{cccc}{\mathrm{L}}_1& 0& 0& 0\\ 0& {\mathrm{L}}_2& 0& 0\\ 0& 0& {\mathrm{L}}_3& 0\\ 0& 0& 0& {\mathrm{L}}_4\end{array}\right)}}\underset{\Theta }{\underbrace{\left(\begin{array}{c}\varphi \\ {\widehat{h}}_1\\ {\widehat{h}}_2\\ h_3\end{array}\right)}}=k_z^2\left(\begin{array}{c}\varphi \\ {\widehat{h}}_1\\ {\widehat{h}}_2\\ h_3\end{array}\right)$$

(45)

where $\mathrm{L}$ is a $4N\times 4N$ differential matrix and expressions for ${\mathrm{L}}_i\left(i=1,2,3,4\right)$ are

$${\mathrm{L}}_1=D^{\left(2\right)}+\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-n^2\cdot \mathrm{diag}\left(\frac{1}{r^2}\right)+\left(\frac{{\omega }^2}{c_1^2}\right)\mathrm{I}$$

(46)

$${\mathrm{L}}_2=D^{\left(2\right)}+\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-{\left(n-1\right)}^2\cdot \mathrm{diag}\left(\frac{1}{r^2}\right)+\left(\frac{{\omega }^2}{c_2^2}\right)\mathrm{I}$$

(47)

$${\mathrm{L}}_3=D^{\left(2\right)}+\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-{\left(n+1\right)}^2\cdot \mathrm{diag}\left(\frac{1}{r^2}\right)+\left(\frac{{\omega }^2}{c_2^2}\right)\mathrm{I}$$

(48)

$${\mathrm{L}}_4=D^{\left(2\right)}+\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-n^2\cdot \mathrm{diag}\left(\frac{1}{r^2}\right)+\left(\frac{{\omega }^2}{c_2^2}\right)\mathrm{I}$$

(49)

In Equations (46)–(49), $\mathrm{diag}$ denotes a diagonal matrix and $\mathrm{I}$ denotes an $N\times N$ unit matrix.

The wave Equation (14) in the fluid is expressed in the form of the differential matrix as

$${\mathrm{L}}_f=D^{\left(2\right)}+\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-n^2\mathrm{diag}\left(\frac{1}{r^2}\right)+\left(\frac{{\omega }^2}{c_{fl}^2}\right)\mathrm{I}$$

(50)

Thus, the differential matrix expression of equations for wave propagation in a fluid-filled elastic tube is as follows:

$$\underset{\mathrm{L}}{\underbrace{\left(\begin{array}{ccccc}{\mathrm{L}}_f& & & & \\ & {\mathrm{L}}_1& & & \\ & & {\mathrm{L}}_2& & \\ & & & {\mathrm{L}}_3& \\ & & & & {\mathrm{L}}_4\end{array}\right)}}\underset{\Theta }{\underbrace{\left(\begin{array}{c}{\varphi }_{fl}\\ \varphi \\ {\widehat{h}}_1\\ {\widehat{h}}_2\\ h_3\end{array}\right)}}=k_z^2\left(\begin{array}{c}{\varphi }_{fl}\\ \varphi \\ {\widehat{h}}_1\\ {\widehat{h}}_2\\ h_3\end{array}\right)$$

(51)

Similarly, the differential matrix expressions for the displacement vector components in the elastic pipe wall are

$$\left(\begin{array}{c}{u}_r\\ u_{\theta }\\ {\widehat{u}}_z\end{array}\right)=\left(\begin{array}{cccc}\mathrm{urs}1& \mathrm{urs}2& \mathrm{urs}3& \mathrm{urs}4\\ \mathrm{uts}1& \mathrm{uts}2& \mathrm{uts}3& \mathrm{uts}4\\ \mathrm{uzs}1& \mathrm{uzs}2& \mathrm{uzs}3& \mathrm{uzs}4\end{array}\right)\left(\begin{array}{c}\varphi \\ {\widehat{h}}_1\\ {\widehat{h}}_2\\ h_3\end{array}\right)$$

(52)

The corresponding coefficients are

$$\begin{array}{l}\mathrm{urs}1=D^{\left(1\right)}\\ \mathrm{urs}2=-\frac{1}{2}\mathrm{I}\\ \mathrm{urs}3=\frac{1}{2}\mathrm{I}\\ \mathrm{urs}4=n\cdot \mathrm{diag}\left(\frac{1}{r}\right)\end{array}$$

(53)

$$\begin{array}{l}\mathrm{uts}1=-n\cdot \mathrm{diag}\left(\frac{1}{r}\right)\\ \mathrm{uts}2=\frac{1}{2}\mathrm{I}\\ \mathrm{uts}3=\frac{1}{2}\mathrm{I}\\ \mathrm{uts}4=-D^{\left(1\right)}\end{array}$$

(54)

$$\begin{array}{l}\mathrm{uzs}1=-D^{\left(2\right)}-\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}+n^2\cdot \mathrm{diag}\left(\frac{1}{r^2}\right)-\left(\frac{{\omega }^2}{c_1^2}\right)\mathrm{I}\\ \mathrm{uzs}2=\frac{1}{2}\left[D^{\left(1\right)}+\left(1-n\right)\cdot \mathrm{diag}\left(\frac{1}{r}\right)\right]\\ \mathrm{uzs}3=-\frac{1}{2}\left[D^{\left(1\right)}+\left(1+n\right)\cdot \mathrm{diag}\left(\frac{1}{r}\right)\right]\\ \mathrm{uzs}4=0\cdot \mathrm{I}\end{array}$$

(55)

The differential matrix expressions for the stress tensor components in the elastic tube wall are

$$\left(\begin{array}{c}{\tau }_{rr}\\ {\tau }_{r\theta }\\ {\widehat{\tau }}_{rz}\end{array}\right)=\left(\begin{array}{cccc}\mathrm{rrs}1& \mathrm{rrs}2& \mathrm{rrs}3& \mathrm{rrs}4\\ \mathrm{rts}1& \mathrm{rts}2& \mathrm{rts}3& \mathrm{rts}4\\ \mathrm{rzs}1& \mathrm{rzs}2& \mathrm{rzs}3& \mathrm{rzs}4\end{array}\right)\left(\begin{array}{c}\varphi \\ {\widehat{h}}_1\\ {\widehat{h}}_2\\ h_3\end{array}\right)$$

(56)

The corresponding coefficients are

$$\begin{array}{l}\mathrm{rrs}1=-\lambda \frac{{\omega }^2}{c_1^2}\mathrm{I}+2\mu D^{\left(2\right)}\\ \mathrm{rrs}2=-\mu D^{\left(1\right)}\\ \mathrm{rrs}3=\mu D^{\left(1\right)}\\ \mathrm{rrs}4=2\mu n\left[\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-\mathrm{diag}\left(\frac{1}{r^2}\right)\right]\end{array}$$

(57)

$$\begin{array}{l}\mathrm{rts}1=-2\mu n\left[\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-\mathrm{diag}\left(\frac{1}{r^2}\right)\right]\\ \mathrm{rts}2=\frac{\mu }{2}\left[D^{\left(1\right)}+\left(n-1\right)\mathrm{diag}\left(\frac{1}{r}\right)\right]\\ \mathrm{rts}3=\frac{\mu }{2}\left[D^{\left(1\right)}-\left(n+1\right)\mathrm{diag}\left(\frac{1}{r}\right)\right]\\ \mathrm{rts}4=\mu \left[-D^{\left(2\right)}+\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-n^2\mathrm{diag}\left(\frac{1}{r^2}\right)\right]\end{array}$$

(58)

$$\begin{array}{l}\mathrm{rzs}1=-2\mu \left[\begin{array}{l}{D}^{\left(3\right)}+\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(2\right)}-\left(n^2+1\right)\mathrm{diag}\left(\frac{1}{r^2}\right)D^{\left(1\right)}\\ +2n^2\mathrm{diag}\left(\frac{1}{r^3}\right)+\frac{{\omega }^2}{c_1^2}{D}^{\left(1\right)}\end{array}\right]\\ \mathrm{rzs}2=\frac{\mu }{2}\left[D^{\left(2\right)}-\left(n-1\right)\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}+\left(n-1\right)\mathrm{diag}\left(\frac{1}{r^2}\right)+{\mathrm{L}}_2\right]\\ \mathrm{rzs}3=-\frac{\mu }{2}\left[D^{\left(2\right)}+\left(n+1\right)\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}-\left(n+1\right)\mathrm{diag}\left(\frac{1}{r^2}\right)+{\mathrm{L}}_3\right]\\ \mathrm{rzs}4=-\mu n\left[\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(2\right)}+\mathrm{diag}\left(\frac{1}{r^2}\right)D^{\left(1\right)}-n^2\mathrm{diag}\left(\frac{1}{r^3}\right)+\frac{{\omega }^2}{c_2^2}\mathrm{diag}\left(\frac{1}{r}\right)\right]\end{array}$$

(59)

The differential matrix expression of Equation (34) is

$$i{k}_z\nabla \cdot \overrightarrow{\Psi }=\left(\begin{array}{cccc}\mathrm{gs}1& \mathrm{gs}2& \mathrm{gs}3& \mathrm{gs}4\end{array}\right)\left(\begin{array}{c}\varphi \\ {\widehat{h}}_1\\ {\widehat{h}}_2\\ h_3\end{array}\right)$$

(60)

The coefficients can be given as

$$\begin{array}{l}\mathrm{gs}1=0\cdot \mathrm{I}\\ \mathrm{gs}2=\frac{1}{2}\left[D^{\left(1\right)}+\left(1-n\right)\mathrm{diag}\left(\frac{1}{r}\right)\right]\\ \mathrm{gs}3=\frac{1}{2}\left[D^{\left(1\right)}+\left(1+n\right)\mathrm{diag}\left(\frac{1}{r}\right)\right]\\ \mathrm{gs}4=-D^{\left(2\right)}-\mathrm{diag}\left(\frac{1}{r}\right)D^{\left(1\right)}+n^2\mathrm{diag}\left(\frac{1}{r^2}\right)-\left(\frac{{\omega }^2}{c_2^2}\right)\mathrm{I}\end{array}$$

(61)

For the fluid, the differential matrix expressions of the displacement vector components are as follows:

$$\left(\begin{array}{c}{u}_{rf}\\ u_{\theta f}\\ {\widehat{u}}_{zf}\end{array}\right)=\left(\begin{array}{c}\mathrm{ufd}1\\ \mathrm{ufd}2\\ \mathrm{ufd}3\end{array}\right){\varphi }_{fl}$$

(62)

$$\mathrm{ufd}1=D^{\left(1\right)}$$

(63)

$$\mathrm{ufd}2=-n\mathrm{diag}\left(\frac{1}{r}\right)$$

(64)

$$\mathrm{ufd}3=-{\mathrm{L}}_f$$

(65)

The differential matrix expression of the stress tensor component ${\tau }_{rrf}$ of the fluid is

$${\tau }_{rrf}=(\mathrm{rs}1){\varphi }_{fl}$$

(66)

$$\mathrm{rs}1=-\lambda \frac{{\omega }^2}{c_{fl}^2}\cdot \mathrm{I}$$

(67)

Taking the fluid-filled elastic tube with a single-layer wall as an example, Figure 2 shows the appropriate position of the boundary conditions in the $\mathrm{L}$ matrix.

This can be performed by introducing a matrix $Q$ on the right-hand side of Equation (51), as shown in Equation (68):

$$\tilde{L}\Theta =k_z^{2}Q\Theta$$

(68)

where $\Theta$ is a matrix consisting of potential functions:

$$\Theta ={\left(\begin{array}{ccccc}{\varphi }_{fl}& \varphi & {\widehat{h}}_1& {\widehat{h}}_2& h_3\end{array}\right)}^T$$

(69)

$Q$ is a $5N\times 5N$ matrix and can be defined as follows:

$$Q=\left(\begin{array}{ccccc}M& & & & \\ & M& & & \\ & & M& & \\ & & & M& \\ & & & & M\end{array}\right)$$

(70)

where $M$ is an $N\times N$ diagonal matrix with the following form:

$$M=\left(\begin{array}{ccccc}0& & & & \\ & 1& & & \\ & & \ddots & & \\ & & & 1& \\ & & & & 0\end{array}\right)$$

(71)

Equation (68) is a generalized eigenvalue problem that can be solved using the MATLAB routine $\mathrm{EIG}\left(\tilde{L},Q\right)$ [24].

By changing the value of the frequency $\omega$, a different dispersion relation $\left(\omega ,k_z\right)$ can be obtained from Equation (68). For a given relationship between $\left(\omega ,k_z\right)$, the dispersion curves of the modes can be drawn. At the same time, the potential functions in Equation (69) are substituted into Equation (52), Equation (56), Equation (62), and Equation (66), which are the displacement vector components and stress tensor components of the modes in the tube wall and in the fluid. One can then obtain the displacement and stress distribution of the axial propagation modes in the fluid-filled tube with frequency $\omega$ with respect to the radial direction.

The flowchart of the proposed spectral method algorithm is shown in Figure 3.

4. Discussion

In this section, the dispersion characteristics of axial propagation modes in a water-filled PMMA tube, whose parameters are from Baik [26], were calculated using the algorithm described in Section 3. The results of the algorithm were verified using a COMSOL simulation. In addition, differences between dispersion curves obtained from the algorithm and Baik [26] were analyzed. Finally, the dispersion curves of a water-filled multi-layer tube were given.

4.1. Comparison between the Spectral Method Algorithm and the COMSOL Simulation

The parameters of the water-filled PMMA tube adopted by Baik [26] are shown in

Table 1. Parameters of the water-filled PMMA tube [26].

$\mathit{b}$ (cm)	$\mathit{h}$ (cm)	${\mathit{v}}_{\mathit{p}}$ (m/s)	${\mathit{v}}_{\mathit{s}}$ (m/s)	${\mathit{c}}_1$ (m/s)	${\mathit{\rho }}_{\mathit{s}}$(kg/m3)	${\mathit{\rho }}_1$ (kg/m3)
4.445	0.5	2690	1340	1479	1190	1000

, where $b$ is the inner radius of the tube, $h$ is the wall thickness, $v_p$, $v_s$ are the longitudinal and shear sound speed of PMMA, respectively, $c_1$ is the sound speed in water, ${\rho }_s$ is the density of PMMA, and ${\rho }_1$ is the density of water.

The two-dimensional model of the water-filled PMMA tube with the above parameters was established in COMSOL, as shown in Figure 4.

The inner part was water, and its physical field was set to be a “Pressure Acoustics” field. The wave equation for time-harmonic acoustic waves in 2D geometry is Equation (72),

$$\nabla \cdot \left(-\frac{1}{{\rho }_c}\left(\nabla p_t-{\overrightarrow{q}}_d\right)\right)-\frac{k_{eq}^2{p}_t}{{\rho }_c}=Q_m$$

(72)

$$p=p\left(x,y\right)e^{-i{k}_{z}z}$$

(73)

Where $p_t$ means total pressure, which is the sum of $p$ and the possible background pressure field $p_b$, as shown in Equation (74):

$$p_t=p+p_b$$

(74)

$k_{eq}^2$ is shown as Equation (75), where $k_z$ is the out-of-plane wave number:

$$k_{eq}^2={\left(\frac{\omega }{c_c}\right)}^2-k_z^2$$

(75)

The subscript “c” in ${\rho }_c$ and $c_c$ means that they may be complex-valued. Lossy media, like porous materials or highly viscous fluids, can be modeled by using the complex-valued speed of sound and density.

In Equation (72), ${\overrightarrow{q}}_d$ is the dipole domain source, which corresponds to a domain force source on the right-hand side of the momentum equation, and $Q_m$ is the monopole domain source, which corresponds to a mass source on the right-hand side of the continuity equation.

In the mode analysis study of COMSOL, $-i{k}_z$ is used as the eigenvalue.

The outer part is the PMMA tube wall, and its physical field, is set to be a “Solid Mechanics” field. The governing equation used here is

$$-\rho {\omega }^2\overrightarrow{u}=\nabla \cdot \left(C:\nabla \cdot \overrightarrow{u}\right)$$

(76)

In Equation (76), $\overrightarrow{u}$ is the structural displacement vector, and $C$ is the elasticity tensor.

The outer tube wall is the free boundary, which satisfies the stress-free boundary condition, and the inner tube wall is set to be the “Acoustic-Structure Boundary”, which is used to couple a Pressure Acoustics model to any structural component. The coupling includes the fluid load on the structure and the structural acceleration as experienced by the fluid.

The condition on exterior boundaries is shown as Equations (77) and (78):

$$-\overrightarrow{n}\cdot \left(-\frac{1}{{\rho }_c}\left(\nabla p_t-{\overrightarrow{q}}_d\right)\right)=-\overrightarrow{n}\cdot u_{tt}$$

(77)

$$F_A=p_t\overrightarrow{n}$$

(78)

where $u_{tt}$ is the structural acceleration, $\overrightarrow{n}$ is the surface normal, $p_t$ is the total acoustic pressure and $F_A$ is the load (force per unit area) experienced by the structure.

The condition on interior boundaries is shown as follows:

$$-\overrightarrow{n}\cdot {\left(-\frac{1}{{\rho }_c}\left(\nabla p_t-{\overrightarrow{q}}_d\right)\right)}_{\mathrm{up}}=-\overrightarrow{n}\cdot u_{tt}$$

(79)

$$-\overrightarrow{n}\cdot {\left(-\frac{1}{{\rho }_c}\left(\nabla p_t-{\overrightarrow{q}}_d\right)\right)}_{\mathrm{down}}=-\overrightarrow{n}\cdot u_{tt}$$

(80)

$$F_A=\left(p_{t,\mathrm{down}}-p_{t,\mathrm{up}}\right)\overrightarrow{n}$$

(81)

such that the acoustic load is given by the pressure drop across the thin structure. The up and down subscripts refer to the two sides of the interior boundary.

In order to find out the influence of mesh grids on the COMSOL simulation results, five different meshes were used to calculate the mode wavenumber, as shown in

Table 2. COMSOL mesh file and their numbers of elements.

Mesh	Mesh 1	Mesh 2	Mesh 3	Mesh 4	Mesh 5
Element Number	2214	3146	4644	8698	27,842

.

The wave number $k_z$ of the axial propagation modes in the water-filled PMMA tube at $f=20$ kHz was calculated using COMSOL with five different mesh files, as shown in in Table 2. The differences in $k_z$ at each mode between mesh 1–4 and mesh 5 are plotted in Figure 5. The horizontal axis is the mode index, meaning that there exist 22 modes at $f=20$ kHz. The vertical axis is the absolute value of the difference between mesh 1–4 and mesh 5 at each mode. Figure 5 shows that the simulation result $k_z$ of each mode tends to be consistent with the increase in the number of grid elements. In order to achieve the most accurate calculation results, the following simulation used mesh 5.

The mode analysis in COMSOL can calculate the wave number and mode shape of a propagation mode at a specific frequency. Combined with sweep analysis, dispersion curves of the axial propagation mode of the water-filled PMMA tube can be obtained. The frequency range calculated here is 0–90 kHz and the frequency interval is 200 Hz.

Figure 6 shows the comparison of the dispersion curves obtained from the COMSOL simulation and the spectral method algorithm.

In Figure 6, the horizontal axis is the dimensionless parameter $k_{1}b$, as shown in Equation (82), and the vertical axis is the normalized phase velocity $c_p/c_1$. The vertical axis uses logarithmic coordinates.

$$k_{1}b=\frac{\omega }{c_1}b$$

(82)

The curves in Figure 6 represent the calculation results from the spectral method algorithm and the dots represent the COMSOL simulation results. It can be seen from Figure 6 that the dispersion curves of the axial propagation modes in the water-filled PMMA tube obtained from the spectral method algorithm were completely consistent with the COMSOL simulation results.

Table 3. The calculation time of the COMSOL simulation and the spectral method algorithm.

Calculation Methods	Calculation Time
COMSOL Multiphysics	Approximately 10 h
The spectral method algorithm	A few tens of seconds

shows the time consumed by the COMSOL Multiphysics 5.6 software and the spectral method algorithm implemented in MATLAB R2022b to obtain the dispersion curves in Figure 6. The computer used here to run the COMSOL simulation and the MATLAB program of the proposed algorithm was a Lenovo desktop workstation, which has 2 Intel^® Xeon^® Silver 4110 CPU, and 128 GB RAM. The spectral method algorithm not only provides an accurate calculation but requires much less time and CPU resources to do so than the COMSOL simulation software.

Table 4. $k_z$ calculated using the COMSOL simulation and the spectral method algorithm at $f=8$ kHz.

	$\mathbf{Axial} \mathbf{Wavenumber} {\mathit{k}}_{\mathit{z}}$
COMSOL Multiphysics	4.3629	21.9621	30.1711	35.6921	37.5116	84.0532
	107.1383	121.4523	130.5509	136.0542	138.9479	139.8424
Spectral Method	4.3678	21.9621	30.1706	35.6932	37.5116	84.0534
	107.1380	121.4523	130.5509	136.0542	138.9479	139.8424
Error	1.12 × 10−3	0	1.66 × 10−5	3.08 × 10−5	0	2.38 × 10−6
	2.80 × 10−6	0	0	0	0	0

compares the wavenumbers $k_z$ of the axial propagation modes in the water-filled PMMA tube at $f=8$ kHz provided by the COMSOL simulation and the spectral method algorithm. The error in the table was calculated using Equation (83), where $k_{z(s)}$ represents the axial wavenumber calculated using the spectral method and $k_{z(c)}$ represents the axial wavenumber calculated using the COMSOL simulation.

$$E=\frac{\left|k_{z(s)}-k_{z(c)}\right|}{k_{z(c)}}$$

(83)

It can be seen from Table 4 that both the COMSOL simulation and the spectral method algorithm can find 12 different modes at 8 kHz, and the axial wavenumber corresponding to the same mode is identical under two different calculation ways.

Figure 7 shows the radial distribution of the displacement components $u_r$, $u_{\theta }$, $u_z$ of some non-axisymmetric modes in the tube wall at $f=8$ kHz. The modes are indexed as $(n,m)$, where $n$ is the circumferential order. The transverse axis $r$ represents the radial component in cylindrical coordinates, and the longitudinal axis is the normalized value of each displacement component. It can be seen from the figure that the distribution of the displacement components of the non-axisymmetric modes ($n\ne 0$) in the tube wall obtained by the spectral method algorithm was consistent with the COMSOL simulation results.

In conclusion, the dispersion curves of axial propagation modes in the water-filled tube obtained using the spectral method algorithm are consistent with the results from the COMSOL simulation. In addition, the algorithm can find all the axial propagation modes at any frequency correctly, and obtain the precise mode displacement distribution in the tube wall. The spectral method algorithm consumes much less time and computation resources to acquire such results.

4.2. Differences between the Results from the Spectral Method Algorithm and Baik [26]

Figure 8 shows the phase velocity curves of the axisymmetric modes ($n=0$) and the non-axisymmetric modes ($n=1,2,3,...$) of a water-filled PMMA tube propagating along the axis, as obtained from the spectral method algorithm. Figure 9 shows the group velocity curves of the corresponding modes. In the graphs, the horizontal axis is the dimensionless parameter $k_{1}b$ and the vertical axis is the phase velocity or group velocity normalized by the sound velocity in water as $c_p/c_1$ or $c_g/c_1$. All figures use $n$ to represent the circumferential order of the modes. The different color of lines in Figure 8 and Figure 9 denote a different group number of the modes with the same $n$. Figure 8 and Figure 9 are compared with Figure 2 in Baik [26], which shows the dispersion curves of phase speed and group speed of modes for the water-filled PMMA tube obtained.

Table 5. Modes calculated using the spectral method versus the same modes from Baik [26].

Circumferential Order	Spectral Method Results	Results from Baik [26]
$n=0$	$(0,1)$  $, (0,3)$  $, (0,4)$  $, (0,5)$	$(0,0)$  $, (0,1)$  $, (0,2)$  $, (0,3)$
$n=1$	$(1,3)$  $, (1,4)$  $, (1,5)$	$(1,1)$  $, (1,2)$  $, (1,3)$
$n=2$	$(2,3)$  $, (2,4)$  $, (2,5)$	$(2,1)$  $, (2,2)$  $, (2,3)$
$n=3$	$(3,3)$  $, (3,4)$  $, (3,5)$	$(3,1)$  $, (3,2)$  $, (3,3)$

shows the modal correspondence between the dispersion curves obtained using the spectral method algorithm and results from Baik [26]. For example, the $(0,1)$ mode in Figure 8 corresponds to the $(0,0)$ mode in Figure 2 of Baik [26], et cetera. This work only analyzed the modes listed in Table 5. The following analysis was based on the mode indices in Figure 8.

For the axisymmetric modes ($n=0$), Baik [26] did not find the $(0,2)$ mode in Figure 8, obtained with the spectral method algorithm. For the non-axisymmetric modes ($n\ne 0$), Baik [26] lacked the modes $(1,1)$, $(1,2)$, $(2,1)$, $(2,2)$, $(3,1)$,$(3,2)$ in Figure 8. The accuracy of the spectral method algorithm was verified using the results from the COMSOL simulation in the previous section.

It can be seen from Figure 8 that the phase velocity of mode $(0,2)$ did not change with the frequency, and was equal to the bulk shear wave velocity of the PMMA tube. Its group velocity was also equal to the bulk shear wave velocity, as shown in Figure 9. This is consistent with the description of the dispersion characteristics of the torsional wave provided by Gazis [3].

The phase velocities of modes $(0,1)$, $(1,1)$, $(2,1)$, and $(3,1)$ in Figure 8 were lower than the sound velocity in water, and the longitudinal and shear bulk wave velocities of the PMMA tube. Figure 10 and Figure 11 show the radial distribution of the displacement components and the stress tensor components of the above four modes at $f=8$ kHz, obtained using the spectral method algorithm.

In Figure 10, the solid line represents $u_r$, the dotted line represents $u_{\theta }$, and the dashed line represents $u_z$. In Figure 11, the solid line represents ${\tau }_{rr}$, the dotted line represents ${\tau }_{r\theta }$, and the dashed line represents ${\tau }_{rz}$. It can be seen from Figure 10 that the maximum amplitude of the axial displacement component $u_z$ of these four modes in the water was larger than that in the tube wall, and the energy carried in the water was higher than that in the tube wall. These phenomena are similar to a Scholte wave propagating at the solid–fluid interface. Baik [26] only obtained the $(0,1)$ mode, which is labeled as the $(0,0)$ mode in Figure 10.

The modes $(1,2)$, $(2,2)$, and $(3,2)$ in Figure 8 are flexural waves propagating along the axial direction.

In conclusion, dispersion curves of the water-filled PMMA tube provided by Baik [26] missed the torsional wave mode $(0,2)$, Scholte wave modes $(1,1)$, $(2,1)$, and $(3,1)$, and flexural wave modes $(1,2)$, $(2,2)$, and $(3,2)$. In contrast, the spectral method algorithm can find all the axial propagation modes precisely without missing any low-order modes.

4.3. A Water-Filled Elastic Multi-Layer Tube Model

The calculation model used in the above has only one layer in the tube wall. Here, the spectral method algorithm was used to obtain the dispersion curves of a water-filled multi-layer elastic tube. The cross section of the tube is shown in Figure 12. The materials of Layer 1, Layer 2 and Layer 3 were PVC, iron and PVC, respectively.

The parameters of the tube are shown in

Table 6. The parameters of the water-filled multi-layer elastic tube.

Dimension	$b$ (cm)	$h_1$ (cm)	$h_2$ (cm)	$h_3$ (cm)
	50	5	5	5
PVC	$v_p$ (m/s)	$v_s$ (m/s)	${\rho }_s$ (kg/m3)
	2388	1060	1380
Iron	$v_p$ (m/s)	$v_s$ (m/s)	${\rho }_s$ (kg/m3)
	5893	3230	7800

, where $b$ is the inner radius of the tube, $h_1$, $h_2$,$h_3$ are thickness of Layer 1, Layer 2 and Layer 3, respectively, $v_p$, $v_s$ are the longitudinal and shear sound speed of the solid material, respectively, ${\rho }_s$ is the density of the solid material, $c_1=1479$ m/s is the sound speed in water, and ${\rho }_1=1000$ kg/m³ is the density of water.

The phase velocity curves of the water-filled multi-layer elastic tube are shown in Figure 13; the horizontal axis is the frequency (kHz) and the vertical axis is $c_p/c_1$. The frequency range in Figure 13 is 0–2 kHz. In Figure 13, the axial propagation modes are plotted in different subgraphs according to the circumferential order $n$. It shows that the spectral method algorithm worked well for the large-dimension, thick-wall, multi-layer tube model.

5. Conclusions

We have shown that the algorithm based on the spectral method can obtain the dispersion characteristics of non-axisymmetric modes in the fluid-filled tube accurately and efficiently. Compared to the commercial software COMSOL, it consumes far fewer resources to obtain the same results. The algorithm is suitable for dealing with a multi-layer tube model, and it has no demands on the dimension of the tube and wall thickness. In addition, the algorithm can find all the axial propagation modes, while some of the low-order modes have been found to be missing in the relevant literature.

The limitations of this paper are that the spectral method algorithm proposed here worked under the assumption that the tube wall material was homogeneous and isotropic. For future work, the algorithm should be adjusted to deal with damped, inhomogeneous and anisotropic problems. In addition, the experiment of wave propagation in a fluid-filled elastic tube should be carried out to validate the numerical algorithm.