Analytical Model and Gas Leak Source Localization Based on Acoustic Emission for Cylindrical Storage

1. Introduction

Cylindrical vessels are more commonly used as pressure vessels than spherical vessels due to their ease of fabrication, space utilization, pressure resistance and ease of maintenance. In a confined space, cylindrical tanks are very effective for the storage and use of fuel gases due to their ease of transport, high storage capacity and pressure regulation. The need for reliable maintenance of high-pressure gas cylinders and tanks, especially those containing hazardous gases, has increased exponentially over the past decade due to the high potential for safety hazards from leaks caused by undetected vessel damage. The most important way to prevent disasters caused by gas leaks is to detect the leak and its location and repair it accordingly. Several leak detection methods have been proposed for the early detection of leaks in pipelines [1]. These include acoustic emission (AE) detection, flow rate or pressure monitoring, or infrared detection [2,3,4,5,6,7,8]. Among the various methods, AE methods have proven to be very effective in detecting and locating leaks by eliminating interfering signals through time and frequency domain analysis and characterization of AE parameters [9,10,11,12].

In general, the most notable methods for locating an AE source are the time domain technique and the zone location technique [13]. The former relies on the time difference of arrival (TDOA) between two or more AE sensors to derive the source location [14]. However, the time delay can also result from different wave paths, different wave modes due to the dispersion of the AE wave for a long-distance pipeline and the velocity variation due to the environment. Therefore, the TDOA technique requires high detection accuracy. However, because the TDOA method is applied to burst-shaped waveforms, it is not suitable for locating leak sources in the case of leaks that produce continuous AE waveforms. To overcome this, the cross-correlation algorithm is sometimes used to obtain the TDOA from continuous waveforms [15,16]. Due to background noise and other uncertainties during the measurement, achieving a good level of correlation coefficient is limited [10,17]. In practical situations, the AE signal generated by a leak has an aperiodic waveform structure and the exact location of the leak source cannot be determined by the cross-correlation based methods. The second one of the basic AE localization techniques is based on the attenuation of the AE amplitude with the distance between the leak source and the detection sensor, taking advantage of the fact that the sensor with the largest amplitude is closest to the leak source, assuming that the sensitivity of the sensors is the same. However, this method has limitations in pinpointing the location of the leak in three dimensions because it ignores the influence of the AE signal due to the structural characteristics of the vessel. To overcome this drawback, AE signal analysis is required using governing equations that take into account the structural characteristics of the storage vessel and the propagation dynamics of AE caused by leaks on the vessel surface. In addition, this method can only be used if the structure has a moderate attenuation of the AE signal. This technique can provide a simple solution where the exact location of the leak source is not critical.

Over the past few decades, numerous studies have been conducted to solve the problems of pipeline leak localization using acoustic-based methods [10,12]. However, there are still limitations and gaps in fundamental research dealing with signal modelling, nonlinear effects and environmental uncertainty in leak AE signals [10]. Previously, a mathematical formula was formulated for AE generated by a point source (PS) as an internal defect in cylindrical structures [18]. The PS was represented by a concentrated force (CF) with both spatial and temporal harmonic characteristics. In cylindrical geometry, the CF generates three potential functions representing a compression wave (P) and two shear waves (horizontal and vertical polarization, SH and SV, respectively), which are referred to as the concentrated force-incorporated potential (CFIP). The radial, tangential and axial displacements were solved by introducing the CFIPs into the Navier-Lamé (NL) equation, based on the model proposed by Morse and Feshbach [19]. Previously, the CFIPs for the pinhole leakage were derived as an excitation AE source, associated with a fluctuating Reynolds stress (FRS) [20]. As the FRS acts on the pinhole wall, the CFIP due to the FRS results only in the axial and tangential displacements. A limited investigation of the angular and axial dependence of the AE signals for the pinhole leakage source has been carried out. The characteristics of the observed AE signals due to the pinhole leakage were well simulated by the proposed theoretical formula. The aim of this study is to investigate the feasibility of applying the proposed mathematical model to leakage source localization in cylindrical geometries. Detailed measurements of the AE signals generated by leakages through pinholes of different diameters in gas cylinders were performed. For application to leak source (LS) localization, the theoretical formula was modified into two semi-empirical equations to describe the angular and axial dependence of the observed AE amplitude, respectively. Finally, the simulation and experimental results were compared to verify the accuracy of the leak source localization.

2. Mathematical Solution

2.1. Concentrated Forces

In general, the CF as the PS excitation of the AE is described as

$$\mathit{f}=\mathit{P}\delta (\mathit{x}-{\mathit{x}}_0)e^{-i\omega t},$$

(1)

where P is the force vector acting at ${\mathit{x}}_0$ and $\omega \left(=2\pi \nu \right)$ is the angular frequency that transfers the AE energy from the PS to a solid space. In Eq. (1), the delta function provides the spatial distribution of the CF at a given time as Green’s function $g\left(\mathit{x};{\mathit{x}}_0\right),$ defined as ${\nabla }^{2}g\left(\mathit{x};{\mathit{x}}_0\right)=\delta \left(\mathit{x}-{\mathit{x}}_0\right).$ The CF generated by the PS can be rewritten in terms of vectors as

$$\mathit{f}=\mathit{P}{\nabla }^{2}g\left(\mathit{x};{\mathit{x}}_0\right)e^{-i\omega t},$$

(2a)

or

$$\mathit{f}=\nabla \left[\nabla ·\mathit{P}g\left(\mathit{x};{\mathit{x}}_0\right)\right]-\nabla \times \left[\nabla \times \mathit{P}g\left(\mathit{x};{\mathit{x}}_0\right)\right]e^{-i\omega t}.$$

(2b)

The derivation of the Green’s function for the cylindrical cell structure [18,20,21] is summarized in Appendix A. Applying the periodicity of 2π ($v$= 0) to Eqs. (A17)–(A19) gives

$$g\left(\xi ,\vartheta \eta ;0\right)=G_0\left(\eta \right)J_0\left({\kappa }_{\eta }\xi \right),$$

(3)

$$G_0=-{\displaystyle \frac{1}{2{\kappa }_z}}{A}_0{J}_0\left({\kappa }_{\eta }{\xi }_0\right)\times \left\{\begin{array}{c}{e}^{{\kappa }_{\eta }\eta }\left(0\le \eta \le l-z_0\right)\\ e^{-{\kappa }_{\eta }\eta }\left(-z_0\le \eta \le 0\right)\end{array}\right.,$$

(4)

$$A_0={\displaystyle \frac{1}{\pi \left[e^{k_{\eta }{z}_0}+e^{-k_{\eta }\left(l-z_0\right)}-2\right]}}\times {\displaystyle \frac{2}{\left\{a^2{\left[J_1\left({\kappa }_{\eta }a\right)\right]}^2-b^2{\left[J_1\left({\kappa }_{\eta }b\right)\right]}^2\right\}}},$$

(5)

where the PS is located at $\xi =0,\vartheta =0$ and $\eta =0$. In Eq. (3), the value of $\xi$ is the shortest distance ${\xi }_P$ between the PS and the detection point on the cylinder surface (Figure A1). If the thickness is much shorter than the diameter of the outer circle, the value of $\xi$ can be given as the length of the arc around the outer circle,

$$\xi \cong {\xi }_P={\displaystyle \frac{\pi \theta }{180}}a={\displaystyle \frac{\pi \left(\vartheta -90\right)}{90}}a,$$

(6)

where the PS and the detector are located on the outer surface of the cylinder (Notes. Corrects a typographical error in ${\xi }_P$ inserted in Figure 2b of Ref. 20). In Refs. 18 and 21, the value of ${\kappa }_z$ was determined by applying boundary conditions to the Bessel function $J_0\left({\kappa }_{z}r\right)$ in Eq. (A6). For the outer surface ${\left.J_0\left({\kappa }_{z}r\right)\right|}_{r=a-r_0}=0$ gives ${\kappa }_z=r_{0n}/\left(a-r_0\right)$. On the other hand, in Ref. 20, the value of ${\kappa }_z$ was determined experimentally from the angular dependence of the AE amplitude. The subscripts $v$ and n of $A_{vn}$ and $G_{vn}$ represent the azimuthal constant and the nth-order root in the boundary conditions, respectively. In this study, the value of ${\kappa }_{\eta }$ was determined using the external excitation method: therefore, the subscript n was omitted in Eqs. (3)–(5).

The FRS is the most effective excitation source for AE due to fluid leakage through a pinhole. For turbulent flow with a radial turbulent velocity of $v_r^{\prime }$ through a circular hole in the cylindrical structure, the CF is exerted on the wall of the pinhole by the FRS. Denoting the fluctuating velocities along the axial (z) and tangential (θ) directions as $v_z^{\prime }$ and $v_{\theta }^{\prime }$, respectively, the FRS at low Mach numbers is defined as

$$T_{rz}^{\prime }=-\rho v_r^{\prime }{v}_z^{\prime } \mathrm{and} T_{r\theta }^{\prime }=-\rho v_r^{\prime }{v}_t^{\prime },$$

(7)

where ρ is the density of the charged gas. From Eq. (7), the force vectors acting on the pinhole wall are given by [22]

$$P_z={\left.-\rho {\displaystyle \frac{\partial \left(\overline{v_r^{\prime }{v}_z^{\prime }}\right)}{\partial x_z}}\right|}_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}\mathrm{and} P_{\theta }={\left.-\rho {\displaystyle \frac{\partial \left(\overline{v_r^{\prime }{v}_t^{\prime }}\right)}{\partial x_{\theta }}}\right|}_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}},$$

(8)

where the upper bar symbol represents the value averaged over a period of time. Since $v_z^{\prime }=v_{\theta }^{\prime }$, $P_w\equiv P_z=P_{\theta }$. The CF vector can be written down as

$$\mathit{P}=P_w\left({\widehat{\mathit{e}}}_z+{\widehat{\mathit{e}}}_{\theta }\right),$$

(9)

From the reported data [23], the CF can be obtained as

$$P_w=6.62\times {10}^{-4}\rho U^2,$$

(10)

where U is the mean velocity of $v_r^{\prime }$ and is given as Eq. (B5).

2.2. Displacement Fields

The displacement fields u generated by the CF in the cylindrical geometry can be derived by using the NL equation [24],

$$\left(\lambda +2\mu \right)\nabla \left(\nabla ·\mathit{u}\right)-\mu \nabla \times \nabla \times \mathit{u}+\mathit{f}=\rho {\displaystyle \frac{d^2\mathit{u}}{d{t}^2}},$$

(11)

and the Morse and Freshbach’s model [19,25],

$$\mathit{u}=\nabla \Phi +\nabla \times \left({\rm X}{\widehat{\mathit{e}}}_z\right)+a\nabla \times \nabla \times \left(\Psi {\widehat{\mathit{e}}}_z\right),$$

(12)

where λ and μ are Lamé constants and ρ is the density of the media. In Eq. (12), Φ is the scalar potential for the P wave, ${\rm X}{\widehat{\mathit{e}}}_z$ is vector potential for the SH wave, $\Psi {\widehat{\mathit{e}}}_z$ is vector potential for the SV wave, and a is the outer radius of the cylinder. Note that the velocities of the P and S waves are given as $c_P=\sqrt{(\lambda +2\mu )/\rho }$ and $c_S=\sqrt{\mu /\rho }$, respectively.

For the PS excitation, the three potentials can be specified as the CFIPs. In $\left(\xi ,\vartheta \eta \right)$ coordinates, the three CFIPs are expressed as

$$\Phi =\left(\nabla ·\mathit{P}\varphi \right)=P_{\theta }{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial \varphi }{\partial \vartheta }}+P_{\eta }{\displaystyle \frac{\partial \varphi }{\partial \eta }},$$

(13)

$${\rm X}{\widehat{\mathit{e}}}_{\eta }=-\left(\nabla \times \mathit{P}\mathit{\chi }\right)=-{\displaystyle \frac{P_{\theta }}{\mathit{\xi }}}{\displaystyle \frac{\partial \left(\xi \chi \right)}{\partial \xi }}{\widehat{\mathit{e}}}_{\eta },$$

(14)

$$\Psi {\widehat{\mathit{e}}}_z=\left(\nabla \times \mathit{P}\psi \right)={\displaystyle \frac{P_{\theta }}{\mathit{\xi }}}{\displaystyle \frac{\partial \left(\xi \psi \right)}{\partial \xi }}{\widehat{\mathit{e}}}_{\eta },$$

(15)

where $\varphi , \chi$ and $\psi$ are scalar functions for $\mathsf{\Phi },\mathsf{{\rm X}}$ and $\mathsf{\Psi }$, respectively. (Notes. Corrects the erratum in the sign of $\Psi {\widehat{\mathit{e}}}_z$ in Refs. 18 and 20.) The combination of Eqs. (2) and (13)–(15) gives a second-order partial differential equation (PDE) for each scalar function. The radial PDEs of $\varphi$ and $\chi$ are nonhomogeneous, and their solutions are a linear combination of the homogeneous and the particular solutions. These three PDEs were solved by introducing $k_{\eta }^2$ $\left(k_{\eta }\ne {\kappa }_{\eta }\right)$ for the imaginary roots of the axial component and $m$ $\left(=0, \pm 1, \pm 2,\cdots \right)$ as an azimuthal constant for the angular component. The solutions of the $\varphi , \chi$ and $\psi$ functions are given by Eqs. (C1)–(C3). By substituting each scale function into the corresponding potential, the CFIPs can be obtained as expressed in Eqs. (C4)–(C9).

When the displacement is produced by the force component $P_f$, u in Eq. (12) can be rewritten as the cylindrical components of the displacement

$$u_f=u_{\xi f}+u_{\vartheta f}+u_{\eta f},$$

(16)

where

$$u_{\xi f}={\displaystyle \frac{\partial {\Phi }_f}{\partial \xi }}+{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial {{\rm X}}_f}{\partial \vartheta }}+a{\displaystyle \frac{{\partial }^2{\Psi }_f}{\partial \xi \partial \eta }},$$

(17)

$$u_{\vartheta f}={\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial {\Phi }_f}{\partial \vartheta }}-{\displaystyle \frac{\partial {{\rm X}}_f}{\partial \xi }}+a{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{{\partial }^2{\Psi }_f}{\partial \vartheta \partial \eta }},$$

(18)

$$u_{\eta f}={\displaystyle \frac{\partial {\Phi }_f}{\partial \eta }}-a\left({\displaystyle \frac{{\partial }^2{\Psi }_f}{\partial {\xi }^2}}+{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial {\Psi }_f}{\partial \xi }}+{\displaystyle \frac{1}{{\xi }^2}}{\displaystyle \frac{{\partial }^2{\Psi }_f}{\partial {\vartheta }^2}}\right).$$

(19)

These three equations can be written as

$$u_{df}=P_f\left(A_{mf}{F}_{df}^1+B_{mf}{F}_{df}^2+C_{mf}{F}_{df}^3+F_{df}^4\right)e^{-i\omega t},$$

(20)

where $d=\xi ,\vartheta$ or $\eta$ and $f=z$ or $\theta$ for the pinhole leakage. In Eq. (20), the $F_{df}^n$ terms are given in Eqs. (C10)–(C15).

Since the Bessell functions $J_m\left(\alpha \xi \right)$ and $J_m\left(\beta \xi \right)$ in the CFIPs become the unity as $\xi$ approaches to the PS, then m = 0. For m = 0, the $F_{df}^n$ terms are as follows: for $P=P_z$ (radial component),

$$F_{\xi z}^1=\left(-i{k}_{\eta }\right){\displaystyle \frac{\partial J_0\left(\alpha \xi \right)}{\partial \xi }}{e}^{-i{k}_{\eta }\eta },$$

(21)

$$F_{\xi z}^2=F_{\xi z}^3=0,$$

(22)

$$F_{\xi z}^4=-{\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{1}{{\alpha }^2-{\kappa }_{\eta }^2}}{\displaystyle \frac{\partial J_0\left({\kappa }_{\eta }\xi \right)}{\partial \xi }}{\displaystyle \frac{\partial G_0}{\partial \eta }},$$

(23)

(tangential component)

$$F_{\vartheta z}^1=F_{\vartheta z}^2=F_{\vartheta z}^3=F_{\vartheta z}^4=0,$$

(24)

(axial component)

$$F_{\eta z}^1=-k_{\eta }^2{J}_0\left(\alpha \xi \right)e^{-i{k}_{\eta }\eta },$$

(25)

$$F_{\eta z}^2=F_{\eta z}^3=0,$$

(26)

$$F_{\eta z}^4=-{\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{J_0\left({\kappa }_{\eta }\xi \right)}{{\alpha }^2-{\kappa }_{\eta }^2}}{\displaystyle \frac{{\partial }^2{G}_0}{{\partial \eta }^2}},$$

(27)

for $P=P_{\theta }$,

(radial component)

$$F_{\xi \theta }^1=F_{\xi \theta }^2=0,$$

(28)

$$F_{\xi \theta }^3=\left(i{k}_{\eta }\right)a\left[{\displaystyle \frac{1}{{\xi }^2}}{J}_0\left(\beta \xi \right)-{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial J_0\left(\beta \xi \right)}{\partial \xi }}-{\displaystyle \frac{{\partial }^2{J}_0\left(\beta \xi \right)}{{\partial \xi }^2}}\right]e^{-i{k}_{\eta }\eta },$$

(29)

$$F_{\xi \theta }^4=0,$$

(30)

(tangential component)

$$F_{\vartheta \theta }^1=F_{\vartheta \theta }^2=F_{\vartheta \theta }^3=0,$$

(31)

$$F_{\vartheta \theta }^4={\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{1}{{\beta }^2-{\kappa }_{\eta }^2}}{G}_0\left\{{\displaystyle \frac{1}{{\xi }^2}}\left[J_0\left({\kappa }_{\eta }\xi \right)+{\displaystyle \frac{\partial J_0\left({\kappa }_{\eta }\xi \right)}{\partial \xi }}\right]-{\displaystyle \frac{1}{\xi }}\left[2{\displaystyle \frac{\partial J_0\left({\kappa }_{\eta }\xi \right)}{\partial \xi }}+{\displaystyle \frac{{\partial }^2{J}_0\left({\kappa }_{\eta }\xi \right)}{{\partial }^2\xi }}\right]\right\},$$

(32)

(axial component)

$$F_{\eta \theta }^1=F_{\eta \theta }^2=0,$$

(33)

$$F_{\eta \theta }^3=a\left[{\displaystyle \frac{1}{{\xi }^3}}{J}_0\left(\beta \xi \right)-{\displaystyle \frac{1}{{\xi }^2}}{\displaystyle \frac{\partial J_0\left(\beta \xi \right)}{\partial \xi }}+{\displaystyle \frac{2}{\xi }}{\displaystyle \frac{{\partial }^2{J}_0\left(\beta \xi \right)}{{\partial \xi }^2}}+{\displaystyle \frac{{\partial }^3{J}_0\left(\beta \xi \right)}{{\partial \xi }^3}}\right]e^{-i{k}_{\eta }\eta },$$

(34)

$$F_{\eta \theta }^4=0.$$

(35)

Substituting the non-zero $F_{dz}^n$ terms into Eq. (20) gives the displacement components.

For $P=P_z$,

$$u_{\xi z}=P_z\left(A_{0z}{F}_{\xi z}^1+F_{\xi z}^4\right)e^{-i\omega t},$$

(36)

$$u_{\vartheta z}=0,$$

(37)

$$u_{\eta z}=P_z\left(A_{0z}{F}_{\eta z}^1+F_{\eta z}^4\right)e^{-i\omega t}.$$

(38)

Eqs. (36)–(38) show that the $P_z$ excitation produces only the P wave AE.

For $P=P_{\theta }$, the following displacements can be obtained as

$$u_{\xi \theta }=P_{\theta }\left(C_{0\theta }{F}_{\xi \theta }^3\right)e^{-i\omega t},$$

(39)

$$u_{\vartheta \theta }=P_{\theta }{F}_{\vartheta \theta }^4,$$

(40)

$$u_{\eta \theta }=P_{\theta }\left(C_{0\theta }{F}_{\eta \theta }^3\right)e^{-i\omega t}.$$

(41)

The $P_{\theta }$ excitation does not produce the displacements because there are no $F_{d\theta }^4$ terms in Eqs. (39) and (41) and no coupling constant in Eq. (40).

To evaluate the coupling constant $A_{0z}$, let us apply the stress-free boundary conditions to the outer surface. Since the cylinder is free from external stresses at its surface, the boundary conditions can be written as

$${\sigma }_{\xi \xi }={\sigma }_{\xi \eta }=0\text{ (at any point on the cylinder surface).}$$

(42)

Expanded expressions for the stress and the strain displacement at any point on the cylinder surface are given as Eq. (C16). Substituting Eqs. (39) and (41) into Eq. (C16) with m = 0 gives

$$a_{11z}=\left({ik}_{\eta }\right)\left[\left(c_{13}{k}_{\eta }^2\right)J_0\left(\alpha \xi \right)-c_{12}{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial J_0\left(\alpha \xi \right)}{\partial \xi }}{-c}_{11}{\displaystyle \frac{{\partial }^2{J}_0\left(\alpha \xi \right)}{{\partial \xi }^2}}\right]e^{{-ik}_{\eta }\eta },$$

(43)

$$a_{13z}=0,$$

(44)

$$a_{31z}=-2k_{\eta }^2{c}_{44}{\displaystyle \frac{\partial J_0\left(\alpha \xi \right)}{\partial \xi }}{e}^{-i{k}_{\eta }\eta },$$

(45)

$$a_{33z}=0,$$

(46)

$$b_{1z}=-\left({\displaystyle \frac{k_p^2}{\rho {\omega }^2}}\right){\displaystyle \frac{1}{{\alpha }^2-{\kappa }_{\eta }^2}}\left[c_{11}{\displaystyle \frac{\partial G_0}{\partial \eta }}{\displaystyle \frac{{\partial }^2{J}_0\left({\kappa }_z\xi \right)}{\partial {\xi }^2}}+c_{12}{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial G_0}{\partial \eta }}{\displaystyle \frac{\partial J_0\left({\kappa }_z\xi \right)}{\partial \xi }}+c_{13}{\displaystyle \frac{{\partial }^3{G}_0}{\partial {\eta }^3}}{J}_0\left({\kappa }_z\xi \right)\right],$$

(47)

$$b_{3z}=-c_{44}\left({\displaystyle \frac{k_p^2}{\rho {\omega }^2}}\right){\displaystyle \frac{2}{{\alpha }^2-{\kappa }_{\eta }^2}}{\displaystyle \frac{{\partial }^2{G}_1}{{\partial \eta }^2}}{\displaystyle \frac{\partial J_0\left({\kappa }_z\xi \right)}{\partial \xi }}.$$

(48)

From these equations, the solution of $A_{0z}$ can be obtained as

$$A_{0z}={\displaystyle \frac{b_{1z}}{a_{11z}}}={\displaystyle \frac{b_{3z}}{a_{31z}}}.$$

(49)

2.3. Analytical Model

The arrival time τ of the signal at position $\left(\xi ,\vartheta ,\eta \right)$ on the outer surface must be introduced to the displacement to analyze the observed AE signal. The arrival time of the P wave propagating with velocities $c_P$ is given as

$${\tau }_P={\displaystyle \frac{\sqrt{{\xi }^2+{\eta }^2}}{c_P}},$$

where $\xi$ is given by Eq. (6). The dominant displacement generated by the gas leak at a given time t can be given by

$$u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}{=P}_z\left(t-{\tau }_P\right)\left(A_{0z}{F}_{dz}^1+F_{dz}^4\right)e^{-i\omega t}\left(d=\xi ,\eta \right).$$

(50)

The displacement $u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}$ in Eq. (50) is associated with a single frequency. Considering that $A_{0z}$ and $F_{dz}^4$ are proportional to $1/{\omega }^2$, Eq. (50) can be rewritten for a multi-frequency wave as

$$u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}\left({\omega }_1,\cdots ,{\omega }_n\right)=\sum _{i=1}^n{w}_{\mathrm{e}\mathrm{m}\mathrm{p},i}{\displaystyle \frac{{\omega }_i^2}{\sum _{j=1}^n{\omega }_j^2}}{u}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}\left({\omega }_i\right),$$

(51)

where $w_{\mathrm{e}\mathrm{m}\mathrm{p},i}$ is the observed weight of component i. In this study, the relative value of the displacement is used to validate the amplitude of the observed AE using Eq. (50).

3. Experiment

The model cylinder system for the laboratory tests is a seamless N₂ gas cylinder (40 L, Mn steel) with a main section (1013 mm long, 232 mm outer diameter and 4.8 mm thick), a spherical shoulder with a neck at one end and a concave bottom at the other (Figure 1). Leakage was simulated by inserting screws with a hole into a borehole in the wall of the main section, located at 0.368 m from the front of the main section. The diameters of the screw holes were approximately 0.20, 0.30, 0.50, 0.80 and 1.20 mm. Note that the experimental test is represented as $\mathrm{D}\mathrm{a}\mathrm{P}\mathrm{b}\mathsf{\eta }\mathrm{c}\mathsf{\vartheta }\mathrm{d}$, where a is the hole diameter (mm), b is the internal pressure (bar), and c and d are the distance $\eta$ (cm) and tangential angle $\vartheta$ (degree) from the leak source (LS), respectively. The data used in the simulation are listed in Table 1.

An 8-channel acquisition system (IDK-AET/E08, Republic of Korea), including a preamplifier, 8-channel DAQ board and computer data storage, was used to analyze the AE waveform and record data. The IDK system runs the built-in AE Studio software to generate acoustic data by sequentially amplifying and Fast Fourier transforming the detected electrical signals. No filters were used and AE hits were recorded over the range of ν = 0–500 kHz. Each hit was recorded for 1.0 ms by a timer controller on the DAQ board. For each hit, the maximum amplitude was selected, which we call the observed amplitude $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$. The signal threshold was 40 dB.

Several broadband AE sensors (IDK-AES-H150, IDK) with a resonant frequency of 150 kHz and ϕ = 16.0 mm were mounted on the cylinder surface using magnetic holders and high vacuum silicone grease. After mounting the sensors on the cylinder surface, it was necessary to calibrate each sensor, since the response of each sensor depends on internal and external influences (its own sensitivity, and contact area, grease thickness, surface roughness, etc.). The reference signal was applied to the actuator (IDK-AES-H150) by a KEYSIGHT 33600A waveform generator. The reference signal was tuned at a frequency of 150.00 kHz with an amplitude of 1.000 V_PP and a square pulse width of 992 ns. In this study, each sensor was calibrated with the reference signal generated by the waveform generator placed horizontally at 100 mm. The actuator was mounted and dismounted at least five times to get the maximum amplitude A_cal for each sensor. For a given sensor, the measured amplitude $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ was divided by its A_cal.

4. Results and Discussion

4.1. Determination of ${\kappa }_{\eta }$

and $k_{\eta }$

The parameters of ${\kappa }_{\eta }$ and $k_{\eta }$ are involved in the Green’s and the scalar functions, respectively. The relationship between these parameters can be found from Eq. (49). From Eq. (49), we can obtain the following equation

$$f\left({\kappa }_{\eta },k_{\eta }\right)={\displaystyle \frac{b_{1z}}{a_{11z}}}-{\displaystyle \frac{b_{3z}}{a_{31z}}}=0.$$

(52)

From the first root of the solution of Eq. (52), we obtained the following relationship [20]

$$k_{\eta }=\left\{\begin{array}{c}1.555-0.076 {\kappa }_{\eta }+7.09\times {10}^{-2}{\kappa }_{\eta }^2-4.20\times {10}^{-3}{\kappa }_{\eta }^3 \left({\kappa }_{\eta }\le 15\right) \\ 0.652-9.95\times {10}^{-3} {\kappa }_{\eta } \left({\kappa }_{\eta }>15\right)\end{array}\right..$$

In Reference 20, the value of ${\kappa }_{\eta }$ was determined by analyzing the angular dependence of the AE amplitude generated by the leak. However, it is desirable to predetermine the value of ${\kappa }_{\eta }$ before the leak occurs in order to apply it to the actual leak source localization. In this study, the value of ${\kappa }_{\eta }$ was determined by applying an external excitation to the surface of the cylinder. Twelve sensors were mounted on the surface of the cylinder at 30-degree intervals from 0 to 270 degrees at a point 10 cm away from the actuator position. Figure 2a shows the mean amplitude with σ value of each sensor signal produced by an external excitation (1.000 V_pp and $\nu =150$ kHz). It can be seen that the angular dependence of the observed signals on the external excitation is clear. Two minima appeared close to 90 and 270 degrees. The dependence shows π symmetry. If the external excitation is treated as the Green’s function, the angular dependence of the signal amplitude can be related to the Bessel function $J_0\left({\kappa }_{\eta }\xi \right)$ in Eq. (3). From Eq. (50), the angular dependence of $u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}$ was simulated for different values of ${\kappa }_{\eta }$ and $\nu =150$ kHz. As shown in Figure 2a, the minimum point at which the Bessel function goes to zero depends on the value of ${\kappa }_{\eta }$. In the range of 0–180°, one minimum point appears at $\vartheta =117$° for ${\kappa }_{\eta }=6$. As the value of ${\kappa }_{\eta }$ increases, the angle of the minimum point decreases: $\vartheta =100°$ for ${\kappa }_{\eta }=7$, $\vartheta =93°$ for ${\kappa }_{\eta }=7.5$ and $\vartheta =88°$ for ${\kappa }_{\eta }=8$. For ${\kappa }_{\eta }=9$, there are two minima at $\vartheta =78°$ and 164°. Considering that the angular separation between two neighboring sensors is 30°, we concluded that the ${\kappa }_{\eta }$ values in the range 7–8 are acceptable. In this study, we set ${\kappa }_{\eta }=7.5$, resulting in the maximum of $u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}$ at $\vartheta =16°$. To verify the determined ${\kappa }_{\eta }$ value, we also measured the axial dependence of the signal generated by the external excitation. We mounted 10 sensors along the axial direction relative to the position of the actuator: 9 sensors in the downward direction and 1 sensor in the upward direction. As shown in Figure 2b, the axial dependence of the observed signal amplitude is in good agreement with the simulated $u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}$ values at ${\kappa }_{\eta }=7.5$ and $\nu =150$ kHz.

4.2. Verification of Angular and Axial Dependences

The ${\kappa }_{\eta }$ parameter is crucial in determining the angular and axial dependencies by associating the coupling constant with the Green’s function, as implied in Eq. (50). In addition, the verification of the dependencies is necessary to apply the proposed mathematical model to the localization of the leak source in cylindrical storage. In this study, the verification has been extended to the whole region of the cylindrical storage.

First, the angular dependence of the AE amplitude was studied as a function of the tangential angle at a fixed η distance from the LS. Similarly to the case of the external excitation experiment, twelve sensors were mounted on the cylinder surface at 30° intervals from 0 to 270° at η = 10 cm. Among the obtained AE parameters, the maximum amplitude $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ of each AE hit was selected to analyze the angular dependence. For a given sensor, the mean amplitude was calculated from over 6 k hits. Figure 3 shows the $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ values observed at η = 10 from five different LSs with D0.2–D1.2 and three different internal pressures. Within the standard deviation, all results show that the minimum of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ occurs close to $\vartheta =90°$ and $\vartheta =270°$, and the maximum of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ occurs close to $\vartheta =30°$ and $\vartheta =330°$. Note that the calculated maximum and minimum values occur at $\vartheta =$16° and $\vartheta =$93° in the first half of the circle, respectively. To simulate the $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ values, it is necessary to convert the calculated $u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}$ value in meters to $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ in millivolts. The analysis of the angular dependence of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ can be done by normalizing the calculated displacement values and then matching the relative amplitude values observed in the 2π region. The simulated angular dependence of the amplitude is formulated as

$$A_{\mathrm{s}\mathrm{i}\mathrm{m}}\left(\vartheta \right)=u_{\mathrm{N}}\left(\vartheta \right)\left(A_{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}-A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}\right)+A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}.$$

(53)

In Eq. (53), $u_{\mathrm{N}}\left(\vartheta \right)$ is the normalized displacement, defined as $u_{\mathrm{N}}\left(\vartheta \right)=u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}\left(\vartheta \right)/u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k},\mathrm{m}\mathrm{a}\mathrm{x}}$, which is independent of the input parameters (such as FRS and frequency) in Eq. (51), and $A_{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}$ and $A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}$ are the maximum and minimum values estimated from $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ at $\vartheta =16°$ and 93°, respectively. As shown in Figure 3a-e, experimental values for D0.2–D0.8 agree well with the simulated values within the standard deviation, but for D1.2P3 and D1.2P4 there is a rather large deviation between $A_{\mathrm{s}\mathrm{i}\mathrm{m}}$ and $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$. The angular dependence illustrates the symmetry property of π and the minimum amplitudes close to 90° and 270°. These results prove that the ${\kappa }_{\eta }$ of 7.5, which was determined by the external excitation, works very well to characterize the AE signal observed due to leakage in the cylindrical geometry. The most striking feature of the observed AE amplitudes is that the value of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ is strongly influenced by the hole size. As expressed in Eq. (10), the CF due to the leakage is only proportional to the mean velocity U. As expressed in Eqs. (B3) and (B4) [26,27,28,29], substituting Q into U cancels the hole diameter factor D in Eq. (B5). Accordingly, the CF is only proportional to $P_0^2$. For a more practical application, we consider the dependence of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ on Reynolds number, expressed as

$$Re={\displaystyle \frac{C_D}{\mu }}D{P}_0\sqrt{{\displaystyle \frac{2\gamma }{\gamma -1}}{\displaystyle \frac{M}{R{T}_1}}\left[{\left({\displaystyle \frac{P_a}{P_0}}\right)}^{{\displaystyle \frac{2}{\gamma }}}-{\left({\displaystyle \frac{P_a}{P_0}}\right)}^{{\displaystyle \frac{\gamma +1}{\gamma }}}\right]{\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma +1}{\gamma -1}}}}\hspace{1em}\hspace{1em}\left(P_0<P_{cr}\right),$$

(54)

$$Re={\displaystyle \frac{C_D}{\mu }}D{P}_0\sqrt{{\displaystyle \frac{\gamma M}{R{T}_1}}{\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma +1}{\gamma -1}}}} \left(P_0>P_{cr}\right).$$

(55)

Figure 3f shows the plot of the $A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}$ values observed at 90° versus Reynolds number (Re). With increasing Re up to 3 × 10⁴, the value of $A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}$ increases very gradually. Above Re = 3 × 10⁴, the value of $A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}$ increases rapidly with increasing Re. The observed $A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}$ values are fitted by an exponential growth $A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}\left(Re\right)=A_{\mathrm{o}}^{Re}{e}^{Re/\tau }+A_0$. As shown in Figure 3f, the observed values are in good agreement with the values calculated by the equation with $A_{\mathrm{o}}^{Re}$ = 5.776×10^-4, τ = 5277 and $A_0=0.1216$ (R² = 0.99553).

We then expanded the angular dependence of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ to four different η positions from 0° to 180° for D0.2. As can be seen in Figure 4, all $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ values agree well with their corresponding $A_{\mathrm{s}\mathrm{i}\mathrm{m}}$ values.

Next, the axial dependence of the AE amplitude was studied as a function of η at $\vartheta =180°$ using five different LSs. Ten sensors were mounted on the cylinder surface at $\vartheta =180°$: eight sensors were mounted downwards and two sensors were mounted upwards relative to the position of the LS. Figure 5 shows the $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ values observed at $\vartheta =180°$ from five different LSs with D0.2–D1.2 and three different internal pressures. Within the standard deviation, all results show that the $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ values decrease exponentially with increasing η. As in the case of the angular dependence, simulations of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ can be performed by normalizing the calculated displacement values and then matching them the observed relative amplitude values in the region of $-0.2 \le \eta /\mathrm{m}\le 0.8$. The simulated angular dependence of the amplitude is formulated as

$$A_{\mathrm{s}\mathrm{i}\mathrm{m}}\left(\eta \right)=u_{\mathrm{N}}\left(\eta \right)\times \left(A_{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}-A_{\mathrm{o},\mathrm{a}\mathrm{x}\mathrm{i}}\right)+A_{\mathrm{o},\mathrm{a}\mathrm{x}\mathrm{i}},$$

(56)

where $A_{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}$ and $A_{\mathrm{o},\mathrm{a}\mathrm{x}\mathrm{i}}$ are the estimated values of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ from extrapolation to η = 0 and η/m = 0.8, respectively. It should be noted that the normalized displacement, $u_{\mathrm{N}}\left(\eta \right)=u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}\left(\eta \right)/u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}\left(0\right)$, is invariant to FRS and frequency. As shown in Figure 5a-e, except the case for D1.2P3 and D1.2P4, all experimental values are in good agreement with the simulated values within the standard deviation. For D1.2P3 and D1.2P4 there is a rather large deviation between $A_{\mathrm{s}\mathrm{i}\mathrm{m}}$ and $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$. Figure 5f shows the plot of the $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$ values versus Re. Similar to the case of the angular dependence, the value of $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$ increases very gradually with increasing Re up to 3 × 10⁴. Above Re = 3 × 10⁴, the value of $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$ increases rapidly with increasing Re. The $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$ values are fitted by an exponential growth $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}\left(Re\right)=A_{\mathrm{o}}^{Re}{e}^{Re/\tau }+A_0$. As shown in Figure 5f, the observed values are in good agreement with the values calculated by the equation with $A_{\mathrm{o}}^{Re}$ = 1.07×10^-3, τ = 5658 and $A_0=0.199$ (R² = 0.972).

For a given value of $\vartheta$, the $u_{leak}$ displacement in Eq. (50) and (51) is exponentially attenuated by the axial distance η via $F_{\xi z}^4$ or $F_{\eta z}^4$ derived from the Green’s function. Note that from Eqs. (4), (23) and (27), $F_{\xi z}^4\left(\eta \right)\propto {\kappa }_{\eta }{e}^{-{\kappa }_{\eta }\left|\eta \right|}$ and $F_{\eta z}^4\left(\eta \right)\propto {\kappa }_{\eta }^2{e}^{-{\kappa }_{\eta }\left|\eta \right|}$. This leads us to propose a semi-empirical equation to fit the observed axial dependence of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ as follows

$$A_{\mathrm{s}\mathrm{i}\mathrm{m}}\left(\eta \right)=A{e}^{-{\kappa }_z\left|\eta \right|}+A_{0,\mathrm{a}\mathrm{x}\mathrm{i}},$$

(57)

where A is the pre-exponential factor, given as $\left(A_{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}-A_{\mathrm{o},\mathrm{a}\mathrm{x}\mathrm{i}}\right)$. This equation reproduces exactly the same values as in the case of Eq. (56); henceforth, we will use Eq. (57) instead of Eq. (56) to analyze the axial dependence of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$.

4.3. Leak Source Localization

A cruciform array of twelve sensors was used to locate the LS. Centered on one sensor, five sensors were mounted at 10 cm intervals along the axial (z) axis and six sensors were mounted at 30° intervals along the tangential (θ) axis (θ ranges from 0 to 180 degrees). Two regions of the cylinder surface were selected for the location of the sensor array, namely the lower and upper regions shown in the first column of Figure 6 (hereafter, referred to as Zone1 and Zone2, respectively). To locate the angular position of the LS, we first estimated the minimum and its position, then shifted Δθ to match the angular dependence of the $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ value with the $A_{\mathrm{s}\mathrm{i}\mathrm{m}}\left(\vartheta \right)$ values calculated by Eq. (53) (Figure 6a and 6c). To locate the axial position, we estimated $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$ by applying the exponential decay to the observed values of $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ and then shifted $\mathsf{\Delta }z$ to match the axial dependence of the $A_{\mathrm{o}\mathrm{b}\mathrm{s}}$ value with the $A_{\mathrm{s}\mathrm{i}\mathrm{m}}\left(\eta \right)$ values calculated by Eq. (57) (Figure 6b and 6d). The matching results for finding the position of the LS location are listed in the second and third columns of Table 2. As listed in Table 2, the error of the LS location is less than 2% for the angular location and less than 1.0% for the axial location. The theoretical model proposed in this study provides a very accurate solution for the LS location. Furthermore, the $A_{0,\mathrm{a}\mathrm{n}\mathrm{g}}$ and $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$ obtained from the localization processes are listed in Table 2. For a given set of experimental conditions, the two values ($A_{0,\mathrm{a}\mathrm{n}\mathrm{g}}$, $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$) agree well within an error of 10%. The Reynolds number was estimated by comparing the $A_{0,\mathrm{a}\mathrm{n}\mathrm{g}}$ and $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$ values obtained from the localization processes with the experimental values shown in Figure 3f and Figure 5f, respectively. As listed in Table 2, the estimated Reynolds number is very close to the theoretical value calculated from the leak hole size and the initial pressure in the cylinder given by Eq. (54). The results suggest that the amplitudes observed from the localization process can provide quantitative information about the hole diameter of the LS if the initial pressure is known.

5. Conclusions

Techniques based on the time difference of arrival (TOD) between two or more AE sensors with the cross-correlation algorithm (CCA) have been widely applied to locate the leak source using the continuous AE waveform data. In practical situations, the AE signal generated by a leak has an aperiodic waveform structure and the exact location of the leak source cannot be determined by the CCA-based methods. A mathematical model applicable to AE due to a pinhole leakage in a cylindrical vessel was summarized, in which the concentrated force-incorporated potentials responsible for the fluctuating Reynolds stress were introduced into the Navier–Lamé equation. In this study, extensive experiments were carried out under different operating conditions to investigate the characteristics of the AE signals due to leakage in a gas cylinder. By validating the angular and the axial dependencies of the AE amplitudes the proposed mathematical model was modified into two semi-empirical equations responsible for the angular and the axial properties of the AE amplitudes, making them well suited for practical applications. Application of the two equations to the data collected from two sets of the cruciform array of 12 sensors mounted on the arbitrary surface of the cylinder demonstrated that the proposed equations were able to accurately localize the leak source in the cylinder within 1% error.

Appendix C.1. Concentrated Force Incorporated Potential

The solutions of the scalar functions, $\varphi , \chi$ and $\psi$ are expressed as

$$\varphi \left(\xi ,\vartheta ,\eta \right)=A_m{J}_m\left(\alpha \xi \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\vartheta \right)e^{{-ik}_{\eta }\eta }-{\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{G_0}{\left({\alpha }^2-{\kappa }_{\eta }^2\right)}}{J}_m\left({\kappa }_{\eta }\xi \right),$$

(C1)

$$\chi \left(\xi ,\vartheta ,\eta \right)=B_m{J}_m\left(\beta \xi \right)\sin \left(\mathit{m\vartheta }\right)e^{{-ik}_{\eta }\eta }-{\displaystyle \frac{k_s^2}{\rho {\omega }^2}}{\displaystyle \frac{G_0}{\left({\beta }^2-{\kappa }_{\eta }^2\right)}}{J}_m\left({\kappa }_{\eta }\xi \right),$$

(C2)

$$\psi \left(\xi ,\theta ,\eta \right)={C_{m}J}_m\left(\beta \xi \right)\cos \left(\mathit{m\vartheta }\right)e^{{-ik}_{\eta }\eta },$$

(C3)

where $k_p={\displaystyle \frac{\omega }{c_p}}$, ${\alpha }^2=k_p^2-k_{\eta }^2$ and $k_s={\displaystyle \frac{\omega }{c_s}}$, ${{\beta }^2=k}_s^2{-k}_{\eta }^2$. In Eqs. (C1)–(C2), A_m, B_m and C_m are coupling constants for P, SH and SV potentials, respectively. The CFIPs can be obtained by substituting each of the scale functions into the corresponding potential.

For Φ

$${\Phi }_{\eta }=P_w\left[\left(-i{k}_{\eta }\right)A_m{J}_m\left(\alpha \xi \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta }-{\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{J_m\left({\kappa }_{\eta }\xi \right)}{{\alpha }^2-{\kappa }_{\eta }^2}}{\displaystyle \frac{\partial G_0}{\partial \eta }}\right],$$

(C4)

$${\Phi }_{\vartheta }=-P_w{\displaystyle \frac{1}{\xi }}\left[A_{m}m{J}_m\left(\alpha \xi \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta }\right].$$

(C5)

For X

$${{\rm X}}_{\xi }=0,$$

(C6)

$$\begin{gathered} {{\rm X}}_{\vartheta }=-P_w{\displaystyle \frac{1}{\xi }}\left\{B_m\left[J_m\left(\beta \xi \right)+\xi {\displaystyle \frac{\partial J_m\left(\beta \xi \right)}{\partial \xi }}\right]\mathrm{s}\mathrm{i}\mathrm{n}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta }\right. \\ \left.-{\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{G_0}{{\beta }^2-{\kappa }_{\eta }^2}}\left[J_m\left({\kappa }_{\eta }\xi \right)+\xi {\displaystyle \frac{\partial J_m\left({\kappa }_{\eta }\xi \right)}{\partial \xi }}\right]\right\}. \end{gathered}$$

(C7)

For Ψ

$${\Psi }_{\xi }=0,$$

(C8)

$${\Psi }_{\vartheta }=P_w{C}_m\left[{\displaystyle \frac{1}{\xi }}{J}_m\left(\beta \xi \right)+{\displaystyle \frac{\partial J_m\left(\beta \xi \right)}{\partial \xi }}\right]\cos \left(\mathit{m\vartheta }\right)e^{-i{k}_{\eta }\eta }.$$

(C9)

For $P=P_z$ (radial component)

$$F_{\xi z}^1=\left(-i{k}_{\eta }\right){\displaystyle \frac{\partial J_m\left(\alpha \xi \right)}{\partial \xi }}\mathrm{c}\mathrm{o}\mathrm{s}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta },$$

(C10a)

$$F_{\xi z}^2=0,$$

(C10b)

$$F_{\xi z}^3=0,$$

(C10c)

$$F_{\xi z}^4=-{\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{1}{{\alpha }^2-{\kappa }_{\eta }^2}}{\displaystyle \frac{\partial J_m\left({\kappa }_{\eta }\xi \right)}{\partial \xi }}{\displaystyle \frac{\partial G_0}{\partial \eta }},$$

(C10d)

(tangential component)

$$F_{\vartheta z}^1=\left(i{k}_{\eta }\right){\displaystyle \frac{m}{\xi }}{J}_m\left(\alpha \xi \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta },$$

(C11a)

$$F_{\vartheta z}^2=0,$$

(C11b)

$$F_{\vartheta z}^3=0,$$

(C11c)

$$F_{\vartheta z}^4=0,$$

(C11d)

(axial component)

$$F_{\eta z}^1=-k_{\eta }^2{J}_m\left(\alpha \xi \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta },$$

(C12a)

$$F_{\eta z}^2=0,$$

(C12b)

$$F_{\eta z}^3=0,$$

(C12c)

$$F_{\eta z}^4=-{\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{J_m\left({\kappa }_{\eta }\xi \right)}{{\alpha }^2-{\kappa }_{\eta }^2}}{\displaystyle \frac{{\partial }^2{G}_0}{{\partial \eta }^2}}.$$

(C12d)

For $P=P_{\theta }$ (radial component)

$$F_{\xi \theta }^1=m\left[{\displaystyle \frac{1}{{\xi }^2}}{J}_m\left(\alpha \xi \right)-{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial J_m\left(\alpha \xi \right)}{\partial \xi }}\right]\mathrm{s}\mathrm{i}\mathrm{n}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta },$$

(C13a)

$$F_{\xi \theta }^2=-m{B}_m\left[{\displaystyle \frac{1}{{\xi }^2}}{J}_m\left(\beta \xi \right)+{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial J_m\left(\beta \xi \right)}{\partial \xi }}\right]\mathrm{c}\mathrm{o}\mathrm{s}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta },$$

(C13b)

$$F_{\xi \theta }^3=-\left(i{k}_{\eta }\right)a\left[{\displaystyle \frac{1}{{\xi }^2}}{J}_m\left(\beta \xi \right)-{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial J_m\left(\beta \xi \right)}{\partial \xi }}-{\displaystyle \frac{{\partial }^2{J}_m\left(\beta \xi \right)}{{\partial \xi }^2}}\right]\cos \left(\mathit{m\vartheta }\right)e^{-i{k}_{\eta }\eta },$$

(C13c)

$$F_{\xi \theta }^4=0,$$

(C13d)

(tangential component)

$$F_{\vartheta \theta }^1=-{\displaystyle \frac{m^2}{{\xi }^2}}{J}_m\left(\alpha \xi \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta },$$

(C14a)

$$F_{\vartheta \theta }^2=\left[-{\displaystyle \frac{1}{{\xi }^2}}{J}_m\left(\beta \xi \right)+{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial J_m\left(\beta \xi \right)}{\partial \xi }}+{\displaystyle \frac{{\partial }^2{J}_m\left(\beta \xi \right)}{{\partial \xi }^2}}\right]\mathrm{s}\mathrm{i}\mathrm{n}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta },$$

(C14b)

$$F_{\vartheta \theta }^3=-{\displaystyle \frac{a}{\xi }}\left(i{k}_{\eta }\right)m\left[{\displaystyle \frac{1}{\xi }}{J}_m\left(\beta \xi \right)+{\displaystyle \frac{\partial J_m\left(\beta \xi \right)}{\partial \xi }}\right]\sin \left(\mathit{m\vartheta }\right)e^{-i{k}_{\eta }\eta },$$

(C14c)

$$F_{\vartheta \theta }^4={\displaystyle \frac{k_p^2}{\rho {\omega }^2}}{\displaystyle \frac{1}{{\beta }^2-{\kappa }_{\eta }^2}}{G}_0\left\{{\displaystyle \frac{1}{{\xi }^2}}\left[J_m\left({\kappa }_{\eta }\xi \right)+{\displaystyle \frac{\partial J_m\left({\kappa }_{\eta }\xi \right)}{\partial \xi }}\right]-{\displaystyle \frac{1}{\xi }}\left[2{\displaystyle \frac{\partial J_m\left({\kappa }_{\eta }\xi \right)}{\partial \xi }}+{\displaystyle \frac{{\partial }^2{J}_m\left({\kappa }_{\eta }\xi \right)}{{\partial }^2\xi }}\right]\right\},$$

(C14d)

(axial component)

$$F_{\eta \theta }^1={\displaystyle \frac{i{k}_{\eta }}{\xi }}m{J}_m\left(\alpha \xi \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\vartheta \right)e^{-i{k}_{\eta }\eta },$$

(C15a)

$$F_{\eta \theta }^2=0,$$

(C15b)

$$\begin{gathered} F_{\eta \theta }^3=a\left[\left({\displaystyle \frac{1-m^2}{{\xi }^3}}\right)J_m\left(\beta \xi \right)-\left({\displaystyle \frac{1+m^2}{{\xi }^2}}\right){\displaystyle \frac{\partial J_m\left(\beta \xi \right)}{\partial \xi }}+{\displaystyle \frac{2}{\xi }}{\displaystyle \frac{{\partial }^2{J}_m\left(\beta \xi \right)}{{\partial \xi }^2}}+{\displaystyle \frac{{\partial }^3{J}_m\left(\beta \xi \right)}{{\partial \xi }^3}}\right] \\ \times \mathrm{c}\mathrm{o}\mathrm{s}\left(m\vartheta \right) e^{{-ik}_{\eta }\eta }, \end{gathered}$$

(C15c)

$$F_{\eta \theta }^4=0.$$

(C15d)

Appendix C.2. Stress-Strain Displacement Relations

The relations between the stress and the strain displacement are given in the $\left(\xi ,\vartheta ,\eta \right)$ coordinates as

$$\begin{gathered} {\sigma }_{\xi \xi f}=c_{11}\left({\displaystyle \frac{\partial u_{\xi f}}{\partial \xi }}\right)+c_{12}\left({\displaystyle \frac{u_{\xi f}}{\xi }}+{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\partial u_{\vartheta f}}{\partial \vartheta }}\right)+c_{13}\left({\displaystyle \frac{\partial u_{\eta f}}{\partial \eta }}\right) \\ =P_f\left(a_{11f}{A}_{mf}+a_{12f}{B}_{mf}+a_{13f}{C}_{mf}+b_{1f}\right)e^{-i\omega t} \\ =0, \end{gathered}$$

(C16a)

$$\begin{gathered} {\sigma }_{\xi \eta f}=c_{44}\left({\displaystyle \frac{\partial u_{\eta f}}{\partial \xi }}+{\displaystyle \frac{\partial u_{\xi f}}{\partial \eta }}\right) \\ =P_f\left(a_{31f}{A}_{mf}+a_{32f}{B}_{mf}+a_{33f}{C}_{mf}+b_{3f}\right)e^{-i\omega t} \\ =0, \end{gathered}$$

(C16b)

where $f=z$ or θ. For m = 0, $B_{0f}=0$. For $f=z$, substituting Eqs. (36) and (38) into Eqs. (C16a) and (C16b), respectively, and for $f=\theta$, substituting Eqs. (39) and (41) into Eqs. (C16a) and (C16b), respectively, we obtain the following two linear algebraic equations.

$$\left[\begin{array}{cc}{a}_{11f}& a_{13f}\\ a_{31f}& a_{33f}\end{array}\right]\left[\begin{array}{c}{A}_{0f}\\ C_{0f}\end{array}\right]=\left[\begin{array}{c}{b}_{1f}\\ b_{3f}\end{array}\right]$$

(C17)