Numerical analysis of linear and circular wedge acoustic waveguides

Abstract

This study employs Hamilton’s principle and a bi-dimensional finite element method (Bi-d FEM) to investigate the dispersion characteristics and resonant modes of linear and circular wedge acoustic waveguides. The ANSYS three-dimensional FEM is commonly used to analyze waveguides and calculate the resonance frequency of vibration modes. However, the accuracy of this approach is constrained by how finely the elements are divided, limiting its applicability to the analysis of lower-order resonant modes. Therefore, to extend the analysis to high-frequency and high-order resonance modes, the Bi-d FEM is used. Based on the separation of variables method, this approach discretizes the time-harmonic propagation factor and the cross-sectional vibration of guided waves into two-dimensional finite elements. The number of divided elements is reduced, resulting in considerably improved calculation accuracy. The Bi-d FEM results indicate that the antisymmetric wave mode of the linear wedge acoustic waveguide decreases with vertex angle and that the wave velocity is lower than the Rayleigh wave velocity, which is consistent with the results obtained using Lagasse’s empirical formula. In both linear and circular wedge acoustic waveguides, the energy of the guided wave is concentrated at the wedge tip; therefore, it exhibit the characteristic behavior of an antisymmetric flexural wave.

Introduction

Research on wedge acoustic waveguides began in the 1970s. Since Lagasse’s discovery of guided wave behavior at the wedge tip in 1972, ¹ scholars have paid increasing attention to guided wave propagation in wedge bodies. Several studies have revealed that the majority of the energy in antisymmetric flexural (ASF) waves propagating along the tip of the wedge is localized at the apex itself. Moreover, the propagation velocity of the waves is slower than the Rayleigh wave velocity of the material, and they exhibit no dispersion characteristics. To describe the physical behavior of the propagation of guided waves in wedge structures, several studies have established mathematical models that accord well with experiment. However, owing to the complex geometric boundary value problems associated with the propagation behavior of wedge guided waves, a precise analytical solution remains elusive. Therefore, only approximate solutions are obtainable through numerical analysis.

In 1972, Lagasse ¹ and Maradudin et al. ² employed numerical calculations to verify the waveguide effect of a wedge with a vertex angle within a specific angular range. They demonstrated that the antisymmetric guided wave propagating along the tip of the wedge has most of its energy confined to a region within approximately one wavelength from the tip of the wedge. This guided wave does not exhibit dispersion when the tip of the wedge is not truncated.

In 1973, Lagasse et al. ³ proposed a simple empirical formula to summarize the numerical results. For integer multiples of the vertex angle less than or equal to 90° (i.e. $n\theta \le {90}^{°}$ ), the wave propagation velocity of each antisymmetric guided wave can be expressed as $V=V_R\sin (n\theta )$ , where $V_R$ represents the Rayleigh wave velocity, $\theta$ represents the vertex angle of the wedge, and $n$ represents the mode of the guided wave.

In 1973, Lagasse ⁴ employed the finite element method with the second-order triangular interpolation function to divide acoustic waveguides of different cross-sectional shapes. This approach investigated the wave propagation behavior of infinitely long linear wedge waveguides and analyzed the dispersion phenomena of guided waves in different geometric cross-sections. This study combined the variational principle and finite element analysis to develop a powerful method for studying wave propagation in acoustic waveguides.

In 1976, Mckenna et al. ⁵ approximated an infinitely long wedge with a small vertex angle as a thin plate. Using plate and shell theory, they derived an expression for the relationship between the phase velocity and the wavenumber of the guided wave for both truncated and untruncated wedges. In 1996, Hladky-Hennion ⁶ conducted numerical analyses, including the FEM to investigate the propagation behavior of linear and circular wedge guided waves. The effects of curvature radius, vertex angle, and vertex angle truncation on the phase velocity of wedge guided waves were studied.

In 1994, Krylov ⁷ adopted geometrical acoustics theory to calculate the propagation velocity of guided waves in a wedge in water. Their results indicated that the wave velocity decreases under the influence of water coupling and that this velocity reduction has a non-linear relationship with the vertex angle of the wedge. Moreover, it was predicted that a propeller exhibiting a motion pattern similar to that of a sting ray could be developed using this low-velocity wedge guided wave. In 1999, Krylov ⁸ derived a simpler and faster analytical solution on the basis of geometric acoustics theory to describe the wave propagation behavior of conical wedges and conical wedge-like guided waves, demonstrating good agreement with the finite element analysis results of Hladky-Hennion. ⁶ In 2007, Krylov et al.^9,10 experimentally validated the concept of wave-like aquatic propulsion for autonomous and crewed vessels; this concept was proposed a decade earlier. The first working prototype of a small catamaran using wave-like propulsion through an attached rubber keel was built and tested at Loughborough University.

In 1992 and 1993, Jia et al.^11,12 applied laser ultrasonic technology to measure the propagation characteristics of wedge guided waves. Employing Nd: YAG pulse laser as the excitation source and an optical heterodyne interferometer for signal reception, they analyzed measured waveform signals to verify the empirical formula proposed by Lagasse. ³ They found that the number of anti-symmetric guided wave modes on the wedge decreased with the vertex angle. In 2016 and 2017, Jia et al.^13,14 investigated the influence of truncation on the dispersion of wedge guided waves propagating along circular wedges using the laser ultrasound technique. A comparison of the dispersion curves revealed that the truncation of the vertex angle led to a dispersion effect in the guided wave, which became more pronounced at larger truncations.

In 1996, de Billy ¹⁵ used a contact transducer on each of the two sections of a linear wedge, an emitter, and a receiver to measure the velocity of guided waves. The results demonstrated the importance of minimizing the contact area between the emitter and sensor on the wedge and aligning the polarization direction of the emitter with the direction of the ASF guided wave vibration to obtain accurate measurements.

In 1998, Yang and Liaw ¹⁶ combined three-dimensional finite element analysis and laser ultrasonic experimental techniques to analyze and measure linear and disk-type wedges. Using a two-dimensional fast Fourier transform (2D-FFT), they mapped line scan wave propagation signal from the time domain to the frequency and wavenumber domains to obtain the phase velocity dispersion curve of the guided wave. This study demonstrated that the propagation behavior of disk-type wedge guided waves is influenced by curvature, resulting in dispersion. Moreover, the phase velocity dispersion decreases with frequency.

In 2003, Ibanescu et al. ¹⁷ analyzed the electromagnetic mode structure of an OmniGuide fiber: a hollow dielectric waveguide in which light is confined by a large index-contrast omnidirectional dielectric mirror. In 2006, Yang and Tsen ¹⁸ investigated the dispersion behavior of ASF modes propagating along the tip of a bilinear wedge using laser ultrasonic measurements and finite element numerical simulations.

In 2011, Mayer et al. ¹⁹ concluded that surface and interface acoustic waves are two-dimensional guided waves because their displacement field exhibit plane-wave-like behavior in spatial coordinates parallel to the guiding plane and decay exponentially along the axis normal to that plane. In the same year, Yu et al. ²⁰ developed a modal separation model to enhance the performance of wedge-wave ultrasonic motors using the Bi-d FEM. In 2015, Yu ²¹ investigated the propagation of flexural waves along the outer edge of a circular wedge, considering phase velocities and corresponding mode displacements. Furthermore, in 2021, Yu ²² investigated the dispersive properties and resonant modes of ridge waves that travel circumferentially around piezoelectric circular ridge waveguides.

In 2013, Sokolova et al. ²³ identified conditions under which acoustic waves guided at the apex of an ideal infinite elastic wedge exhibit weak dispersion. These conditions include (i) the coating of one or both surfaces of the infinite wedge, (ii) the truncation of the wedge apex or the use of a different material for the tip, and (iii) slight changes to material constants in an extended spatial region near the wedge tip.

In 2014, Dundu ²⁴ studied the behavior of circular hollow members with flattened edges in compression. Variables in the experiments included the diameter, thickness, and length of the sections, and the number of bolts in the connection. In 2019, Serey et al. ²⁵ introduced a general methodology for mode-selective generation in finite-cross-section waveguides using multiple transducers.

In 2021, Wang et al. ²⁶ explored to explore the effect of the micro-structure interaction, the distribution of physical properties, and the structure size on the thermo-mechanical response of a typical functionally graded hollow cylindrical structure whose inner boundary is subjected to a sudden temperature rise. They obtained and validated layered-formed solutions.

In 2022, Spada et al. ²⁷ presented a three-dimensional global–local formulation for predicting guided wave scattering patterns induced by discontinuities (e.g. defects). In 2023, Niu et al. ²⁸ explored techniques to minimize the number of incident/scattered directions required for accurate simulation of the scattering matrix of an arbitrary defect. In 2024, Wei et al. ²⁹ investigated the deformation of square plates with circular holes subjected to blast loads through experiment and simulation, and conducted explosion experiments and 2D to 3D mapping simulations.

This study proceeded in two stages. First, it used ANSYS for the modal analysis of a three-dimensional circular wedge model, simulated the wave propagation behavior of the circular wedge, established the dispersion curves, and investigated the convergence of finite elements. Second, it used the Bi-d FEM to calculate the dispersion curves and modes of linear and circular wedge acoustic waveguides with varying wedge angles and examined the influence of boundary conditions and curvature on the propagation of wedge guided waves.

Basic theory

Relationship between stress, strain, and displacement

For a linear wedge acoustic waveguide, in terms of the Cartesian coordinate system shown in Figure 1(a), the displacement vector $\mathbf{u}$ , strain vector $\mathbf{S}$ , and stress vector $\mathbf{T}$ are expressed as follows ³⁰ :

$$\mathbf{u}=[\begin{array}{ccc}{u}_x\hfill & u_y\hfill & u_z\hfill \end{array}{]}^T$$

(1)

$$\mathbf{S}=[\begin{array}{cccccc}{\epsilon }_{\mathrm{xx}}\hfill & {\epsilon }_{\mathrm{yy}}\hfill & {\epsilon }_{\mathrm{zz}}\hfill & 2{\epsilon }_{\mathrm{yz}}\hfill & 2{\epsilon }_{\mathrm{xz}}\hfill & 2{\epsilon }_{\mathrm{xy}}\hfill \end{array}{]}^T$$

(2)

$$\mathbf{T}=[\begin{array}{cccccc}{\sigma }_{\mathrm{xx}}\hfill & {\sigma }_{\mathrm{yy}}\hfill & {\sigma }_{\mathrm{zz}}\hfill & {\sigma }_{\mathrm{yz}}\hfill & {\sigma }_{\mathrm{xz}}\hfill & {\sigma }_{\mathrm{xy}}\hfill \end{array}{]}^T$$

(3)

Figure 1.

Linear and circular wedge acoustic waveguides: (a) linear wedge and (b) circular wedge.

If the acoustic waveguide is considered as an orthotropic material, the relationship between stress and strain can be obtained using the generalized Hooke’s law as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}{\sigma }_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\sigma }_{\mathrm{yy}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{yz}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xy}}}\end{array}\right\}=\left[\begin{array}{cccccc}{C}_{11}\hfill & C_{12}\hfill & C_{13}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{C_{12}}& C_{22}\hfill & C_{23}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{C_{13}}& C_{23}\hfill & C_{33}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & C_{44}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & C_{55}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & C_{66}\hfill \end{array}\right]\left\{\begin{array}{c}{\epsilon }_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{yy}}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{2{\epsilon }_{\mathrm{yz}}}\\ \multicolumn{1}{c}{2{\epsilon }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{2{\epsilon }_{\mathrm{xy}}}\end{array}\right\}}\end{array}$$

(4)

The relationship between strain and displacement is

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}{S}_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{S_{\mathrm{yy}}}\\ \multicolumn{1}{c}{S_{\mathrm{zz}}}\\ \multicolumn{1}{c}{S_{\mathrm{yz}}}\\ \multicolumn{1}{c}{S_{\mathrm{xz}}}\\ \multicolumn{1}{c}{S_{\mathrm{xy}}}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{yy}}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{2{\epsilon }_{\mathrm{yz}}}\\ \multicolumn{1}{c}{2{\epsilon }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{2{\epsilon }_{\mathrm{xy}}}\end{array}\right\}=\left\{\begin{array}{c}\partial u_x/\partial x\hfill \\ \multicolumn{1}{c}{\partial u_y/\partial y}\\ \multicolumn{1}{c}{\partial u_z/\partial z}\\ \multicolumn{1}{c}{\partial u_y/\partial z+\partial u_z/\partial y}\\ \multicolumn{1}{c}{\partial u_x/\partial z+\partial u_z/\partial x}\\ \multicolumn{1}{c}{\partial u_x/\partial y+\partial u_y/\partial x}\end{array}\right\}}\end{array}$$

(5)

For a circular wedge acoustic waveguide, in terms of the polar coordinate system shown in Figure 1(b), the displacement vector $\mathbf{u}$ , strain vector $\mathbf{S}$ , and stress vector $\mathbf{T}$ are expressed as follows:

$$\mathbf{u}=[\begin{array}{ccc}{u}_r\hfill & u_x\hfill & u_{\theta }\hfill \end{array}{]}^T$$

(6)

$$\mathbf{S}=[\begin{array}{cccccc}{\epsilon }_{\mathrm{rr}}\hfill & {\epsilon }_{\mathrm{xx}}\hfill & {\epsilon }_{\theta \theta }\hfill & 2{\epsilon }_{x\theta }\hfill & 2{\epsilon }_{r\theta }\hfill & 2{\epsilon }_{\mathrm{rx}}\hfill \end{array}{]}^T$$

(7)

$$\mathbf{T}=[\begin{array}{cccccc}{\sigma }_{\mathrm{rr}}\hfill & {\sigma }_{\mathrm{xx}}\hfill & {\sigma }_{\theta \theta }\hfill & {\sigma }_{x\theta }\hfill & {\sigma }_{r\theta }\hfill & {\sigma }_{\mathrm{rx}}\hfill \end{array}{]}^T$$

(8)

Similarly, if the acoustic waveguide is considered an orthotropic material, the relationship between stress and strain is expressed as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}{\sigma }_{\mathrm{rr}}\hfill \\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xx}}}\\ \multicolumn{1}{c}{{\sigma }_{\theta \theta }}\\ \multicolumn{1}{c}{{\sigma }_{x\theta }}\\ \multicolumn{1}{c}{{\sigma }_{r\theta }}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{rx}}}\end{array}\right\}=\left[\begin{array}{cccccc}{C}_{11}\hfill & C_{12}\hfill & C_{13}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{C_{12}}& C_{22}\hfill & C_{23}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{C_{13}}& C_{23}\hfill & C_{33}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & C_{44}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & C_{55}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & C_{66}\hfill \end{array}\right]\left\{\begin{array}{c}{\epsilon }_{\mathrm{rr}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{xx}}}\\ \multicolumn{1}{c}{{\epsilon }_{\theta \theta }}\\ \multicolumn{1}{c}{2{\epsilon }_{x\theta }}\\ \multicolumn{1}{c}{2{\epsilon }_{r\theta }}\\ \multicolumn{1}{c}{2{\epsilon }_{\mathrm{rx}}}\end{array}\right\}}\end{array}$$

(9)

The relationship between strain and displacement is:

$$\left\{\begin{array}{c}{S}_{\mathrm{rr}}\hfill \\ \multicolumn{1}{c}{S_{\mathrm{xx}}}\\ \multicolumn{1}{c}{S_{\theta \theta }}\\ \multicolumn{1}{c}{S_{x\theta }}\\ \multicolumn{1}{c}{S_{r\theta }}\\ \multicolumn{1}{c}{S_{\mathrm{rx}}}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_{\mathrm{rr}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{xx}}}\\ \multicolumn{1}{c}{{\epsilon }_{\theta \theta }}\\ \multicolumn{1}{c}{2{\epsilon }_{x\theta }}\\ \multicolumn{1}{c}{2{\epsilon }_{r\theta }}\\ \multicolumn{1}{c}{2{\epsilon }_{\mathrm{rx}}}\end{array}\right\}=\left\{\begin{array}{c}\partial u_r/\partial r\hfill \\ \multicolumn{1}{c}{\partial u_x/\partial x}\\ \multicolumn{1}{c}{u_r/r+(\partial u_{\theta }/\partial \theta )/r}\\ \multicolumn{1}{c}{(1/r)(\partial u_x/\partial \theta )+\partial u_{\theta }/\partial x}\\ \multicolumn{1}{c}{(1/r)(\partial u_r/\partial \theta )+\partial u_{\theta }/\partial r-u_{\theta }/r}\\ \multicolumn{1}{c}{\partial u_r/\partial x+\partial u_x/\partial r}\end{array}\right\}$$

(10)

Hamilton’s principle

If zero surface traction is assumed for an infinitely long linear wedge acoustic waveguide, the waveguide’s strain energy $\mathrm{U}$ and kinetic energy $T$ are expressed in matrix form as follows:

$$U=\frac{1}{2}{\int }_{\mathrm{vol}}{\mathit{T}}^H\mathit{S}\mathrm{dV}=\frac{1}{2}{\int }_z{\int }_x{\int }_y{\mathit{S}}^H\mathit{CS}\mathrm{dxdydz}$$

(11a)

$$T=\frac{1}{2}{\int }_{\mathrm{vol}}{\stackrel{·}{\mathit{u}}}^H\mathbf{\rho }\stackrel{·}{\mathit{u}}\mathrm{dV}=\frac{1}{2}{\int }_z{\int }_x{\int }_y{\stackrel{·}{\mathit{u}}}^H\mathbf{\rho }\stackrel{·}{\mathit{u}}\mathrm{dxdydz}$$

(11b)

Following the same approach, considering the case where the surface traction of a circular acoustic waveguide is zero, its strain energy $\mathrm{U}$ and kinetic energy $T$ are expressed in matrix form as follows:

$$U=\frac{1}{2}{\int }_{\mathrm{vol}}{\mathit{T}}^H\mathit{S}\mathrm{dV}=\frac{1}{2}{\int }_{\theta }{\int }_r{\int }_x{\mathit{S}}^H\mathit{CS}r\mathrm{drdxd}\theta$$

(12a)

$$T=\frac{1}{2}{\int }_{\mathrm{vol}}{\stackrel{·}{\mathit{u}}}^H\mathbf{\rho }\stackrel{·}{\mathit{u}}\mathrm{dV}=\frac{1}{2}{\int }_{\theta }{\int }_r{\int }_x{\stackrel{·}{\mathit{u}}}^H\mathbf{\rho }\stackrel{·}{\mathit{u}}r\mathrm{drdxd}\theta$$

(12b)

where $\mathbf{C}$ denotes the material stiffness matrix and $\mathbf{\rho }$ denotes a $3\times 3$ diagonal matrix, which represents the density of the elastic acoustic waveguide. $H$ denotes the Hermitian transpose. The dot denotes the partial derivative with respect to time.

If $\mathbf{C}$ and $\mathbf{\rho }$ are symmetric real matrices, the strain energy and kinetic energy are positive definite. Therefore, Lagrangian functions $L$ of linear and circular acoustic waveguide systems are expressed as follows:

$$L=T-U=\frac{1}{2}{\int }_z{\int }_x{\int }_y({\stackrel{·}{\mathit{u}}}^H\mathbf{\rho }\stackrel{·}{\mathit{u}}-{\mathit{S}}^H\mathit{C}S)\mathrm{dxdydz}$$

(13a)

$$L=T-U=\frac{1}{2}{\int }_{\theta }{\int }_r{\int }_x({\stackrel{·}{\mathit{u}}}^H\mathbf{\rho }\stackrel{·}{\mathit{u}}-{\mathit{S}}^H\mathit{CS})\mathrm{rdrdxd}\theta$$

(13b)

Applying Hamilton’s principle, the first variation of the Lagrangian time integral is set to zero:

$$\delta \underset{{\mathrm{t}}_1}{\overset{{\mathrm{t}}_2}{\int }}\mathrm{L}\mathrm{dt}=0$$

(14)

Subsequently, the equations of motion for linear and circular acoustic waveguides can be derived.

Bi-dimensional finite element analysis

In equation (13a), the separation of variables method can be used to extract the Lagrangian of a linear acoustic waveguide from the integral over z. Similarly, in equation (13b), the Lagrangian of a circular acoustic waveguide can be extracted from the integral along the circumferential $\theta$ direction.

The Bi-d FEM is then used to discretize the cross-sections of the linear and circular acoustic waveguides into small elements (Figure 2). The displacement of any point within an element is expressed as the product of the element node displacement and the interpolation function. The element used in this study is a two-dimensional four-node isoparametric element (Q4).

Figure 2.

Mesh map of linear and circular wedge acoustic waveguide: (a) linear wedge and (b) circular wedge.

A linear acoustic waveguide is considered in which the wave propagates along the positive z-direction, with the section perpendicular to the z-axis maintaining the resonant mode.

Therefore, referring to the coordinate system shown in Figure 1(a), the displacement vector $\mathbf{u}$ is assumed to be

$$\mathbf{u}=\left\{\begin{array}{c}{u}_x\hfill \\ \multicolumn{1}{c}{u_y}\\ \multicolumn{1}{c}{u_z}\end{array}\right\}=\left\{\begin{array}{c}{u}_x(x,y)\hfill \\ \multicolumn{1}{c}{u_y(x,y)}\\ \multicolumn{1}{c}{u_z(x,y)}\end{array}\right\}{e}^{\mathrm{jkz}}$$

(15)

where $e^{\mathrm{jkz}}$ is the wave propagation factor, $k$ is the wavenumber, and $j$ is the imaginary unit $\sqrt{-1}$ .

The displacement components of each element inside the acoustic waveguide are expressed as node displacements and interpolation functions.

$$u_x=(N_1{u}_x^1+N_2{u}_x^2+N_3{u}_x^3+N_4{u}_x^4)e^{\mathrm{jkz}}$$

(16a)

$$u_y=(N_1{u}_y^1+N_2{u}_y^2+N_3{u}_y^3+N_4{u}_y^4)e^{\mathrm{jkz}}$$

(16b)

$$u_z=j(N_1{u}_z^1+N_2{u}_z^2+N_3{u}_z^3+N_4{u}_z^4)e^{\mathrm{jkz}}$$

(16c)

where $u_x$ , $u_y$ , and $u_z$ represent real-valued functions; however, to satisfy the equation of motion in the frequency domain, $u_z$ should be multiplied by the complex number $j$ . In the physical sense, $u_x$ , $u_y$ , and $u_z$ are ${90}^{°}$ out of phase.

Equation (16) can be expressed as follows (in matrix form):

$$\mathbf{u}=\mathbf{Nd}$$

(17a)

$$\mathbf{u}=[\begin{array}{ccc}{u}_x\hfill & u_y\hfill & u_z\hfill \end{array}{]}^T$$

(17b)

$$\mathbf{N}=[\begin{array}{cccc}{\mathbf{N}}_1\hfill & {\mathbf{N}}_2\hfill & {\mathbf{N}}_3\hfill & {\mathbf{N}}_4\hfill \end{array}]$$

(17c)

$${\mathbf{N}}_i=\left[\begin{array}{ccc}{N}_i\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_i\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & j{N}_i\hfill \end{array}\right]e^{\mathrm{jkz}},N_i=\frac{1}{4}(1+\xi {\xi }_i)(1+\eta {\eta }_i)$$

(17d)

$$\mathbf{d}=[\begin{array}{cccc}{\mathbf{d}}_1\hfill & {\mathbf{d}}_2\hfill & {\mathbf{d}}_3\hfill & {\mathbf{d}}_4\hfill \end{array}{]}^T,{\mathbf{d}}_i=[\begin{array}{ccc}{u}_x^i\hfill & u_y^i\hfill & u_z^i\hfill \end{array}{]}^T$$

(17e)

where $\mathbf{N}$ denotes the interpolation function matrix, $\mathbf{d}$ denotes the node displacement vector, and $N_i$ denotes the interpolation function that maps real coordinates to natural coordinates ( $i=1,2,3,\mathrm{or}4$ ). Substituting equation (16) into the strain-displacement relationship of equation (5) yields the following expression for the strain vector $\mathbf{S}$ :

$$\mathbf{S}=\mathbf{Bd}$$

(18a)

$$\mathbf{B}=\left[\begin{array}{ccc}\begin{array}{cc}{\mathbf{B}}_1\hfill & {\mathbf{B}}_2\hfill \end{array}\hfill & {\mathbf{B}}_3\hfill & {\mathbf{B}}_4\hfill \end{array}\right]$$

(18b)

$${\mathbf{B}}_i=\left[\begin{array}{ccc}\partial N_i/\partial x\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \partial N_i/\partial y\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -k{N}_i\hfill \\ \multicolumn{1}{c}{0}& \mathrm{jk}{N}_i\hfill & j\partial N_i/\partial y)\hfill \\ \multicolumn{1}{c}{\mathrm{jk}{N}_i}& 0\hfill & j(\partial N_i/\partial x)\hfill \\ \multicolumn{1}{c}{\partial N_i/\partial y}& \partial N_i/\partial x\hfill & 0\hfill \end{array}\right]e^{\mathrm{jkz}}$$

(18c)

Considering that the guided wave of a circular acoustic waveguide propagates along the circumferential direction, it maintains a resonance pattern in the cross-section perpendicular to the circumferential direction. Therefore, referring to the coordinate system of Figure 1(b), the displacement vector is assumed to be

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{u}=\left\{\begin{array}{c}{u}_r\hfill \\ \multicolumn{1}{c}{u_x}\\ \multicolumn{1}{c}{u_{\theta }}\end{array}\right\}=\left\{\begin{array}{c}{u}_r(r,x)\hfill \\ \multicolumn{1}{c}{u_x(r,x)}\\ \multicolumn{1}{c}{u_{\theta }(r,x)}\end{array}\right\}{e}^{\mathrm{jn}\theta }}\end{array}$$

(19)

where $2\pi R=n\lambda =n(2\pi /k)$ , $n=\mathrm{kR}$ , and $n$ in the wave propagation factor $e^{\mathrm{jn}\theta }$ represents the number of circumferential waves propagating along the circumferential direction.

The displacement components of each element inside the acoustic waveguide are expressed as node displacements and interpolation functions.

$$u_r=(N_1{u}_r^1+N_2{u}_r^2+N_3{u}_r^3+N_4{u}_r^4)e^{\mathrm{jn}\theta }$$

(20a)

$$u_x=(N_1{u}_x^1+N_2{u}_x^2+N_3{u}_x^3+N_4{u}_x^4)e^{\mathrm{jn}\theta }$$

(20b)

$$u_{\theta }=j(N_1{u}_{\theta }^1+N_2{u}_{\theta }^2+N_3{u}_{\theta }^3+N_4{u}_{\theta }^4)e^{\mathrm{jn}\theta }$$

(20c)

where $u_{\theta }$ is multiplied by $j$ to satisfy the equation of motion in the frequency domain.

Equation (20) can be expressed in the following matrix form:

$$\mathbf{u}=\mathbf{N}\mathbf{d}$$

(21a)

where

$$\mathbf{u}=[\begin{array}{ccc}{u}_r\hfill & u_x\hfill & u_{\theta }\hfill \end{array}{]}^T$$

(21b)

$$\mathbf{N}=[\begin{array}{cccc}{\mathbf{N}}_1\hfill & {\mathbf{N}}_2\hfill & {\mathbf{N}}_3\hfill & {\mathbf{N}}_4\hfill \end{array}]$$

(21c)

$${\mathbf{N}}_i=\left[\begin{array}{ccc}{N}_i\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_i\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & j{N}_i\hfill \end{array}\right]e^{\mathrm{jn}\theta },N_i=\frac{1}{4}(1+\xi {\xi }_i)(1+\eta {\eta }_i)$$

(21d)

$$\mathbf{d}=[\begin{array}{cccc}{\mathbf{d}}_1\hfill & {\mathbf{d}}_2\hfill & {\mathbf{d}}_3\hfill & {\mathbf{d}}_4\hfill \end{array}{]}^T,{\mathbf{d}}_i=[\begin{array}{ccc}{u}_r^i\hfill & u_x^i\hfill & u_{\theta }^i\hfill \end{array}]$$

(21e)

Substituting equation (20) into equation (10) for the polar coordinate system and expressing the strain vector in a matrix form result in

$$\mathbf{S}=\mathbf{Bd}$$

(22a)

where

$$\mathbf{B}=[\begin{array}{ccc}\begin{array}{cc}{\mathbf{B}}_1\hfill & {\mathbf{B}}_2\hfill \end{array}\hfill & {\mathbf{B}}_3\hfill & {\mathbf{B}}_4\hfill \end{array}]$$

(22b)

$${\mathbf{B}}_i=\left[\begin{array}{ccc}\partial N_i/\partial r\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \partial N_i/\partial x\hfill & 0\hfill \\ \multicolumn{1}{c}{N_i/r}& 0\hfill & -n(N_i/r)\hfill \\ \multicolumn{1}{c}{0}& \mathrm{jn}(N_i/r)\hfill & j(\partial N_i/\partial x)\hfill \\ \multicolumn{1}{c}{\mathrm{jn}(N_i/r)}& 0\hfill & j(\partial N_i/\partial r-N_i/r)\hfill \\ \multicolumn{1}{c}{\partial N_i/\partial x}& \partial N_i/\partial r\hfill & 0\hfill \end{array}\right]e^{\mathrm{jn}\theta }$$

(22c)

Dispersion equations

Dispersion equation for a linear acoustic waveguide

By substituting the expressions of the displacement $\mathbf{u}$ , strain $\mathbf{S}$ , and interpolation function at any node of an element from equations (21a) and (22a) into equations (11) and (13a), the strain energy $U$ , kinetic energy $T$ , and Lagrangian $L$ of each element are obtained as follows:

$$U=\frac{1}{2}\sum {\int }_z{\int }_x{\int }_y{\mathit{d}}^H{\mathit{B}}^H\mathit{CBd}\mathrm{dxdydz}=\frac{1}{2}\sum {\int }_z{\mathit{d}}^H\mathit{kd}\mathrm{dz}$$

(23a)

$$T=\frac{1}{2}\sum {\int }_z{\int }_x{\int }_y{\stackrel{·}{\mathit{d}}}^H{\mathit{N}}^H\mathit{\rho }\mathit{N}\stackrel{·}{\mathit{d}}\mathrm{dxdydz}=\frac{1}{2}\sum {\int }_z{\stackrel{·}{\mathit{d}}}^H\mathit{m}\stackrel{·}{\mathit{d}}\mathrm{dz}$$

(23b)

$$L=T-U=\frac{1}{2}\sum {\int }_z({\stackrel{·}{\mathit{d}}}^H\mathit{m}\stackrel{·}{\mathit{d}}-{\mathit{d}}^H\mathit{kd})\mathrm{dz}$$

(23c)

where

$$\mathbf{m}=\underset{x}{\int }\underset{y}{\int }{\mathbf{N}}^H\mathit{\rho }\mathbf{N}\mathrm{dxdy}$$

(24a)

$$\mathbf{k}=\underset{x}{\int }\underset{y}{\int }{\mathbf{B}}^H\mathbf{CB}\mathrm{dxdy}$$

(24b)

Here, $\mathbf{m}$ and $\mathbf{k}$ are the mass and stiffness matrices of the element, respectively. By substituting equation (23) into equation (14), we find that

$$\delta \sum {\int }_{t_1}^{t_2}\{\frac{1}{2}{\int }_z({\stackrel{·}{\mathit{d}}}^H\mathit{m}\stackrel{·}{\mathit{d}}-{\mathit{d}}^H\mathit{kd})\mathrm{dz}\}\mathrm{dt}=0$$

(25)

By performing variational operations on equation (25), we obtain

$$\begin{array}{c}\frac{1}{2}\sum {\int }_z\hfill \\ \{{\int }_{t_1}^{t_2}\delta {\stackrel{·}{\mathit{d}}}^H\mathit{m}\stackrel{·}{\mathit{d}}\mathrm{dt}+{\int }_{t_1}^{t_2}{\stackrel{·}{\mathit{d}}}^H\mathit{m}\delta \stackrel{·}{\mathit{d}}\mathrm{dt}-{\int }_{t_1}^{t_2}\delta {\mathit{d}}^H\mathit{kd}\mathrm{dt}-{\int }_{t_1}^{t_2}{\mathit{d}}^H\mathit{k}\delta \mathit{d}\mathrm{dt}\}\mathrm{dz}=0\hfill \end{array}$$

(26)

Then, using integration by parts to expand equation (26) yields

$$\begin{array}{l}\frac{1}{2}\sum {\int }_z\{{\delta {\mathit{d}}^H\mathit{m}\stackrel{.}{\mathit{d}}|}_{t_1}^{t_2}-{\int }_{t_1}^{t_2}\delta {\mathit{d}}^H\mathit{m}\stackrel{··}{\mathit{d}}dt+{{\stackrel{.}{\mathit{d}}}^H\mathit{m}\delta \mathit{d}|}_{t_1}^{t_2}\\ -{\int }_{t_1}^{t_2}{\stackrel{··}{\mathit{d}}}^H\mathit{m}\delta \mathit{d}dt-{\int }_{t_1}^{t_2}\delta {\mathit{d}}^H\mathit{k}\mathit{d}dt-{\int }_{t_1}^{t_2}{\mathit{d}}^H\mathit{k}\delta \mathit{d}dt\}dz=0\end{array}$$

(27)

Considering that constraints $\mathrm{\delta }{\mathbf{d}}^H=\mathrm{\delta }\mathbf{d}=0$ at $t_1$ and $t_2$ , by rearranging equation (27) and extracting $\mathrm{\delta }{\mathbf{d}}^H$ and $\mathrm{\delta }\mathbf{d}$ , we obtain

$$\frac{1}{2}\sum \{\delta {\mathit{d}}^H{\int }_{t_1}^{t_2}(\mathit{m}\stackrel{··}{\mathit{d}}+\mathit{kd})\mathrm{dt}+{\int }_{t_1}^{t_2}({\stackrel{··}{\mathit{d}}}^H\mathit{m}+{\mathit{d}}^H\mathit{k})\mathrm{dt}\delta \mathit{d}\}=0$$

(28)

Therefore, the equation of motion of the linear acoustic waveguide can be expressed as follows:

$$\sum \{\mathit{m}\stackrel{··}{\mathit{d}}+\mathit{kd}\}=0$$

(29)

$$\sum \{{\stackrel{··}{\mathit{d}}}^H\mathit{m}+{\mathit{d}}^H\mathit{k}\}=0$$

(30)

where $\mathbf{m}$ and $\mathbf{k}$ denote real symmetric matrices; therefore, equations (29) and (30) are identical.

Combining the discrete element matrices into a global matrix, the equation of motion of the system can be obtained as follows:

$$\mathit{KD}+\mathit{M}\stackrel{··}{\mathit{D}}=0$$

(31)

If the node displacement has a time harmonic factor $e^{-i\omega t}$ , that is, $\mathit{D}=\overline{\mathit{D}}{e}^{-i\omega t}$ and $\stackrel{··}{\mathit{D}}=-{\omega }^2\overline{\mathit{D}}{e}^{-i\omega t}$ , and these are substituted into equation (31), the dispersion equation of the relationship between the wavenumber $k$ and the angular frequency $\omega$ is expressed as follows:

$$(\mathit{K}-{\omega }^2\mathit{M})\overline{\mathit{D}}=0$$

(32)

Equation (32) is an eigenvalue problem, and non-zero solution exists if and only if the determinant of the global matrix is zero.

$$\det [\mathit{K}-{\omega }^2\mathit{M}]=0$$

(33)

Dispersion equation for a circular acoustic waveguide

By substituting the expressions for the displacement $\mathbf{u}$ , strain $\mathbf{S}$ , and interpolation function of any node of a single element from equations (21a) and (22a) into equations (12) and (13b), the strain energy $U$ , kinetic energy $T$ , and Lagrangian $L$ of each element are obtained as follows:

$$U=\frac{1}{2}\sum {\int }_{\theta }{\int }_x{\int }_r{\mathit{d}}^H{\mathit{B}}^H\mathit{CBd}r\mathrm{drdx}d\theta =\frac{1}{2}\sum {\int }_{\theta }{\mathit{d}}^H\mathit{kd}d\theta$$

(34a)

$$T=\frac{1}{2}\sum {\int }_{\theta }{\int }_x{\int }_r{\stackrel{·}{\mathit{d}}}^H{\mathit{N}}^H\mathbf{\rho }\mathit{N}\stackrel{·}{\mathit{d}}\mathrm{rdrdxd}\theta =\frac{1}{2}\sum {\int }_{\theta }{\stackrel{·}{\mathit{d}}}^H\mathit{m}\stackrel{·}{\mathit{d}}d\theta$$

(34b)

$$L=T-U=\frac{1}{2}\sum {\int }_{\theta }({\stackrel{·}{\mathit{d}}}^H\mathit{m}\stackrel{·}{\mathit{d}}-{\mathit{d}}^H\mathit{kd})d\theta$$

(34c)

where

$$\mathbf{m}=\underset{x}{\int }\underset{r}{\int }{\mathbf{N}}^H\mathbf{\rho }\mathbf{N}\mathrm{rdrdx}$$

(35a)

$$\mathbf{k}=\underset{x}{\int }\underset{r}{\int }{\mathbf{B}}^H\mathbf{CB}\mathrm{rdrdx}$$

(35b)

Moreover, $\mathbf{m}$ and $\mathbf{k}$ denote the mass and stiffness matrices of the element, respectively.

Referring to equations (25)–(30), the equation of motion for the circular acoustic waveguide without external force is expressed as follows:

$$\mathit{K}\mathit{D}+\mathit{M}\stackrel{··}{\mathit{D}}=0$$

(36)

If the node displacement has a time harmonic factor $e^{-i\omega t}$ , that is, $\mathbf{D}=\overline{\mathbf{D}}{e}^{-i\omega t}$ , the dispersion equation of the relationship between the wavenumber $k$ and the angular frequency $\omega$ is established from equation (36) as follows:

$$(\mathit{K}-{\omega }^2\mathit{M})\overline{\mathit{D}}=0$$

(37)

Equation (37) is an eigenvalue problem, and a non-zero solution exists if and only if the determinant of the global matrix is zero.

$$\det [\mathit{K}-{\omega }^2\mathit{M}]=0$$

(38)

Transformation of element matrices

The interpolation function of equation (17d) is expressed in natural coordinates, and equation (18c) includes the derivatives of the interpolation function with respect to the actual coordinates $x$ and $y$ . To perform area or volume integrals of elements in natural coordinates, we must establish transformation relationships for the differential operator and the area of an element between these coordinate systems (Figure 3).

Figure 3.

Conversion maps for global and natural coordinates.

According to the chain rule for differentiating two-variable functions,^31,32 the partial derivatives of the interpolation function $N_i$ with respect to the natural coordinates $\xi$ and $\eta$ can be expressed as

$$\frac{\partial N_i}{\partial \xi }=\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \xi }+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \xi }$$

(39)

$$\frac{\partial N_i}{\partial \eta }=\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \eta }+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \eta }$$

(40)

These equations can be written in matrix form as

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}\frac{\partial N_i}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial \eta }}\end{array}\right\}=\mathbf{J}\left\{\begin{array}{c}\frac{\partial N_i}{\partial x}\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial y}}\end{array}\right\}}\end{array}$$

(41)

where $\mathbf{J}$ represents the Jacobian matrix, which is expressed as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{J}=\left[\begin{array}{cc}\frac{\partial x}{\partial \xi }\hfill & \frac{\partial y}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial x}{\partial \eta }}& \frac{\partial y}{\partial \eta }\hfill \end{array}\right]}\end{array}$$

(42)

Therefore, the partial derivatives of $x$ and $y$ can be expressed in terms of natural coordinates.

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}\frac{\partial N_i}{\partial x}\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial y}}\end{array}\right\}={\mathbf{J}}^{-1}\left\{\begin{array}{c}\frac{\partial N_i}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial \eta }}\end{array}\right\}=\left[\begin{array}{cc}{\mathrm{\Gamma }}_{11}\hfill & {\mathrm{\Gamma }}_{12}\hfill \\ \multicolumn{1}{c}{{\mathrm{\Gamma }}_{21}}& {\mathrm{\Gamma }}_{22}\hfill \end{array}\right]\left[\begin{array}{c}\frac{\partial N_i}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial N_i}{\partial \eta }}\end{array}\right]}\end{array}$$

(43)

where

$${\mathrm{\Gamma }}_{11}=\frac{J_{22}}{det(J)},{\mathrm{\Gamma }}_{12}=-\frac{J_{12}}{det(J)}$$

(44a)

$${\mathrm{\Gamma }}_{21}=-\frac{J_{21}}{det(J)},{\mathrm{\Gamma }}_{22}=\frac{J_{11}}{det(J)}$$

(44b)

The unit area $\mathrm{dA}=\mathrm{dxdy}$ of the actual coordinate system can be expressed in the natural coordinate system as

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{dA}=\left|\begin{array}{cc}\frac{\partial x}{\partial \xi }\hfill & \frac{\partial y}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial x}{\partial \eta }}& \frac{\partial y}{\partial \eta }\hfill \end{array}\right|d\xi d\eta =det(J)d\xi d\eta }\end{array}$$

(45)

According to the aforementioned relationships, the element mass matrix (24a) of the linear acoustic waveguide can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{m}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathbf{N}}^H(\xi ,\eta )\mathbf{\rho }\mathbf{N}(\xi ,\eta )det(J)d\xi d\eta }\\ =\left[\begin{array}{cccc}{\mathbf{m}}_{11}\hfill & {\mathbf{m}}_{12}\hfill & {\mathbf{m}}_{13}\hfill & {\mathbf{m}}_{14}\hfill \\ \multicolumn{1}{c}{{\mathbf{m}}_{21}}& {\mathbf{m}}_{22}\hfill & {\mathbf{m}}_{23}\hfill & {\mathbf{m}}_{24}\hfill \\ \multicolumn{1}{c}{{\mathbf{m}}_{31}}& {\mathbf{m}}_{32}\hfill & {\mathbf{m}}_{33}\hfill & {\mathbf{m}}_{34}\hfill \\ \multicolumn{1}{c}{{\mathbf{m}}_{41}}& {\mathbf{m}}_{42}\hfill & {\mathbf{m}}_{43}\hfill & {\mathbf{m}}_{44}\hfill \end{array}\right]\hfill \end{array}$$

(46)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{m}}_{\mathrm{ij}}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathbf{N}}_i^H\mathbf{\rho }{\mathbf{N}}_{j}det(J)d\xi d\eta }\\ =\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\rho \left[\begin{array}{ccc}{N}_i{N}_j\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_i{N}_j\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & N_i{N}_j\hfill \end{array}\right]det(J)d\xi d\eta \hfill \end{array}$$

(47)

Here, $i$ and $j$ denote integers from 1 to 4. Therefore, the element stiffness matrix (24b) of the linear acoustic waveguide can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{k}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathbf{B}}^H(\xi ,\eta )\mathbf{CB}(\xi ,\eta )det(J)d\xi d\eta }\\ =\left[\begin{array}{cccc}{\mathbf{k}}_{11}\hfill & {\mathbf{k}}_{12}\hfill & {\mathbf{k}}_{13}\hfill & {\mathbf{k}}_{14}\hfill \\ \multicolumn{1}{c}{{\mathbf{k}}_{21}}& {\mathbf{k}}_{22}\hfill & {\mathbf{k}}_{23}\hfill & {\mathbf{k}}_{24}\hfill \\ \multicolumn{1}{c}{{\mathbf{k}}_{31}}& {\mathbf{k}}_{32}\hfill & {\mathbf{k}}_{33}\hfill & {\mathbf{k}}_{34}\hfill \\ \multicolumn{1}{c}{{\mathbf{k}}_{41}}& {\mathbf{k}}_{42}\hfill & {\mathbf{k}}_{43}\hfill & {\mathbf{k}}_{44}\hfill \end{array}\right]\hfill \end{array}$$

(48)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{k}}_{\mathrm{ij}}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathbf{B}}_i^H(\xi ,\eta )\mathbf{C}{\mathbf{B}}_j(\xi ,\eta )det(J)d\xi d\eta }\\ =\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\left[\begin{array}{ccc}A\hfill & B\hfill & C\hfill \\ \multicolumn{1}{c}{\overline{B}}& D\hfill & E\hfill \\ \multicolumn{1}{c}{\overline{C}}& \overline{E}\hfill & F\hfill \end{array}\right]det(J)d\xi d\eta \hfill \end{array}$$

(49)

The matrices $\mathbf{m}$ and $\mathbf{k}$ of the linear acoustic waveguide expressed in equations (47) and (49) are real and symmetric.

Similarly, the element mass matrix (35a) of the circular acoustic waveguide can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{m}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathbf{N}}^H(\xi ,\eta )\mathbf{\rho }\mathbf{N}(\xi ,\eta )rdet(J)d\xi d\eta }\\ =\left[\begin{array}{cccc}{\mathbf{m}}_{11}\hfill & {\mathbf{m}}_{12}\hfill & {\mathbf{m}}_{13}\hfill & {\mathbf{m}}_{14}\hfill \\ \multicolumn{1}{c}{{\mathbf{m}}_{21}}& {\mathbf{m}}_{22}\hfill & {\mathbf{m}}_{23}\hfill & {\mathbf{m}}_{24}\hfill \\ \multicolumn{1}{c}{{\mathbf{m}}_{31}}& {\mathbf{m}}_{32}\hfill & {\mathbf{m}}_{33}\hfill & {\mathbf{m}}_{34}\hfill \\ \multicolumn{1}{c}{{\mathbf{m}}_{41}}& {\mathbf{m}}_{42}\hfill & {\mathbf{m}}_{43}\hfill & {\mathbf{m}}_{44}\hfill \end{array}\right]\hfill \end{array}$$

(50)

where

$$\begin{array}{c}{\mathbf{m}}_{\mathrm{ij}}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathbf{N}}_i^H(\xi ,\eta )\mathbf{\rho }{\mathbf{N}}_j(\xi ,\eta )rdet(J)d\xi d\eta \hfill \\ =\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\rho \left[\begin{array}{ccc}{N}_i{N}_j\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_i{N}_j\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & N_i{N}_j\hfill \end{array}\right]rdet(J)d\xi d\eta \hfill \end{array}$$

(51)

Here, $i$ and $j$ are represent integers from 1 to 4. The Jacobian matrix of the circular acoustic waveguide is expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{J}=\left[\begin{array}{cc}\frac{\partial r}{\partial \xi }\hfill & \frac{\partial x}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial r}{\partial \eta }}& \frac{\partial x}{\partial \eta }\hfill \end{array}\right]}\end{array}$$

(52)

The element stiffness matrix (35b) of the circular acoustic waveguide is expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{k}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathbf{B}}^H(\xi ,\eta )\mathbf{CB}(\xi ,\eta )rdet(J)d\xi d\eta }\\ =\left[\begin{array}{cccc}{\mathbf{k}}_{11}\hfill & {\mathbf{k}}_{12}\hfill & {\mathbf{k}}_{13}\hfill & {\mathbf{k}}_{14}\hfill \\ \multicolumn{1}{c}{{\mathbf{k}}_{21}}& {\mathbf{k}}_{22}\hfill & {\mathbf{k}}_{23}\hfill & {\mathbf{k}}_{24}\hfill \\ \multicolumn{1}{c}{{\mathbf{k}}_{31}}& {\mathbf{k}}_{32}\hfill & {\mathbf{k}}_{33}\hfill & {\mathbf{k}}_{34}\hfill \\ \multicolumn{1}{c}{{\mathbf{k}}_{41}}& {\mathbf{k}}_{42}\hfill & {\mathbf{k}}_{43}\hfill & {\mathbf{k}}_{44}\hfill \end{array}\right]\hfill \end{array}$$

(53)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{k}}_{\mathrm{ij}}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathbf{B}}_i^H(\xi ,\eta )\mathbf{C}{\mathbf{B}}_j(\xi ,\eta )rdet(J)d\xi d\eta }\\ =\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\left[\begin{array}{ccc}A\hfill & B\hfill & C\hfill \\ \multicolumn{1}{c}{\overline{B}}& D\hfill & E\hfill \\ \multicolumn{1}{c}{\overline{C}}& \overline{E}\hfill & F\hfill \end{array}\right]rdet(J)d\xi d\eta \hfill \end{array}$$

(54)

Moreover, the matrices $\mathbf{m}$ and $\mathbf{k}$ of the circular acoustic waveguide expressed in equations (51) and (54) are real and symmetric.

Simulation results

Three-dimensional analysis of resonance modes

Finite element dynamic analysis

Modal analysis was performed using ANSYS to simulate the propagation behavior of guided waves in a circular wedge with a vertex wedge angle of 15°.^33–35 Dynamic ANSYS finite element analysis proceeds in three stages: (i) pre-processing, (ii) solution, and (iii) post-processing.

In the pre-processing stage, a geometric model of the circular wedge was established with the following dimensions: vertex angle $\theta$  = 15°, height H = 18.67 mm, outer radius 15 mm, and inner radius 10 mm. Isotropic materials and titanium (Ti) were considered in this analysis. The parameters of Ti are listed in Table 1. ³⁶ The geometric model was discretized using SOLID45 elements. These are 3D elements with eight nodes, each having three degrees of freedom: $x$ , $y$ , and $z$ . After completing the above steps, a 3D mesh map for finite element analysis was obtained (Figure 4).

Table 1.

Parameters of Ti. ³⁶

Density $\rho (g/c{m}^3)$	Young’s modulus $E(\mathrm{GPa})$	Poisson’s ratio $\upsilon$
4.85	102	0.3

Figure 4.

3D ANSYS mesh map of a circular wedge.

In the solution stage, the modal analysis was performed, and the degrees of freedom of each node at the bottom of the circular wedge were assumed to be zero. Numerical calculations were then carried out.

In the post-processing stage, all resonant frequencies and corresponding resonant modes after modal analysis were determined. Antisymmetric wave modes and their associated resonant frequencies were subsequently selected for analysis. The phase velocity of the guided wave is calculated using the following formula:

$$C=\lambda f=\frac{2\pi R}{n}f$$

(55)

where $C$ is the phase velocity, $\lambda$ is the wavelength, $f$ is the resonance frequency, and $n$ is the number of guided waves on the circumference. The dispersion curve of the relationship between the wavenumber and phase velocity was established on the basis of the phase velocity of the guided waves.

Convergence analysis

The number of elements employed in finite element analysis influences the solution convergence. Divergence is likely when only a few elements are used, but the use of too many elements is computationally wasteful.

Thus, convergence analysis is conducted to determine the number of elements of mesh elements that strike the best trade-off between convergence and computational cost. In this study, the element mesh adheres to the following empirically derived guidelines for finite element wave propagation problems^30,31:

$$\frac{h}{\lambda }\le \frac{1}{16}$$

(56)

where $h$ represents the minimum distance between element nodes and $\lambda$ represents the wavelength; the minimum distance between two adjacent nodes must be less than 1/16 of the overall wavelength.

$h$ and $\lambda$ are expressed as follows:

$$h=\frac{2\pi R}{N}$$

(57a)

$$\lambda =\frac{2\pi R}{n}$$

(57b)

where $N$ is the number of circumferential elements, $R$ is the outer radius, and $n$ is the number of circumferential waves (i.e. the number of standing waves). By substituting (57a) and (57b) into equation (56), the following condition can be obtained:

$$\frac{n}{N}\le \frac{1}{16}$$

(58)

According to the modal analysis results of this study, when the frequency is f = 90,575 Hz, the number of standing waves $n$ on the circular wedge is 10.

On the basis of equation (58), the number of circumferential elements $N$ must be greater than or equal to 160 to satisfy the convergence limit of equation (55). In this study, the number of circumferential mesh elements was 160 (Table 2), which meets the convergence requirements.

Table 2.

Circular wedge element mesh.

Total number of elements	Number of circumferential elements $N$	Number of axial elements
12,160	160	20

Resonant modes and dispersion curve.

After the 3D ANSYS method was used to analyze the resonant modes of the circular wedge acoustic waveguide, the natural frequencies of the circular wedge and the corresponding resonant modes did not wholly correspond to sequence ASF waves and indirectly combined with other vibration modes; thus, the appropriate frequency and mode must be selected to obtain the dispersion curve. During the selection process, natural frequencies may get duplicated and resonant modes may remain similar, differing only by a certain angle of rotation with respect to the z-axis.

Table 3 lists the frequencies at which antisymmetric modes occurred. Figure 5 shows the corresponding resonant mode diagram at each frequency, with the number of circumferential antisymmetric waves n on the circular wedge ranging from 1 to 10. Substituting the values of $n$ and the natural frequency $f$ into equation (55) yielded the phase velocity $C$ and wavenumber $k$ of guided wave propagation.

Table 3.

Frequencies and wave velocities of a Ti circular wedge.

$n$	Frequency $f$ (Hz)	Wavenumber $k$ (1/mm)	Phase velocity $C$ (mm/µs)
1	21,492	0.067	2.02557
2	18,120	0.133	0.85388
3	24,582	0.200	0.77227
4	32,653	0.267	0.76937
5	41,015	0.333	0.77311
6	49,768	0.400	0.78175
7	59,040	0.467	0.79491
8	68,906	0.533	0.81178
9	79,409	0.600	0.83157
10	90,575	0.667	0.85365

Figure 5.

Resonant mode shapes obtained using 3D ANSYS: (a) n = 1, (b) n = 2, (c) n = 3, (d) n = 4, (e) n = 5, (f) n = 6, (g) n = 7, (h) n = 8, (i) n = 9, and (j) n = 10.

A comparison of these results with those obtained via Bi-d FEM analysis (Figure 6) revealed minor discrepancies in the dispersion curves found using the two approaches. When the wavenumber $k$ was low, the phase velocity $C$ of the dispersion curve was affected by the bottom boundary conditions and increased rapidly.

Figure 6.

Comparison between the 3D ANSYS and Bi-d FEM results for the first mode.

Bi-d FEM analysis results

The results of the 3D ANSYS analysis of the resonant modes indicated that the method faced limitations related to element discretization, restricting its effectiveness to low-frequency modes. In this study, the Bi-d FEM was employed to analyze guided wave propagation behavior. The separation of variables method was used to decouple the time-harmonic propagation factor from the cross-sectional vibrational modes. This enabled an analysis of the high frequency modes, providing insights into guided wave propagation and motion with respect to cross-sectional resonant modes.

Equations (32) and (37) represent the dispersion equations for linear and circular acoustic waveguides, respectively. Solving for eigenvalues and corresponding eigenvectors of these characteristic equations yields dispersion curves and resonant modes, elucidating wave propagation characteristics of linear and circular guided waves. Numerical solutions and simulation programs for the Bi-d FEM were implemented using FORTRAN.

The method used to solve the eigenvalue problems involves fixing the wavenumber $k$ , employing the half-interval method to determine a phase velocity satisfying the dispersion equation within an acceptable error tolerance, and subsequently plotting the dispersion curve with wavenumber $k$ as the abscissa and phase velocity as the ordinate.

The resonant mode is calculated on the basis of the direction of the dispersion curve, thereby fixing the wavenumber $k$ , and sequentially determining the resonant mode displacement corresponding to different phase velocities $C$ .

The following calculations are based on the specific cases of linear and circular acoustic waveguides with different cross-sectional shapes, in which Ti metal is used as the material. The Bi-d FEM was used to establish dispersion curves and resonant mode shapes, and the effects of boundary conditions and curvature on wedge guided wave propagation were investigated.

Dispersion curve and mode shapes of a linear wedge

This study considered a linear wedge with a vertex angle of 30° and a cross-sectional height H = 18.67 mm (Figure 2(a)). The cross-section was divided into 680 elements with 757 nodes, and the number of degrees of freedom of each node at the bottom was set to zero. The results of the analysis are shown in Figure 7. For each resonant mode, the solid lines represent dispersion curves, and dotted lines represent phase velocities obtained using Lagasse’s empirical formulas $V=V_R\sin (n\theta )$ and $n\theta \le {90}^{°}$ . ³ The Rayleigh wave velocity V _R  = 2.634 mm/µs of Ti metal was calculated using Achenbach’s method. ³³

Figure 7.

Dispersion curves of a linear wedge with a vertex angle of ${30}^{°}$ .

As shown in Figure 7, a linear wedge with a vertex angle of ${30}^{°}$ had three antisymmetric guided waves, all with phase velocities below the Rayleigh wave velocity. Comparing these phase velocities with the results of Lagasse’s empirical formula, ³ we observed that the first resonant mode’s phase velocity exceeded its predicted counterpart and that those of the second and third resonant modes were lower than the predicted values.

Figure 7 indicates a rapid increase in phase velocities for each mode’s dispersion curve at low wavenumbers. This is affected by the constraints of the degrees of freedom of each node at the bottom. However, beyond a certain wavenumber, the dispersion curves become nearly horizontal, indicating minimal dispersion effects during guided wave propagation.

Figure 8 shows 3D mesh maps of linear wedges with wedge angles of 30° and 15°. Figures 9 to 11 illustrate resonant mode shapes of linear wedges with a vertex angle of 30° at various wavenumbers; the resonant mode corresponding to each dispersion curve changed substantially as the wavenumber increased. At low wavenumber $k$ , the guided wave mode was influenced by the bottom boundary conditions, making it impossible for the guided wave to be concentrated at the wedge tip. Above a certain wavenumber $k$ , the dispersion curve flattened. At this point, most of the energy of the guided wave became concentrated at the wedge tip, and the wave pattern was antisymmetric.

Figure 8.

3D models of linear wedges: (a) 30° wedge and (b) 15° wedge.

Figure 9.

First-mode shapes of a linear wedge with a vertex angle of ${30}^{°}$ : (a) k = 0.1, (b) k = 0.5, (c) k = 1.0, (d) k = 2.0, (e) k = 2.5, (f) k = 3.1, (g) k = 3.5, and (h) k = 4.0.

Figure 10.

Second-mode shapes of a linear wedge with a vertex angle of ${30}^{°}$ : (a) k = 0.1, (b) k = 0.5, (c) k = 1.0, (d) k = 2.0, (e) k = 2.5, (f) k = 3.1, (g) k = 3.5, and (h) k = 4.0.

Figure 11.

Third-mode shapes of a linear wedge with a vertex angle of ${30}^{°}$ : (a) k = 0.1, (b) k = 0.5, (c) k = 1.0, (d) k = 2.0, (e) k = 2.5, (f) k = 3.1, (g) k = 3.5, and (h) k = 4.0.

Figure 12 compares the dispersion curves of a linear wedge with a vertex angle of ${30}^{°}$ under different bottom boundary conditions. The solid line represents a clamped bottom boundary, and the dotted line represents a simply supported bottom boundary. The relationship between the wavenumber $k$ and the wavelength $\lambda$ is $k=2\pi /\lambda$ . At lower wavenumber $k$ (longer $\lambda$ ), the influence of the bottom boundary conditions became more pronounced, leading to differences in the dispersion characteristics of the guided wave. Conversely, at higher $k$ (shorter $\lambda$ ), the wave energy was concentrated at the tip, reducing the influence of the bottom boundary conditions. Therefore, at higher values of the wavenumber $k$ , the dispersion curves overlap.

Figure 12.

Comparison of the dispersion curves of a linear wedge with a vertex angle of 30° under different bottom boundary conditions.

Subsequently, we considered a linear wedge acoustic waveguide with a vertex angle of ${15}^{°}$ . The outer boundary of the wedge was perpendicular to the bottom boundary (Figure 8(b)). The cross-section was divided into 254 elements with 284 nodes. The number of degrees of freedom of each node at the bottom boundary was set to zero. The dispersion curves obtained using the Bi-d FEM analysis are shown in Figure 13. For each resonant mode, the solid lines represent the dispersion curves, and the dotted lines represent the phase velocities estimated using Lagasse’s empirical formula. ³ At ${15}^{°}$ , the dispersion curves indicated six antisymmetric wave modes having phase velocities below the Rayleigh wave velocity. The phase velocities of the first three modes exceeded the estimated velocities and those of the remaining three higher modes were less than the estimated velocities. The linear wedge with a vertex angle of ${15}^{°}$ had three more modes than the linear wedge with a vertex angle of ${30}^{°}$ corroborating Lagasse’s findings ³ ; that mode count increases with the vertex angle of a linear wedge.

Figure 13.

Comparison of the dispersion curves of linear wedges with wedge angles of 30° and 15°.

Figure 13 shows that in the dispersion curve of each mode at low wavenumbers, the phase velocity is affected by the constraints of the number of degrees of freedom of each node at the bottom boundary and increases rapidly. When the wavenumber was greater than a certain value, the dispersion curves of each mode almost flattened, indicating that the numerical analysis results are consistent with Lagasse’s conclusion ³ that no dispersion phenomenon exists in the guided waves of ideal linear wedges.

Dispersion curve and mode shapes of a circular wedge

The cross-sectional shape of the circular wedge acoustic waveguide analyzed in this study is shown in Figure 2(b), and its corresponding 3D mesh map analyzed using ANSYS is shown in Figure 4.

The wedge had a vertex angle of ${15}^{°}$ , a cross-sectional aspect ratio of $W/H=\tan ({15}^{°})$ , and an outer radius R = 15 mm. The cross-section was divided into 254 elements with 284 nodes, with the bottom degrees of freedom of each node set to zero. Figure 14 shows the dispersion curves of the guided wave, with a total of 5 resonant modes.

Figure 14.

Dispersion curve of a circular wedge with a vertex angle of 15°.

Figures 15 to 19 illustrate the resonant mode shapes at different wavenumbers. As $k$ increased, the resonant mode deformation for each dispersion curve became increasingly concentrated at the wedge tip. At a low wavenumber $k$ , the guided waves in each mode were influenced by the bottom boundary conditions, preventing energy concentration at the wedge tip. As the wavenumber $k$ increased, the energy of the wave gradually became more concentrated at the wedge tip, and the wave pattern of the guided wave began to exhibit antisymmetric behavior.

Figure 15.

First-mode shapes of the 15° circular wedge: (a) k = 0.1, (b) k = 0.5, (c) k = 1.0, (d) k = 1.5, (e) k = 2.0, (f) k = 3.0, (g) k = 4.0, and (h) k = 5.0.

Figure 16.

Second-mode shapes of the 15° circular wedge: (a) k = 0.1, (b) k = 0.5, (c) k = 1.0, (d) k = 1.5, (e) k = 2.0, (f) k = 3.0, (g) k = 4.0, and (h) k = 5.0.

Figure 17.

Third-mode shapes of the 15° circular wedge: (a) k = 0.1, (b) k = 0.5, (c) k = 1.0, (d) k = 1.5, (e) k = 2.0, (f) k = 3.0, (g) k = 4.0, and (h) k = 5.0.

Figure 18.

Fourth-mode shapes of the 15° circular wedge: (a) k = 0.1, (b) k = 0.5, (c) k = 1.0, (d) k = 1.5, (e) k = 2.0, (f) k = 3.0, (g) k = 4.0, and (h) k = 5.0.

Figure 19.

Fifth-mode shapes of the 15° circular wedge: (a) k = 0.1, (b) k = 0.5, (c) k = 1.0, (d) k = 1.5, (e) k = 2.0, (f) k = 3.0, (g) k = 4.0, and (h) k = 5.0.

Figure 20 shows the comparison of the dispersion curves for circular (R = 15 mm) and linear (R = ∞) wedges acoustic waveguides, both with a vertex angle of ${15}^{°}$ under fixed bottom boundary conditions. For the first mode When the wavenumber ( $k$  > 0.6) and second mode ( $k$  > 1.4), the dispersion curves of the linear and circular wedge acoustic waveguides were nearly identical with each other, with a relative error below 2%. In other ranges, the difference between the dispersion curves of the linear and circular acoustic waveguides increased with the order of the resonant modes.

Figure 20.

Comparison between the linear wedge (R = ∞) and circular wedge (R = 15 mm).

We conclude that for low-order modes, the influence of the curvature of the acoustic waveguide is minimal, resulting in similar phase velocity dispersion curves for linear and circular wedges. However, for higher-order modes, the influence of the curvature on the phase velocity becomes more pronounced.

Conclusion and discussion

3D ANSYS resonant mode analysis yields resonant frequencies and modes in ascending order; however, not all are antisymmetric. To obtain a dispersion curve, appropriate frequencies and modes must be selected. During the selection process, resonant modes may coincide with subtly distinct ones differing by a rotational angle. This arises owing to the irrational nature of a circle’s circumference and the integer discretization of elements, resulting in an error in the number of significant digits and thus causing differences in numerical calculations. The accuracy of ANSYS analysis of circular wedges with 3D finite elements is limited by the number of element divisions. Therefore, to precisely investigate the wave propagation behavior of antisymmetric guided waves in acoustic waveguides, this study uses the Bi-d FEM instead of resonance analysis with 3D finite elements.

Our Bi-d FEM calculation confirm that the number of antisymmetric wave modes in a linear wedge acoustic waveguide decreases as the vertex angle increases, with the wave velocity of each mode being less than the Rayleigh wave velocity. Assuming a phase velocity below the Rayleigh wave velocity (V _R  = 2.634 mm/µs), linear acoustic waveguide with vertex angle of ${30}^{°}$ and ${15}^{°}$ have three and six antisymmetric wave modes, respectively, which agree with the results of Lagasse’s empirical formulas, $V=V_R\sin (n\theta )$ and $n\theta \le {90}^{°}$ . ³

Owing to the boundary conditions at the bottom of the dispersion curves of the linear and circular wedge acoustic waveguides, the phase velocity increases considerably at low frequencies. This indicates that the boundary conditions only influence wedge acoustic waveguides in low-frequency bands.

When the wavenumber $k$ decreases, the wavelength $\lambda$ increases and the bottom boundary conditions considerably affect the phase velocity of the dispersion curves. However, when $k$ increases, $\lambda$ decreases, the energy of the wave is concentrated at the tip of the wedge, and the dispersion curves gradually disappear owing to the changes in the bottom boundary conditions and tend to overlap.

The influence of a structural curvature on the guided wave phase velocity of a wedge becomes more pronounced for higher modes. The guided wave phase velocity in a circular wedge increases relative to its linear counterpart as the resonant mode order increases.

The obtained resonant mode shapes confirm that the wedge acoustic waveguide exhibits the largest displacement at its tip. This indicates that the energy of the wave is concentrated at the tip, and the guided wave is an ASF wave. The Bi-d FEM proves to be a powerful analysis method for the study of wave propagation in circular and linear acoustic waveguides.

This study methodology can be extended in future research to investigate linear and circular wedge acoustic waveguides with different cross-sectional geometries, wedge angles, truncations, and coatings. Such investigations can lead to the development of ultrasonic motors and related applications.