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

Abstract

A theoretical model is presented for the accurate detection of a gas leak source through a pinhole in a cylindrical storage vessel using the acoustic emission (AE) technique. Pinholes of various diameters ranging from 0.20 to 1.2 mm were installed as leak sources, and safe N₂ was used as a filler gas. AE signals were measured and analyzed in terms of AE parameters (such as frequency, amplitude and RMS) as a function of angle and axial distance. Among them, the amplitude characteristic was the most important parameter to determine the leakage dynamics of AE with a continuous waveform. The simulation of AE amplitude was performed using the theoretical model for AE. For practical applications, the theoretical formula was modified into two semi-empirical equations by introducing the normalization method to fit the angular and axial characteristics of the observed AE amplitude, respectively. The main finding of this study is that the semi-empirical equations provide an accurate solution for leak source localization in the cylindrical vessel. As a priori knowledge, the value of ${\kappa }_{\eta }$ in Green’s function, which determines the angular and axial dependence of the AE amplitude, was determined by applying external excitation to the cylinder surface. The proposed formulas provide a suitable approach for practical application in the localization of leak sources in cylindrical storage tanks.

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. However, high-pressure gas cylinders and tanks are plagued by the risk of gas leakage due to container damage. The need for reliable maintenance has increased dramatically in recent years, because the most important way to prevent disasters caused by gas leaks is to detect and repair leaks at an early stage. 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 [2,9,10].

In general, the most notable methods for locating an AE source are the time domain technique and the zone location technique [11]. The former relies on the time difference of arrival (TDOA) between two or more AE sensors to derive the source location [12]. 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 [12,13,14]. Due to background noise and other uncertainties during the measurement, achieving a good level of correlation coefficient is limited [9,15]. 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 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. Basically, a TDOA-based source location is very difficult to apply to the case of low signal-to-noise ratio (SNR) AE signals. Recently, machine learning (ML) and deep learning (DL) techniques have been adapted to solve this problem. Banjarra et al. used support vector machine (SVM) learning and relevance vector machine (RVM) pattern recognition algorithms to classify the healthy and leakage AE data for identification and localization of leakage in pipelines [16]. Zhang et al. modified the original residual neural network (ResNet) DL architecture to overcome the limitations of conventional TDOA-based source localization techniques [17]. Saleem et al. presented a novel acoustic emission (AE)-based pipeline monitoring approach with a customized one-dimensional dense convolutional network (DenseNet) architecture to achieve accurate leak detection and size classification with an average accuracy of 99.7% in energy pipeline networks [18]. They also introduced a hybrid DL technique by integrating the convolutional neural network (CNN) and the long short-term memory (LSTM) models, which transforms time domain AE signals into scalogram images. The technique resulted in accuracy rates of 99.69% [19]. However, the TDOA-based source localization methods have 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 [9,10]. However, there are still limitations and gaps in fundamental research dealing with signal modeling, nonlinear effects and environmental uncertainty in AE signals from leaks [9]. Previously, a mathematical formula was formulated for AE generated by a point source (PS) as an internal defect in cylindrical structures [20]. The PS was treated as a concentrated time harmonic force (CF). In cylindrical geometry, the CF produces three potential functions representing a compression wave (P) and two split shear waves, known as a vertically polarized (SV) wave and a horizontally polarized (SH) wave, which are referred to as the concentrated force-incorporated potential (CFIP). The displacements responsible for the three waves were developed by introducing the CFIPs into the Navier–Lamé (NL) equation. Previously, the CFIPs for the pinhole leakage were derived as an excitation AE source, associated with a fluctuating Reynolds stress (FRS) [21]. As the FRS acts on the pinhole wall, the CFIP only results in axial and tangential displacements due to the FRS. 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 proposed two semi-empirical equations 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 follows:

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

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 Equation (1), the delta function is associated with Green’s function $g\left(\mathit{x};{\mathit{x}}_0\right)$, providing the spatial distribution of the CF at a given time. Since ${\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 the following:

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

or

(2b)$\mathit{f}=\left\{\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]\right\}{e}^{-i\omega t}.$

The derivation of the Green’s function for the cylindrical cell structure [22,23,24] is summarized in Appendix A. Applying the periodicity of 2π ($v$ = 0) to Equations (A18)–(A20) gives the following:

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

(4)$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.,$

(5)$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\}}},$

where the PS is located at $\xi =0,\vartheta =0$ and $\eta =0$. In Equation (3), the variable $\xi$ corresponds to the distance between the PS and a given point. Since the thickness is much shorter than the diameter of the inner circle, the value of $\xi$ can be given as the length of the arc (${\xi }_P)$ around the outer circle (Figure A1):

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

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. [22].) In Refs. [20,22], the value of ${\kappa }_z$, defined as Equation (A4), was determined by applying boundary conditions to the Bessel function $J_0\left({\kappa }_{z}r\right)$ in Equation (A7). 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. [22], 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 Equations (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 follows:

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

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

(8)$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}},$

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 follows:

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

From the reported data [24], the CF can be obtained as follows:

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

where U is the mean velocity of $v_r^{\prime }$ and is given as Equation (A25) in Appendix B.

2.2. Displacement Fields

The displacement fields u that resulted from the CF in the cylindrical geometry can be derived by using the NL equation [25],

(11)$\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}},$

and the Morse and Feshbach’s model [26,27],

(12)$\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),$

where λ and μ are Lamé constants, and ρ is the density of the media. In Equation (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 follows:

(13)$\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 }},$

(14)${\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 },$

(15)$\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 },$

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. [20,23].) The combination of Equations (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 Equations (A26)–(A28) in Appendix C. By substituting each scale function into the corresponding potential, the CFIPs can be obtained as expressed in Equations (A29)–(A34).

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

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

where

(17)$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 }},$

(18)$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 }},$

(19)$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).$

These three equations can be written as follows:

(20)$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},$

where $d=\xi ,\vartheta$ or $\eta$ and $f=z$ or $\theta$ for the pinhole leakage. In Equation (20), the $F_{df}^n$ terms are given in Equations (A35)–(A58) in Appendix C.

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

for $P=P_z$

(radial component),

(21)$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 },$

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

(23)$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 }},$

(tangential component)

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

(axial component)

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

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

(27)$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}},$

for $P=P_{\theta }$,

(radial component)

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

(29)$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 },$

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

(tangential component)

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

(32)$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\},$

(axial component)

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

(34)$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 },$

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

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

For $P=P_z$,

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

(37)$u_{\vartheta z}=0,$

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

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

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

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

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

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

The $P_{\theta }$ excitation does not produce displacements because there are no $F_{d\theta }^4$ terms in Equations (39) and (41) and no coupling constant in Equation (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 follows:

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

Expanded expressions for stress and strain displacement at any point on the cylinder surface are given as Equations (A59) and (A60). Substituting Equations (39) and (41) into Equations (A59) and (A60) with m = 0 gives the following:

(43)$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 },$

(44)$a_{13z}=0,$

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

(46)$a_{33z}=0,$

(47)$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],$

(48)$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 }}.$

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

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

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 follows:

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

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

(50)$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).$

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

(51)$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),$

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 Equation (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 eight-channel acquisition system (IDK-AET/E08, Daejeon, Republic of Korea), including a preamplifier, eight-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 (v. 3) 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, Daejeon, Republic of Korea) 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 (Santa Rose, CA, USA). 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, which was placed horizontally at 100 mm. The actuator was mounted and dismounted at least five times to obtain 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 in Equation (49). From Equation (49), we can obtain the following equation:

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

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

$k_{\eta }=\left\{\begin{array}{ll}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 [22], 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. For the determination of the ${\kappa }_{\eta }$ value 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 the σ 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 Equation (3). From Equation (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 becomes 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 of 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: nine sensors in the downward direction and one 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 Equation (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 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 performed 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 follows:

(53)$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}}.$

In Equation (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 Equation (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 symmetrical 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 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 Equation (10), the CF due to the leakage is only proportional to the mean velocity U. As expressed in Equations (A23) and (A24) in Appendix B [28,29,30,31], substituting Q into U cancels the hole diameter factor D in Equation (A25). Accordingly, 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 follows:

(54)$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}\left(P_0<P_{cr}\right),$

(55)$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}}}}\hspace{1em}\left(P_0>P_{cr}\right).$

Figure 3f shows the plot of the $A_{\mathrm{o},\mathrm{a}\mathrm{n}\mathrm{g}}$ values observed at 90° versus the 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⁻⁴, τ = 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 to 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 follows:

(56)$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}},$

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⁻³, τ = 5658 and $A_0=0.199$ (R² = 0.972).

For a given value of $\vartheta$, the $u_{leak}$ displacement in Equations (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 Equations (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:

(57)$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}},$

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 Equation (56); henceforth, we will use Equation (57) instead of Equation (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 and 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 Equation (53) (Figure 6a,c). 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 $\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 Equation (57) (Figure 6b,d). 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 Equation (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 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 a 1% error. The research highlights include the following: (1) the predetermination of the Green’s function parameter ${\kappa }_{\eta }$ through external stimuli to enhance the universality of the model; (2) the analysis of the angular/axial attenuation law of AE signals by establishing normalized simulation equations; and (3) the three-dimensional localization of the leak accurately achieved using a cruciform sensor array. The model provides a new approach to the safety monitoring of industrial storage tanks as well as gas transmission pipelines installed in the interstitial spaces of structures, especially for early warning of hazardous gas leaks, which has important practical implications.

Author Contributions

Conceptualization, data curation, methodology, investigation and writing—original draft, J.-G.K.; software, validation and formal analysis, K.B.K.; resources, data curation and funding acquisition, K.H.K.; writing—review and editing and visualization, K.B.K., K.H.K. and B.K.K.; supervision and project administration, B.K.K. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

All authors were employed by the company Integrity Diagnostics Korea (IDK). All authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Schematic diagram of the main parts of the leak test apparatus (h: pinhole, p: AE sensor (IDK-AES-H150), PA: preamplifier, PC/DAQ (IDK-AET/E08): PC with 8-channel DAQ).

Figure 2 (a) Angular dependence of the mean amplitudes (circles) and σs (bars) of signal obtained by an external excitation (1.000 V_pp and $\nu =150$ kHz) and the simulated $u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}$ displacements (lines) at different values of ${\kappa }_{\eta }$ (1: 6.0, 2: 7.0, 3: 7.5, 4: 8.0, 5: 9.0); (b) axial dependence of the mean amplitudes of signal and the simulated $u_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}$ displacement (line) with ${\kappa }_{\eta }$ = 7.5.

Figure 3 (a–e) Observed mean maximum amplitudes (circles) and σs (bars) obtained from five different hole diameters: (a) D0.2η10, (b) D0.3η10, (c) D0.5η10, (d) D0.8η10 and (e) D1.2η10 at $P_0$ = 2 bar (black), 3 bar (red) and 4 bar (green), simulated amplitudes (lines) and (f) plot of mean amplitudes (circle) at $\vartheta =90$° vs. Reynolds number and the fit (line) by one-phase exponential growth.

Figure 4 Observed mean maximum amplitudes (circles) and σs (bars) obtained from four different axial positions: (a) D0.2η3, (b) D0.2η30, (c) D0.2η50 and (d) D0.2η70 at $P_0$ = 2 bar (black), 3 bar (red) and 4 bar (green) and simulated amplitudes (lines).

Figure 5 (a–e) Observed mean maximum amplitudes (circles) and σs (bars) obtained from five different hole diameters: (a) D0.2$\vartheta$180, (b) D0.3$\vartheta$180, (c) D0.5$\vartheta$180, (d) D0.8$\vartheta$180 and (e) D1.2$\vartheta$180 at P = 2 bar (black), 3 bar (red) and 4 bar (green), simulated amplitudes (lines) and (f) plot of $A_{0,\mathrm{a}\mathrm{x}\mathrm{i}}$ values (circles) vs. Reynolds NO and the fit (line) by one-phase exponential growth.

Figure 6 Fitting the angular (a,c) and the axial (b,d) dependences of the observed mean maximum amplitudes (circles) and σs (bars) obtained from the sensor array in Zone1 (a,b) and Zone2 (c,d) to the simulated values (solid lines) (three different gas pressures used).

Table 1

Material and geometric properties of the leak test cylinder, and molecular and thermal properties of the charged gas.

Steel Cylinder
$\mathrm{Outer}\mathrm{radius},a=0.116$ m,$\mathrm{Density},\rho =7.8\times {10}^3$ kg m−3, $\mathrm{P}\mathrm{wave}\mathrm{velocity},c_P=5.98\times {10}^3$ m s−1,			$\mathrm{Inner}\mathrm{radius},b=0.1112$ m$\mathrm{Length},L=1.013$ m$\mathrm{S}\mathrm{wave}\mathrm{velocity},c_S=3.30\times {10}^3$ m s−1
Stiffness parameters (kg m−1 ms−2)
c11 = 208.6 × 109,c12 = 146.5 × 109,c44 = 126.9 × 109,		c22 = 208.6 × 109,c13 = 146.5 × 109,c55 = 126.9 × 109,		c33 = 208.6 × 109,c23 = 146.5 × 109,c66 = 126.9 × 109
Nitrogen gas (25 °C, 1 bar) *
Molecular weight, M = 0.028 kg mol−1$\mathrm{Dynamic}\mathrm{viscosity},\mu =1.782\times {10}^{-5}$ Pa s $\mathrm{Kinematic}\mathrm{density},\nu =1.573\times {10}^{-5}$ m2 s−1			$\mathrm{Density},\rho =1.13$ kg m−3,$\mathrm{Isentropic}\mathrm{index},\gamma =1.40$,

* https://www.engineeringtoolbox.com/nitrogen-N2-dynamic-kinematic-viscosity-temperature-pressure-d_2067.html (accessed on 3 April 2025).

Table 2

Results for the localization of the leak source determined by fitting the observed signal amplitudes to the simulated values. The leak source is located at ϑ = 0° and η = 0.

ArrayPosition	P0(bar)	$\mathit{\vartheta }$(error%)	$\mathit{\eta }$/mm(error%)	${\mathit{A}}_{\mathbf{o},\mathbf{a}\mathbf{n}\mathbf{g}}$ $/{\mathit{A}}_{\mathbf{o},\mathbf{a}\mathbf{x}\mathbf{i}}$	Reynolds NO/kTheory/Estimation
Zone1 a	2.07	−6.1° (1.7)	5.0 (0.5)	0.05/0.05	5.3/5.5–5.9
	3.06	−6.1° (1.7)	5.0 (0.5)	0.09/0.08	7.2/7.0–8.4
	3.97	−6.1° (1.7)	5.0 (0.5)	0.13/0.12	8.5/8.5–11.7
Zone2 b	2.01	3.6° (1.0)	7.0 (0.7)	0.12/0.11	13.3/8.0–11.0
	3.01	3.6° (1.0)	7.0 (0.7)	0.26/0.25	18.1/18.1–23.0
	4.01	3.6° (1.0)	7.0 (0.7)	0.40/0.38	21.5/28.9

^a,b Two different hole sizes were used: a, D0.2; b, D0.5

Appendix A. Green’s Function

In cylindrical $\left(r,\theta ,z\right)$ coordinates, the relationship between the Green’s and the delta functions is expressed as follows:(A1)$${\nabla }^{2}g\left(\mathit{r},\mathit{\theta },\mathit{z};{\mathit{r}}_0,{\mathit{\theta }}_0,{\mathit{z}}_0\right)={\displaystyle \frac{\delta \left(r-r_0\right)\delta \left(\theta -{\theta }_0\right)\delta \left(z-z_0\right)}{r}},$$ where the Green’s function is given by the following:(A2)$$g\left(\mathit{r},\mathit{\theta },\mathit{z};{\mathit{r}}_0,{\mathit{\theta }}_0,{\mathit{z}}_0\right)=A_v{g}_r\left(r;r_0\right)g_{\theta }\left(\theta ;{\theta }_0\right)g_z\left(z;z_0\right),$$ where $A_v$ is a coupling constant, which will be specified later by an azimuthal constant $v\left(=0,\pm 1,\pm 2,\cdots \right)$.

At any point in the cylindrical domain other than the PS locating point, the values of the delta functions are zero:$${\nabla }^2\left[g_r\left(r;r_0\right)g_{\theta }\left(\theta ;{\theta }_0\right)g_z\left(z;z_0\right)\right]=0,$$ and Equation (A1) can be rewritten as follows:(A3)$${\displaystyle \frac{r^2}{g_r}}{\displaystyle \frac{{\partial }^2{g}_r}{{\partial r}^2}}+{\displaystyle \frac{r}{g_r}}{\displaystyle \frac{\partial g_r}{\partial r}}+{\displaystyle \frac{1}{g_{\theta }}}{\displaystyle \frac{{\partial }^2{g}_{\theta }}{{\partial \theta }^2}}+{\displaystyle \frac{r^2}{g_z}}{\displaystyle \frac{{\partial }^2{g}_z}{\partial z^2}}=0.$$

Introducing ${\kappa }_z^2$ for the real roots of $g_z$ and $v$ as an azimuthal constant for $g_{\theta }$ into Equation (A3) gives three individual differential equations for $g_z$, $g_{\theta }$ and $g_r$.(A4)$${\displaystyle \frac{1}{g_z}}{\displaystyle \frac{{\partial }^2{g}_z}{{\partial z}^2}}={\kappa }_z^2,$$$${\displaystyle \frac{1}{g_{\theta }}}{\displaystyle \frac{{\partial }^2{g}_{\theta }}{{\partial \theta }^2}}=-v^2,$$$$r^2{\displaystyle \frac{{\partial }^2{g}_r}{{\partial r}^2}}+r{\displaystyle \frac{\partial g_r}{\partial r}}+\left({\kappa }_z^2{r}^2-v^2\right)g_r=0.$$

The solutions of the three functions can be obtained by applying the conditions for the delta function, such as the continuity, the discontinuity and/or the symmetry principles around the PS.(A5)$$g_z\left(z;z_0\right)=\left\{\begin{array}{ll}-{\displaystyle \frac{1}{2{\kappa }_z}}{e}^{-{\kappa }_z(z_0-z)}& 0<z<z_0<l\\ -{\displaystyle \frac{1}{2{\kappa }_z}}{e}^{-{\kappa }_z(z-z_0)}& 0<z_0<z<l\end{array}\right.,$$(A6)$$g_{\theta }\left(\theta ;{\theta }_0\right)=\cos \left[v\left(\theta -{\theta }_0\right)\right](v=0,\pm 1,\pm 2,\cdots ),$$

(A7) $$g_r\left(r;r_0\right)=J_v\left({\kappa }_z{r}_0\right)J_v\left({\kappa }_{z}r\right)\left(b<r<a\right).$$

The coefficient $A_v$ can be obtained by substituting Equations (A5)–(A7) into Equation (A1) and integrating over a given cylindrical domain $\left(0\le z\le l,0\le \theta \le 2\pi ,b\le r\le a\right)$. Since ${\int }_0^l\delta \left(z-z_0\right)dz=1$, and(A8)$${\left.{\displaystyle \frac{\partial g_z}{\partial z}}\right|}_0^{z=l}={\left.-{\displaystyle \frac{e^{-{\kappa }_z{z}_0}}{2}}{e}^{kz}\right|}_0^{z_0}+{\left.{\displaystyle \frac{e^{{{\kappa }_{z}z}_0}}{2}}{e}^{-{\kappa }_{z}z}\right|}_{z_0}^l={\displaystyle \frac{1}{2}}\left[e^{-{\kappa }_z{z}_0}+e^{-{\kappa }_z\left(l-z_0\right)}-2\right],$$ then(A9)$$\begin{gathered} A_v{\int }_0^l\left\{{\displaystyle \frac{g_{\theta }{g}_z}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial g_r}{\partial r}}\right)+{\displaystyle \frac{g_r{g}_z}{r^2}}{\displaystyle \frac{{\partial }^2{g}_{\theta }}{{\partial \theta }^2}}+{g_{r}g}_{\theta }{\displaystyle \frac{{\partial }^2{g}_z}{\partial z^2}}\right\}dz \\ =A_v{g_{r}g}_{\theta }{\int }_0^l{\displaystyle \frac{{\partial }^2{g}_z}{\partial z^2}}dz \\ ={\displaystyle \frac{A_v}{2}}\left[e^{-{\kappa }_z{z}_0}+e^{-{\kappa }_z\left(l-z_0\right)}-2\right]{g_{r}g}_{\theta } \end{gathered}$$

Similarly, since ${\int }_0^{2\pi }\delta \left(\theta -{\theta }_0\right)d\theta =1$,(A10)$${\int }_0^{2\pi }{g}_{\theta }d\theta ={\int }_0^{2\pi }\cos \left[v\left(\theta -{\theta }_0\right)\right]d\theta =\left\{\begin{array}{c}2\pi \hspace{1em}\mathrm{i}\mathrm{f} \theta ={\theta }_0\\ 0\hspace{1em}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right..$$

From Equations (A1), (A9) and (A10), we obtained the following:(A11)$$A_v{\pi g}_r\left[e^{-{\kappa }_z{z}_0}+e^{-{\kappa }_z\left(l-z_0\right)}-2\right]={\displaystyle \frac{\delta \left(r-r_0\right)}{r}}.$$

Substituting Equation (A7) into Equation (A11) gives the following:(A12)$$\pi A_v{J}_v\left({\kappa }_z{r}_0\right)J_v\left({\kappa }_{z}r\right)\left[e^{-{\kappa }_z{z}_0}+e^{-{\kappa }_z\left(l-z_0\right)}-2\right]={\displaystyle \frac{\delta \left(r-r_0\right)}{r}}.$$

Multiplying both sides of Equation (A12) by ${rJ}_v\left({\kappa }_{z}r\right)$ and integrating over (b, a) gives the solution for $\delta \left(r-r_0\right)$ as follows:(A13)$$\begin{gathered} \pi \left[e^{-{\kappa }_z{z}_0}+e^{-{\kappa }_z\left(l-z_0\right)}-2\right]J_v\left({\kappa }_z{r}_0\right)A_v{\int }_b^a{J}_v\left({\kappa }_{z}r\right)J_v\left({\kappa }_{z}r\right)rdr \\ ={\int }_b^a\delta \left(r-r_0\right)J_p\left({\kappa }_{z}r\right)dr \\ =J_v\left(k_z{r}_0\right). \end{gathered}$$

Applying the normalization and the orthogonality of Bessel function,$${\int }_0^c{J}_p\left(\alpha {\displaystyle \frac{r}{c}}\right)J_q\left(\alpha {\displaystyle \frac{r}{c}}\right)rdr={\displaystyle \frac{c^2}{2}}{\left[J_{p+1}\left(\alpha \right)\right]}^2{\delta }_{pq},$$ to the integration of the Bessel function in Equation (A13) gives the following:(A14)$$\begin{gathered} {\int }_b^a{J}_v\left({\kappa }_{z}r\right)J_v\left({\kappa }_{z}r\right)rdr \\ ={\int }_0^a{J}_v\left({\kappa }_{z}r\right)J_v\left({\kappa }_{z}r\right)rdr-{\int }_0^b{J}_v\left({\kappa }_{z}r\right)J_v\left({\kappa }_{z}r\right)rdr \\ ={\displaystyle \frac{a^2}{2}}{\left[J_{v+1}\left({\kappa }_{z}a\right)\right]}^2-{\displaystyle \frac{b^2}{2}}{\left[J_{v+1}\left({\kappa }_{z}b\right)\right]}^2. \end{gathered}$$

Substituting Equation (A14) into Equation (A12) gives the following:$$A_v={\displaystyle \frac{1}{\pi \left[e^{-k_z{z}_0}+e^{-k_z\left(l-z_0\right)}-2\right]}}\times {\displaystyle \frac{2}{\left\{a^2{\left[J_{v+1}\left({\kappa }_{z}a\right)\right]}^2-b^2{\left[J_{v+1}\left({\kappa }_{z}b\right)\right]}^2\right\}}}.$$

Finally, the Green’s function can be obtained as follows:(A15)$$g\left(\mathit{r},\mathit{\theta },\mathit{z};{\mathit{r}}_0,{\mathit{\theta }}_0,{\mathit{z}}_0\right)=G_v{J}_v\left({\kappa }_{z}r\right)\cos \left[v\left(\theta -{\theta }_0\right)\right],$$(A16)$$G_v=-{\displaystyle \frac{1}{2{\kappa }_z}}{A}_v{J}_v\left({\kappa }_z{r}_0\right)\times \left\{\begin{array}{c}{e}^{-{\kappa }_z\left(z_0-z\right)}\hspace{1em}0<z<z_0<l\\ e^{-{\kappa }_z\left(z-z_0\right)}\hspace{1em}0<z_0<z<l\end{array}\right.,$$(A17)$$A_v={\displaystyle \frac{1}{\pi \left[e^{{-k}_z{z}_0}+e^{-k_z\left(l-z_0\right)}-2\right]}}\times {\displaystyle \frac{2}{\left\{a^2{\left[J_{v+1}\left({\kappa }_{z}a\right)\right]}^2-b^2{\left[J_{v+1}\left({\kappa }_{z}b\right)\right]}^2\right\}}}.$$

In addition, replacing the $(r,\theta ,z)$ coordinates by relative coordinates, defined as $\xi =r-r_0$, $\vartheta =\theta -{\theta }_0$ and $\eta =z-z_0$, sets the location of the PS as the new origin at ${\xi }_0=0,{\vartheta }_0=0$ and ${\eta }_0=0$ (Figure A1). In the $\left(\xi ,\vartheta ,\eta ;0\right)$ coordinates, Equations (A15)–(A17) are rewritten as follows:(A18)$$g\left(\mathit{\xi },\vartheta ,\eta ;0\right)=G_v\left(\eta \right)J_v\left({\kappa }_{\eta }\xi \right)\cos v\vartheta ,$$(A19)$$G_v=-{\displaystyle \frac{1}{2{\kappa }_z}}{A}_v{J}_v\left({\kappa }_{\eta }{\xi }_0\right)\times \left\{\begin{array}{c}{e}^{{\kappa }_{\eta }\eta }\hspace{1em}\left(0\le \eta \le l-z_0\right)\\ e^{-{\kappa }_{\eta }\eta }\hspace{1em}\left(-z_0\le \eta \le 0\right)\end{array}\right.,$$(A20)$$A_v={\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_{v+1}\left({\kappa }_{\eta }a\right)\right]}^2-b^2{\left[J_{v+1}\left({\kappa }_{\eta }b\right)\right]}^2\right\}}}.$$

Figure A1 A radial cross-section containing PS and a given point S. The red line represents the arc connecting PS and S.

Appendix B. Leak Characteristics

When the gas is leaking from a pressurized container into atmospheric air, the Reynolds number (Re) for the leakage flow will be much higher than for the atmospheric air leaking into an evacuated vessel and will reach into the turbulent region. For the turbulent flow through a small hole, the flow characteristic must be averaged over a period of time. The Reynolds number for the small hole model can be described as follows:(A21)$$Re={\displaystyle \frac{\rho UD}{\mu }},$$ where $U$ is the mean leakage velocity, D is the diameter of a hole and μ is the dynamic viscosity of the fluid.

The mass flow rate Q of the gas at the leak hole is divided into sonic flow and subsonic flow on the basis of the critical pressure ratio (CPR):(A22)$$\mathrm{C}\mathrm{P}\mathrm{R}={\displaystyle \frac{P_a}{P_{cr}}}={\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma }{\gamma -1}}},$$(A23)$$Q=C_D{A}_{lh}{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}{P}_0<P_{cr},$$(A24)$$Q=C_D{A}_{lh}{P}_0\sqrt{{\displaystyle \frac{\gamma M}{R{T}_1}}{\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma +1}{\gamma -1}}}}\hspace{1em}{P}_0\ge P_{cr},$$ where $P_a$ is the ambient atmospheric pressure, $P_{cr}$ is the critical pressure, $P_0$ is the charged gas pressure, $A_{lh}$ is the leakage hole area, $C_D$ is the flow correction factor of the leakage hole (0.6~1.0), γ is the isentropic index and M is the molar mass of the gas. For nitrogen gas, CPR = 0.528. When $p_a=1.01$ bar, $p_{cr}=1.91$ bar.

Since the mass flow rate of the gas at the leakage hole is given as $Q=\rho A_{lh}U$, the mean flow velocity can be given by the following:(A25)$$U={\displaystyle \frac{Q}{\rho A_{lh}}}={\displaystyle \frac{4Q}{\rho \pi D^2}}.$$

Appendix C

Appendix C.1. Concentrated Force Incorporated Potential

The solutions of the scalar functions, $\varphi ,\chi$ and $\psi$ are expressed as follows:(A26)$$\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),$$(A27)$$\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),$$(A28)$$\psi \left(\xi ,\theta ,\eta \right)={C_{m}J}_m\left(\beta \xi \right)\cos \left(\mathit{m\vartheta }\right)e^{{-ik}_{\eta }\eta },$$ 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 Equations (A26) and (A27), 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 Φ,(A29)$${\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],$$(A30)$${\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].$$

For X,(A31)$${{\rm X}}_{\xi }=0,$$(A32)$$\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}$$

For Ψ,(A33)$${\Psi }_{\xi }=0,$$(A34)$${\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 }.$$

For $P=P_z$,

(radial component)(A35)$$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 },$$(A36)$$F_{\xi z}^2=0,$$(A37)$$F_{\xi z}^3=0,$$(A38)$$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 }},$$

(tangential component)(A39)$$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 },$$(A40)$$F_{\vartheta z}^2=0,$$(A41)$$F_{\vartheta z}^3=0,$$(A42)$$F_{\vartheta z}^4=0,$$

(axial component)(A43)$$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 },$$(A44)$$F_{\eta z}^2=0,$$(A45)$$F_{\eta z}^3=0,$$(A46)$$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}}.$$

For $P=P_{\theta }$,

(radial component)(A47)$$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 },$$(A48)$$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 },$$(A49)$$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 },$$(A50)$$F_{\xi \theta }^4=0,$$

(tangential component)(A51)$$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 },$$(A52)$$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 },$$(A53)$$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 },$$(A54)$$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\},$$

(axial component)(A55)$$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 },$$(A56)$$F_{\eta \theta }^2=0,$$(A57)$$\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}$$(A58)$$F_{\eta \theta }^4=0.$$

Appendix C.2. Stress–Strain Displacement Relations

The relations between stress and strain displacement are given in the $\left(\xi ,\vartheta ,\eta \right)$ coordinates as follows:(A59)$$\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}$$(A60)$$\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}$$ where $f=z$ or θ. For m = 0, $B_{0f}=0$. For $f=z$, substituting Equations (36) and (38) into Equations (A59) and (A60), respectively, and for $f=\theta$, substituting Equations (39) and (41) into Equations (A59) and (A60), respectively, we obtain the following two linear algebraic equations:(A61)$$\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]$$

References

1. Korlapati, N.V.S.; Khan, F.; Noora, Q.; Mirza, S.; Vaddiraju, S. Review and Analysis of Pipeline Leak Detection Methods. J. Pipeline Sci. Eng.; 2022; 2, 100074. [DOI: https://dx.doi.org/10.1016/j.jpse.2022.100074]

2. Meng, L.; Li, Y.; Wang, W.; Fu, J. Experimental Study on Leak Detection and Location for Gas Pipeline Based on Acoustic Method. J. Loss Prev. Process Ind.; 2012; 25, pp. 90-102. [DOI: https://dx.doi.org/10.1016/j.jlp.2011.07.001]

3. Hao, Y.; Du, Z.; Jiang, J.; Xing, Z.; Yan, X.; Wang, S.; Rao, Y. Research on Multipoint Leak Location of Gas Pipeline Based on Variational Mode Decomposition and Relative Entropy. Shock Vib.; 2020; 2020, 8868963. [DOI: https://dx.doi.org/10.1155/2020/8868963]

4. Quy, T.B.; Kim, J.-M. Pipeline Leak Detection Using Acoustic Emission and State Estimate in Feature Space. IEEE Trans. Instrum. Meas.; 2022; 71, 2518709. [DOI: https://dx.doi.org/10.1109/TIM.2022.3206833]

5. Nguyen, D.-T.; Nguyen, T.K.; Ahmad, Z.; Kim, J.-M. A Reliable Pipeline Leak Detection Method Using Acoustic Emission with Time Difference of Arrival and Kolmogorov–Smirnov Test. Sensors; 2023; 23, 9296. [DOI: https://dx.doi.org/10.3390/s23239296]

6. Xiao, R.; Joseph, P.; Li, J. The Leak Noise Spectrum in Gas Pipeline Systems: Theoretical and Experimental Investigation. J. Sound Vib.; 2020; 488, 115646. [DOI: https://dx.doi.org/10.1016/j.jsv.2020.115646]

7. Tian, X.; Jiao, W.; Liu, T. Intelligent Leak Detection Method for Low-Pressure Gas Pipeline inside Buildings Based on Pressure Fluctuation Identification. J. Civ. Struct Health Moni.; 2022; 12, pp. 1191-1208. [DOI: https://dx.doi.org/10.1007/s13349-022-00607-y]

8. Meribout, M. Gas Leak-detection and Measurement Systems: Prospects and Future trends. IEEE Trans. Inst. Meas.; 2021; 70, 4505813. [DOI: https://dx.doi.org/10.1109/TIM.2021.3096596]

9. Hu, Z.; Tariq, S.; Zayed, T. A Comprehensive Review of Acoustic Based Leak Localization Method in Pressurized Pipelines. Mech. Syst. Signal Process.; 2021; 161, 107994. [DOI: https://dx.doi.org/10.1016/j.ymssp.2021.107994]

10. Liu, Y.; Habibi, D.; Chai, D.; Wang, X.; Chen, H.; Gao, Y.; Li, S. A Comprehensive Review of Acoustic Methods for Locating Underground Pipelines. Appl. Sci.; 2020; 10, 1031. [DOI: https://dx.doi.org/10.3390/app10031031]

11. Miller, R.K.; Carlos, M.F.; Finlayson, R.D.; Godinez-Azcuaga, V.; Rhodes, M.R.; Shu, F.; Wang, W.D. Acoustic Emission Source Location. Nondestructive Testing Handbook; 3rd ed. Miller, R.K.; Hill, E.K.; Moore, P.O. American Society for Nondestructive Testing: Columbus, OH, USA, 2005; pp. 121-130.

12. Kund, T. Mechanics of Elastic Waves and Ultrasonic Nondestructive Evaluation; CRC Press: Boca Raton, FL, USA, 2019; pp. 317-369.

13. Gao, Y.; Brennan, M.J.; Joseph, P.F.; Muggleton, J.M.; Hunaidi, O. A Model of the Correlation Function of Leak Noise in Buried Plastic Pipes. J. Sound Vib.; 2004; 277, pp. 133-148. [DOI: https://dx.doi.org/10.1016/j.jsv.2003.08.045]

14. Kafle, M.D.; Fong, S.; Narasimhan, S. Active Acoustic Leak Detection and Localization in a Plastic Pipe Using Time Delay Estimation. Appl. Acoustics; 2022; 187, 108482. [DOI: https://dx.doi.org/10.1016/j.apacoust.2021.108482]

15. Xiao, R.; Joseph, P.F.; Muggleton, J.M.; Li, J. Limits for Leak Noise Detection in Gas Pipes Using Cross Correlation. J. Sound Vib.; 2022; 520, 116639. [DOI: https://dx.doi.org/10.1016/j.jsv.2021.116639]

16. Banjara, N.K.; Sasmal, S.; Voggu, S. Machine Learning Supported Acoustic Emission Technique for Leakage Detection in Pipelines. Int. J. Press. Vessels Pip.; 2020; 188, 104243. [DOI: https://dx.doi.org/10.1016/j.ijpvp.2020.104243]

17. Zhang, T.; Xu, C.; Ozevin, D. Acoustic Emission Source Localization in Complex Pipe Structure Using Multi-task Deep Learning Models. Adv. Struct. Eng.; 2025; 28, pp. 23-37. [DOI: https://dx.doi.org/10.1177/13694332241269250]

18. Saleem, F.; Ahmad, Z.; Siddique, M.F.; Umar, M.; Kim, J.-M. Acoustic Emission-Based Pipeline Leak Detection and Size Identification Using a Customized One-Dimensional DenseNet. Sensors; 2025; 25, 1112. [DOI: https://dx.doi.org/10.3390/s25041112]

19. Ullah, N.; Siddique, M.F.; Ullah, S.; Ahmad, Z.; Kim, J.-M. Pipeline Leak Detection System for a Smart City: Leveraging Acoustic Emission Sensing and Sequential Deep Learning. Smart Cities; 2024; 7, pp. 2318-2338. [DOI: https://dx.doi.org/10.3390/smartcities7040091]

20. Kim, K.B.; Kim, B.K.; Kang, J.-G. Modeling Acoustic Emission Due to an Internal Point Source in Circular Cylindrical Structures. Appl. Sci.; 2022; 12, 12032. [DOI: https://dx.doi.org/10.3390/app122312032]

21. Morse, P.M.; Feshbach, H. Methods of Theoretical Physics; McGraw-Hill: New York, NY, USA, 1953; pp. 1764-1767.

22. Kim, K.B.; Kim, J.-H.; Jin, J.-E.; Kim, H.-J.; Kim, C.-I.; Kim, B.K.; Kang, J.-G. The Characteristics of Acoustic Emissions due to Gas Leaks in Circular Cylinders: A Theoretical and Experimental Investigation. Appl. Sci.; 2023; 13, 9814. [DOI: https://dx.doi.org/10.3390/app13179814]

23. Kim, K.B.; Kim, B.K.; Lee, S.G.; Kang, J.-G. Analytical Modeling of Acoustic Emission due to an Internal Point Source in a Transversely Isotropic Cylinder. Appl. Sci.; 2022; 12, 5272. [DOI: https://dx.doi.org/10.3390/app12105272]

24. Lighthill, M.J. On Sound Generated Aerodynamics. I. General Theory. Proc. R. Soc. A; 1952; 211, pp. 564-587.

25. Brennen, C.E. An Internet Book on Fluid Dynamics. 2006; Available online: http://brennen.caltech.edu/fluidbook/basicfluiddynamics/turbulence.htm (accessed on 20 August 2023).

26. Eringen, A.C.; Suhubi, E.S. Fundamentals of Linear Elastodynamics. Elastodynamics; Academic Press: New York, NY, USA, 1975; Volume II, pp. 343-366.

27. Honarvar, F.; Enjilel, E.; Sinclair, A.N.; Mirnezami, S.A. Wave Propagation in Transversely Isotropic Cylinders. Int. J. Solids Struct.; 2007; 44, pp. 5236-5246. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2006.12.029]

28. Bomelburg, H.J. Estimation of Gas Leak Rates Through Very Small Orifices and Channels; Report for Nuclear Regulatory Commission. BNWL-2223-77 Pacific Northwest Laboratory: Richland, WA, USA, 1977.

29. Levenspiel, O. Engineering Flow and Heat Exchange; 3rd ed. Springer Science: New York, NY, USA, 2014; pp. 52-54.

30. Bariha, N.; Mishra, I.M.; Srivastava, V.C. Hazard Analysis of Failure of Natural Gas and Petroleum Gas Pipelines. J. Loss Prev. Process. Ind.; 2016; 40, pp. 217-226. [DOI: https://dx.doi.org/10.1016/j.jlp.2015.12.025]

31. Li, Y.; Zhou, P.; Zhuang, Y.; Wu, X.; Liu, Y.; Han, X.; Chen, G. An Improved Gas Leakage Model and Research on the Leakage Field Strength Characteristics of R290 in Limited Space. Appl. Sci.; 2022; 12, 5657. [DOI: https://dx.doi.org/10.3390/app12115657]