**Prediction of HIFU Propagation in a Dispersive Medium via Khokhlov–Zabolotskaya–Kuznetsov Model Combined with a Fractional Order Derivative**

**Shilei Liu <sup>1</sup> ID , Yanye Yang 1, Chenghai Li 1, Xiasheng Guo 1,\*, Juan Tu <sup>1</sup> and Dong Zhang 1,2,\* ID**


Received: 17 March 2018; Accepted: 10 April 2018; Published: 12 April 2018

**Abstract:** High intensity focused ultrasound (HIFU) has been proven to be promising in non-invasive therapies, in which precise prediction of the focused ultrasound field is crucial for its accurate and safe application. Although the Khokhlov–Zabolotskaya–Kuznetsov (KZK) equation has been widely used in the calculation of the nonlinear acoustic field of HIFU, some deviations still exist when it comes to dispersive medium. This problem also exists as an obstacle to the Westervelt model and the Spherical Beam Equation. Considering that the KZK equation is the most prevalent model in HIFU applications due to its accurate and simple simulation algorithms, there is an urgent need to improve its performance in dispersive medium. In this work, a modified KZK (mKZK) equation derived from a fractional order derivative is proposed to calculate the nonlinear acoustic field in a dispersive medium. By correcting the power index in the attenuation term, this model is capable of providing improved prediction accuracy, especially in the axial position of the focal area. Simulation results using the obtained model were further compared with the experimental results from a gel phantom. Good agreements were found, indicating the applicability of the proposed model. The findings of this work will be helpful in making more accurate treatment plans for HIFU therapies, as well as facilitating the application of ultrasound in acoustic hyperthermia therapy.

**Keywords:** KZK equation; fractional order derivative; ultrasound hyperthermia; HIFU; acoustic simulation; Kramers–Kronig relation

#### **1. Introduction**

Although pioneering clinical studies of focused ultrasound were carried out as early as the 1940s [1,2], it did not attract intensive research interest until the end of the 20th century and the beginning of 21st century, during which several theoretical models were developed, improved and then broadly accepted [3–9]. In the past decades, high intensity focused ultrasound (HIFU) has played an increasingly significant role in the study of non-invasive therapies by demonstrating unique advantages in safety, effectiveness and high efficiency [10–14]. However, the applications of HIFU are still limited, and clinical treatments are only available for limited sites [11,14–16].

One of the challenges confronting HIFU treatments is the spatial precision of tissue ablation. Several techniques such as ultrasound B-Scan and Magnetic Resonance Imaging (MRI) have been combined with HIFU to achieve real-time monitoring of focal areas [17–20]. With these methods, the actual focal profiles were usually found to deviate from those predicted through theoretical models [21,22]. As was indicated by Petrusca et al., the shift of the focal point away from the prescribed

position caused by acoustic aberrations and non-linear wave propagating effects make it mandatory to evaluate the spatial accuracy of HIFU ablation [21]. Li et al. observed a focal shift of 1–2 mm and ascribed it to the layered distribution of tissues. Several different mechanisms might contribute to these deviations, such as the thermos-lensing effect [23,24], bubble formation [25–27], acoustic radiation force [27] and the nonlinear nature of acoustic waves [21,22]. Connor et al. pointed out that the positioning error of the focal spots could be mainly related to the thermos-lensing effect and nonlinear propagation of ultrasonic waves [23]. When accounting for the complexity of wave propagation, e.g., ribs, abdomen tissues, blood vessels and other celiac organs between the transducers and the targeted area, these could constitute sources of acoustic scattering, diffraction, attenuation and dispersion etc. [28,29], and negatively affect the precision of HIFU treatments. Therefore, it is very important to take into account the complexities of the acoustic paths by developing theoretical models for higher accuracy.

Using the state of the art methods, the spatial distribution of HIFU field can be simulated with the well-known Khokhlov–Zabolotskaya–Kuznetsov (KZK) equation [5,6], the Westervelt model [22] or the Spherical Beam Equation (SBE) [9]. Because of its reasonable parabolic approximation, the KZK equation has been widely used for describing the propagation of finite amplitude acoustic beams emitted by focused transducers, with the only restriction being that the half angle of divergence of the transducer does not exceed 16◦ [7]. Existing algorithms to solve the KZK equation are usually balanced between accuracy and simplicity, making it possible to calculate and adjust HIFU fields in real time during treatments. In contrast, the heavy calculation burden of the Westervelt equation and the SBE limits their application. Therefore, it is not unexpected that the KZK equation has been preferred against the other two in both clinical situations and industry. In the KZK equation, the effects of diffraction, nonlinearity and attenuation have been taken into consideration and are described by separate terms. Therefore, it has been recognized to be more effective than the Rayleigh integral in simulating HIFU propagation [8]. Meanwhile, the parabolic approximation is also accurate enough in practice, by which the near-axis HIFU field can be calculated relatively accurately [7,30,31]. However as Meaney et al. have pointed out, it is hard to explain the shifts of focal positions with the traditional KZK equation [24]. Therefore, several questions are still open to discussion. Firstly, the Kramers–Kronig dispersion relations indicate a power law between the attenuation and the working frequency, i.e., *<sup>α</sup>* <sup>∝</sup> *<sup>ω</sup>y*, *<sup>y</sup>* = <sup>1</sup> − 2. For most tissues, the attenuation factor *<sup>y</sup>* sits in the range from 1 to 1.7 [32,33]. However, in the KZK model, the propagation loss caused by viscosity and thermal conduction is considered to be proportional to the square of the working frequency (*y* = 2), which is actually only valid for fresh water [32,33]. Furthermore, *y* = 2 for fresh water causes a third derivative term to appear in the KZK equation. The third derivative as well as the derivative of acceleration have not yet been well clarified in the physical overview. In the context of Newton's law of motion, acceleration is directly affected by the force, and the derivative of acceleration has no obvious physical meaning. Secondly, sound velocity, which is treated as a constant in the KZK equation, usually changes with frequency in biological tissues. In HIFU fields, propagation nonlinearity could be non-negligible due to the high acoustic pressure, and the orders of harmonic waves could be rather high in many cases [5–7,9,29]. Therefore, the descriptions of both attenuation and sound velocity should be modified when biological tissues are present on the wave path.

To overcome these shortcomings in the existing models, efforts need to be made to account for the case of a non-integer power index in the Kramers–Kronig relationship. In the frequency domain, with the help of the classic Laplace Transform, the acoustic field could be easily accessible, but only valid for linear models. The Rayleigh model proposed by Wojcik et al. in 1995 was also not applicable in cases using wide-bandwidth acoustic pulses [34]. The fractional order derivative method described by Makris and Constantinou [35] was believed to be an effective approach to solve this problem. However, the complexity in mathematics has made numerical analysis hard to achieve, hence restrained its application in developing HIFU theories. After that, Szabo et al. utilized a convolution integral and further developed the theory of the fractional order derivative, in which Fourier Transform was

adopted to define a derivative of non-integer order, while the basic principles of the Kramers–Kronig relationship was still preserved [36–38]. Some other calculation algorithms have also been proposed. Treeby et al. developed a popular, open source nonlinear simulation tool named the k-wave toolbox to simulate nonlinear wave propagation [39]. Prieur et al. proposed time-fractional acoustic wave equations [40]. Inspired by these algorithms, simulation of the KZK equation has made definite progress in the past few years.

In this paper, a fractional order derivative was introduced to modify the KZK model to address the power-law relationship (non-integer power index), and a frequency-dependent sound velocity was employed to account for the dispersive behaviors of media. Numerical and experimental verifications were carried out to demonstrate the different behaviors between the modified model and the original KZK model. Results showed that the obtained modified KZK (mKZK) equation is in better agreement with the experimental results. The findings in this paper will promote not only the prediction and design of focused sound fields and acoustic transducers, but also the therapeutic applications of ultrasonic treatments.

#### **2. Theory and Experiments**

#### *2.1. The KZK Equation*

The KZK equation is an extended form of Burgers model, and an approximation of the Westervelt equation, written as [41]

$$\frac{\partial^2 p}{\partial z \partial \tau} = \frac{c\_0}{2} \Delta\_\perp p + \frac{\delta}{2c\_0^3} \frac{\partial^3 p}{\partial \tau^3} + \frac{\beta}{2\rho\_0 c\_0^3} \frac{\partial^2 p^2}{\partial \tau^2},\tag{1}$$

in which *p* is the acoustic pressure, *c*<sup>0</sup> is the sound velocity, *δ* is the sound diffusivity, *β* is the nonlinearity coefficient, *ρ*<sup>0</sup> is the ambient density of the medium, Δ<sup>⊥</sup> is the transverse Laplace operator (defined as <sup>Δ</sup><sup>⊥</sup> <sup>=</sup> <sup>1</sup> *r ∂ ∂r* , *r ∂ ∂r* - + 1 *r*2 *∂*2 *∂θ*<sup>2</sup> in cylindrical coordinates), and *<sup>τ</sup>* <sup>=</sup> *<sup>t</sup>* <sup>−</sup> *<sup>z</sup>*/*c*<sup>0</sup> is the time delay at the axial distance of *z* with *t* being the time. The three terms on the right side of Equation (1) represent the diffraction, attenuation and nonlinearity, respectively. The attenuation caused by viscosity and thermal conduction is hence proportional to the square of the angular frequency *ω*, i.e., *α*(*ω*) = *ω*2*δ*/(2*c*<sup>3</sup> 0). However, in biological tissues, the attenuation factor is usually a non-integer less than 2, the attenuation term is hence only valid for describing wave propagation in fresh water (attenuation factor *y* = 2). If the attenuation parameter is directly set for fresh water rather than considering the actual properties of the media, the accuracy of calculation would certainly be undermined. Furthermore, the sound velocity *c*<sup>0</sup> here is regarded as a constant for all harmonic components, which is in conflict with the inherent Kramers–Kronig dispersion relations where the phase velocity varies with increasing frequency.

#### *2.2. The Modified KZK Model*

By introducing Fourier Transform, the *n*th time derivative of time-dependent acoustic pressure *p*(*t*) could be written as

$$\frac{d^n p(t)}{dt^n} = F^- \left\{ (i\omega)^n F^+ \left[ p(t) \right] \right\}. \tag{2}$$

Here *F*<sup>+</sup> and *F*<sup>−</sup> represents operators of the Fourier transform and its inversion, respectively. As the order number *y* is a non-integer, the fractional order derivative is then defined as a convolution [38],

$$\frac{d^y p(t)}{dt^y} = \frac{1}{\Gamma(-y)} \int\_{-\infty}^{t} \frac{p(t')}{(t - t')^{(1+y)}} dt',\tag{3}$$

where <sup>Γ</sup>(·) is the Gamma Function, i.e., <sup>Γ</sup>(*x*) <sup>=</sup> ; <sup>∞</sup> <sup>0</sup> *<sup>ξ</sup>x*−1*e*−*ξdξ*. The definition in Equation (3) is hence the Riemann–Liouville fractional derivative [42]. Therefore, the fractional order derivative is not only determined by the pressure value at the time point *t*, but is also related to its history range from −∞ to *t*.

Since the attenuation factor *y* is a non-integer in lossy media, the wavenumber 6*k* could then be considered as a complex, while its squared form in the low-frequency approximation is

$$
\widehat{k}^2 \approx \frac{\omega^2}{c\_0^2} + 2i \frac{\omega}{c\_0} a\_0 \omega^y = \frac{\omega^2}{c\_0^2} - 2 \frac{a\_0}{c\_0} \frac{(-i\omega)^{y+1}}{(-i)^y}.\tag{4}
$$

in which *i* is the imaginary unit, and (−*i*) *<sup>y</sup>* can be expressed with trigonometric functions as (−*i*) *<sup>y</sup>*=cos(*yπ*/2) − *<sup>i</sup>*sin(*yπ*/2). With the wave number <sup>6</sup>*<sup>k</sup>* = *<sup>ω</sup>*/*c*<sup>0</sup> + *<sup>i</sup>α*0(−*iω*) *y* /[cos(*yπ*/2)], the phase velocity can be calculated as

$$\frac{1}{\mathcal{L}(\omega)} = \frac{\text{Re}\left(\stackrel{\sim}{k}\right)}{\omega} = \frac{1}{c\_0} + \alpha\_0 \tan(y\pi/2) |\omega|^{y-1},\tag{5}$$

where Re(·) represents the real part of the complex value.

Considering the non-integer attenuation factor *y* and the dispersion of phase velocity, the classical KZK equation could be modified as

$$\frac{\partial^2 p}{\partial z \partial \tau} = \frac{c(\omega)}{2} \nabla^2\_\perp p + \frac{\delta}{2c(\omega)^3} \frac{\partial^{y+1} p}{\partial \tau^{y+1}} + \frac{\beta}{2\rho\_0 c(\omega)^3} \frac{\partial^2 p^2}{\partial \tau^2},\tag{6}$$

where the attenuation term and the phase velocity *c*(*ω*) could be calculated according to Equations (4) and (5), respectively. Consistent conclusions can be found between Equation (6) and the work of Zhao et al. [43], which is an extension of an earlier model for ultrasound propagation in power-law media proposed by Kelly et al. [44]. In addition, Equation (6) will not hold if *y* = 1. However, in the framework of discussing HIFU propagation problems, the *y* = 1 case is of no importance [37,38] and is not of interest here.

#### *2.3. The Numerical Algorithm*

In the numerical analysis, both the KZK model and its modified form were solved with the finite difference time domain (FDTD) method for the acoustic field emitted from a single-element self-focusing transducer. For the traditional model, coordinate transformation *Z* = *z*/*F*, *R* = *r*/*a*, *T* = *ωt*, *P* = *p*/*P*<sup>0</sup> were introduced, with *F*, *r*, *a*, *P*<sup>0</sup> being the geometrical focal length, the radial coordinate, the aperture radius of the transducer, and the surficial acoustic pressure, respectively. The following assumptions were then made,

$$G = \frac{ka^2}{2F} \tag{7}$$

$$A = \frac{\omega^2 \delta}{2\rho\_0 c\_0^3} F = \mathfrak{a}' F\_\prime \tag{8}$$

$$N = \frac{F}{\rho\_0 c\_0^3 / (P\_0 \beta \omega)} = \frac{F}{l\_d} \,\tag{9}$$

in order that the KZK equation could be normalized to the following form,

$$\frac{\partial^2 P}{\partial T \partial Z} = \frac{1}{4G} \triangle\_\perp P + A \frac{\partial^3 P}{\partial T^3} + \frac{N}{2} \frac{\partial^2 P^2}{\partial T^2}. \tag{10}$$

The values of the physical constants used for acoustic modeling were *ρ<sup>0</sup>* = 1000 kg/m3, *c*<sup>0</sup> = 1486 m/s, *β* = 3.5, *α* = 0.025 Np/m at 1 MHz, and *μ* = 2 for fresh water [45]. The FDTD algorithm adopted here was generally the same as that used in a previous study [45,46], in which Equation (10) was decomposed into three independent equations, accounting for the diffraction, attenuation and

nonlinearity, respectively. Based on an orthogonal spatial grid, the discretized forms of these equations were expressed as,

$$\frac{P\_{i,j+1}^{n} - P\_{i,j}^{n}}{dZ} = \frac{1}{4G} \int\_{T\_{\text{min}}}^{T} \left( \frac{P\_{i+1,j+1}^{n} - 2P\_{i,j+1}^{n} + P\_{i-1,j+1}^{n}}{dR^{2}} + \frac{P\_{i+1,j+1}^{n} - P\_{i-1,j+1}^{n}}{2kdR^{2}} \right) dt,\tag{11}$$

$$\frac{P\_{i,j+1}^u - P\_{i,j}^u}{dZ} = A \frac{P\_{i,j+1}^{u+1} - 2P\_{i,j+1}^u + P\_{i,j+1}^{u-1}}{dT^2},\tag{12}$$

And

$$P\_{j+1}^{n} = \begin{cases} P\_j^n \left( 1 - N \frac{P\_j^{n+1} - P\_j^n}{dT} dZ \right)^{-1} & P\_j^n \ge 0 \\\ P\_j^n \left( 1 - N \frac{P\_j^n - P\_j^{n-1}}{dT} dZ \right)^{-1} & P\_j^n < 0 \end{cases} \tag{13}$$

Here *i* and *j* were the spatial coordinate indexes in the radial and axial directions, respectively, and *n* was the time step, i.e., *P<sup>n</sup> <sup>i</sup>*,*<sup>j</sup>* = *P*(*i* · *dr*, *j* · *dz*; *n* · *dt*).

In the simulations for the mKZK equation, the major difference was the differential form of the fractional derivative in Equation (6),

$$\frac{\partial^{y+1}P}{\partial \tau^{y+1}} = \frac{A}{\Delta \tau^y} \times \left[\Delta \tau^{-1} \sum\_{r=0}^{n} \omega\_r^{(y)} \left(P\_{i,j}^{n-r+1} + P\_{i,j}^{n-r}\right) + \omega\_{n+1}^{(y)} P'(0)\_{i,j}^k\right],\tag{14}$$

where

$$A = 2\Gamma(-y)\Gamma(y+1)\cos[(y+1)\pi/2]/\pi,\tag{15}$$

$$
\omega\_r^{(y)} = (-1)^r y (y - 1) \Gamma(y - r + 1) / r! \,. \tag{16}
$$

The FDTD simulations were then accomplished through a self-developed FORTRAN code package running on a ×64 PC platform. In the calculation, the boundary condition was symmetrical and had equal amplitude acoustic pressure driving conditions, which was also the boundary condition generally used to calculate HIFU [29,41].

#### *2.4. Experimental Methods*

#### 2.4.1. Phantom Preparation

A tissue phantom was prepared based on the recipe of polyacrylamide electrophoresis gel [47], in which micron-sized polystyrene microspheres were added to adjust its attenuation and phase velocity dispersion. The formula of phantom contained 100 mL degassed water, 10 g acrylamide (A9099, Sigma-Aldrich, St Louis, MO, USA), 0.05 g ammonium persulfate (A9164, Sigma-Aldrich), 0.3 g methylene double acrylamide (146072, Sigma-Aldrich), 0.2 mL TEMED (411019, Sigma-Aldrich), and 4 mL 10-micron microsphere solution (P107798, Aladdin, Shanghai, China, original concentration 5% *w*/*v*).

#### 2.4.2. Experimental Setup

The thickness of each phantom sample was *L*<sup>s</sup> = 42.3 mm, with a density of *ρ*<sup>1</sup> = 1000 kg/m3 so that it could stably suspend in the water. The distance between transducer and the phantom was 37.7 mm. Following the same protocol used in our previous work [45], the nonlinearity parameter was measured as *β* = 4.2. Since the only difference in gel recipe between the two works is the introduction of amino polystyrene microspheres in this paper, the same *β* value indicates that microspheres did not influence the nonlinear propagation of waves. To confirm this, we measured the ratio between the second harmonic and fundamental components, and found it was identical to that in [45] under the same sonication conditions. The attenuation coefficient and sound velocity of phantom, as functions of frequency, were measured through a broadband spectrum method [48,49]. In the measurement, two planar piston transducers (Immersion, Unfocused, Panametrics, Waltham, MA, USA) calibrated with a needle hydrophone (HNC-1000, ONDA Corp., Sunnyvale, CA, USA), were placed on the opposite sides of the phantom, with one of them driven by a broadband pulse generator (5900PR, Panametrics). The reflected and transmitted acoustic signals were then acquired by the transducers and digitalized with a digital oscilloscope (54830B, Agilent, Santa Clara, CA, USA). (Device connection was similar to that in [49]). In the measurements, 8 continuous pulse sequences were acquired and averaged to reduce the signal-to-noise ratio (SNR).

Prior to the HIFU experiments, a low-level driving voltage was used to drive a customized HIFU transducer (Chongqing Haifu Med. Tech. Co., Ltd., Chongqing, China) to emit a linear sound field. Simultaneously, effective parameters of the transducer, such as effective radius, radiation profile, and angle of divergence were obtained by adjusting the transducer parameters in the KZK calculations, such that the linear field predicted via KZK was consistent with that measured [45]. The effective parameters were then used in the simulations of both the KZK and mKZK models. As a result, the HIFU transducer (working frequency 1.12 MHz) had an effective aperture radius of 48.6 mm and a geometrical focal length of 101.5 mm. As illustrated in Figure 1, the transducer was immersed in water and driven with signals from a signal generator (33250A, Agilent, Santa Clara, CA, USA) amplified by a broadband power amplifier (2200L, E&I, Rochester, NY, USA). The input voltage was set as 465 mV (20 cycles; burst period, 10 ms; duty cycle, 0.18%). Another needle hydrophone (HNA-0400, ONDA Corp., Sunnyvale, CA, USA) was mounted on a customized three-dimensional (3D) scanning system (Controller Model: XPS-C8, Newport, CA, USA) to scan the HIFU field. To suppress possible acoustic cavitation in surrounding liquid, the water was processed with a self-developed water degassing and deionizing system. The temporal and spatial scanning procedure was controlled via the GPIB interface (National Instruments, Austin, TX, USA).

**Figure 1.** The experimental setup.

#### **3. Results and Discussions**

#### *3.1. Non-Dispersive Water*

In order to verify the validity of the modified model, the mKZK equation was used to predict the sound field distributions generated from the transducer. In this case, the media in the direction of propagation was degassed water, and the surface pressure of the transducer was set to be 0.4 MPa. Then, the results were compared with those obtained from experimental measurements as well as from the original KZK model. The results are presented in Figure 2 for comparison. Due to the axial-symmetry of the sound field, only the sound field distribution in the axial direction was studied. Also, since the major concern of this paper is to investigate how the focal-shift could be accurately predicted, the pressure profiles are all presented in a normalized way, so that the focal-shift effect is more intuitive and easy to observe. For the total pressure distributions presented in Figure 2a, the axial distribution of acoustic pressure seems identical for the modified and original KZK models, which provides reliable proof that the current theoretical modification did not compromise the accuracy of the sound field prediction in non-dispersive media. With the help of fast Fourier transform (FFT) algorithm, further analysis was then carried out by decomposing the total sound pressure into the superposition of linear and nonlinear components. In Figure 2b–d, all the components exhibited good agreement between the results calculated from the modified and the original KZK models, showing that both models are applicable for sound field prediction under the experimental conditions mentioned above. It should also be noted that, in Figure 2 the locations of the pressure peaks from the two models were exactly the same for all components, although the measured axial beam-widths seem a bit narrower than both theoretical predictions, especially for harmonic components. Since the linear fields were found to be almost identical for the measured and predicted results, we speculate that some far-field attenuation factors such as dissolved oxygen in water, or bubbles might exist. However, this does not affect the conclusion on the location of the maximum pressure.

**Figure 2.** The normalized acoustic pressure distributions along the axis of the HIFU transducer without the phantom: (**a**) the overall pressure; (**b**) the fundamental component; (**c**) the second harmonic; (**d**) the third harmonic.

#### *3.2. Dispersive Phantom*

To incorporate the phantom model into the study, the acoustic parameters of the phantom was firstly characterized according to procedures described in earlier studies [48,49]. In brief, the acoustic phase velocity *c*(*f*) inside the phantom material was determined by [48]

$$\mathcal{L}(f) = c\_{\rm w} \left[ 1 + 2 \frac{\theta\_{\rm w}(f) - \theta\_{\rm s}(f)}{\theta\_2(f) - \theta\_1(f)} \right],\tag{17}$$

where the sound velocity in water *c*<sup>w</sup> was considered as 1500 m/s. When the sound velocity inside the phantom was measured, the acoustic signals *p*<sup>s</sup> and *p*<sup>f</sup> were the acoustic pressure acquired by a transducer before and after the phantom was inserted into the acoustic path. While *p*<sup>1</sup> and *p*<sup>2</sup> were the pressure of the reflected signals from the first and second water/phantom interfaces, respectively, *θ*w(*f*), *θ*s(*f*), *θ*1(*f*) and *θ*2(*f*) were the corresponding phase spectra. The frequency dependence of the attenuation coefficient *α*(*f*) was calculated according to [48,49]

$$\ln(f) = \frac{1}{L\_s} \left[ \ln \left( \frac{A\_1}{A\_2} \right) - \ln \left( \frac{A\_w}{A\_s} \right) \right],\tag{18}$$

where *A*w, *A*s, *A*<sup>1</sup> and *A*<sup>2</sup> were the amplitude spectra corresponding to the above-mentioned phase responses. Figure 3 plots the measured attenuation coefficient and acoustic velocity as a function of frequency. In the frequency range 1.5–3.1 MHz, the sound velocity increased by about 26 m/s. The acoustic attenuation coefficient increased from 0.59 dB/cm at 1.5 MHz to 1.79 dB/cm at 3.1 MHz, indicating the attenuation factor being *y* = 1.83 for the phantom. The attenuation factor was obtained from the exponential fitting using a curve fitting toolbox in Matlab. It should be noted that sound velocity could fluctuate due to the thermal effect of focused ultrasound. However, in this work the duty cycle was as low as 0.18%, and a medium-level surface pressure of up to 0.4 MPa was chosen for the transducer. Thus, no significant temperature elevation was observed during the experiments, and the thermal-induced change in sound velocity could be neglected [50]. It should also be mentioned, that to account for the thermal-effect, an appropriate bio-heat transferring equation should be incorporated with the current model. The main concern then lies in the dispersion due to microsphere scattering. In clinical applications, the major dispersion originates from tissues like fat, whose attenuation coefficient is generally larger than that of body tissue. Therefore, in comparison with other research in which the parameters of the phantom were nearly the same as body tissue (e.g., liver and spleen), the choice of phantoms with larger attenuation coefficients in the present experiments might provide results in better agreement with the actual situation. Meanwhile, since the measured physical parameters are close to those of dense fat, this setup could be regarded as a simplified mimic of the abdomen in HIFU therapies.

**Figure 3.** Measured frequency-dependent acoustic velocity and attenuation coefficient of the phantom sample.

As is illustrated in Figure 4, the acoustic field distribution along the axial direction is examined by sitting the phantom on the acoustic path of the transducer. The results show the comparison between measured data and simulated results obtained from both the modified and traditional KZK models. Figure 4a describes the total pressure and Figure 4b–d shows the fundamental, second harmonic and third harmonic components, respectively. For the total pressure distribution, it is clearly observed that the peak-pressure location predicted through the mKZK model agrees well with that acquired in experiments, while the data calculated from the traditional KZK model show a deviation from the previous two groups. Note that the axial peak-pressure deviation is quite small for the fundamental components, but it gradually becomes significant when more harmonic components appear. Thus, the overall peak-pressure deviation observed in Figure 4a is mainly caused by higher-order harmonic components, indicating that the significance of the theoretical modification relies on high nonlinearity, such as in HIFU. This phenomenon could be addressed with the results shown in Figure 3, where higher-order harmonic waves that occupy higher frequency bands exhibit more notable sound velocity/attenuation dispersions, thus play a more dominant role in the modification of the KZK model. Meanwhile, it can be seen from Figure 4 that the modified equation has larger variation than the KZK equation. The effect of the mKZK equation on the simulation results can be thus clarified as providing more precise prediction for experimental results. The results in Figure 4 give persuasive proof that theoretical modifications made previously are necessary and valid. The beam narrowing effect caused by data normalization still exists. However, in actual treatment more attention is paid to the location of the focus point, because the location determines the heat distribution area, which significantly alters the biological properties of the treatment area.

**Figure 4.** The normalized acoustic pressure distributions along the axis of the HIFU transducer, with the phantom placed at 37.7 mm away from the transducer. (**a**) the overall pressure; (**b**) the fundamental component; (**c**) the second harmonic; (**d**) the third harmonic.

Among existing studies, most researchers have focused on the thermos-lensing effect [23] and bubble formation induced focal region distortion [25,26], while some also mentioned the acoustic radiation force induced tissue displacement [27]. With the results presenting the dispersion of sound velocity and attenuation in tissues, the deviation of focal spots can be well explained in combination with the strong nonlinearity nature of HIFU.

#### *3.3. Dispersion-Induced Focus Shift*

It is of great importance to evaluate how the mKZK model demonstrates its significance in HIFU applications. In Figure 5, comparisons are carried out to display how the wave distribution behaves differently before and after inserting the phantom sample into the wave propagation path. It can be observed that, for either data from experiments or from the mKZK model, although only a slight difference is seen in the axial wave profile when examining the fundamental components, an axial focus shift is more evident for the overall acoustic pressure since it is highly affected by the harmonic components. It is explained here that, due to the dispersive nature of the phantom, higher-order harmonic components require larger values of both sound velocity and attenuation coefficient in the KZK model, urging the overall wave profile to move forward to the transducer.

**Figure 5.** Comparison of the axial distributed acoustic pressure between the cases of with and without the phantom sample: (**a**) the overall pressure; (**b**) the fundamental component; (**c**) the second harmonic; (**d**) the third harmonic.

In the present study, the focus shift distance was quantified for the overall axial acoustic pressure distribution and its decomposed components, and this is listed in Table 1.

**Table 1.** Focus shift induced by acoustic-dispersive phantom sample (in mm).


The focus shift, which is found to be of millimeter magnitude for the studied condition, cannot be ignored, especially for clinical HIFU studies. On the one hand, in typical HIFU applications, the much higher surface pressure would induce even stronger acoustic nonlinearity in the focus area, and up to tens of orders of harmonic components might contribute to the overall acoustic responses. In that case, dispersion in acoustic velocity as well as the attenuation coefficient could cause larger focus shifts. On the other hand, even for the millimeter-level focus shift exhibited with the current setup, an impressive amount of acoustic energy would be deposited outside the designated focus area. As demonstrated in Figure 6, about 25% of the −3 dB focal region would fall outside of the one predicted by the traditional KZK equation, when the phantom sample is introduced into the transducer axis. In addition, the absence of shockwaves in the experiment should be noted, and this indicates that the influence of shockwaves could be ignored in the theoretical model under such experimental parameters.

**Figure 6.** Illustration of the focus shift induced by the presence of phantom sample: (**a**) the theoretical −3 dB area; (**b**) the measured −3 dB area.

#### *3.4. Discussion*

In previous studies, shifts between the predicted and observed focal regions have been frequently reported in HIFU-related studies. Usually, the shift was measured to be around 1–3 mm for 1-MHz HIFU transducers [24], 4–5 mm for ~2.2-MHz excitation [51], and an empirical formula was also used to calculate the focus shift [22]. Although different possible mechanisms have been proposed and a series of correlation studies have been carried out [21–27], detailed theoretical proof that could quantitatively explain the inherent mechanisms of the observed focal shifts is still lacking. The proposed model modification here clarifies how the acoustic dispersion played a role in the complicated physics of this problem.

The experiments carried out here have demonstrated that the observed focal shift should come from the dispersive behavior of sound velocity and attenuation in media. For different harmonic components in the HIFU beam, their sound velocities could be different at the gel/water interface. Speculating from Snell's law, these different components could actually propagate along slightly different paths in the phantom, causing the focal point to shift its location. That is also why the main difference between the two models is observed for the harmonics in Figures 4 and 5—because the sound velocity of the harmonics deviated further away as their frequencies were higher. During this process, the attenuation should also have contributed in a dispersive way. However, the influence of dispersive attenuation could not be separated from that of sound velocity.

Despite the general recognition that the Westervelt equation and SBE provide higher precision in non-dispersive sound field prediction than KZK, difficulties in predicting accurate sound field distribution in strong dispersion media still remains a great challenge [22]. Moreover, due to the existing complexity of the Westervelt equation and SBE, adding more modifications to these two equations could be over-whelming and/or time-consuming for clinical applications that require real-time monitoring. Considering that KZK shows a balance between accuracy and computation burden, this paper chose to further modify the KZK equation as a straightforward way to achieve improved simulation of HIFU propagation, especially to work out how the focal-shift could be predicted accurately. In this work, for consistency of the experiment, the mKZK and KZK in ordinary media was first confirmed. Then, when measuring the parameters in the strong dispersion medium, mKZK gave a more accurate prediction of the focal shift. Thus, by using the mKZK equation we can quickly and accurately predict the HIFU field in complex media.

However, further theoretical studies are still needed since the current model is unable to eradicate all the possible unfavorable factors in predicting the characteristics of HIFU. For example, defocusing effects or shifts in focal position might be caused by the layered tissue effect, where sound speed/attenuation may vary in different tissues layers. The thermal lesion effect, in which the tissue properties change due to the heating of HIFU, could be more difficult to include in the modeling. A full solution could be even more challenging if other possible mechanisms, including acoustic radiation force, acoustic cavitation and inconsistent thermal deposition [27] are also considered.

The focal shifts observed in phased-array-based HIFU devices are also notable. For instance, a focal shift of about 2-mm was observed along the transducer axis in multiple-layered soft tissues sonicated with a 65-element phased array transducer [22]. To overcome this problem in phased-array HIFU, possible solutions could be obtained by drawing lessons from the underlying mechanisms discussed above.

#### **4. Conclusions**

Although HIFU technology based on the KZK equation calculation has been widely accepted and used in the clinical setting and transducer designs, the absence of an accurate theory to predict the sound field inevitably limits the application of the ultrasound focusing. In this work, a mKZK equation is proposed to predict the HIFU field established with a spherical focusing transducer. Meanwhile, the accuracy of the methods is verified through experiment. This method could improve the computational accuracy of the KZK equation in dispersive media, which is similar to human tissue. Simulation and experimental results show that the focus area will shift towards the transducer and the offset increases as the nonlinearity becomes higher. Therefore, in the process of HIFU transducer design, the impact of dispersion on the results need to be taken into account, in order that accurate sonication can be achieved. By modifying the KZK equation, the findings will help with transducer design and the application of the HIFU. This will also help to ensure the stability and safety of HIFU and further accelerate its clinical applications.

**Acknowledgments:** This work was partially supported by the National Natural Science Foundation of China (Grant Nos. 81627802, 11774166, 11474161, 11474001 and 11674173), QingLan Project and Nanjing University Innovation and Creative Program for Ph.D. candidate, No. CXCY17-13.

**Author Contributions:** Xiasheng Guo and Dong Zhang conceived and designed the experiments. Shilei Liu and Yanye Yang performed the experiments. Shilei Liu, Xiasheng Guo and Juan Tu analyzed the data. Shilei Liu and Chenghai Li performed the derivation theory. Shilei Liu and Xiasheng Guo contributed reagents/materials/analysis tools. Shilei Liu, Xiasheng Guo and Dong Zhang wrote the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Lynn, J.G.; Zwemer, R.L.; Chick, A.J.; Miller, A.E. A new method for the generation and use of focused ultrasound in experimental biology. *J. Gen. Physiol.* **1942**, *26*, 179. [CrossRef] [PubMed]


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Application of Elastic Wave Velocity for Estimation of Soil Depth**

#### **Hyunwook Choo 1, Hwandon Jun <sup>2</sup> and Hyung-Koo Yoon 3,\***


Received: 21 March 2018; Accepted: 5 April 2018; Published: 11 April 2018

**Abstract:** Because soil depth is a crucial factor for predicting the stability at landslide and debris flow sites, various techniques have been developed to determine soil depth. The objective of this study is to suggest the graphical bilinear method to estimate soil depth through seismic wave velocity. Seismic wave velocity rapidly changes at the interface of two different layers due to the change in material type, packing type, and contact force of particles and thus, it is possible to pick the soil depth based on seismic wave velocity. An area, which is susceptible to debris flow, was selected, and an aerial survey was performed to obtain a topographic map and digital elevation model. In addition, a seismic survey and a dynamic cone penetration test were performed in this study. The comparison between the soil depth based on dynamic cone tests and the graphical bilinear method shows good agreement, indicating that the newly suggested soil depth estimating method may be usefully applied to predict soil depth.

**Keywords:** graphical bilinear method; seismic survey; dynamic cone penetration test; soil depth; time-distance curve

#### **1. Introduction**

In assessing the stability of landslide or debris flow areas, both hydrological and geotechnical properties are the key parameters [1–6] among the various geotechnical properties such as soil strength, hydraulic conductivity, and friction angle, the soil depth is the most important parameter because the capacity for inflow and outflow of water is related to the soil thickness. Note that soil depth can be defined as the thickness of the soil from ground surface to consolidated medium [7,8].

The most reliable method to estimate soil depth at a given location is the test pit method, which involves direct excavation of a testing site in a square shape. However, excavation of multiple pits for the estimation of soil depth is very expensive and time-consuming. To overcome these limitations, Ref. [9] used a cone-tipped metal probe with a diameter of 18 mm for estimating local soil depth based on probe penetration resistance. Note that the methods mentioned above can only be performed at the selected local points; therefore, the soil depth of unmeasured areas is assumed to be the same as the measured depth of nearby points. However, soil depth shows substantial spatial variation; thus, the reliability of the above methods may be low, leading to the development of several models that consider the spatial distribution of soil depth. Ref. [10] developed the regolith-mantle slope method to predict the thickness of original parent material on a slope. Ref. [11] estimated soil depth based on the mass balance between soil production from underlying bedrock and divergence of diffusive soil transport. Even though these models are advantageous for obtaining soil depth across the whole area, the use of these models may not be easy because they require detailed information

regarding hydrological, geotechnical, and geochemical parameters. Therefore, the geophysical method has been widely applied to estimate soil depth because it can quickly and cost-effectively provide various soil characteristics across an entire area. Electromagnetic waves have been used to measure electrical resistivity [12] and electrical conductivity [2], and the measured values can be converted into soil depths. Refs. [13,14] used seismic survey to predict soil depth and however, the studies merely showed the distributions of results without methodological content for picking soil depth. In addition, the applied previous geophysical methods require special reference values to determine soil depth and additional invasive experiments to enhance the reliability of the estimated soil depth.

This study proposes a graphical bilinear method to estimate soil depth based only on seismic survey. A seismic survey was performed, and the dynamic cone penetration test (DCPT) was also used to verify the soil depths deduced by the suggested technique. The measured information, including geological map, location, particle distribution, and distribution of seismic wave velocity were introduced first. Then the graphical bilinear method for determining soil depth, is demonstrated. The estimated soil depths using the suggested method are compared with values deduced by the dynamic cone penetration tests and their reliability is assessed.

#### **2. Testing Site**

#### *2.1. Site Description*

The selected testing site experienced debris flow a few years ago and is still susceptible to additional debris flow because of several geological characteristics of the site, including steep slope angle (over 32◦) and saturated soil condition on slope. According to the geological map issued by the Korea Institute of Geology, Analysis, Mining (KIGAM), the testing area mainly consists of gneissose granite terrane. The testing site belongs to Mt. Geohwa, South Korea, where the altitude and area are approximately 200 m and 1.2 km2, respectively, and the latitude and longitude of the top of the mountain are N 36◦29 15.7313 and E 127◦18 40.5986. A drone aerial survey was performed to obtain a topographic map and digital elevation model. Figure 1 shows the topographic map: the area consists of one main stream and several branch streams. The vertical length and area of the main stream are approximately 206 m and 14,675 m2, respectively. Figure 1 also indicates the presence of the debris barrier and check dam at the site to prevent additional debris flow (N 36◦29 04.3147, E 127◦18 44.7017). The main stream was divided at four points, selected by considering slope characteristics: point A (N 36◦29 10.8530, E 127◦18 44.1891), point B (N 36◦29 12.045, E 127◦18 39.5912), point C (N 36◦29 14.8329, E 127◦18 30.7079), and point D (N 36◦29 18.2046, E 127◦18 23.6158). Point A, which is located at the bottom of the main stream, shows various fallen trees and weeds undergoing decay and rot, suggesting that considerable time has passed since the occurrence of the last debris flow. The slope of point A was measured to be approximately 10◦, which is low compared to the slopes at points B and C. Thus, point A is covered with various flowed materials from the streams above due to debris flow. Point B shows a dramatically steeper slope of approximately 27◦ and the area mainly consists of sedimentary basin. Point C has the largest rapid slope (~32◦) and width (maximum 105 m) among the selected points. Note that point C is located at a large catchment basin. Finally, point D indicates the initial zone where the debris flow occurred as ascertained by the presence of upright trees and a stable subsurface without collapsed conditions.

**Figure 1.** Aerial photography of the testing area. Points A and D denote the bottom and top streams, respectively. Points B and C indicate the middle stream. The field tests were performed from points A to D.

#### *2.2. Soil Classification*

A hand auger with an outer diameter of 20 cm was used to excavate the subsurface and a soil sampler was used to gather the specimen at the wall of a borehole. The length of the hand auger was limited to 2 m due to a shortage of engine output, and thus the maximum depth of extracted soil in this study was 2 m. Even though the hand auger can potentially excavate soil down to a maximum of 2 m, the actual excavated depth was reduced due to the presence of gravel and weathered rock. Therefore, the actual maximum extracted depths were 1, 1, and 2 m at locations near points A, B, and C, respectively. The results from the excavations at points A and B show various deposited materials, including gravel and boulder (Figure 2). At point C near the initial zone, a larger extraction depth (~2 m) could be achieved due to the presence of a relatively deep soil layer (Figure 2). Soil samples were collected at 4–5 different depths: 10, 40, 60, and 100 cm for points A and B; and 10, 40, 90, 140, and 200 cm for point C. Sieve tests were performed using the extracted specimens according to [15].

The results of sieve tests are shown in Figure 2. The diameters at passing percentages of 10, 30 and 60% are calculated for every specimen, to classify the specimen based on the unified soil classification system. The calculated coefficients of uniformity are 10, 5.92, and 9.33 for extracted soils at points A, B, and C, respectively. The coefficients of curvature at points A, B, and C are determined to be 1.34, 1.01, and 1.51, respectively. The measured coefficients of uniformity show values almost greater than 6 and the coefficients of curvature are in the range of 1~3; hence, the specimens at testing sites are classified as SW (well-graded sand). Additionally, Figure 2 indicates that the grain size distributions of extracted samples at different depths are very similar.

**Figure 2.** Sieve test results: (**a**) bottom stream (point A in Figure 1); (**b**) middle stream (point B); (**c**) top stream (point C). D10, D30, and D60 denote the diameters at passing percentages of 10, 30, and 60%, respectively. Cu and Cg are coefficients of uniformity and curvature, respectively.

#### **3. Methodology**

A seismic survey and dynamic cone penetrometer test were performed to obtain primary wave velocity and dynamic cone penetration index (DPI), respectively, and the detailed descriptions are as follows.

#### *3.1. Seismic Survey*

A seismic survey has the advantage of assessing a whole area without altering the fabric state; thus, the measured values can greatly reflect the state and behavior of soils. A seismic wave propagates in solid media by travelling through particle connections, and thus the wave velocity increases with an increase in the particle contact area (or applied confining stress or soil depth). The seismic refraction method was applied in this study to obtain a profile of compressional (or primary) wave velocity. Four transects were determined with consideration for spatial variability in the geological characterization of each stream. Figure 3 shows schematic drawings of the seismic transections. The bottom-top stream (BTS) line reflects the main stream from south to north (from points A to D in Figure 1), whereas the bottom stream (BS), middle stream (MS), and top stream (TS) lines are horizontally set to reflect bottom, middle, and top valleys, respectively. The center of the MS line is located between points B and C in Figure 1. The lengths of the transection lines are 90, 20, 20, and 20 m for BTS, BS, MS, and TS, respectively. A geophone was installed as a sensor for gathering seismic waves every 2 m; therefore, the installed geophone numbers in the profiles of BTS, BS, MS and TS are 45, 10, 10, and 10, respectively. A sledgehammer was used to generate vibrations, and the impactions were performed at start, middle, and end positions to enhance signal-to-noise ratios. Impactions were performed five times in the same position for each test to minimize random noise.

**Figure 3.** Locations of the field investigations, including seismic survey and dynamic cone penetrometer (DCP). Note: the numbers with red circles denote the DCP testing sites (BS: 3 holes, MS: 3 holes, TS: 3 holes, and BTS: 10 holes); and the distance between each number in the figure is 10 m.

#### *3.2. Dynamic Cone Penetrometer Test*

A dynamic cone penetrometer (DCP) test was selected as an invasive method to verify the soil depths estimated by seismic waves. Note that the DCP method can readily detect the thickness of soil layer because the continuous penetration of DCP enables the easy recognition of the presence of different layers [16]: the penetration depth through stiff material is relatively small. The DCP technique has been widely applied to investigate soil properties in geotechnical engineering, especially for railways and roads [16–18]. The test records the penetration depth of DCP when a hammer drops from a fixed height according to [19]. The hammer weight and drop height were fixed to 78.8 N and 575 mm, respectively. The tip diameter of the DCP was 20 mm with an apex angle of 60◦, and driving energy was 45 J. The penetration depth was converted into dynamic penetration index (DPI), which represents (mm/blow) as follows.

$$\text{DPI} = \frac{\mathbf{P}\_{i+1} - \mathbf{P}\_i}{\mathbf{B}\_{i+1} - \mathbf{B}\_i} \tag{1}$$

where P = penetration depth (mm); B = blow count; and *i* and *i* + 1 = experimental number.

#### **4. Results and Analysis**

The first arrival time of the measured seismic wave signals was determined through PickWin software (OYO corporation, Tokyo, Japan), and the time-distance curve was obtained by using the SeisImager-Poltrefa program (OYO corporation, Tokyo, Japan). An interactive process was performed to increase resolution when carrying out the inversion process. Two groups, which are based on the picked first-arrival times and the results analyzed by the simultaneous iterative method, are plotted on a travel-time curve in Figure 4. Figure 4 shows that the two groups are nearly identical, reflecting that the quality of the measured data is high. The distribution of seismic wave velocity is plotted in Figure 5, and the soil layers are divided into various different layers according to the reference P-wave velocities suggested by [20]: 0–0.7 km/s (landfill and alluvial soil), 0.7–1.2 km/s(weathered soil), 1.2–1.9 km/s (weathered rock), and over 1.9 km/s (soft rock). The BTS, BS, MS, and TS lines were found to comprise four, four, three, and three layers, respectively based on Reynolds (2003). Note that 0.7 km/s of primary wave velocity is an accepted criterion in South Korea for distinguishing soil layers from other geomaterials [21,22].

**Figure 4.** Travel time-distance curves through seismic survey: (**a**) BTS; (**b**) BS; (**c**) MS; and (**d**) TS. Note 10, 20, 30, and 100 m distances indicate points 1, 2, 3, and 10, respectively, in Figure 3.

**Figure 5.** Converted primary wave velocity profiles: (**a**) BTS; (**b**) BS; (**c**) MS; and (**d**) TS.

The calculated DPI values from dynamic cone penetration tests are shown in Figure 6 for each stream. The BTS line shows high variation of DPI values according to the testing locations, which are described in Figure 3, because the line covers the bottom, middle, and top streams. Note that the high variation of DPI values means that the soil depth changes along the stream in the BTS line. Number 4, 8 and 9 holes in Figure 6a show high DPI values at shallow depth, reflecting the presence of weak soil layers. Being different from other testing holes, the DPI was measured at deeper depths at No. 7 and 9 holes, demonstrating that the soil depth near the top stream is deeper than that of other locations. Note that Figure 5a and the maximum extracted soil depth using a hand auger also support the presence of deep soil layer near the top stream in the BTS line. In the case of the BS line, high DPI values at the initial stage and high penetrated depth were observed at the right side of the stream. A similar DPI trend was observed in the MS line, and the distributions of soil depth in the No. 1, 2, and 3 holes look alike, reflecting similar soil depths along the MS line. Even though the DPI of the TS line shows a similar penetrated depth near the surface, the DPI was continuously measured at the only No. 2 hole. Therefore, it is predicted that the soil depth at the No. 2 hole, which is the center of the TS line, is high. The depth where the DPI value is close to zero is referred to as the soil depth of the testing site in this study, and the average soil depths based on DPI are calculated to be approximately 0.92, 0.75, 0.86 and 0.85 m for the BTS, BS, MS, and TS lines, respectively.

Comparison between Figures 5 and 6 reveals that the soil depth of the testing site based on dynamic cone penetration test is generally shallower than 1 m; while, the depth based on seismic wave velocity (0.7 km/s) is around 2 m or more (Table 1), reflecting that the subsurface classification based on previous reference values of seismic wave velocity cannot precisely estimate the soil depth. Therefore, in this study, the graphical bilinear method is newly suggested for determining soil depth based on seismic wave velocity.

**Figure 6.** Measured DPI values of each stream: (**a**) BTS; (**b**) BS; (**c**) MS; and (**d**) TS. Note the numbers in each figure indicate the location of dynamic cone penetration test described in Figure 3.


**Table 1.** Comparison between Measured and Estimated Soil Depths.


**Table 1.** *Cont.*

Note: error ratio = (estimated value − measured value)/measured value; the numbers in Position column indicate the testing locations in Figure 3.

#### **5. Discussion**

The seismic refraction method measures the travel time of the waves refracted at the interface between different sublayers with different wave velocities or impedances. Thus, the thickness of each layer can be calculated using the wave velocities of two consecutive layers and the time intercept in the travel time-distance curve as:

$$Z\_i = \frac{t\_i \cdot V\_i \cdot V\_{i+1}}{2\sqrt{(V\_{i+1}^2 - V\_i^2)}}\tag{2}$$

where *Z* = thickness of layer; *t* = time intercept; *V* = wave velocity; and *i* and *i* + 1 = different layers. The thickness of the soil layer (or the first layer) can be determined by the velocities of the first layer and the second layer (*V*<sup>1</sup> and *V*2), and the time intercept (*t*1). However, the clear determination of *V*1, *V*2, and t1 is not easy because wave velocity in a given layer is not constant, resulting from the natural soil deposits rarely being homogeneous and wave velocity slightly increasing with an increase in depth at a given layer. Note that seismic waves propagate in a medium through connected particles, and seismic wave velocity depends on soil structure and stress condition. Seismic wave propagation is primarily affected by the stiffness of the fabric in a particular material, and thus seismic wave velocity and effective stress (*σ* ) have a certain relationship with experimentally determined coefficients (*α* and *β*).

$$\text{Seismic wave velocity} = a(\sigma')^{\beta} \tag{3}$$

The *α* and *β* coefficients are dependent on the packing type (porosity and coordination number) and the contact force (Hertzian contact and Coulombic force).

The seismic wave velocity slightly increases with depth even in one layer because of the increase in effective stress (Equation (3)), which is the function of mass density and depth. However, the wave velocity changes rapidly at layer boundaries due to the change in material type, mass density, particle size, and others. Hence, the relation between seismic wave velocity and depth can produce different slopes along the depth. Figure 7 shows the seismic wave velocity with the depth at which the geophones are installed. It can be observed in Figure 7 that there are four different lines "a", "b", "c", and "d": Line "a" indicates the initial slope between wave velocity and depth; thus, it is related to the first soil layer; Line "b" indicates the second slope between wave velocity and depth; thus, it is related to the second layer. Therefore, the intersection point between lines "a" and "b" may correspond to the thickness of the first soil layer, which is point "e". Figure 8 shows the calculated soil depths based on the suggested graphical bilinear method at the selected testing sites. It is shown in Figure 8 that

various soil depths can be estimated by the proposed method. The detailed soil depths based on the suggested graphical bilinear method are summarized in Table 1.

**Figure 7.** Method for selecting soil depth through measured elastic wave velocity. Line "a" denotes the first slope between wave velocity and depth. Lines "b", "c", and "d" represent the second, third and fourth slopes, respectively. Point "e" is the intersection point between lines "a" and "b". Note the figure is the relation between wave velocity and depth of point 1 (or 10 m distance) in the BTS line (Figure 3).

**Figure 8.** Estimated soil depth profile based on graphical bilinear method: (**a**) 30 m of BTS; (**b**) 80 m of BTS; (**c**) 100 m of BTS; (**d**) 20 m of BS; (**e**) 30 m of MS; and (**f**) 20 m of TS.

The measured DPI values at the selected testing sites are plotted in Figure 9 to determine soil depth based on the impaction study. Soil depth is estimated to the maximum penetrated depth at which the DPI value is nearly zero. Even though the DPI shows variation with depth within a borehole, the final penetrated depth can be easily calculated based on the depth where the value of DPI is close to zero. Figure 9c shows that DPI does not converge close to zero until a depth of 2 m. The maximum possible experiment depth is fixed to 2 m due to the rod length and energy transfer during impaction, thus the soil depth is expected to be greater than 2 m. This deeper expected soil depth can be anticipated because the soil thickness was estimated to be approximately 3.0 m based on the graphical bilinear method shown in Figure 8c for the same position as in Figure 9c. Table 1 summarizes the soil depth calculated by dynamic cone penetration tests.

**Figure 9.** Soil depth based on dynamic cone penetration test: (**a**) 30 m of BTS; (**b**) 80 m of BTS; (**c**) 100 m of BTS; (**d**) 20 m of BS; (**e**) 30 m of MS; and (**f**) 20 m of TS.

Figure 10 shows the comparison of soil depths predicted by the reference P-wave velocity (0.7 km/s) and the graphical bilinear method with the measured depths using dynamic cone penetration test. Note that the depth based on the results of the dynamic cone penetration test (DPI values) has relatively high reliability because the data is gathered through direct penetration of the probe into soil; therefore, the soil depth measured by DPI can be regarded as the real soil depth. The soil depth predicted based on the reference P-wave velocity shows high variation and high soil depth compared with those based on DPI values. In contrast, the soil depth based on the graphical bilinear method shows small variation and the estimated values are comparable with the measured values by dynamic cone test, reflecting the enhanced reliability of the estimated soil depth by using the suggested method. Furthermore, the deduced soil depths based on DPI, P-wave velocity and graphical bilinear method are compared through box and whisker plot as shown in Figure 11. The median values of DPI and graphical bilinear method show 0.35 m and 0.3 m, and however the median soil depth based on P-wave velocity is highly deduced to 0.5 m. The first and third quartiles of DPI and

graphical bilinear method show also similar ranges of 0.15–0.55 m and 0.2–0.6 m, respectively and it shows the difference is just 0.5 m. However, the soil depth estimated by existing method shows unreasonable ranges of 1.05–1.8 m. The ranges of minimum and maximum values based on DPI (0.2–2 m) and suggested method (0.2–3 m) are also demonstrated to similar trends while those deduced by P-wave velocity exhibited huge variations (1.4–5.6 m). And thus, the Figure 11 shows that the suggested method can provide the reliable soil depth under ≈3 m with consideration of 2 m intervals of geophone.

Table 1 shows the calculated error ratios (error ratio = (estimated value − measured value)/ measured value). The estimated soil depths based on reference P-wave velocity show high error ratios, ranging from 55% to 1000%, and the average value shows approximately over 100%, reflecting the estimation of soil depth using the classification of soil layers based on reference wave velocity is not reliable. In contrast, the graphical bilinear method shows relatively small average error ratios of 28.2, 6.6, 10.5, and 21.9% for BTS, BS, MS, and TS, respectively, demonstrating that the graphical bilinear method suggested in this study provides reliable estimation of soil depth.

**Figure 10.** Comparison of soil depths estimated by seismic wave velocity and measured by dynamic cone penetration test.

**Figure 11.** Comparison of soil depths through box plot.

#### **6. Conclusions**

This paper suggests methods for the determination of soil depth through seismic wave velocity measurements. The graphical bilinear method is newly introduced, and the soil depth estimated with the suggested method is compared with the result of dynamic cone penetration tests. The results of this study demonstrate the following:


**Acknowledgments:** This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (NRF-2017R1A2B4008157).

**Author Contributions:** Hyunwook Choo, Hwandon Jun and Hyung-Koo Yoon performed field experiments and wrote paper together.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **A CFD Results-Based Approach to Investigating Acoustic Attenuation Performance and Pressure Loss of Car Perforated Tube Silencers**

#### **Hao Zhang, Wei Fan and Li-Xin Guo \***

School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China; peterat2011@foxmail.com (H.Z.); 1410110@stu.neu.edu.cn (W.F.)

**\*** Correspondence: lxguo@mail.neu.edu.cn

Received: 1 March 2018; Accepted: 29 March 2018; Published: 2 April 2018

**Abstract:** This paper proposes an approach to investigating the effect of different temperatures and flow velocities on the acoustic performance of silencers in a more accurate and meticulous fashion, based on steady computational results of the flow field inside the silencer using computational fluid dynamics (CFD). This approach can transfer the CFD results—including temperature and flow velocity distribution—to acoustic meshes by mesh mapping. A numerical simulation on the sound field inside the silencer is then performed, using the CFD results as a boundary condition. This approach facilitates the analysis of complex silencer designs such as perforated tube silencers, and the numerical predictions are verified by a comparison with available experimental data. In the case of the three-pass perforated tube silencer of a car, the proposed approach is implemented to calculate the transmission loss (*TL*) of the silencer at different temperatures and flow velocities. We found that increasing the air temperature shifts the TL curve to a higher frequency and reduces the acoustic attenuation at most frequencies. As the air flow increases, the curve moves to a slightly lower frequency and the acoustic attenuation increases slightly. Additionally, the pressure loss of perforated tube silencers could be calculated according to the total pressure distribution of their inlet and outlet from the steady computational results using CFD.

**Keywords:** perforate tube silencer; transmission loss (*TL*); pressure loss; computational fluid dynamics (CFD); temperature; air flow velocity

#### **1. Introduction**

The silencer is the main device used for suppressing automobile noise. The most important goal for a high-performance silencer is to simultaneously have good aerodynamic and acoustic attenuation characteristics. However, these often contradict each other. Perforated tube silencers are well-balanced in terms of both characteristics. Thus, they are applied extensively to the exhaust systems of automobiles.

Pressure loss is the main indicator used for evaluating the aerodynamic performance of a silencer. It has a negative influence on the efficiency of an engine. If the loss exceeds the backpressure limit, it will result in reduced engine power and an increase in fuel consumption [1]. Currently, the pressure loss is predicted by performing a 3D steady computation using computational fluid dynamics (CFD). This method has been widely accepted among researchers because of its high level of accuracy and adaptability [2–4]. Middelberg et al. [5] computed the pressure loss of a simple expansion chamber muffler using a CFD simulation. Lee et al. [6] predicted the pressure loss of concentric tube silencers with five different patterns of perforated elements using CFD analysis. Ren et al. [7] employed the CFD approach to predict the pressure loss of an exhaust muffler, which was influenced by the insert duct, the position of the baffle, and the inlet air velocity.

Transmission loss (*TL*) is one of the most important indicators used for evaluating the acoustic attenuation performance of a silencer. The gas from an automotive engine has a high temperature and speed, which will have a strong influence on the acoustic attenuation performance. High temperature changes, gas density, sound velocity and acoustic impedance [8], and the effects of flow velocity embody two sides: one is to affect the sound wave propagation in a medium, and the other is to reduce aerodynamic noise resulting from turbulence [9]. Various theoretical and experimental studies were conducted in the past in order to investigate the acoustic attenuation performance of a silencer under the influence of different temperatures and air flow [10–12]. Kim et al. [13] investigated the acoustic characteristics of an expansion chamber with a constant mass flow and a steady temperature gradient. Tsuji et al. [14] applied finite element and boundary element methods to evaluate the acoustic wave transmission characteristics in a medium with uniform and steady mean flow. Kirby [15] compared two numerical methods and two analytical methods of modelling automotive dissipative silencers with a uniform mean gas flow of Mach number (M). The comparison indicated a close similarity to the *TL* predictions that were obtained for the silencers that were examined. With the rapid development of high-performance computers, the 3D numerical simulation method plays an increasingly significant role in sound field analysis. Broatch et al. [16] proposed a 3D time-domain technique based on the CFD approach to calculate the *TL* of exhaust mufflers with different chambers. Sánchez-Orgaz et al. [17] proposed a hybrid finite element approach—combining an acoustic velocity potential formulation in the central airway with a pressure-based wave equation in the outer chamber—to study the *TL* of perforated dissipative silencers with heterogeneous properties in the presence of mean flow. Dong et al. [18] employed the CFD approach to perform a 3D steady computation and obtain a temperature distribution inside a two-pass perforated tube silencer, and then defined the elements with a 5 ◦C difference in temperature as a collection. The air parameters of the corresponding temperatures were assigned to the defined collections using SYSNOISE software, in order to calculate the *TL*.

However, most of the present works associated with predicting the *TL* of perforated tube silencers assume the mean flow to be uniform. Similarly, the temperature fields are simplified as constant or linear fields, which may result in inaccurate analytical results—especially in the case of silencers with complicated temperature and flow velocity distributions. In light of the aforementioned disadvantages, this paper proposes an approach based on CFD results to investigate the acoustic attenuation performance of silencers. This approach comprises the following steps: (1) the flow field inside the silencer is calculated by performing a 3D steady computation using CFD; (2) the CFD results—including temperature and air flow velocity—are transferred to acoustic meshes by mesh mapping, so that the results can be used as a boundary condition of sound field analysis; (3) the *TL* of the silencer is calculated at a high temperature and air flow. Compared to previous works, the proposed approach avoids the use of simplified temperature and flow velocity distributions. Consequently, the effects of temperature variation and air flow on the acoustic attenuation performance of silencers can be calculated and observed in greater detail.

This paper proceeds as follows: Following the introduction, Section 2 describes the proposed computational approach and verifies its accuracy. Section 3 builds an internal fluid model of a perforated tube silencer of a car, and generates CFD meshes and acoustic meshes of the model, respectively. Following this, a 3D steady computation using CFD is performed in order to calculate the pressure loss of the silencer and obtain its temperature and flow velocity distributions. Section 4 employs the proposed approach to investigate the acoustic attenuation performance of the silencer under the influence of different temperatures and air flow. Section 5 concludes the study.

#### **2. Methods**

#### *2.1. Mesh Mapping*

The most critical step in the proposed approach to calculating the acoustic attenuation of silencers is to set a mesh mapping for transferring data. The mesh mapping is used to transfer the node data of the CFD mesh (source mesh) into the acoustic mesh (target mesh). However, there is not usually a one-to-one correspondence between the nodes of different meshes. Thus, an appropriate mapping algorithm should be employed.

A maximum distance algorithm is often applied in Virtual.lab acoustics software to set a mesh mapping between different meshes with the same geometrical shape, but a different density of nodes. The algorithm includes the following two necessary parameters [19]:


The *N* closest nodes to a given source node are used to transfer data to the target node. The data value assigned to the target node is a weighted average of the values at the *N* source nodes. The weights are:

$$\mathcal{W}\_{i} = \frac{1}{d\_{i}} \Big/ \sum\_{i=1}^{N} \frac{1}{d\_{i}}.\tag{1}$$

The transferred value of the target node is then:

$$P\_{\text{Target}} = \sum\_{i=1}^{N} \frac{P\_i^{\text{Source}}}{d\_i} \bigg/ \sum\_{i=1}^{N} \frac{1}{d\_i} \,\prime \tag{2}$$

where *di* is the distance between the source node and the target node, and *P*Source *<sup>i</sup>* is the value of the source node.

For example, when the target node is defined as a center, there are three source nodes (*N* = 3) lying inside a sphere with *R* = 100 mm. Thus, the mapping data transfer relation is depicted in Figure 1. The value at *A* is given by:

$$P\_{\text{Target}} = \frac{\frac{1}{d\_1} P\_1^{\text{Source}} + \frac{1}{d\_2} P\_2^{\text{Source}} + \frac{1}{d\_3} P\_3^{\text{Source}}}{\frac{1}{d\_1} + \frac{1}{d\_2} + \frac{1}{d\_3}}. \tag{3}$$

After transferring data by the mesh mapping method, Virtual.lab software is used to conduct a numerical simulation on the internal sound field of the silencer, and the *TL* is then determined by

$$TL = 10\log\_{10}\left(\frac{W\_1}{W\_2}\right),\tag{4}$$

with

$$\mathcal{W}\_1 = |p\_1|^2 A\_1 / Z\_1 \,\mathcal{W}\_2 = |p\_2|^2 A\_2 / Z\_2. \tag{5}$$

Substituting (5) into (4) yields

$$TL = 10\log\_{10}\left(\frac{p\_1\overline{p\_1}Z\_2}{p\_2\overline{p\_2}Z\_1}\frac{A\_1}{A\_2}\right),\tag{6}$$

where *W*<sup>1</sup> and *W*<sup>2</sup> are the acoustic power of the inlet and outlet of the silencer, respectively; *p*<sup>1</sup> and *p*<sup>2</sup> are the sound pressure of the incident and the transmitted waves, respectively; *p*<sup>1</sup> and *p*<sup>2</sup> are conjugate complex numbers of *p*<sup>1</sup> and *p*2, respectively; *Z*<sup>1</sup> and *Z*<sup>2</sup> are the acoustic impedance of the inlet and outlet of the silencer, respectively; and *A*<sup>1</sup> and *A*<sup>2</sup> are the cross-sectional areas of the inlet and outlet of the silencer, respectively.

**Figure 1.** Sketch of data transfer (Number of nodes *N* = 3, Maximum distance *R* = 100 mm).

#### *2.2. Method Validation*

In order to verify the accuracy of the proposed approach, a straight-through perforated tube silencer with two different perforated patterns that was proposed by the authors of [20,21] and a cross-flow perforated tube silencer that was proposed by the authors of [10] was considered. Previously, Liu et al. [22] employed the 3D time-domain CFD approach to calculate the *TL* of two perforated tube silencers at different flow velocities and temperatures, and their predictions were verified by experimental data. Moreover, they found that the distribution of flow velocity inside the silencer was anisotropic and nonhomogeneous. Therefore, it is not sufficiently accurate to consider the actual gas flow as the mean flow. The approach proposed in this paper avoids this problem effectively.

Figure 2 presents the straight-through perforated tube silencer. The diameters of the inner and outer cavities are *d* = 32 mm and *D* = 110 mm, respectively; the length of the silencer is *l* = 200 mm; the thickness of the wall is 2 mm; and the diameter of the hole and the porosityare *dh* = 4 mm and *σ* = 4.7%, respectively, for Pattern 1, and *dh* = 8 mm and *σ* = 14.7%, respectively, for Pattern 2. The cross-flow perforated silencer is illustrated in Figure 3. The diameters of the inner and outer cavities are *d* = 49.3 mm and *D* = 101.6 mm, respectively; the lengths of tubes on both sides of the baffle are *l*<sup>1</sup> = *l*<sup>2</sup> = 128.6 mm; each tube is perforated with 160 orifices, with a porosity of 3.9%, a diameter of 2.49 mm; and a wall thickness of 0.81 mm.

**Figure 2.** Straight-through perforated tube silencer.

**Figure 3.** Cross-flow perforated tube silencer.

The proposed approach is applied to calculate the *TL* for the silencers under different boundary conditions, and the predictions were compared with the available measurement results. Figure 4a and 4b compare the *TL* curves of the straight-through perforated tube silencer for Pattern 1 and Pattern 2, respectively. Considering the predictions and measurements at an air flow of *M* = 0.1 and a temperature of *T* = 288 K, it is evident that the predictions are in line with the experimental data. The predicted and measured TL curves for Pattern 1, with *M* = 0.2 and *T* = 288 K, are presented in Figure 5, and it is evident that these predictions are also in line with the experimental data. Figure 6 presents a comparison of the numerical results and the measurements for the cross-flow silencer with an air flow velocity of *v* = 17 m/s and *T* = 347 K, which again presents an excellent agreement.

**Figure 4.** Comparison of the predicted and measured transmission loss (*TL*) for the straight-through perforated tube silencer (Mach number *M* = 0.1, Temperature *T* = 288 K): (**a**) Pattern 1; (**b**) Pattern 2.

Based on the above comparisons, it is evident that the proposed approach displayed a high degree of accuracy in investigating the acoustic attenuation performance of perforated tube silencers. This high degree of accuracy should be attributed to the avoidance of using simplified temperature fields and air flow in the calculation process. In doing so, the actual working condition of the silencer can be better reflected. Additionally, it must be stated that the computations may be highly time-consuming because there are too many orifices in the silencers, which increases the amount of mesh.

**Figure 5.** Comparison of the predicted and measured *TL* for Pattern 1 (*M* = 0.2, *T* = 288 K).

**Figure 6.** Comparison of the calculated and measured *TL* for the cross-flow silencer (*v* = 17 m/s, *T* = 374 K).

#### **3. Modeling and Steady Computation Using CFD**

#### *3.1. Modeling*

The three-pass perforated tube silencer of a car that is considered in this paper is presented in Figure 7. The silencer is divided into three chambers by two baffles. The inlet tube and the middle buffer tube are respectively perforated by 40 and 48 orifices with 4 mm diameter, and the porosity of the tubes are 2.7% and 5.8%, respectively. There is no perforation in the outlet tube. Additionally, the axes of the three tubes are not on the same plane. Consequently, it is difficult to build a fluid model directly inside the silencer because of its complicated structure. Therefore, firstly CATIA (version V5R20, Dassault Systemes, Paris, France, 2010) software was used to build its structural model, and then it was filled using ANSYS Workbench (version 14.5, ANSYS Inc., Canonsburg, PA, USA, 2012) software to obtain the fluid model. Following this, the fluid model was split into several parts to generate mesh individually. Tetrahedral mesh was applied to the perforation area and the transition tube, and hexahedral mesh was applied to the remaining parts.

Figure 8a,b illustrate the CFD mesh model and the acoustic mesh model of the fluid model, which will be used for CFD steady computation and acoustic computation, respectively. There is a different node density between the two models. The CFD mesh is refined in the perforation area in order to obtain more accurate results. However, the accuracy of the acoustic computation is determined by the overall mesh. Local mesh refinement can not improve its accuracy, and thus the size of the acoustic mesh element should be kept as uniform as possible.

**Figure 7.** 3D geometrical model for the three-pass perforated tube silencer: (**a**) structural model; (**b**) flow field model.

**Figure 8.** Mesh model for the three-pass perforated tube silencer: (**a**) computational fluid dynamics (CFD) mesh model; (**b**) acoustic mesh model.

The steady flow computation is carried out using the CFD mesh model. The governing equations for pressure velocity coupling—based on a finite volume method—are solved by the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm. Turbulence is examined using the standard k-epsilon model [23]. The fluid material is air—with the density conforming to the ideal gas law—and is considered as incompressible. The boundary conditions are concretely set as follows: (1) the velocity inlet boundary condition is defined at the inlet of the silencer. According to the information provided by the manufacturer, when the engine runs at 5500 r/min, the velocity and temperature of the inlet are 55 m/s and 760 K, respectively; (2) the pressure outlet boundary condition is defined at the outlet of the silencer, and the gauge pressure is 0 Pa—which is relative to one standard atmospheric pressure; and, (3) the walls are assumed to be stationary—with no slip condition—and adiabatic.

#### *3.2. The CFD Results*

Figure 9 shows the temperature distribution in the axial cross sections of the three tubes of the silencer. Overall, the temperature value decreases along the direction of the air flow. There is little difference in temperature between the first and second chambers, where the value is approximately 574 K. However, in the third chamber, the value increases to approximately 670 K. The highest temperature—which is approximately 751 K—occurs in the inlet tube. The temperature in the outlet tube is close to the value of the first chamber, and it is distributed evenly. A greater temperature gradient exists in the perforation area. Additionally, the temperature of the edge chamber is the lowest, which results from the large difference in temperature between the inner and outer walls.

Figure 10 depicts the velocity distribution in the axial cross sections of the three tubes of the silencer. When air flows through the perforation area of the inlet tube and the middle buffer tube, some of the air will enter into the second chamber through the orifices. At this moment, an eddy forms in the perforation area (as shown in Figure 11) because of the disturbance of the airflows with different velocities and flow directions, which may induce a greater velocity gradient in the perforation area. Additionally, there exists a greater velocity gradient next to the nozzle of each tube, resulting from the sudden change of the flow sections.

**Figure 9.** Contour of temperature in the cross section of the silencer.

**Figure 10.** Contour of velocity in the cross section of the silencer.

**Figure 11.** Contour of turbulence kinetic energy in the cross section of the silencer.

Figure 12 illustrates the total pressure distribution in the axial cross sections of the three tubes of the silencer. The highest pressure occurs in the inlet tube—especially in the perforation area—where the value is approximately 23,660 Pa. In the second and third chambers, the pressure is reduced to approximately 19,560 Pa. The middle buffer tube plays the role of the transition region with a greater pressure gradient. Gradually, however, the pressure distribution tends to be uniform in the first chamber. The pressure is distributed more evenly in the outlet tube than in the other regions, and the value is approximately −2313 Pa (the pressure, which is lower than the atmospheric pressure, may be a result of the Venturi effect).

**Figure 12.** Contour of total pressure in the cross section of the silencer.

Based on the CFD results, it is evident that the distributions of temperature and flow velocity are non-uniform. Therefore, it is not accurate enough to consider them as uniform distributions.

#### *3.3. Pressure Loss Calculation*

The pressure loss equals the total pressure difference between the inlet and the outlet of a silencer. Following the method of calculating the pressure loss that was proposed by the authors of [24]—whose accuracy has been validated by experiments—nine evenly distributed points were selected in the cross sections of the inlet and the outlet, respectively. The average total pressure of the nine points was considered to be the total pressure of the corresponding cross section. Therefore, the pressure loss could be expressed as

$$
\Delta p = \overline{p\_{t\_1}} - \overline{p\_{t\_2}} \tag{7}
$$

where *pt*<sup>1</sup> is the total pressure of the inlet, and *pt*<sup>2</sup> is the total pressure of the outlet.

Figure 13 presents the cross section of the silencer in question for the pressure loss calculation. The total pressure distribution—predicted by the steady computation using CFD—is applied to obtain the pressure loss of the silencer. The value of the pressure loss is approximately 37,710 Pa. The following explanation is a plausible explanation for the pressure loss that occurred: When air flow moves from the chamber to the tube—or vice versa—the mechanical energy of the air flow is reduced dramatically because of the sudden change of the flow sections. Moreover, there are large-scale eddy zones present in the perforation area that restrict air flow and lead to energy dissipation, which results from the increased turbulence kinetic energy (as shown in Figure 11).

**Figure 13.** The selected points for the total pressure calculation in the cross sections of the inlet and outlet of the silencer.

#### **4. Sound Field**

#### *4.1. Check the Maximum Frequency Value of the Acoustic Mesh*

When there is mean flow, the maximum frequency value that is allowed by the acoustic mesh is defined by

$$f\_{\text{max}} = \frac{c(1 - M)}{6L},\tag{8}$$

where *c* is the local sound velocity, *M* is the mean flow Mach number, and *L* is the mesh element size [19]. Therefore, the calculated frequency must satisfy the condition *f* ≤ *f*max in order to guarantee accuracy during the acoustic calculation. This means that the maximum size of the acoustic mesh should be small enough to allow at least six elements to fit into the wavelength of the calculated frequency.

The CFD results that were acquired in Section 3.2—including the temperature and flow velocity—are transferred to the acoustic mesh (presented in Figure 8b) by mesh mapping. Following the data transfer, the maximum frequency value of the investigated acoustic mesh was obtained through checking the maximum frequency report in Virtual.Lab. The observed value is 3676.1 Hz. The noise frequency that needs to be calculated for the silencer is usually below 3000 Hz. Therefore, the investigated acoustic mesh could guarantee computational accuracy.

#### *4.2. Effects of Temperature Changes on the TL of the Silencer*

Prior to performing the acoustic computation, the inlet of the silencer is supplied with a unit vibration velocity source. Additionally, an anechoic end duct property is defined at the end of the outlet to simulate the anechoic termination condition. It is known that the existing acoustic mesh is valid up to a frequency of 3676.1 Hz. Therefore, the analysis is conducted from 20 Hz to 2000 Hz, in linear steps of 10 Hz. Additionally, an input point and an output point are chosen from the inlet and the outlet of the silencer, respectively. Figure 14 presents the sound pressure level of the inlet and outlet of the silencer, calculated by the proposed approach. Figure 14a–c present the sound pressure level of the inlet and outlet of the silencer with the same flow velocity but different temperatures.

Based on the sound pressure of the inlet and the outlet, the *TL* of the silencer is calculated using Equation (6). Figure 15 indicates the *TL* predictions with the same flow velocity (55 m/s) and different temperatures. It may be noted that, as the temperature rises, the *TL* curves move to a higher frequency and the acoustic attenuation is reduced at most frequencies. Considering the *TL* curves from Figure 15a as an example, a valley at 850 Hz is moved to 1160 Hz with an increase in temperature from 560 K to 760 K. Similarly, a peak moves from 1240 Hz to 1560 Hz, while their corresponding TL values are reduced by 10 dB and 20 dB, respectively.

**Figure 14.** *Cont.*

**Figure 14.** The sound pressure level of the inlet and outlet of the silencer (*v* = 55 m/s): (**a**) *T* = 760 K; (**b**) *T* = 560 K; (**c**) *T* = 960 K.

**Figure 15.** Effects of temperature on the TL of the silencer (*v* = 55 m/s) : (**a**) 760 K versus 560 K; (**b**) 760 K versus 960 K.

#### *4.3. Effects of Flow Velocity Changes on the TL of the Silencer*

It is evident from the flow velocity distribution in Figure 10 that the velocity inside the silencer is lower. Therefore, it is assumed that the air flow mainly influences the sound wave propagation, and the turbulence noise is neglected [25]. Following the proposed approach, the sound pressure level curves of the inlet and outlet of the silencer—with the same temperature but different flow velocities—are presented in Figure 16.

**Figure 16.** The sound pressure level of the inlet and outlet of the silencer (*T* = 760 K): (**a**) *v* = 20 m/s; (**b**) *v* = 35 m/s.

Figure 17 indicates the effects of increasing the flow velocity on the *TL* predictions at the same temperature (760 K). It can be observed that, as the velocity goes up, the *TL* curves move to a slightly lower frequency and the acoustic attenuation is increased slightly at most frequencies. The variation of the TL curve may be attributed to the fact that the higher velocity increases the acoustic resistance of the perforations and reduces the effective flow area of the orifices [26]. This causes the frequency response of the silencer to differ from a static state. Generally, however, the propagation of the sound wave is only changed dramatically if the air flow Mach number is greater than 0.3 [27]. Therefore, the variation of the curves in Figure 17 is small.

**Figure 17.** Effects of flow velocity on the *TL* of the silencer (*T* = 760 K).

#### **5. Conclusions**

This paper proposes a new approach—based on CFD results—to investigating the effects of temperature and flow velocity on the acoustic attenuation performance of a perforated tube silencer of a car. The validation results suggest that the proposed approach can make accurate predictions. In actual fact, the proposed approach is distinguished from the existing approaches by defining the fluid properties of every acoustic mesh element under the CFD results (including temperature and flow velocity). In doing so, the simplification of the temperature and flow velocity distribution is avoided, thus the effect can be predicted in greater detail. The results indicate that air flow velocity and temperature variation may have an effect on the investigated three-pass perforated tube silencer. An increase in temperature will shift the TL curve to a higher frequency and reduce the acoustic attenuation at most frequencies. As flow velocity increases and the curve is moved to a slightly lower frequency, the acoustic attenuation is enhanced, although the tendency is not particularly obvious. Additionally, the pressure loss of the perforated tube silencer can be calculated through the total pressure distribution from the CFD results. When the engine runs at 5500 r/min, the calculated value is approximately 37,710 Pa.

**Acknowledgments:** This work was supported by the National Natural Science Foundation of China (51275082, 11272273).

**Author Contributions:** Li-Xin Guo conceived and designed the research idea and the framework; Wei Fan performed the simulations; Hao Zhang analyzed the data and wrote the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **A Pseudo-3D Model for Electromagnetic Acoustic Transducers (EMATs)**

**Wuliang Yin 1,\*, Yuedong Xie 2,\*, Zhigang Qu <sup>1</sup> and Zenghua Liu <sup>3</sup>**


Received: 31 January 2018; Accepted: 9 March 2018; Published: 15 March 2018

**Abstract:** Previous methods for modelling Rayleigh waves produced by a meander-line-coil electromagnetic acoustic transducer (EMAT) consisted mostly of two-dimensional (2D) simulations that focussed on the vertical plane of the material. This paper presents a pseudo-three-dimensional (3D) model that extends the simulation space to both vertical and horizontal planes. For the vertical plane, we combines analytical and finite-difference time-domain (FDTD) methods to model Rayleigh waves' propagation within an aluminium plate and their scattering behaviours by cracks. For the horizontal surface plane, we employ an analytical method to investigate the radiation pattern of Rayleigh waves at various depths. The experimental results suggest that the models and the modelling techniques are valid.

**Keywords:** Analytical solutions; FDTD; EMATs; beam directivity

#### **1. Introduction**

A wide group of non-destructive testing (NDT) techniques are commonly used in biomedical industries, such as ultrasonic techniques, electromagnetic techniques, and laser testing [1–4]. Due to the non-contact nature, more and more attention has been paid to the NDT technique with electromagnetic acoustic transducers (EMATs), and EMATs have gradually been used in industrial applications, such as thickness gauging and defect detection [5–8].

A classic EMAT sensor is made of a meander-line-coil and a permanent magnet (Figure 1). There are two major coupling principles for EMATs: magnetostriction is for ferromagnetic metallic materials, and the Lorentz force mechanism is for conductive and ferromagnetic materials [9]. This work focussed on only the Lorentz force mechanism performing on an aluminium plate. The Lorentz force mechanism is: the meander-line-coil placed above the sample generates eddy currents **J** within the sample. A permanent magnet placed above the coil generates a static magnetic field **B** to the sample. The interaction between **J** and **B** produces Lorentz force density **F,** as shown in Equation (1):

$$\mathbf{F} = \mathbf{J} \times \mathbf{B} \tag{1}$$

Substantial works have been reported on EMATs modelling, which comprises an electromagnetic (EM) model and an ultrasonic (US) model [10–12]. The EM model was accomplished by the finite element method (FEM) and the analytical method, while the US model was accomplished with FEM, the finite-difference time-domain method (FDTD), and the analytical solutions. Some of the previous

work modelled EMATs by combining FEM and analytical solutions, i.e., FEM for EM modelling, and analytical solutions for US modelling [10–12]. On the other hand, some of the previous work utilised FEM for both EM and US modelling, i.e., COMSOL (a commercial EM simulation package) for EM modelling, and Abaqus for US modelling [13,14]. Authors have proposed several methods to model EMATs, including a method combining FEM and FDTD, a method combining analytical solutions and FDTD, and a wholly analytical method [15–18].

**Figure 1.** The configuration of a typical meander-line-coil electromagnetic acoustic transducer (EMAT). Reproduced with permission from [15], IEEE, 2016.

Most of EMATs' simulations were two-dimensional (2D), and can only focus on the specific plane. This article is attempting to build a pseudo-three-dimensional (3D) model in order to further study EMATs by combining a surface plane 2D model (the *x-y* cross-section shown in Figure 1) and a vertical plane 2D model (the *y-z* cross-section shown in Figure 1) together. More specifically, the Lorentz force density obtained from the vertical plane of the sample is imported to the surface plane of the sample as the driving source to generate Rayleigh waves. Previously, only the beam directivity at the surface of the sample (*z* = 0 in Figure 1) was investigated [18,19]. However, Rayleigh waves not only distribute along the surface (*z* = 0), they also distribute within a depth of the Rayleigh waves' wavelength. Some industrial defects are within the sample instead of on the surface of the sample; in order to use Rayleigh waves to detect such defects, the beam directivity of the Rayleigh waves at various depths is worth investigating. In this article, the equations to study the beam directivity are more complete compared to the approximate equations presented in [18] by Xie et al., and the beam directivity at various depths are investigated. This study lays a solid industrial foundation for near-surface defects detection using Rayleigh waves, and can be a starting point to build an advanced 3D EMAT model in the future. Except for near-surface detection, the proposed strategy can be used to perform body detection, i.e., to model bulk waves, including longitudinal waves and shear waves. In addition, the proposed 3D EMAT model can be used to characterise other EMAT structures to generate surface waves, such as unidirectional Rayleigh waves EMATs and multiple directional Rayleigh waves EMATs. The 2D simulation on the vertical plane utilizes the analytical method and FDTD; the 2D simulation on the vertical plane and experimental validations are introduced in Section 2. The pseudo-3D model is presented in Section 3 to investigate the radiation pattern of Rayleigh waves at various depths by utilizing a wholly analytical solution. This work is an extension of the work published in [15,18] by Xie et al.

#### **2. Vertical Plane Modelling**

Previously, the authors have conducted EMAT modelling for Rayleigh waves focussing on the vertical plane of the material. This model combines the analytical method and FDTD to model EMATs [15]. The dimension and material of the test piece, the coil, and the permanent magnet are the same as the ones used in [15] by Xie et al. Based on such design, the working frequency used to form the interference of Rayleigh waves is 483 kHz.

#### *2.1. EMAT-EM Model*

This section introduces the EMAT-EM model to analyse the distribution of **F** (Lorentz force density). Firstly, the classic Dodd and Deeds solution [20] to the vector potential is described, and the strategy of adapting the circular analytical solutions for a straight wire is introduced (Section 2.1.1). The FEM is employed to validate the adapted solution (Section 2.1.2). Finally, based on the adapted analytical solutions, the distribution of Lorentz force density is presented (Section 2.1.3).

#### 2.1.1. Adapted Analytical Solutions to the Vector Potential for a Straight Wire

The governing equations for calculating eddy currents are described in Equations (2)–(4):

$$
\nabla^2 \mathbf{A} = -\mu \mathbf{I} + \mu \sigma \frac{\partial \mathbf{A}}{\partial t} + \mu \varepsilon \frac{\partial^2 \mathbf{A}}{\partial t^2} + \mu \nabla(\frac{1}{\mu}) \times (\nabla \times \mathbf{A}) \tag{2}
$$

$$\mathbf{E} = -j\omega\mathbf{A} \tag{3}$$

$$
\mathbf{J} = \sigma \mathbf{E} \tag{4}
$$

where **A** is the vector potential generated by **I**, *ω* and **I** are the angular frequency and the density of the applied alternating current (AC), respectively, *ε*, *μ* and *σ* are the permittivity, permeability and conductivity of the test piece respectively, and **E** and **J** are the induced electric field and eddy current density, respectively [20].

As described in [15] by Xie et al., for a small radius circular coil, the distribution of the vector potential **A** at *z* = 0 (surface of the sample) is not symmetrical with the radius due to the bent wire.

The coil used in this work consists of straight wires; thus, the analytical solution for a straight wire is needed. The adapted solution for a straight wire has been presented in [15] by Xie et al. Here is a brief introduction. A hypothesis is proposed: if the radius of the circular coil, compared with its width, is very large, a bent wire can be viewed as a straight wire, and the distribution of **A** should be symmetrical. To validate such a hypothesis, a model is built with a large-radius circular coil above the aluminium plate. The aluminium sample used has a dimension of 80 mm × 30 mm, and the inner radius and the outer radius of the circular coil are set to 5.0395 m and 5.0405 m, respectively. At kHz, the current density applied to the circular coil is 1 A/m2, and the lift-off of the coil is 1 mm. The permeability and the conductivity of the aluminium plate are 1.257 × <sup>10</sup>−<sup>6</sup> H/m and 3.8 × 107 Siemens/m, respectively.

The distribution of the magnitude of **A** based on the adapted solution is shown in Figure 2. The red marker denotes the maximum magnitude of the vector potential. The distribution of the magnitude of **A** is symmetrical with the radius of 5.04 m, where the coil is located. The result verifies the hypothesis that, when the radius of the circular coil is very large, a bent wire serves as a straight wire.

**Figure 2.** The magnitude distribution of A under a large-radius circular coil. Reproduced with permission from [15], IEEE, 2016.

#### 2.1.2. Comparison between the Adapted Solution and FEM

In order to compare the adapted solutions to FEM, Maxwell Ansoft, which is a FE solver, is utilised. The FEM model has a rectangular cross-sectional coil located above the cross-sectional aluminium plate, and is surrounded by a vacuum region that is four times larger than the sample. The FEM subdivides the large model to smaller elements, and this FEM solves the calculation by minimising the energy error. In this work, when the elements number is beyond 20,000, the energy error is as low as 0.05%, which is sufficiently accurate for the FEM computation. In this work, the mesh number used is 20,395, and the boundary used is a balloon boundary to simulate an infinite space. The distribution of **A** at the sample's surface (*z* = 0) is presented in Figure 3. The analytical solution and FEM present a good agreement at an operational frequency of 1 kHz. However, at a working frequency of 1 MHz, the distribution of **A** from the FEM is not smooth compared to that from the analytical solution; the reason is that the FEM is affected by the elements density and numerical approximation is unavoidable, etc. Therefore, the adapted analytical method presents a more accurate result compared to FEM, especially for a high working frequency.

**Figure 3.** The distribution of A from the adapted analytical solution and the finite element method (FEM). The left curves are the results at 1 kHz, while the right curves are the results at 1 MHz. The red curve is the result from the FEM, and the blue curve is the result from the analytical solutions. Reproduced with permission from [15], IEEE, 2016.

#### 2.1.3. EMAT-Lorentz Force Calculation

The analytical solution to a straight wire has been described in Section 2.1.2. **A** (vector potential) generated by a meander-line-coil is the addition of **A** generated by every single straight wire. The zone where the meander-line-coil mainly operates on is selected to model EM simulation to increase the modelling effectiveness. The distribution of **A** and **F** at *z* = 0 are shown in Figures 4 and 5. The generated periodic fields have different directions for any neighbouring wires, since their applied ACs are opposite, and therefore, the periodic distribution of **A** and **F** has six positive values and six negative values. The outermost **A** and outermost **F** are largest, because **A** is under the outermost wires, and thus is only determined by the fields on one side.

**Figure 4.** For a meander-line-coil, the distribution of the vector potential **A** at *z* = 0. Reproduced with permission from [15], IEEE, 2016.

**Figure 5.** For a meander-line-coil, the distribution of the Lorentz force density **F** at *z* = 0. Reproduced with permission from [15], IEEE, 2016.

#### *2.2. EMAT-US Simulation*

#### 2.2.1. Elastodynamic Equations

Elastodynamic equations (Equations (5) and (6)) are a group of partial differential equations to model the wave propagation:

$$\rho(\mathbf{x})\frac{\partial v\_i}{\partial \mathbf{t}}(\mathbf{x},t) = \sum\_{j=1}^d \frac{\partial T\_{ij}}{\partial \mathbf{x}\_j}(\mathbf{x},t) + f\_l(\mathbf{x},t) \tag{5}$$

$$\frac{\partial T\_{\bar{l}}}{\partial \mathbf{t}} = \sum\_{j=1}^{d} \sum\_{i=1}^{d} c\_{i\bar{j}kl}(\mathbf{x}) \frac{\partial \upsilon\_{k}}{\partial \mathbf{x}\_{l}} + \theta\_{l\bar{j}}(\mathbf{x}, t) \tag{6}$$

where *ρ* and *cijkl* are the density and the fourth stiffness tensor of the material, and *fi* and *θij* are the force source and strain tensor rate source, respectively.

#### 2.2.2. Combination of EMAT-EM and EMAT-US Models

As described in Section 2.1.3, **F** is obtained from the EMAT-EM calculation. In this section, **F**, which is used as the force source, is imported to the EMAT-US model to produce ultrasound (Figure 6). Since **F** is calculated in the frequency domain and FDTD is a time-domain solver, the excitation signal for the EMAT-US model is a time sequence signal with the peak equalling the peak values of **F**. The excitation signal used is a Gaussian-modulated sinusoidal with a fractional width of 0.18. A crack and a receiver R are located within the sample, as shown in Figure 6. Regarding the FDTD setup in the ultrasonic model, there are two main parameters: the spatial step, and the time step. The spatial step used is 0.2 mm, which approximately equals to 1/30th of the wavelength. The time step is 0.0222 μs, which is calculated based on the Courant–Friedrichs–Lewy (CFL) condition. Free surface conditions are utilised on the surface of the sample. Perfectly-matched layers (PML) with a thickness of 16 mm are utilised to absorb ultrasound.

**Figure 6.** On the vertical plane of the material, the combination of the EMAT-electromagnetic (EM) and EMAT-ultrasonic (US) models. Reproduced with permission from [15], IEEE, 2016.

#### 2.2.3. Wave Propagations

Based on FDTD calculations, the velocity fields of ultrasound waves are obtained. The velocity fields at 27 μs and 83 μs are shown in Figure 7, which describes the Rayleigh waves' propagation and their scattering behaviours, respectively. The white arrows denote the propagation path.

**Figure 7.** Wave propagations. (**a**,**b**) denote the velocity fields at 27 μs and 83 μs, respectively. Reproduced with permission from [15], IEEE, 2016.

#### *2.3. EMAT-Reception Simulation*

The EMAT reception process has been reported quite a lot [21]. The received signal from simulations is shown in Figure 8, where DRW and RRW denote the directly transmitted Rayleigh waves and the reflected Rayleigh waves, respectively. The propagation distance of DRW and RRW is 100 mm and 300 mm, respectively, and the velocity of Rayleigh waves is 2.93 mm/μs; hence, the theoretically arrival time of DRW and RRW is 34 μs and 102.4 μs, respectively. Figure 8 shows the numerically arrival time of DRW and RRW; these numerical and experimental results present a good agreement.

**Figure 8.** The received signal from simulations. Reproduced with permission from [15], IEEE, 2016.

#### *2.4. Experimental Validations*

Experiments were carried out to validate the proposed modelling method. The setup of the experiments is the same as that of our previous work [15]. The received signals from experiments are shown in Figure 9, where three signals are observed. The "Main bang" is the interference signal due to a high power excitation, arriving before DRW and RRW. The red curve and the blue curve are the envelope and the time series signal, respectively.

**Figure 9.** The received signal from experiments. The blue curve denotes the induced voltage in the received coil, and the red curve denotes the envelope of the received signal. Reproduced with permission from [15], IEEE, 2016.

Figure 10 shows the envelope of DRW and RRW from experiments and simulations. The experimental arrival times of DRW is 34 μs, which is consistent with the numerical results. However, the experimental arrival time of RRW is slightly different from the simulations. This is because of the approximated model used and the inevitable experimental noise.

**Figure 10.** The envelope of the received signals. The blue curve and the red curve are the envelope of the received signal from simulations and experiments, respectively.

#### **3. Horizontal Surface Plane Modelling—Directivity Analysis of Rayleigh Waves**

The main mode of propagation on the horizontal surface of the material is the Rayleigh wave. The radiation pattern and the beam directivity of Rayleigh waves at *z* = 0 (Figure 1) were reported in [18] by Xie et al. Rayleigh waves not only distribute along the surface (*z* = 0), they also distribute within a depth of the Rayleigh waves' wavelengths. However, Rayleigh waves not only concentrate on the surface, they also concentrate within a depth of one wavelength. In this section, the beam directivity of Rayleigh waves at various depths are investigated utilising an analytical method.

#### *3.1. The Analytical Solution to the Displacement of Rayleigh Waves*

N. A. Haskell proposed an analytical solution to Rayleigh waves [22,23]. Ref. [24] by Love introduced these solutions in detail, and investigated the beam directivity of Rayleigh waves at *z* = 0 with approximated equations. However, Rayleigh waves not only concentrate on the surface; they also concentrate within a depth of one wavelength. The complete equations of the Rayleigh waves' displacement at various depths are:

$$u\_{\mathbf{r}} = A(\mathbf{x}, r, z)e^{-\frac{\eta i}{4}} \frac{2\kappa(\gamma - 1)}{v\_{\beta}} \mathbf{F}(\frac{2}{\gamma} - 1 + \frac{\gamma - 1}{\gamma} e^{-z(v\_{\mathbf{a}} - v\_{\beta})}) \tag{7}$$

$$\mu\_z = \frac{i\gamma\upsilon\_a\mu\_r}{\kappa(\gamma - 1)}\tag{8}$$

where:

$$A(\mathbf{x}, r, z) = \frac{\mathbf{x}^2 \gamma v\_{\beta}}{4\rho(\frac{2\gamma^2 v\_{\mathbf{z}} v\_{\beta}}{\mathbf{K}^3})} \sqrt{\frac{2}{\tau \mathbf{K} r}} e^{(-i\mathbf{K} r - z v\_{\beta})} \tag{9}$$

$$\gamma = \cos\left(\theta\right) \tag{10}$$

$$v\_{\rm tr} = \begin{cases} \frac{\sqrt{\kappa^2 - \left(\omega/c\_L\right)^2}}{i\sqrt{\left(\omega/c\_L\right)^2 - \kappa^2}} & \mathbf{x} > \omega/c\_L\\ \ i\sqrt{\left(\omega/c\_L\right)^2 - \kappa^2} & \mathbf{x} < \omega/c\_L \end{cases} \tag{11}$$

$$w\_{\beta} = \left\{ \begin{array}{c} \sqrt{\mathbf{\kappa}^2 - \left(\omega/c\_S\right)^2} \\ i\sqrt{\left(\omega/c\_S\right)^2 - \mathbf{\kappa}^2} \end{array} \quad \begin{array}{c} \mathbf{\kappa} > \omega/c\_{\mathbf{s}} \\ \mathbf{\kappa} < \omega/c\_{\mathbf{s}} \end{array} \tag{12}$$

$$\mathbf{x} = \frac{\omega}{\mathcal{L}\_R} \tag{13}$$

where *ur* and *uz* are the in-plane and the out-of-plane displacement to be calculated, respectively; *r* is the distance between the source point and the field point (Figure 11); *F* is the excitation source; *ρ* is the density of the material; *θ* is the angle between the force vector and the in-plane displacement vector; *ω* is the operational angular frequency; z is the depth; and *cL*, *cS*, *cR*, and Ë are the velocity of the longitudinal waves, shear waves, Rayleigh waves, and the wave number, respectively.

**Figure 11.** The model used to simulate Rayleigh waves. Reproduced with permission from [18], Elsevier, 2017.

#### *3.2. Linking EMAT-EM and EMAT-US Models*

On the surface plane of the material, the combination between the EMAT-EM model and the EMAT-US model is described in Figure 12. The calculated *F* values, acting as the excitation source, are imported to each surface layer at various depths. The Rayleigh waves' distribution is the addition of the Rayleigh waves generated by each point source. Table 1 illustrates the parameters used for the EMAT-US model.

**Figure 12.** On the surface plane of the material, the transformation between the EM and US models. Reproduced with permission from [18], Elsevier, 2017.



#### *3.3. Analysis of the Beam Directivity of Rayleigh Waves*

Figure 13 shows the calculated Rayleigh waves' radiation pattern at *z* = 0, which is symmetrical with the center of the coil. Rayleigh waves are mainly generated along the y direction (main lobe) and some undeniable directions (side lobe). The characteristics of Rayleigh waves are quantitatively investigated by means of beam directivity, as shown in the red arc in Figure 13 (*r* = 250 mm, *θ*<sup>1</sup> = −70◦, and θ<sup>2</sup> = 70◦).

**Figure 13.** The radiation pattern of Rayleigh waves generated by the meander-line-coil EMAT. Reproduced with permission from [18], Elsevier, 2017.

The beam directivity of Rayleigh waves at *z* = 0 is shown in Figure 14. A main lobe contains a larger magnitude compared to the side lobes, which contain a smaller magnitude. The main lobe is centered at 0◦, and the side lobes are roughly centered at −25.5◦, −18◦, −10◦, 10◦, 18◦, and 25.5◦, respectively. The largest magnitude of the side lobes is 25.87% that of the main lobe. In most applications, side lobes are usually undesirable.

**Figure 14.** The simulated beam directivity of Rayleigh waves at *z* = 0. Reproduced with permission from [15], IEEE, 2016.

The beam directivity of the Rayleigh waves at various depths is shown in Figure 15. The magnitude of the beam directivity is normalised. From this image, the magnitude of the Rayleigh waves decays with the depth, especially for the depth larger than one Rayleigh wavelength. At a depth equalling to one Rayleigh wave's wavelength, *z* = 6 mm, the magnitude of the Rayleigh waves is 34.9% of that at *z* = 0. At a depth of 7 mm, the magnitude of the Rayleigh waves decays to 22.8% of that at the surface of the test piece (*z* = 0). This observation confirms that Rayleigh waves mainly distribute within a depth equal to one Rayleigh wavelength.

**Figure 15.** The simulated beam directivity of Rayleigh waves at various depths.

#### *3.4. Experimental Validations*

In Figure 15, based on the analytical simulations, the distribution of the Rayleigh waves at various depths is presented. In this part, the measured results at *z* = 0 are picked to compare with the simulation results. The experimental setup was the same as the one in Section 2.4. The measured beam directivity at *z* = 0 is obtained by placing the receiver along the scanning path, as shown in Figure 16. Figure 17 shows the beam directivity results at *z* = 0 from simulations and experiments. Thirty-three measuring points on the scanning path were picked with a moving step of 2.5◦. The measured beam directivity at *z* = 0 agrees well with the simulated beam directivity, which validates the proposed method. There are some non-overlapping points between these two curves due to several factors, which include the experimental noise and the tolerance of the receiver's position, etc.

**Figure 16.** The scanning path of the receiver.

**Figure 17.** The comparison between the simulated beam directivity and the experimental beam directivity at *z* = 0. The red curve is the beam directivity from simulations, while the blue curve is the beam directivity from experiments.

#### **4. Conclusions**

A pseudo-3D model for simulating meander-line-coil EMATs was proposed. A method combining the analytical method for the EM model and FDTD for the US model was utilised to simulate the Rayleigh waves' propagation. On the other hand, a wholly analytical method was utilised to simulate the radiation pattern of the Rayleigh waves. For both cases, analytical solutions to the EM model were adapted from the classic Dodd and Deeds solution in order to calculate eddy currents under straight wires. By comparing with the FEM, the analytical solutions are more accurate. Experiments were conducted in order to validate the proposed method, and these showed a good consistency. The beam directivity of Rayleigh waves at various depths were investigated, and results confirmed that Rayleigh waves mainly distribute within a depth of one wavelength of Rayleigh waves. Overall, this pseudo-3D model combines both the vertical plane and surface plane of EMAT models, and provides the beam directivity of Rayleigh waves at various depths, which have not been reported previously. Therefore, this work can be a starting point to build an advanced EMAT 3D model in the future. There are some limitations of the proposed 3D EMAT model. Firstly, it can only be applied in a homogeneous medium. It is worth investigating the 3D EMAT model for a multiple-layer medium, since multiple-layer samples are widely used in applicable industries. Secondly, the EMAT-US model on the vertical plane of the sample uses an approximated model with only point sources (Lorentz force density) to generate Rayleigh waves. A more detailed model with volume sources within the skin depth is worth considering in the future. In addition, it is worth investigating the scattering behaviours of defects in other orientations, as it is a variable in practical application.

**Acknowledgments:** This work was financially supported by Engineering and Physical Sciences Research Council (Grant No. EP/M020835/1).

**Author Contributions:** Wuliang Yin conceived and reviewed the manuscript, did some simulations and designed the experiments. Yuedong Xie performed the experiments and wrote this manuscript. Zhigang Qu contributed some simulation results. Zenghua Liu provided related instruments for carrying out experiments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Calculation of Noise Barrier Insertion Loss Based on Varied Vehicle Frequencies**

**Haibo Wang 1,2, Peng Luo <sup>3</sup> and Ming Cai 1,2,\***


Received: 14 December 2017; Accepted: 5 January 2018; Published: 11 January 2018

**Abstract:** A single frequency of 500 Hz is used as the equivalent frequency for traffic noise to calculate the approximate diffraction in current road barrier designs. However, the noise frequency changes according to the different types of vehicles moving at various speeds. The primary objective of this study is the development of a method of calculating the insertion loss based on frequencies. First, the noise emissions of a large number of vehicles classified by speed and type were measured to obtain data the noise spectrum. The corresponding relation between vehicle type, speed, and noise frequency was obtained. Next, the impact of different frequencies on the insertion loss was analyzed and was verified to be reasonable in experiments with different propagation distances compared to the analysis of a pure 500 Hz sound. In addition, calculations were applied in a case with different traffic flows, and the effect of a road noise barrier with different types of constituents and flow speeds were analyzed. The results show that sound pressure levels behind a barrier of a heavy vehicle flow or with a high speed are notably elevated.

**Keywords:** noise barrier; insertion loss; vehicle frequencies; diffraction; flow speed

#### **1. Introduction**

The rapidly increasing number of vehicles is in the world, especially in many developing countries, has raised the serious problem of traffic noise [1–3]. Traffic noise is disturbing to the daily routine of residents along the roads [4]. Many issues, physiological [5] and psychological [6] health problems, are demonstrated to be related to a noisy environment. To reduce traffic noise, one of the most effective methods is to set noise barriers along the roads [7].

Current research on noise barriers primarily focuses on noise calculation [8,9], the effects of barrier shapes [10,11], and the design of barriers [12–15]. As one of the important aspects, prediction of the insertion loss by the noise barrier has been attracting the attention of scholars, and both mathematical and experimental approaches have been developed. The mathematical formula was first developed by Sommerfeld in 1896 [16] and has been subsequently developed by many other scholars. As such, rigorous diffraction results can be calculated in mathematical methods, e.g., the Boundary Element Method [17], the Finite Element Method (FEM) [18], and the Finite Difference Method (FDM) [19]. In many cases, the typical procedure for the mathematical computation is too complicated to be applicable to the engineering application. Alternatively, comprehensive sets of data have been measured to plot the experimental attenuation of noise barrier since 1940 [20]. The most famous data set is the Maekawa chart [21], and many experimental formulae were developed based on Maekawa's original chart [22,23], including Kurze and Anderson's derived formula [24]. Kurze and Anderson's formula is widely used and is regarded as the simplest experimental practice for application in China and many other countries. Fresnel number, which is the single parameter in Kurze and

Anderson's diffraction formulae, is the ratio between the path length difference and half of the sound wavelength [25,26].

In practical application, the path length difference of a noise barrier is defined by a particular source-barrier-receiver geometry and easily measured. However, because of the multiple frequency components of traffic noise, the computation of diffraction is always time consuming. The noise frequency changes according to the different types of vehicles travelling at various speeds. Hence, the noise protection performance of the barrier depends on the traffic flow states. For simplicity, an equivalent frequency of 550 Hz can be set to represent the total frequency band in the calculation of diffraction suggested in the FHWA traffic noise model [27]. In the 1/3rd-octave-band spectrum, the frequency of 500 Hz is used in the approximate calculation of diffraction and noise reduction effect of a sound barrier [9,28–31]. However, the data used was collected from 1993 and 1995, and the characteristics of vehicle noise have changed since then because of the variety of vehicle types and speeds [32]. Moreover, in the micro traffic noise prediction model, the reduction effect of a noise barrier should consider every vehicle to be point source. Therefore, a single 500 Hz frequency is not adequate in the noise reduction calculation of a sound barrier.

In this paper, data of the noise spectrum of different vehicles in a variety of speed ranges were measured to calculate the equivalent frequency. The corresponding relationships among vehicle type, speed, and noise frequency were obtained. Several path length differences were preset to calculate the 1/3rd-octave-band diffraction and the total diffraction. Next, verification experiments were established to ensure the calculation accuracy of barrier diffraction based on different equivalent frequencies. Last, an application of the approach to different traffic flows is implemented, and the effects of a barrier to traffic noise for different speeds and different types of constituents are analyzed.

#### **2. Methods**

#### *2.1. Measurement Method of Vehicle Noise Spectrum Data*

The measurements were conducted on the roadside of seven typical dry roads in Guangzhou, China, with a surface constructed of asphalt concrete. The parameters of surface could refer to the Chinese standard JTG F40-2004 (Technical Specification for Construction of Asphalt Pavements) [33]. To avoid possible disturbances, the samples were collected in places far away from intersections, rivers, and populated areas. The chosen urban roads have the following speed limits: 2 roads with 70 km/h, 3 roads with 60 km/h, 1 road with 50 km/h, and 1 road with 40 km/h. The spectrum was captured when the A-weighted sound pressure level during a vehicle pass-by was at its maximum, and the sampling frequency of sound level meter was 23 Hz, i.e., 23 spectra per second. Only one vehicle was measured at a time. In Figure 1, the locations of the measurement sites are shown. Both sites were 1.2 m high and 6 m to the roadside with a distance of 50 m to each other. In the experiment, the subject vehicles only moved from left to right alone. As the single vehicle passed by measurement site 1, the noise spectrum data were recorded by a digital noise recorder. The vehicle speed had little variation between points 1 and 2 while traveling on a long-straight road with a very sparse traffic flow. The vehicle speed was detected by a laser speedometer at measurement site 2.

In this paper, the vehicles are classified into three categories: heavy cars, middle cars, and light cars. The 1/3rd-octave-band sound pressure level (SPL) and the total SPL are included in the noise spectral data. All the noise data are A-weighted.

**Figure 1.** Locations of the Measurement Spots.

In this paper, 1351 sets of valid data were collected, among which 973 sets were light car data, 166 sets were middle car data, and 212 sets were heavy car data. The distribution of the samples is shown in Table 1.

**Table 1.** Distribution of the samples of vehicle noise spectra.


#### *2.2. Noise Spectra of Different Vehicle Types and Speeds*

To analyze the frequency characteristics of traffic noise of different types of vehicle with various speed, the 1/3rd octave band spectrum data of vehicles were recorded, including 34 frequencies from 10 Hz to 20,000 Hz. To facilitate the analysis of the frequency characteristics, the following ranges are defined: frequencies from 10 Hz to 315 Hz are low frequency; 400 Hz to 800 Hz are medium-low frequency; 1000 Hz to 2500 Hz are medium-high frequency; 3150 Hz to 20,000 Hz are high frequency. To analyses the frequency characteristics of traffic noise and exclude the effect of SPL, the noise energy percentage of the 1/3rd octave band spectrum was calculated [34], as given by the following formulas (Equations (1)–(3)):

$$E(i,j) = 10^{0.1 \ast L(i,j)} \tag{1}$$

$$E'(i,j) = \frac{E(i,j)}{\sum\_{j} E(i,j)}\tag{2}$$

$$E\_S(j) = \frac{\sum E'(i, j)}{i} \tag{3}$$

where *i* is the number of vehicles, *j* is the number of frequencies; *L*(*i*, *j*) is the SPL (dB), *E*(*i*, *j*) is the relative value of noise energy (dimensionless quantities), *E* (*i*, *j*) is the noise energy percentage (%), and *ES*(*j*) is the mean percentage of noise energy (%).

From the equations, the SPL of each frequency of each vehicle is transformed to the noise energy percentage of each frequency. Figure 2 depicts the average noise energy percentage versus frequency for each type of vehicle with different speeds.

**Figure 2.** Noise energy percentage with the 1/3rd octave band spectrum frequencies at different speeds of three types of vehicles. (**a**) light vehicles; (**b**) middle vehicles; (**c**) heavy vehicles.

As shown in Figure 2, the noise energy percentage of the different speeds of three types of vehicles were performed with the 1/3rd octave band spectrum frequencies in the range of 10–20,000 Hz. The following points can be made:

A. For a light vehicle, the noise energy is concentrated in the range of 500–2500 Hz, with peak frequency at approximately 1250 Hz. For a middle vehicle, the noise energy is concentrated in the range of 315–2500 Hz, with peak frequency at approximately 1600 Hz. For a heavy vehicle, the noise energy is concentrated in the range of 315–2500 Hz, with peak frequency at approximately 500 Hz.


#### *2.3. The Calculation Method of Insertion Loss*

The center frequency with the minimum mean deviation between the 1/3rd-octave-band diffraction and the total diffraction is the equivalent frequency [36]. The equivalent frequency could be used in the approximate calculation of barrier attenuation. The results of the approximate calculation are close to the direct calculation using the whole spectrum as parameters, and the total error is always less than 1 dB. The procedures for calculating the insertion loss in this paper are presented below. Generally, several path length differences were preset first to calculate the 1/3rd-octave-band diffraction and the total diffraction. The detailed calculations are shown as follows:


$$
\Delta L\_{di} = 20 \lg \frac{\sqrt{2\pi N}}{\tanh \sqrt{2\pi N}} + 5 \text{dB}, N = \frac{\delta \cdot f}{170} \tag{4}
$$

where *δ* is the path length difference and *f* is the 1/3rd-octave-band center frequency;

(4) For a certain path length difference, the total SPL *L* with the reduction effect of noise barrier can be calculated as (Equations (5) and (6))

$$L\_i{}^{\prime} = L\_i - \Delta L\_{di} \tag{5}$$

$$L' = 10 \lg(\sum\_{i} 10^{L\_i'/10}) \tag{6}$$

where *Li* is the 1/3rd-octave-band SPL with the reduction effect of barrier;

(5) The following equation is used to yield the total diffraction Δ*Ld*(Equation (7)):

$$
\Delta L\_d = L - L'\tag{7}
$$

(6) After calculating the total diffraction Δ*Ld*(*δj*) and the 1/3rd-octave-band diffraction Δ*Ldi*(*δj*) with seven preset path length differences, the mean deviation Δ*di* between Δ*Ld*(*δj*) and Δ*Ldi*(*δj*) can be obtained (Equation (8)).

$$
\Delta d\_{\bar{i}} = \frac{1}{7} \sum\_{j=1}^{7} \left| \Delta L\_d(\delta\_{\bar{j}}) - \Delta L\_{d\bar{i}}(\delta\_{\bar{j}}) \right| \tag{8}
$$

#### **3. Experimental Verification**

Two verifying experiments were conducted: an insertion loss experiment and a distance attenuation experiment. In the experiments, vehicle noise was recorded next to the roads and subsequently played in the audio amplifier to emulate a point sound source that was used to measure actual insertion loss of noise barrier. For every type of car at different speed range, the synthesized point source consisted of five recorded audio samples. The length of each audio sample is 8 s.

The insertion loss experiment was conducted to collect sound level pressure data in front of and behind the noise barrier. A schematic of the experiment is shown in Figure 3.

**Figure 3.** Insertion loss experiment.

Two noise recorders were installed, one at measurement site 1 and the other at measurement site 2, to collect noise data and compute the noise levels *L*<sup>1</sup> and *L*2, respectively. The insertion loss of barrier is the difference between *L*<sup>1</sup> and *L*2. The computed equivalent noise levels are shown in Table 2.


**Table 2.** Noise levels computed in the insertion loss experiment (dB).

To obtain the actual value of diffraction, noise attenuation as a function of distance between the two measurement sites should be considered. The noise data were measured in an open space near the sound barrier to avoid possible errors introduced by the change of the experimental locations. Because noise attenuation with distance is independent of the frequency of the sound source [20], the sound source audios were merged into three categories in the experiment: heavy car, middle car, and light car. The geometric configuration of measurement is shown in Figure 4.

**Figure 4.** Noise attenuation experiment with distance.

Noise recorders were installed in measurement site 1 and measurement site 2. The data are summarized in Table 3. Noise attenuation with distance for heavy, middle, and light vehicles are 11.16 dB, 11.03 dB, and 11.13 dB, respectively. The average of the three is defined as Δ*L*. In this instance, Δ*L* = 11.11 dB.

**Table 3.** Data collected in the noise attenuation experiment with distance (dB).


From the experimental verification, the actual value of diffraction by barrier was determined as follows (Equation (9)):

$$
\Delta L\_d = (L\_1 - L\_2) - \Delta L \tag{9}
$$

The theoretical value of diffraction is defined as (Equation (10))

$$
\Delta L'\_d = \begin{cases}
20 \text{lg} \frac{\sqrt{2\pi N}}{\tanh\sqrt{2\pi N}} + 5dB & \text{, } N > 0 \\
 & 5dB & \text{, } N = 0
\end{cases}
\tag{10}
$$

where *N* is the Fresnel number, *N* = *<sup>δ</sup>*· *<sup>f</sup>* <sup>170</sup> .

Equations (9) and (10) are mainly applied to near-field conditions. According to reference [20], Equation (11) should be satisfied (Equation (11)).

$$D < 2d^2/\lambda \tag{11}$$

where *D* is the distance from the sound source to the receivers' center, *d* is the step of the receivers, and *λ* is the wave length. The speed of sound is 340 m/s.

In the experiment, *D* is calculated as approximately 5.26 m, *d* is set as 0.5 m in environmental acoustics calculation, and the wave length *λ* is in a range of [2.125, 0.054] with a primary frequency range of 160 Hz to 6300 Hz. Thus, Equation (11) is satisfied. The experimental setup can be considered with the theoretical formulation of the Fresnel number *N*.

Two types of equivalent frequencies were used to yield the theoretical value of diffraction: one was the 500 Hz frequency, and the other was the varied equivalent frequencies computed in this paper. The actual value and the theoretical values are shown in Figure 5. The figure reveals that the error of diffraction calculated with the variety of noise equivalent frequencies is less than that calculated with 500 Hz frequency alone. Under the same experimental conditions, the average error using the 500 Hz frequency is approximately 2.3 dB, whereas the average error with the variety of noise equivalent frequencies is approximately 0.9 dB. With a difference of approximately 1.4 dB, the accuracy of the proposed method in calculating the noise barrier diffraction is undoubtedly more justified.

**Figure 5.** Comparisons between 500 Hz frequency and the variety of equivalent frequencies for diffraction by barrier. (**a**) light vehicle; (**b**) middle vehicle; (**c**) heavy vehicle.

#### **4. Effects of a Road Noise Barrier with Different Flow States**

As shown in Figure 6, there is a two-dimensional sound field involving a sufficiently long enough sound barrier and a line road traffic noise source that is parallel to both the sound barrier and the ground. Three receivers (R1, R2, and R3) are chosen to analyze the effect of a road noise barrier with different type constituents and flow speed; all the receivers are sheltered by the barrier. The insertion losses of all receivers are calculated with different flow speed and the proportion of heavy vehicles. For each point, the SPL is calculated. The insertion loss of barrier can be defined as (Equation (12))

$$L\_D = SPL\_0 - SPL\_b \tag{12}$$

where *SPL*<sup>0</sup> is the SPL at the receiver when the barrier is absent and *SPLb* is the SPL at the same receiver when the barrier is present.

**Figure 6.** A case of noise calculation behind a road sound barrier.

#### *4.1. Effect of Vehicle Type*

Traffic noise caused by heavy vehicles is quite different in frequencies compared to other vehicles based on the results of Section 2.2. Hence, the effects of the proportion of heavy vehicles on insertion loss of barrier were analyzed. The speed of the vehicles in this part is set as 50 km/h, and the flow of 1000 vehicles per hour consists of light and heavy vehicles with 11 different proportions. To ensure that the results are not affected by the different vehicles' sound source strong, the emission intensity of each vehicle source is set as 85 dB. Traffic noise for different heavy vehicle proportions (0, 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 90%, and 100%) on each point was calculated.

With different proportions of heavy vehicles, the insertion losses at three points are shown in Figure 7. *LD* presents a distinct pattern regarding points R1, R2, and R3, and the insertion loss at all of the receivers have relationships with the proportion of heavy vehicles. For point 1, *LD* and heavy vehicles *a* (%) follow the linear relationship of *LD* = −0.0293*a* + 24.661, i.e., 10 percent of heavy vehicles can cause a 0.29 dB decline on insertion loss in this scene. The results show that the SPL behind a barrier of a heavy vehicle flow is approximately 2.3 dB higher than that a light vehicle flow with the same source emission intensity. The same rule is also applicable to the overall shadow area. The heavy vehicles have more acoustical constituents at low frequencies compared to light vehicles, which are more significantly shaded by barriers. The insertion loss of heavy vehicle flow is the highest, followed by the insertion loss of mixed traffic flow, and the insertion loss of light vehicle flow is the lowest; this result indicates that the lower-frequency sound can bypass the barrier more easily. Thus, flow control of heavy vehicles whose sound is concentrated in low frequencies is an effective measure to improve the acoustic environment near roadways.

**Figure 7.** Insertion loss of the road barrier with different proportions of heavy vehicles.

#### *4.2. Effect of Flow Speed*

Vehicle noise characteristics of various speeds are different in frequencies according to the results of Section 2.2. Hence, the effects of different traffic flow speeds on the insertion loss of a noise barrier were analyzed. The flow of 1000 vehicles per hour consists of light and heavy vehicles with seven different flow speeds. The emission intensity of each vehicle source is set as 85 dB. Traffic noise with different traffic flow speeds (30 km/h, 40 km/h, 50 km/h, 60 km/h, 70 km/h, 80 km/h, and 90 km/h) on each point was calculated. Light vehicle flow, middle traffic flow, heavy vehicle flow, and mixed traffic flow (with 30% heavy vehicles) were considered in this section.

With different type constituents and speeds of traffic flow, the insertion losses at three points are shown in Figure 8. From the results of the calculated *LD* values, the analysis is as follows:


**Figure 8.** *Cont*.

**Figure 8.** Insertion loss of barrier with different flow speeds. (**a**) Point R1; (**b**) Point R2; (**c**) Point R3.

#### **5. Conclusions**

A method for insertion loss calculation of a road noise barrier based on frequencies was presented in this paper. The proposed method realizes a more accurate calculation compared to the commonly used method based on a single 500-Hz sound. Several path length differences were preset to calculate the 1/3rd-octave-band diffraction and the total diffraction. Using different weights of spectrum according to the experimental noise data and traffic flow state, the insertion loss can be accurately calculated.

The method was verified by experiments with approximately 1.4 dB improvement in accuracy. The varied equivalent frequencies can be applied in the microscopic simulation of traffic noise, in which all different types of cars are considered to be the point source in the road network.

The noise protection of a road noise barrier with different flow states was analyzed based on a case study. The sound pressure level behind a barrier of a heavy vehicle flow is approximately 2.3 dB higher than that of a light vehicle flow with the same source emission intensity. The sound pressure level behind a barrier of a common mixed traffic flow with high speed is approximately 2 dB lower than the level with a low speed with the same source emission intensity. A road with high speed and high light vehicle proportion flow is recommended to establish a noise barrier.

**Acknowledgments:** This work was supported by the National Natural Science Foundation of China (No. 11574407) and the Fundamental Research Funds for the Central Universities (No. 17lgpy55).

**Author Contributions:** This paper is a result of the full collaboration of all the authors. Ming Cai and Haibo Wang conceived and designed the experiments; Peng Luo performed the experiments; Haibo Wang and Peng Luo analyzed the data; Haibo Wang wrote the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **The Boundary Element Method in Acoustics: A Survey**

#### **Stephen Kirkup**

School of Engineering, University of Central Lancashire, Preston PR1 2HE, UK; smkirkup@uclan.ac.uk; Tel.: +44-779-442-2554

Received: 23 January 2019; Accepted: 9 April 2019; Published: 19 April 2019

**Abstract:** The boundary element method (BEM) in the context of acoustics or Helmholtz problems is reviewed in this paper. The basis of the BEM is initially developed for Laplace's equation. The boundary integral equation formulations for the standard interior and exterior acoustic problems are stated and the boundary element methods are derived through collocation. It is shown how interior modal analysis can be carried out via the boundary element method. Further extensions in the BEM in acoustics are also reviewed, including half-space problems and modelling the acoustic field surrounding thin screens. Current research in linking the boundary element method to other methods in order to solve coupled vibro-acoustic and aero-acoustic problems and methods for solving inverse problems via the BEM are surveyed. Applications of the BEM in each area of acoustics are referenced. The computational complexity of the problem is considered and methods for improving its general efficiency are reviewed. The significant maintenance issues of the standard exterior acoustic solution are considered, in particular the weighting parameter in combined formulations such as Burton and Miller's equation. The commonality of the integral operators across formulations and hence the potential for development of a software library approach is emphasised.

**Keywords:** boundary element method; acoustics; Helmholtz equation

#### **1. Introduction**

Acoustics is the science of sound and has particular applications in sound reproduction, noise and sensing [1,2]. Acoustics may be interpreted as an area of applied mathematics [1,3]. In special cases, the mathematical equations governing acoustic problems can be solved analytically, for example if the equations are linear and the geometry is separable. However, for realistic acoustic problems, numerical methods provide a much more flexible means of solution [4,5]. Numerical methods are only useful in practice when they are implemented on computer. It is over the last sixty years or so that numerical methods have become increasingly developed and computers have become faster, with increasing data storage and more widespread. This brings us to the wide-ranging area of computational acoustics; the solution of acoustic problems on computer [6–11], which has significantly developed in this timescale.

In the context of this work, the acoustic field is assumed to be present within a fluid domain. If there is an obstacle within an existing acoustic field, then the disturbance it causes is termed *scattering*. If an object is vibrating and hence exciting an acoustic response, or contributing to an existing acoustic field, then this is termed *radiation*. If the fluid influences the vibration of the object or structure, and energy is exchanged in both directions between them, then this is termed coupled fluid-structure interaction [12] or, in the acoustic context, vibro-acoustics [13,14]. Alternatively, if the acoustic field is an outcome of a background flow, for example the generation of noise by turbulence, then this is termed aero-acoustics [11,15]. The determination of the properties of a vibrating or scattering object (for example, its shape, the surface impedance or the sound intensity) from acoustic measurements in the field is termed inverse acoustics [16].

Vibration analysis is not the focus of this work, but it clearly cannot be divorced from the study of acoustics. In vibration analysis, the domain of the structure oscillates, perhaps under variable loading or forcing. In the case of vibro-acoustics, the vibration can excite or be excited by the acoustic field and, in any analysis, both need to be modelled [17–19], and the equations coupled. The computational modelling of vibration is normally carried out by the finite element method (FEM) [20,21], and this method is significantly established in this application area. Vibrational modelling may be carried out in the time domain. However, with the nature of structural vibration in that it is normally dominated by modes and their corresponding resonant frequencies, and hence vibration is often close to periodic. Similarly, the excitation forces on a structure, such as the driver on a loudspeaker of the explosions in an engine can similarly be close to periodic in short time scales. It is, therefore, found that vibratory problems are routinely analysed in the frequency domain, both in the practical work and computational modelling.

An acoustic field is most straightforwardly interpreted as a sound pressure field, with the sound pressure varying over the extent of the domain, and with time. The transient acoustic field can be computationally modelled by the finite element method [22,23], the finite difference—time domain method (a particular finite difference method that was originally developed for electromagnetic simulation [24–26]) can also be applied in acoustics [27–30]) and the boundary element method [31–34], but acoustics and vibration are more often analysed and modelled in the frequency domain. In acoustic problems, the most likely fluid domains are air or water and, in many cases in these fluids, the linear wave equation is an acceptable model. By observing one frequency at a time, the wave equation can be simplified as a sequence of Helmholtz equations. Again, there are a variety of methods for the numerical solution of Helmholtz problems; the finite difference method, the finite element method [23,35–43] and, the subject of this paper, the boundary element method.

The boundary element method is one of several established numerical or computational methods for solving boundary-value problems that arise in mathematics, the physical sciences and engineering [44]. Typically, a boundary-value problem consists of a domain within or outside a boundary in which a variable of interest or physical property is governed by an identified partial differential equation (PDE). A computational method for solving a boundary-value problem is tasked with finding an approximation to the unknown variable within the domain. Mostly, this is carried out by domain methods, such as the finite element method [45], in which a mesh is applied to the domain. However, the boundary element method works in an entirely different way in that in the first stage further information is found on the boundary only; the solution at the domain points is found in the second stage by integrating over the known boundary data. The boundary element method requires a mesh of the boundary only, and hence is generally easier to apply than domain methods. The number of elements in the BEM is therefore expected to be much less than in the corresponding finite element method (for the same level of required accuracy or element size), and there is therefore often a potential for significant efficiency savings. The BEM is not as widely applicable as the domain methods; when problems are non-linear, for example, the application of the development of a suitable BEM requires significant further adaption [46–49]. Typically, the boundary element method has found application in sound reproduction modelling, such as loudspeakers [50,51], sonar transducers [52] and in modelling noise from vehicles [53–56] and, more recently, aircraft [57–59].

For the acoustic boundary element method to be accessible, and hence widely used, it has to be implemented in software. This was precisely the rationale behind the development of the software [60] and the monograph [61], the latter also serving as a manual. Several texts on the same theme were published at about the same time [62,63], following on from the earlier collection of works [64]. A chapter of a recent book contains a modern introduction to the acoustic boundary element method [65].

Implementing the acoustic BEM as software continues to be challenging in terms of scoping, the choice of sub-methods and efficiency. The method requires matrices to be formed with each component being the result of an integration, with some of the integrals being singular or hypersingular. The standard BEM requires the solution of a linear system of equations that is formed from the matrices, or, if the BEM is used for modal analysis, a non-linear eigenvalue problem. There have been significant reliability issues with the BEM for exterior problems. These challenges have been met, but much of the research is focused on future-proofing the method, with the increasing expectations in scaling-up and maintaining reasonable processing time. There is a continual desire to progress to higher resolutions of elements, particularly as this is necessary for modelling high frequency problems.

#### **2. Acoustic Model**

In this Section, the underlying acoustic model of the wave equation that governs the sound pressure in the domain is stated. It is shown how the model can be revised into a sequence of Helmholtz problems for periodic signals. The classes of domains that can be solved by boundary element methods are summarised and a generalised boundary condition is adopted. The other acoustic properties, such as sound intensity, radiation efficiency and sound power that that are mostly used in exterior 'noise' problems, are defined. Traditionally, acoustics properties are presented in the decibel scale, and it is shown how they are converted.

#### *2.1. The Wave Equation and the Helmholtz Equation*

The acoustic field is assumed to be present in the domain of a homogeneous isotropic fluid and it is modelled by the linear wave equation,

$$
\nabla^2 \Psi(\mathbf{p}, t) = \frac{1}{c^2} \left. \frac{\partial^2}{\partial t^2} \right| \Psi(\mathbf{p}, t), \tag{1}
$$

where Ψ(**p**, *t*) is the scalar time-dependent velocity potential related to the time-dependent particle velocity *V*(**p**, *t*) by *V*(**p**, *t*) = ∇Ψ(**p**, *t*) and *c* is the propagation velocity (**p** and *t* are the spatial and time variables). The time-dependent sound pressure Q(**p**, *t*) is given in terms of the velocity potential by *<sup>Q</sup>*(**p**, *<sup>t</sup>*) = <sup>−</sup>ρ<sup>∂</sup> <sup>∂</sup>*<sup>t</sup>* Ψ(**p**, *t*) where ρ is the density of the acoustic medium.

The time-dependent velocity potential Ψ(**p**, *t*) can be reduced to a sum of components each of the form

$$\Psi^{\prime}(\mathbf{p},t) = \text{Re}\{\boldsymbol{\varrho}(\mathbf{p})\,\,\boldsymbol{e}^{-i\omega t}\},\tag{2}$$

where ω is the angular frequency (ω = 2πν, where ν is the frequency in hertz) and ϕ(**p**) is the (time-independent) velocity potential. The substitution of the above expression into the wave equation reduces it to the Helmholtz (reduced wave) equation:

$$
\nabla^2 \varphi(\mathbf{p}) + k^2 \varphi(\mathbf{p}) = 0,\tag{3}
$$

where *<sup>k</sup>*<sup>2</sup> <sup>=</sup> <sup>ω</sup><sup>2</sup> *<sup>c</sup>*<sup>2</sup> and *<sup>k</sup>* is the wavenumber. The complex-valued function <sup>ϕ</sup> relates the magnitude and phase of the potential field.

Similarly, the components of the particle velocity have the form *<sup>V</sup>*(**p**, *<sup>t</sup>*) <sup>=</sup> *Re* <sup>∇</sup>ϕ(**p**) *<sup>e</sup>*−*i*ω*<sup>t</sup>* . Often the boundary normal velocity *v*(*p*) is given as a condition or is required and this is defined as follows,

$$
\psi(\mathbf{p}) = \nabla \varphi(\mathbf{p}) \cdot \mathfrak{n}\_{\mathcal{P}} = \frac{\partial \varphi(\mathbf{p})}{\partial \mathfrak{n}\_{\mathcal{P}}},\tag{4}
$$

where *np* is the unit normal to the boundary at *p*.

#### *2.2. Acoustic Properties*

To carry out a complete solution, the wave equation is written as a series of Helmholtz problems, through expressing the boundary condition as a Fourier series with components of the form in Equation (2). For each wavenumber and its associated boundary and other conditions, the Helmholtz

equation is then solved. The time-dependent velocity potential Ψ(**p**, *t*) can then be constituted from the separate solutions, but it is more usual that the results are considered in the frequency domain.

The sound pressure *p*(**p**) at the point **p** in the acoustic domain is one of the most useful acoustic properties, and it is related to the velocity potential by the formula *p*(**p**) = *i*ωρϕ(**p**). In practice, the magnitude of the sound pressure is measured on the decibel scale in which it is evaluated as the sound pressure level as 20 log10 *p*(*p*) √ 2*p*∗ , where *<sup>p</sup>*<sup>∗</sup> is the reference pressure of 2 <sup>×</sup> <sup>10</sup>−<sup>5</sup> Pa.

Particularly for 'noise' problems, the sound power, the time-averaged sound intensity and radiation efficiency are often considered to be useful properties. The normal sound intensity *I*(**p**) at points *p* on a boundary is defined by the formula

$$I(\mathbf{p}) = \frac{1}{2} \text{Re} \langle \overline{p}(\mathbf{p}) v(\mathbf{p}) \rangle,\tag{5}$$

where *p* represents the complex conjugate of *p*. The sound power *W* is an aggregation of the sound intensity to one value by direct integration,

$$\mathcal{W} = \int\_{H} I(\mathbf{q}) dH\_{\mathbf{q}}.\tag{6}$$

where *H* is a boundary. The sound power is also often expressed in decibels as the sound power level as 10 log10 *W W*∗ , where *<sup>W</sup>*<sup>∗</sup> is the reference sound power of 10−<sup>12</sup> watts. The radiation ratio is defined as *<sup>W</sup>* <sup>1</sup> <sup>2</sup> ρ*c <sup>H</sup> <sup>v</sup>*∗(*q*)*v*(*q*)*dHq* .

#### *2.3. The Scope of the Boundary Element Method in the Solution of Acoustic*/*Helmholtz Problems*

In applying the boundary element method to acoustic or Helmholtz problems, the user is in effect adopting a model for the boundary/ies and domain(s), the nature of which determines the integral equation(s) that is/are employed. Traditionally, the BEM has been developed to solve the acoustic or Helmholtz problem interior or exterior to a closed boundary in the standard physical domains are two-dimensional, three-dimensional and axisymmetric three-dimensional space. The basis of the method is the integral equations that arise through applying Green's theorems (direct BEM) or by using layer potentials (indirect BEM). The exterior problem has received much more attention, because of its value in solving over an infinite domain from a surface mesh and because of the difficulty in resolving its reliability problem. Domain methods, such as the finite element method, may also be applied; more elements are required, but the matrices are sparse and structured and the FEM is more established than the BEM in general. If the FEM is applied to exterior problems then techniques such as infinite elements [66] can be used to complete the outer mesh or the perfectly matched layer may be used to absorb outgoing waves [67–69]. The solution of the interior and exterior acoustic problem by the BEM is considered in Sections 4.1 and 4.3.

An outline of the equations that arise in the finite element and related methods is given in Section 5.3.1. In domain methods like the FEM, whether it is applied to a structure or an enclosed fluid, a linear eigenvalue problem is a natural consequence. With the FEM, it is routine to extract the modes and resonant frequencies and use the modal basis to determine the response under excitation. From that perspective, it is natural to include the eigenvalue problem within the boundary element library. However, the application of the BEM leads to a non-linear eigenvalue problem, which requires special solution techniques and hence the construction of solutions through the modal basis is not a natural pathway in the BEM. Acoustic modal analysis via the BEM is considered in Section 4.2.

Although significantly restrictive, one of the simplest acoustic models is the Rayleigh integral [70]. The model is that of a vibrating plate set in an infinite rigid baffle and radiating into a half-space. The Rayleigh integral relates the velocity potential or sound pressure at any point to the velocity map on the plate. As the Rayleigh integral shares its one operator with the integral equation formulations in the boundary element method, the Rayleigh integral method (RIM) [71] can be adopted into the

boundary element method fold. The Rayleigh integral is a special case of half-space problems and these are considered in Section 5.1

An acoustic model that uses the same operators—and hence fits in the BEM context—is that of a shell discontinuity. For example, this can be used to model the acoustic field around a thin screen or shield. The integral equation formulation for the shell are derived from those used in the traditional BEM, but taking the limit as the boundary becomes thinner [72,73], leaving a model that relates a discontinuity in the field across the shell. Shell elements in acoustics are considered in Section 5.2.

The development and application of the BEM from the models described are reviewed in this paper. Although this is a fairly exhaustive list of basic models as things stand, hybrid models can also be developed by using superposition or applying continuity and they will also be considered. The BEM may be applied in vibro-acoustics, aero-acoustics and inverse acoustics and these areas are discussed in Sections 5.3–5.5. Hybrid boundary element methods that include half-space formulation are covered in Section 5.1.2 and shells in Section 5.2.2.

There have been recent developments in the development of the BEM for periodic structures, with noise barriers being the usual application [74–81]. This special case will not be developed further in this paper.

#### *2.4. Boundary Conditions*

To maintain reasonable generality in the software, the author has generally worked with the boundary condition of the following general (Robin) form

$$
\alpha(p)\varphi(p) + \beta(p)\upsilon(p) = f(p),
\tag{7}
$$

with the condition that α(*p*) and β(*p*) cannot both be zero at any value of *p*. This model includes the Dirichlet boundary condition by setting β(*p*) = 0, the Neumann boundary condition by setting α(*p*) = 0 and an impedance condition by setting *f*(*p*) = 0. In the boundary condition (7), *p* is any point on the boundary and α, β and *f* are complex-valued functions that may vary across the boundary. Although the boundary condition for the shell is an adaption of this, this generalised boundary condition model is apparently achievable and seems to cover most current expectations. An explanation of typical boundary conditions that occur in acoustics can be found in Marburg and Anderssohn [82].

For exterior problems, it is also necessary to introduce a condition at infinity, in order to ensure that all scattered and radiated waves are outgoing. This is termed the Sommerfeld radiation condition. In two-dimensional space the condition is

$$\lim\_{r \to \infty} r^{\frac{1}{2}} \left( \frac{\partial \!\!\!\! \!\! \! \! / \partial \!\!\! \! / \partial \!\! / \partial \!\!\! / \partial \!\!\! / \partial \!\!\! / \partial \!\!\! / \partial \!\!\! / \partial \!\!\! / \partial \!\!\! / \partial \!\!\! / \partial \!\!\! / \partial \!\!\!\/) = 0$$

and

$$\lim\_{r \to \infty} r \left( \frac{\partial \varphi}{\partial r} - ik\varphi \right) = 0. \tag{9}$$

in three dimensions. A thorough discussion on the Sommerfeld radiation condition can be found in Ihlenburg [43] (pp. 6–8).

#### **3. The Boundary Element Method and the Laplace Equation Stem**

Laplace's equation is the special case of the Helmholtz Equation (3) with *k* = 0,

$$
\nabla^2 \varphi(\mathbf{p}) = 0.\tag{10}
$$

Although Laplace's equation models many phenomena, it is not directly useful in acoustic modelling. However, many of the issues that have to be tackled in solving the Helmholtz equation by the boundary element method are also found with Laplace's equation. The methods applied in developing the BEM for Laplace's equation also can be developed to be used in the Helmholtz problems. As such, applying the BEM to Laplace's equation is a useful entry point in initiating work on the Helmholtz equation. For example, much of the author's work in solving Helmholtz problems has been underpinned by corresponding work on Laplace's equation [83–87]. Thus, in this section, the BEM is communicated in its simplest, but still realistic and practical context as a stepping-stone on the journey to the full scope of the BEM in acoustics, the subject of this paper. The theoretical basis is set out in Kellogg [88] and development of integral equation methods in Jaswon and Symm [89]. In this Section, the boundary element method is derived for Laplace's equation, and this forms the foundation for its application to acoustic problems. In terms of comunication, it is helpful to sustain a notational conformance, that continues through this paper, and is helpful in relating mathematical expressions and their discrete equivalent as software components. Integration methods for finding the matrix components are considered. Interesting and useful properties of some of the operators and matrices are revealed.

#### *3.1. Elementary Integral Equation Formulation for the Interior Problem*

The boundary element method is not based on the direct solution of the PDE, but rather its reformulation as an integral equation. Historically, the integral equation reformulation of the PDE has followed two distinct routes, termed the direct method and the indirect method. The direct method is based on Green's second identity

$$\int\_{D} \left( \wp(q) \nabla^{2} \psi(q) - \psi(q) \nabla^{2} \wp(q) \right) \, dV\_{q} = \int\_{S} \wp(q) \frac{\partial \psi(q)}{\partial n\_{q}} - \psi(q) \frac{\partial \wp(q)}{\partial n\_{q}} \, dS\_{q} \tag{11}$$

where ϕ and ψ are twice-differentiable scalar function in a domain *D* that is bounded by the closed surface *<sup>S</sup>*. If <sup>ϕ</sup> is a solution of Laplace's equation, <sup>∇</sup>2<sup>ϕ</sup> = 0, then

$$\int\_{D} \varrho(\boldsymbol{q}) \nabla^{2} \psi(\boldsymbol{q}) \, dV\_{\boldsymbol{q}} = \int\_{S} \varphi(\boldsymbol{q}) \frac{\partial \psi(\boldsymbol{q})}{\partial \boldsymbol{n}\_{q}} - \psi(\boldsymbol{q}) \frac{\partial \wp(\boldsymbol{q})}{\partial \boldsymbol{n}\_{q}} \, dS\_{\boldsymbol{q}}.\tag{12}$$

Green's third identity can be derived from the second identity by choosing ψ(*q*) = *G*(*p*, *q*), where *G* is a Green's function. A Green's function is a fundamental solution of the partial differential equation, in this case Laplace's equation, that is the effect or influence a unit source at the point *<sup>p</sup>* has at the point *<sup>q</sup>*, and is defined by <sup>∇</sup>2*G*(*p*, *<sup>q</sup>*) = <sup>−</sup>δ(*<sup>p</sup>* <sup>−</sup> *<sup>q</sup>*). For the two-dimensional Laplace equation *<sup>G</sup>*(*p*, *<sup>q</sup>*) = <sup>−</sup> <sup>1</sup> <sup>2</sup><sup>π</sup> ln(*r*) and for three-dimensional problems *<sup>G</sup>*(*p*, *<sup>q</sup>*) <sup>=</sup> <sup>1</sup> <sup>4</sup>π*<sup>r</sup>* where *<sup>r</sup>* <sup>=</sup> *p* − *q* . The substitution of the Green's function into Equation (12) gives

$$
\int\_{D} \boldsymbol{\wp}(\boldsymbol{q}) \nabla^{2} \boldsymbol{G}(\boldsymbol{p}, \boldsymbol{q}) \, dV\_{\boldsymbol{q}} = \int\_{S} \boldsymbol{\wp}(\boldsymbol{q}) \frac{\partial \boldsymbol{G}(\boldsymbol{p}, \boldsymbol{q})}{\partial \boldsymbol{n}\_{\boldsymbol{q}}} - \boldsymbol{G}(\boldsymbol{p}, \boldsymbol{q}) \frac{\partial \boldsymbol{\wp}(\boldsymbol{q})}{\partial \boldsymbol{n}\_{\boldsymbol{q}}} \, dS\_{\boldsymbol{q}}
$$

$$
$$

or

$$\text{Hence, as a result of the properties of the Dirac delta function.}$$

$$\int\_{S} \varphi(\boldsymbol{q}) \frac{\partial G(\boldsymbol{p}, \boldsymbol{q})}{\partial \boldsymbol{n}\_{\boldsymbol{q}}} - G(\boldsymbol{p}, \boldsymbol{q}) \frac{\partial \varphi(\boldsymbol{q})}{\partial \boldsymbol{n}\_{\boldsymbol{q}}} \, dS\_{\boldsymbol{q}} = \begin{cases} -\varphi(\boldsymbol{p}) \text{ if } \boldsymbol{p} \in D \\\ 0 \text{ if } \boldsymbol{p} \in E \\\ -c(\boldsymbol{p}) \varphi(\boldsymbol{p}) \text{ if } \boldsymbol{p} \in S \end{cases},\tag{14}$$

where *c*(*p*) is the angle (in 2D, divided by 2π) or solid angle (in 3D, divided by 4π) subtended by the interior domain at the boundary point *p*. (Similarly, for exterior problems, *c*(*p*) similarly relates to the exterior angle). If the boundary is smooth at *p* then *c*(*p*) = <sup>1</sup> <sup>2</sup> , and, for simplicity, this is the value that will be used for the remainder of this paper.

In the indirect method, ϕ is presumed to be related by a layer potential σ on the boundary

$$\varphi(p) = \int\_{S} G(p, q) \, \sigma(q) \, dS\_{\mathbb{Q}} \quad (p \in D \cup S). \tag{15}$$

Differentiating Equation (15) with respect to a unit outward normal to a point on the boundary that that passes through *p*, gives the following equation:

$$\frac{\partial \phi(p)}{\partial n\_p} = \int\_S \frac{\partial G(p, q)}{\partial n\_p} \, \sigma(q) \, dS\_q \quad (p \in D). \tag{16}$$

As *p* approaches the boundary, the integral operator has a 'jump' similar to the direct method (14) resulting in the following equation:

$$\frac{\partial \rho(p)}{\partial n\_p} = \int\_S \frac{\partial G(p, q)}{\partial n\_p} \, \sigma(q) \, dS\_q + \frac{1}{2} \sigma(p) \quad (p \in S). \tag{17}$$

The BEM can be derived from the equations in this Section. The equation defined on the boundary (the last one in Equation (14) for the direct method, Equations (15) and (17) for the indirect method (*p* ∈ *S*)), the boundary integral equations, are used to find the unknown functions from the known data on the boundary. The corresponding integrals defined in the domain (the first one in Equation (14) for the direct method and Equation (15) for (*p* ∈ *D*) in the indirect method) return the solution at the chosen domain points.

#### *3.2. Operator Notation and Further Integral Equations for the Laplace Problem*

In this Subsection, operator notation is introduced and further integral equations are introduced in order to illuminate some useful properties. Operator notation provides a shorthand and is an aid to communication. The Laplace integral operators are defined as follows:

$$
\langle L\mu \rangle\_{\Gamma}(p) = \int\_{\Gamma} G(p, q) \,\mu(q) dS\_{q} \,\,\,\,\tag{18}
$$

$$\{\mathcal{M}\mu\}\_{\Gamma}(p) = \int\_{\Gamma} \frac{\partial G(p, q)}{\partial n\_q} \, \mu(q) d\mathcal{S}\_{q} \, \, \, \, \, \tag{19}$$

$$\left\{\mathcal{M}^{t}\mu\right\}\_{\Gamma}(p;\upsilon\_{\mathcal{P}}) = \frac{\partial}{\partial \upsilon\_{\mathcal{P}}} \int\_{\Gamma} G(p,\mathfrak{q})\,\mu(\mathfrak{q})dS\_{\mathfrak{q}}\,\,\,\,\tag{20}$$

where Γ is a boundary (not necessarily closed), *n<sup>q</sup>* is the unique unit normal vector to Γ at *q*, *v<sup>p</sup>* is a unit directional vector passing through *p* and μ*(q)* is a function defined for *q* ∈ Γ. With this notation, the Equations (14), which form the basis of the elementary direct method for the interior problem, can be written as

$$\{M\rho\}\_S(p) - \langle Lx \rangle\_S(p) = \begin{cases} -\varphi(p) \text{ if } p \in D \\ 0 \text{ if } p \in E \\ -\frac{1}{2}\varphi(p) \text{ if } p \in S \end{cases} \tag{21}$$

Similarly, the equations for the indirect method ((15) and (17)) in operator notation for the interior problem are as follows:

$$
\varphi(p) = \langle L\sigma \rangle\_S(p) \quad (p \in D \cup S), \tag{22}
$$

$$w(p) = \begin{cases} M^t \sigma \big|\_{\Gamma} \big( p; n\_p \big) + \frac{1}{2} \sigma(p) \quad (p \in S). \tag{23}$$

Further integral equations for the direct formulation can be obtained by differentiating (21) (as in the derivation of Equations (16), (17)) and for the indirect formulation by introducing a *double* layer potential. Although they are generally unnecessary in solving Laplace problems, the equations resulting from differentiating the elementary integral equations or using double layer potentials are often used in the exterior acoustic problem, considered in the next Section. These equations are also useful in general and they will illuminate some important aspects. For example, differentiating Equation (21) with respect to the outward normal to the boundary,

$$\frac{\partial}{\partial n\_p} \{M\varphi\}\_S(p) - \frac{\partial}{\partial n\_p} \{Lv\}\_S(p) = -\frac{1}{2} \frac{\partial \varphi(p)}{\partial n\_p} (p \in S),\tag{24}$$

which can be written in operator notation,

$$
\langle \{N\boldsymbol{\varphi}\}\_{\mathcal{S}} \{\boldsymbol{p}; \boldsymbol{n}\_{\mathcal{P}}\} - \left\{\boldsymbol{M}^{t}\boldsymbol{\upsilon}\right\}\_{\mathcal{S}} \{\boldsymbol{p}; \boldsymbol{n}\_{\mathcal{P}}\} = -\frac{1}{2}\boldsymbol{\upsilon}(\boldsymbol{p})(\boldsymbol{p}\in\mathcal{S}),\tag{25}
$$

in which a new operator, *N*, has been introduced,

$$
\langle N\mu \rangle\_{\Gamma}(p; v\_{\mathbb{P}}) = \frac{\partial}{\partial v\_{p}} \int\_{\Gamma} \frac{\partial G(p, q)}{\partial n\_{q}} \, \mu(q) dS\_{q}. \tag{26}
$$

Further indirect integral equations may be obtained through presuming the field is the result of a double layer potential

$$
\varphi(p) = \langle M\zeta \rangle\_S(p) \quad (p \in D), \tag{27}
$$

$$w(p) = \langle \mathcal{N}\zeta \rangle\_S(p; n\_p) \quad (p \in S), \tag{28}$$

and, taking into account the jump discontinuity,

$$
\varphi(p) = \langle M\check{\varsigma}\rangle\_S(p) - \frac{1}{2}\check{\varsigma}(p) \quad (p \in S). \tag{29}
$$

The integral equations for the exterior Laplace problem may be derived in the same way. In general, the equations are the same as for the interior problem, except for a change of sign. The equations that make up the direct formulation are

$$\{M\varphi\}\_{\mathbb{S}}(p) - \langle Lv \rangle\_{\mathbb{S}}(p) = \begin{cases} \varphi(p) \text{ if } p \in E \\ 0 \text{ if } p \in D \\ \frac{1}{2}\varphi(p) \text{ if } p \in \mathbb{S} \end{cases} \tag{30}$$

$$\{N\varphi\}\_{\mathcal{S}}\Big(p;n\_{\mathcal{P}}\big) - \left\{M^{t}\upsilon\right\}\_{\mathcal{S}}\Big(p;n\_{\mathcal{P}}\big) = \frac{1}{2}\upsilon(p) \text{ if } p \in \mathcal{S} \text{ .}\tag{31}$$

The equivalent of equations for the indirect method ((15) and (17)) in operator notation for the exterior problem are similar, with just a couple of sign changes to indicate that the jump discontinuity in *M* and *Mt* in the limit from the exterior, rather than the interior:

$$
\varphi(\mathfrak{p}) = \langle L\sigma \rangle\_S(\mathfrak{p}) \quad (\mathfrak{p} \in E \cup S),
\tag{32}
$$

$$\sigma(p) = \begin{cases} M^t \sigma \big|\_{\mathbb{S}} \big( p; n\_p \big) - \frac{1}{2} \sigma(p) \quad (p \in \mathbb{S}). \tag{33} \\ \end{cases} \tag{33}$$

$$
\varphi(\mathfrak{p}) = \langle M\zeta \rangle\_S(\mathfrak{p}) \quad (\mathfrak{p} \in E), \tag{34}
$$

$$
v(p) = \langle \mathcal{N}\check{\mathsf{U}}\rangle\_{\mathcal{S}}(p; n\_p) \quad (p \in \mathcal{S}),\tag{35}$$

$$
\varphi(p) = \langle M\check{\varsigma}\rangle\_S(p) + \frac{1}{2}\check{\varsigma}(p) \quad (p \in S). \tag{36}
$$

#### *3.3. Derivation of the Boundary Element Method*

There are several integral equation methods that can be applied to transform integral equations into boundary element methods, but the most popular and straightforward method is that of collocation [90]. The development of the boundary element method from the selected integral equation requires that the boundary is represented by a set of *panels* or a mesh. An integral equation method is applied to solve the boundary integral equation. The solution in the domain can then be achieved by effecting the appropriate integration over the boundary.

#### 3.3.1. Boundary Element Approximation

In the spirit of the author's previous work [61,91], in order to maintain generality, the discrete forms of the Laplace operators (18)–(20), (26) are sought, in order to effectively become a software component. Let the boundary Γ in Equation (18) be represented by the approximation Γ˜, a set of *n* panels:

$$
\Gamma \approx \tilde{\Gamma} = \sum\_{j=1}^{n} \Delta \tilde{\Gamma}\_{j} \, , \tag{37}
$$

The boundary function μ is replaced by its equivalent μ˜ on the approximate boundary Γ˜:

$$\langle L\mu\rangle\_{\Gamma}(p) \approx \langle L\bar{\mu}\rangle\_{\Gamma}(p) = \int\_{\tilde{\Gamma}} G(p,q) \,\, \bar{\mu}(q) dS\_{\tilde{q}} = \sum\_{j=1}^{n} \int\_{\Delta \tilde{\Gamma}\_{j}} G(p,q) \,\, \bar{\mu}(q) dS\_{\tilde{q}}.\tag{38}$$

In general, the function on the boundary is replaced by a sum of a set of basis functions. The simplest approximation is that of approximating the boundary functions by a constant on each panel:

$$\sum\_{j=1}^{n} \int\_{\Lambda\widetilde{\Gamma}\_{j}} G(p,q) \, \left[ \mathfrak{a}(q) dS\_{q} \approx \sum\_{j=1}^{n} \int\_{\Lambda\widetilde{\Gamma}\_{j}} G(p,q) \, \left[ \mathfrak{a}\_{j} dS\_{q} = \sum\_{j=1}^{n} \, \up|\_{j} \langle \mathrm{L\!} \mathfrak{} \rangle\_{\Lambda\widetilde{\Gamma}\_{j}}(p),\tag{39} \right]$$

where *e*˜ = 1.

For example, for the simple boundary integral Equation (22) with *p* ∈ *S*,

$$\varphi(\boldsymbol{p}) = \langle \boldsymbol{L}\boldsymbol{\sigma} \rangle\_{\mathcal{S}}(\boldsymbol{p}) \approx \sum\_{j=1}^{n} \vec{\sigma}\_{j} \langle \boldsymbol{L}\boldsymbol{\tilde{x}} \rangle\_{\Delta \tilde{S}\_{j}}(\boldsymbol{p})\,. \tag{40}$$

The most common method of solving boundary integral equations is collocation, in which a linear system is formed through setting *p* to take the value of each collocation point in turn:

$$\varphi\_{Si} = \varphi(p\_{Si}) = \langle L\sigma \rangle\_S(p\_{Si}) \approx \sum\_{j=1}^n \sigma\_j \langle L\vec{\sigma} \rangle\_{\Delta \widetilde{S}\_j}(p\_{Si}) \text{ for } i = 1, 2, \dots, n \text{ .} \tag{41}$$

3.3.2. Solution by Collocation

The linear system of approximations (41) may be written in matrix-vector form

$$
\underbrace{\mathcal{Q}}\_{\text{S}} \approx L\_{\text{SS}} \underline{\mathcal{Q}}\_{\text{S}} \, \, \, \, \tag{42}
$$

where <sup>ϕ</sup>*<sup>S</sup>* <sup>=</sup> ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ ϕ*S*<sup>1</sup> ϕ*S*<sup>2</sup> : : ϕ*Sn* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ , σ*<sup>S</sup>* = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ σ*S*<sup>1</sup> σ*S*<sup>2</sup> : : σ*Sn* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ , and [*LSS*]*ij* = {*Le*˜} Δ*S j pSi* . For example, for the Dirichlet

problem in which <sup>ϕ</sup>*<sup>S</sup>* is known, the solution of system (42) (now as an equation relating approximate

values) returns an approximation <sup>σ</sup> *<sup>S</sup>* to <sup>σ</sup>*S*. To find the solution at a set of *<sup>m</sup>* points in the domain the integral (22) is evaluated at the points *pDi* ∈ *D* for *i* = 1, 2, ... *m*;

$$\varphi\_{Di} = \varphi(p\_{Di}) = \langle L\sigma \rangle\_D(p\_{Di}) \approx \sum\_{j=1}^n \tilde{\sigma}\_j \langle L\vec{x} \rangle\_{\Delta \overleftarrow{S}\_j}(p\_{Di}) \text{ for } i = 1, 2, \dots, m \rangle$$

or, more concisely,

$$
\underline{\underline{\underline{\mathcal{O}\_D}}} \approx L\_{DS} \underline{\underline{\mathcal{O}\_S}} \, \, \, \, \, \, \tag{43}
$$

where <sup>ϕ</sup>*<sup>D</sup>* <sup>=</sup> ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ ϕ*D*<sup>1</sup> ϕ*D*<sup>2</sup> : : ϕ*Dn* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ and [*LDS*]*ij* = {*Le*˜} Δ*S j <sup>p</sup>Di* . Hence approximations to the solution within the

domain <sup>ϕ</sup>!*<sup>D</sup>* may be found by the matrix-vector multiplication

$$
\widehat{\underline{\underline{\widehat{\rho\_D}}}} = L\_{DS} \, \widehat{\underline{\underline{\mathcal{G}}}}\,.\tag{44}
$$

#### 3.3.3. The Galerkin Method

Although most of the implementations of the boundary element method are derived through collocation, other techniques can be used, most typically the Galerkin method. The Galerkin method and collocation can both be derived from a more generalised approach that are termed weighted residual methods. For example, for the approximation (40), the residual is the difference between the approximation and the exact solution;

$$R(\underline{\sigma}\_S; p) = \sum\_{j=1}^n \overline{\sigma}\_j \langle L\tilde{\varepsilon} \rangle\_{\Delta \widetilde{S}\_j}(p) - \varphi(p) \ (p \in S).$$

In weighted residual methods, *R* σ*S*; *p* is integrated with *test* basis functions χ*i*(*p*) (*p*·*S*) and the methods arise by setting the result to zero;

$$\int\_{S} \mathcal{R}(\underline{\partial}\_{S}; p\_{i})\_{\chi\_{i}} (p\_{i}) dS = 0$$

at points *pi* on the boundary, for *i* = 1, 2, ... *nS*.

If χ*i*(*p*) = δ *p* − *p<sup>i</sup>* , the Dirac delta function, the collocation method is derived;

$$\int\_{S} \mathcal{R}(\underline{\hat{\mathbf{q}}}\_{S}; p) \delta(p - p\_{i}) \, dS = \, \mathcal{R}(\underline{\hat{\mathbf{q}}}\_{S}; p\_{i}) = 0,$$

which leads to the methods outlined above. In the Galerkin method, the test functions are the same as the original basis functions. For example, for constant elements, the basis and the test functions are defined as

$$\chi\_i(p) = \begin{cases} 1 \text{ or } \overline{e}, \text{ if } p \in \Delta \overline{S}\_i \\\ 0, \text{ otherwise} \end{cases}$$

Substituting this and the definition of the residual equation

$$\int\_{S} R(\underline{\hat{\varphi}}\_{i}; p\_{i}) \chi\_{i}(p) dS = \int\_{\Delta \mathcal{S}\_{i}} \sum\_{j=1}^{n} \partial\_{j} \langle L\bar{x} \rangle\_{\Delta \widetilde{S}\_{j}}(p\_{i}) - \varphi(p\_{i}) \, dS = \begin{array}{c} 0 \ \dots \end{array}$$

The reason for the relative unpopularity of the Galerkin approach in boundary elements is that the matrix is now the result of a double integration, rather than the single integration in the collocation method. However, the matrix is understood to be symmetric.

#### *3.4. Properties and Further Details*

In the boundary element method, the boundary is represented by a mesh of *panels*. The boundary functions are approximated by a linear combination of basis functions on each panel. The boundary element is the combination the panel and the functional representation. In the previous Subsection, constant elements were mentioned as an example. In the finite element method, the degree of the approximation has to be at least half the order of the PDE. Fortunately, this is not the case for the application of the boundary element method, where constant elements are widely used. The panels that make up the boundaries are most simply represented by straight line panels in 2D, triangles in 3D and conic sections for axisymmetric 3D problems.

In this Subsection, important properties and further details of Laplace's equation, the related boundary integral equations and operators are discussed. An overview of methods for carrying out the integrations is included. The issue of non-uniquesness, that is an important feature of the boundary integral equation formulations for the exterior acoustic problem, is first addressed with Laplace's equation and useful outcomes from this are outlined.

#### 3.4.1. Integration

In the previous Subsection, it was shown that up to four operators are involved in the boundary element method for Laplace problems, and these four operators extend to acoustic/Helmholtz problems. For the Laplace problem, it may be possible to develop analytic expressions for the integrals [92], but in general—and particularly for acoustic/Helmholtz problems—numerical integration is necessary. In the main, the integrals are regular, and these are most efficiently evaluated by Gauss-Legendre quadrature [93]. This is straightforward to apply to straight-line panels and in the generator and (typically in composite form) azimuthally. There are also published points and weights for Gaussian quadrature on a triangle [94].

However, the integrals corresponding to the diagonal components of the LSS matrix are weakly singular and the diagonal components of the NSS matrix are hypersingular, that is when the point lies on the panel. Moreover, if the point is close to the panel, for example at the centre of a neighbouring panel then the integrals are said to be nearly singular and may require a more accurate numerical integration rule, or special treatment [95]. Special numerical integration methods can be applied in order to evaluate the singular integrals and expressions for the hypersingular integrals may be found through limiting process. The weakly singular integrals, corresponding to the diagonal components of the LSS matrix have a O(ln r) singularity in 2D and an O 1 r singularity in 3D, where r is the distance from the central collocation point. For the axisymmetric 3D case, the azimuthal integration resolves the O 1 r singularity to an O(ln r) singularity on the generator. For the simple elements discussed, analytic expressions for the diagonal components of the LSS and NSS matrices are available for the 2D and 3D (non-axisymmetic) cases for Laplace's equation. Further work on singular integration can be found in the following references [96–101].

Although for simplicity, the four operators are often lumped together as integral operators, *N* is not an integral operator: *N* is termed a pseudo-differential operator, and it therefore has distinct properties. The full expressions for the straight line and triangular panels are given in the author's book [61] (pages 49–50). For illustration the expressions are given for the straight-line panel of length *h* (2D) and for an equilateral triangular panel (3D) with each side of length *h*:

$$\{N\wp\}\_{\Lambda}\{p;n\_{\mathcal{P}}\} = -\frac{2}{\pi\hbar} \text{ for the straight-line panel of length } h \text{ (2D)}\tag{45}$$

$$(N\wp)\_{\Delta}(\mathfrak{p};\mathfrak{n}\_{\mathfrak{p}}) = -\frac{3}{2\pi\hbar} \text{ for the equilateral triangular panel with side of length } \hbar \text{ (3D)}\tag{46}$$

There is one simple but noteworthy remark to be made about the expressions (45) and (46). Normally, if the domain of integration is reduced in size, the integral is similarly reduced, at least in the limit as the domain size converges to zero. However, with these hypersingular integrals, the opposite is found to be the case! A significant consequence of this property is considered in Section 6.2.

#### 3.4.2. Non-Uniqueness and Its Useful Outcomes

The non-uniqueness of some boundary integral equations in the exterior acoustic problem at a set of frequencies is well-known, and the issue and methods of resolution will be considered in the next Section. Because of the importance of the non-uniqueness of solution within the context of this paper, it is helpful to visit this early, with a simpler equation. Although the term non-uniqueness often prefixes the word problem, it is found that some very useful outcomes arise from this analysis.

The non-uniqueness is found with Laplace's equation itself; if ϕ is a solution of the interior Neumann problem then ϕ + *c* is also solution, where *c* is any constant. The non-uniqueness in the interior Laplace problem with a Neumann boundary condition must be reflected in the boundary integral equations. It therefore follows from Equation (21) that the operator *M* + <sup>1</sup> <sup>2</sup> *I* is degenerate, and similarly for *Mt* + <sup>1</sup> <sup>2</sup> *I* from Equation (23) and *N* from Equations (25) and (28).

As we will see in the next Section on acoustic/Helmholtz problems, the operators from the interior problem equations are shared with the exterior problem, and this is shown for Laplace's equation in Section 3.2. For the exterior problem, the derivative direct formulation (31) is unsuitable for solving exterior Laplace problems as both the *N* and *M<sup>t</sup>* + <sup>1</sup> <sup>2</sup> *I* operators are degenerate. Similarly, the indirect formulations, formed from double layer potentials, (35) and (36) are unsuitable. This correlation between the interior Neumann problem and the exterior derivative (direct) and double-layer potential (indirect) formulations carries through to the Helmholtz equation.

Most simply, if *v*(*p*) = 0 on the boundary then ϕ is any constant (e.g., ϕ = 1) throughout the domain is a solution. Substituting these values into the discrete form of the boundary integral Equation (21), the resulting matrix-vector equation is # *MSS* + <sup>1</sup> 2 *I* \$ 1 ≈ 0, where 0 is a vector of zeros and <sup>1</sup> is a vector of ones; every row of the *MSS* matrix must approximately sum to <sup>−</sup><sup>1</sup> <sup>2</sup> . Similarly, from Equation (25), NSS1 ≈ 0; every row of the NSS matrix must approximately sum to zero. A potentially useful outcome of this is that the hypersingular diagonal components of the NSS matrix can be computed from the others, which are all 'regular' (at least for simple elements). This can be taken a step further and, the panel in question may be linked to a fictitious boundary made up of panels and the value of the hypersingular integral determined by summing the other integrals. The values of the singular and hypersingular integrals in the BEM solution of Laplace's equation may be stored and be used to subtract out the singular and hypersingular components of the same integrals for the Helmholtz equation. The method of inventing a fictitious surface to evaluate the hypersingular integrals is precisely the method used in the author's axisymmetric programs [60], as there were analytic expressions for the hypersingular integrals for the panels used in 2D and general 3D, as discussed in the previous Sub-subection, but no other similar way forward for the axisymmetric panels.

The results in the previous paragraph also provide useful methods of validation. One row of the *MSS* and the *NSS* matrices need to be evaluated that that is when the point *p* is on the boundary. The result of summing the rows of *MSS* can verify that a 'closed' boundary is actually closed, if their values are <sup>−</sup><sup>1</sup> <sup>2</sup> and zero, and, if they are not then it indicated that the boundary may be open or there are errors in the mesh. Finally, substituting these values into Equation (30) gives the following,

$$\langle Me \rangle\_S(p) = \begin{cases} 1 \text{ if } p \in E \\ 0 \text{ if } p \in D \\ \frac{1}{2} \text{ if } p \in S \end{cases},\tag{47}$$

where *e* is the unit function; useful in validating that the solution points are within the domain.

These simple validation methods, arising from potential theory, are applicable to all BEMs and beyond; any computer simulation involving surface meshes may benefit from these simple techniques. The author includes these methods of validation routinely in his boundary element codes.

#### **4. The Standard Boundary Element Method in Acoustics**

In the author's book and software [60,61], three classes of acoustic problem were considered; the determination of the acoustic field within a closed boundary, outside of a closed boundary and the interior modal analysis problem. In this Section, the three methods are revisited and recent developments and applications are included. Recently, the software has been re-written in Python [102].

#### *4.1. The Interior Acoustic Problem*

In this Subsection, the BEM is developed for the solution of acoustic problems in a domain that is interior to a closed boundary [61,103]. The method has been applied to room acoustics [104–108], the interior of a vehicle [109–111], modelling sound in the human lung [112,113] and in biological cells [114].

#### 4.1.1. Integral Equation Formulation

The direct integral equations for the interior acoustic problem, reformulating the Helmholtz Equation (10), but follows the same format as the formulation for Laplace's Equation (21);

$$\langle \!\langle M\_k \!\!\varphi \rangle\_S(\mathfrak{p}) - \langle L\_k \!\!\psi \rangle\_S(\mathfrak{p}) = \begin{cases} -\!\!\varphi(\mathfrak{p}) \text{ if } \mathfrak{p} \in D \\ 0 \text{ if } \mathfrak{p} \in E \\ -\frac{1}{2}\!\!\varphi(\mathfrak{p}) \text{ if } \mathfrak{p} \in S \end{cases}.\tag{48}$$

Equation (48) introduces two Helmholtz integral operators that are analogous to the *L* and *M* operators for Laplace's equation, and are defined similarly;

$$\{L\_k \mu\}\_\Gamma(p) = \int\_\Gamma G\_k(p, q) \,\mu(q) dS\_{\mathbb{T}^\*} \tag{49}$$

$$
\langle M\_k \mu \rangle\_\Gamma(p) = \int\_\Gamma \frac{\partial G\_k(p, q)}{\partial n\_q} \, \mu(q) dS\_{q^\*}, \tag{50}
$$

where the Green's function is defined as follows:

$$G\_k(p,q) = \frac{i}{4} \, H\_0^{(1)}(kr), \text{ for two-dimensional problems, and} \tag{51}$$

$$G\_k(p,q) = \frac{e^{jkr}}{4\pi r'}\text{ for three-dimensional problems.}\tag{52}$$

As with the interior Laplace Equations (22) and (23), the indirect integral equation is derived through presuming that the field is generated by a layer potential, defined on the boundary;

$$
\varphi(\mathfrak{p}) = \langle L\_k \sigma \rangle\_S(\mathfrak{p}) \quad (\mathfrak{p} \in D \cup S), \tag{53}
$$

$$w(p) = \left\{ \left( M\_k^t + \frac{1}{2}I \right) \sigma \right\}\_S \left( p; n\_p \right) (p \in S), \tag{54}$$

where, similarly with Equation (20), the operator *M<sup>t</sup> <sup>k</sup>* is defined as follows:

$$\left\{\mathcal{M}\_k^t \mu\right\}\_\Gamma \left(p; \upsilon\_p\right) = \frac{\partial}{\partial \upsilon\_p} \int\_\Gamma \mathcal{G}\_k(p, q) \left\{\mu(q) dS\_{q\_p}\right\} \tag{55}$$

#### 4.1.2. The Boundary Element Method for the Generalised Boundary Condition

Substituting the expressions (53) and (54) into the boundary condition (7) gives

$$
\alpha(p)\langle L\_k \sigma \rangle\_S(p) + \beta(p)\left\langle \left(M\_k^t + \frac{1}{2}I\right)\sigma \right\rangle\_S(p; n\_p) = f(p). \tag{56}
$$

Following the derivation of the BEM through collocation, this resolves to a linear system of the form

$$[D\_a L\_{SS,k} + D\_\beta \Big(M\_{SS,k}^t + \frac{1}{2}I\Big)]\widehat{\underline{\sigma}}\_S = \underline{f}\_{S'} \tag{57}$$

where *D*<sup>α</sup> and *D*<sup>β</sup> are diagonal matrices, with the values of α and β aligned along the diagonal. Equation (56) is solved in the primary stage of the boundary element method, yielding an approximation <sup>σ</sup> *<sup>S</sup>* to <sup>σ</sup>*S*. In the secondary stage, the discrete equivalent of Equation (53) returns the solution at the interior points:

$$
\overleftarrow{\underline{\underline{\underline{D}}}\_{D}} = L\_{DS,k} \, \overleftarrow{\underline{\underline{\underline{D}}}\_{S}} \,. \tag{58}
$$

For the direct formulation, the discrete equivalent of Equation (48) (*p* ∈ *S*),

$$(M\_{SS,k} + \frac{1}{2}I)\widehat{\underline{\widehat{\varphi\_{S}}}} = L\_{SS,k}\widehat{\underline{\widehat{\varphi\_{S'}}}}\tag{59}$$

which is solved with the discrete equivalent of the boundary condition (7),

$$D\_{\alpha}\widehat{\underline{\boldsymbol{\sigma}}}\_{\mathcal{S}} + D\_{\beta}\widehat{\underline{\boldsymbol{\sigma}}}\_{\mathcal{S}} = \underline{\boldsymbol{f}}\_{\mathcal{S}}\,. \tag{60}$$

Comparing the linear systems for the indirect method (56) with that of the direct method (59) and (60) illustrates an apparent significant advantage in the indirect approach, the number of unknowns in the system corresponding to the direct method is twice that of the indirect method. However, through exchanging columns, the system for the direct method can be reduced to an *n* × *n* system, and this method has been automated [115].

There have been several papers discussing the accuracy of the interior acoustic boundary element method near to the resonance frequencies. It is considered that real eigenvaulues are shifted into the complex plane following discretisation. This phenomenon is termed numerical damping [116–119].

#### 4.1.3. Equations of the First and Second Kind

Integral equations with a fixed region of integration are termed Fredholm integral equations. Fredholm integral equations are categorised as first kind or second kind. Equation (53) is a typical first kind equation, in which we are solving over integral operator(s) alone, in this case the *Lk* operator, to find σ from ϕ. With second kind equations, we are solving not just over integral operators, but also the identity (or diagonal) operator. For example, Equation (54) is a second kind equation, solving over the *M<sup>t</sup> <sup>k</sup>* <sup>+</sup> <sup>1</sup> <sup>2</sup> *I* operator, in order to find σ from *v*.

In general, first kind equations have poor numerical properties and are avoided. Although pure first kind equations only occur in particular Dirichlet cases, they can be avoided, through using the derivatives of the integral equation formulations (direct) or double layer potentials (indirect). However, from experience, first kind equations with singular kernels, as we find with the *L* or *Lk* operator, do not have poor convergence properties. It follows, therefore, that the formulations and methods developed thus far in this Subsection are generally suitable in solving the interior acoustic problem.

4.1.4. Derivative and Double-Layer Potential Integral Equations for the Interior Helmholtz Problem

As with the equations listed in Section 4.1.1, the derivative equations are unnecessary in solving the interior Helmholtz problem. However, the derivative equations provide an alternative formulation and help with the general understanding. The form of the equations is analogous with the equations for the interior Laplace equation, developed in Section 3.2.

For the direct method, the derivative boundary integral equation is analogous to Equation (25):

$$\{\mathcal{N}\_k \boldsymbol{\uprho}\}\_S \{\boldsymbol{p}; \boldsymbol{n}\_p\} - \left\{\mathcal{M}\_k^t \boldsymbol{\uprho}\right\}\_S \{\boldsymbol{p}; \boldsymbol{n}\_p\} = -\frac{1}{2} \boldsymbol{\uprho}(\boldsymbol{p}) \ (\boldsymbol{p} \in \mathcal{S}).\tag{61}$$

where

$$
\langle N\_k \mu \rangle\_\Gamma \Big( p; v\_p \big) = \frac{\partial}{\partial v\_p} \int\_\Gamma \frac{\partial G\_k(p, q)}{\partial n\_q} \, \mu(q) dS\_q. \tag{62}
$$

For the indirect method, the boundary integral equations follow the form of Equations (28) and (29):

$$\left[\!\!\!/p(p) = \!\!/(M\_k - \frac{1}{2}I)\!\!/\!\!/ \!\!/ \!\!/ \!\!/ \!\!/ \!\!/ \!\!/ \!\!/ p\right] \left(p \in S\right), \tag{63}$$

$$
\psi(p) = \langle \mathcal{N}\_k \check{\mathsf{\mathsf{\zeta}}} \rangle\_{\mathcal{S}}(p; u\_p) \ (p \in \mathcal{S}).\tag{64}
$$

#### *4.2. Interior Acoustic Modal Analysis: The Helmholtz Eigenvalue Problem*

An enclosed acoustic domain has resonance frequencies and associated mode shapes. Mathematically, these are the solutions *k*∗ (the eigenvalues—that relate to the resonance frequencies) and ϕ∗ (the eigenfunctions—that relate to the mode shapes) of the Helmholtz equation with the homogeneous form of the boundary condition (7) (i.e., with *f*(*p*) = 0). This problem is more typically solved by the finite element method and a finite element model of an acoustic or structural problem is developed in Section 5.3.1. In this Subsection the modal analysis of an enclosed fluid via the boundary element method is outlined.

#### 4.2.1. The Generalised Non-Linear Eigenvalue Problem from the Boundary Element Method

In the indirect boundary element method, the eigenvalues *k*∗ and the eigenfunctions σ∗ are found through solving

$$
\alpha(p)\langle l\_k \sigma \rangle\_S(p) + \beta(p) \left\{ \left( M\_k^t + \frac{1}{2}I \right) \sigma \right\}\_S(p; n\_p) = 0,\tag{65}
$$

that follows from Equation (56), with the true eigenfunctions ϕ∗ can then be found with Equation (53). Through applying collocation, this is equivalent to solving the non-linear algebraic eigenvalue problem

$$[D\_a L\_{\mathbb{S}\mathbb{S}}{}^{\mathbb{A}} + D\_{\beta}(\mathcal{M}\_{\mathbb{S}\mathbb{S},\mathbb{A}^\*}^t + \frac{1}{2}I)]\widehat{\underline{\sigma}\_{\mathbb{S}}} = \underline{0},\tag{66}$$

which follows from Equation (57).

The eigen-solution of Equation (66) returns the approximations ˆ *k*∗ to the eigenvalues and the approximation <sup>σ</sup> *<sup>S</sup>* to the layer potential eigenfunction. The approximations to the physical eigenfunctions at the chosen domain points can then be found with Equation (58). The method described can also be applied with the direct integral equation formulation but requiring the row-exchanging method mentioned in Section 4.1.2.

Although in this work, the focus is on the generalised boundary condition, most of the examples in the literature consider the Dirichlet and Neumann eigenvalues, and it is revealing to focus on those. The Dirichlet/Neumann interior eigenvalues and eigenfunctions are the solutions of the equations of this Subsection with ϕ(*p*) = 0 / *v*(*p*) = 0 (*p* ∈ *S*). Hence the Dirichlet eigenvalues are the eigenvalues of the *Lk* operator from Equations (48) and (53), of the *Mk* <sup>−</sup> <sup>1</sup> <sup>2</sup> *I* operator from Equation (63) and the *M<sup>t</sup> <sup>k</sup>* <sup>−</sup> <sup>1</sup> <sup>2</sup> *I* operator from Equation (61). Similarly, the Neumann eigenvalues are the eigenvalues of the *Mk* + <sup>1</sup> <sup>2</sup> *<sup>I</sup>* operator from Equation (48), *Mt <sup>k</sup>* <sup>+</sup> <sup>1</sup> <sup>2</sup> *I* from Equation (54) and the *Nk* operator from Equations (61) and (64).

#### 4.2.2. Solving the Non-linear Eigenvalue Problem

The standard algebraic eigenvalue problem has the form

$$A\underline{\mathbf{x}} = \lambda \underline{\mathbf{x}}\_{\prime} \tag{67}$$

where *A* is a square matrix. Because of the analogy between the Helmholtz Equation (3) and the standard algebraic eigenvalue problem (67), the Helmholtz eigenvalues are sometimes termed the eigenvalues of the Laplacian [120–124].

The generalised algebraic eigenvalue problem has the form

$$A\underline{\mathbf{x}} = \lambda B\underline{\mathbf{x}}\tag{68}$$

where *A* and *B* are square matrices. Standard methods exist for solving these problems, termed the *QR* and *QZ* algorithm. As we have noted, the application of the boundary element method in Equation (66) results in a non-linear eigenvalue problem, and this has the form

$$A\_k \underline{\mathfrak{x}} = \underline{\mathfrak{Q}}.\tag{69}$$

Although this is a significantly complicated problem, at least—in the case of the boundary element matrices—the individual components of the components matrix are continuous. Earlier methods tended to find the eigenvalues by finding the zeros of |*Ak*|. Further developments and applications can be found in the following research papers [61,122–148].

#### *4.3. The Exterior Acoustic Problem*

The boundary element solution of the exterior acoustic problem is the most popular area of research in the context of this paper. A method that can solve over a theoretically infinite domain from data on a surface mesh has significant value in the context of acoustics. However, it was found more than half a century ago that the boundary integral equations, derived in the same way as the ones for the interior problem, resulted in unreliable boundary element methods. Much has been achieved on the road to developing a more successful outcome. However, to obtain a reliable exterior acoustic boundary element method remains a tantalising goal.

Two reasonably successful pathways for developing a boundary element solution of the exterior acoustic problem were developed about 50 years ago, and these remain today. In this paper, the first is termed the Schenck method [149] and the second is termed the combined boundary integral equation method (CBIEM). There are several references on a general review and evaluation of the methods [55,150–157] and of software implementations [158–160]. The methods have been used in a range of applications: loudspeakers [38,50,51,161–178], transducers [52,179–183], hearing aid [184], diffusers [185,186], sound within or around the human body [159,187–190], scattering by blood cells [191], engine/machine noise [53,56,192–208], aircraft noise [209–211], rail noise [212], non-destructive defect evaluation [142,213–215], noise barriers [74,75,216,217], environmental noise [218–224], underwater acoustics [225,226], detecting fish in the ocean [227,228] and perforated panels [229].

#### 4.3.1. The Integral Equation Formulations of the Exterior Helmholtz Equation and Their Properties

The direct boundary element reformulation of the Helmholtz equation follows the format for Laplace's Equation (30) and is as follows:

$$\{M\_k \wp\}\_S(p) - \langle L\_k \upsilon \rangle\_S(p) = \begin{cases} 0 \text{ if } p \in D \\ \wp(p) \text{ if } p \in E \\ \frac{1}{2} \wp(p) \text{ if } p \in S \end{cases}.\tag{70}$$

The indirect formulation is

$$
\varphi(\mathfrak{p}) = \langle L\_k \sigma \rangle\_S(\mathfrak{p}) \quad (\mathfrak{p} \in E \cup \mathbb{S}), \tag{71}
$$

$$w(p) = \left\{M\_k^t - \frac{1}{2}I\right\}\_{\Gamma} \sigma(p; \mathfrak{n}\_{\mathbb{P}}) \quad (p \in S). \tag{72}$$

However, the *Lk*, *Mk* <sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>I</sup>* and *Mt <sup>k</sup>* <sup>−</sup> <sup>1</sup> <sup>2</sup> *I* operators are degenerate at the eigenvalues of the interior Dirichlet problem (see Section 4.2.1). The issues in using the boundary integral equations in (70)–(72) as general purpose exterior acoustic/Helmholtz equation solvers has been known for over 50 years. The wavenumbers or frequencies in which these boundary integral equations are unsuitable are often termed the characteristic or irregular wavenumbers or frequencies in the literature. These characteristic wavenumbers are physical in the interior problem, but they are unphysical in the exterior problem in which they are not a property of the Helmholtz equation model but are manifest in the boundary integral equations. In Section 3.4.2 the effects of the non-uniqueness of the solution of the interior Laplace problem with a Neumann boundary condition were discussed. The boundary integral Equations (70)–(72) have a similar property at the characteristic wavenumbers, and this is also often termed the non-uniqueness problem.

Although it is 'unlikely' in practice that the value of wavenumber *k* is equal to a characteristic wavenumber *k*∗ , it is shown in Amini and Kirkup [230] that the numerical error as a consequence of being 'close' to a characteristic wavenumber is *O*( <sup>1</sup> |*k*−*k*∗| ). Hence, one technique of resolving the problem is to use finer meshes in the neighbourhood of the characteristic wavenumber in order to offset this error [231]. However, this strategy is likely to be prohibitive, requiring the overhead of creating a range of meshes and increased computational cost. The values of the characteristic wavenumbers are also generally unknown, although they can be found (as discussed in Section 4.2), but at a substantial computational cost. It is also found that the characteristic wavenumbers cluster more and more as the frequency increases. In conclusion, therefore, the boundary integral equations are—in practice—only useful for frequencies reasonably below a conservatively estimated first characteristic wavenumber, severely restricting the methods to the low frequency range in practice.

#### 4.3.2. The Derivative and Double Layer Potential Integral Equations

For the direct method, the derivative boundary integral equation is analogous to Equation (31):

$$
\langle N\_k \boldsymbol{\uprho} \rangle\_S \begin{pmatrix} \boldsymbol{\uprho} \; \boldsymbol{\uprho} \end{pmatrix}\_S = \left\{ (M\_k^t + \frac{1}{2}I)\boldsymbol{\uprho} \right\}\_S \begin{pmatrix} \boldsymbol{\uprho} \; \boldsymbol{\uprho} \; \boldsymbol{\uprho} \end{pmatrix}.\tag{73}
$$

For the indirect method, the boundary integral equations follow the form of Equations (35) and (36):

$$\left\|\varphi(\mathbf{p}) = \left\|(M\_k + \frac{1}{2}I)\zeta\right\|\_{\mathcal{S}}(\mathbf{p}),\tag{74}$$

$$\upsilon(p) = \langle \mathsf{N}\_k \check{\mathsf{L}} \rangle\_S(p; \mathfrak{n}\_p)\_\prime \tag{75}$$

Again, some of the operators are shared with the interior formulation. The operators *Mk* + <sup>1</sup> <sup>2</sup> , *<sup>M</sup><sup>t</sup> <sup>k</sup>* <sup>+</sup> <sup>1</sup> 2 *I* and *Nk* are degenerate at the eigenvalues of the interior Neumann problem, and hence these equations are also unsuitable as the basis of methods of solution at those frequencies. These equations mirror the properties of the elementary equations discussed earlier and their straightforward solution does not result in an acceptable boundary element method.

The characteristic wavenumber for the derivative and double-layer potential Equations (73)–(75) (the interior Neumann eigenvalues) are generally different from those of the elementary Equations (70)–(72) (the interior Dirichlet eigenvalues), and at least therefore they provide alternative methods. However, a more useful path involves combining the elementary Equations (70)–(72) with these equations and these methods are considered in Section 4.3.4.

#### 4.3.3. The Schenck Method

The Schenck method [149] is often termed the CHIEF method in the literature and it is a development of the standard direct method, based on Equation (70). Given the potential non-uniqueness, discussed in Section 4.3.1, the system of equations that form the discrete equivalent of the boundary integral equations are regarded as potentially underdetermined and they are augmented with equations related to a set of points in the interior D:

$$\left[\begin{array}{c} M\_{\text{SS},k} - \frac{1}{2}I \\ M\_{\text{DS},k} \end{array}\right] \underline{\oplus}\_{\text{S}} = \left[\begin{array}{c} L\_{\text{SS},k} \\ L\_{\text{DS},k} \end{array}\right] \underline{\oplus}\_{\text{S}} \dots$$

By adding more equations, the expectation is to eliminate the non-uniqueness and determine a solution. The equations can be solved by the least-squares method for Dirichlet and Neumann problems, the column exchanging algorithm [115] could be used in the case of the general boundary condition (7). There are several reported implementations and applications of the Schenck method [216,232–234].

Obviously, there are immediate questions about the number and position of the interior CHIEF points. Juhl [235] develops an iterative method for selecting points and halting when the results are unchanged. Equation (47) can verify that CHIEF points are interior, so this could be usefully included in the method. There has been a significant number of implementations and testing of the Schenck method. In general, more and more points are required to offset the non-uniqueness, as *k* increases. This increases the computational overhead with respect to the wavenumber, additional to the potential need for more elements at higher wavenumbers. Several improvements in the method have been put forward and these are summarised in Marburg and Wu [236]. There have been several evaluations of the CHIEF and combined methods and their variants [157,237,238].

#### 4.3.4. The Combined Integral Equation Method

In this method, a boundary integral equation is formed through a linear combination of the initial boundary integral equation and its corresponding derivative equation (for the direct method and double-layer potential for the indirect metod). This concept was initially derived for the indirect method and is attributed to Brakhage and Werner [239], Leis [240], Panich [241] and Kussmaul [242]. The corresponding direct integral equations are attributed to Burton and Miller [243].

The Burton and Miller method is based on a boundary integral equation that is a linear combination of the initial one (70) with the derivative (73),

$$\left\{ (M\_k - \frac{1}{2}I + \mu N\_k)\varphi \right\}\_S \left( \mathfrak{p}; \mathfrak{n}\_{\mathbb{P}} \right) = \left\{ (L\_k + \mu (M\_k^t + \frac{1}{2}I)\upsilon \right\}\_S \left( \mathfrak{p}; \mathfrak{n}\_{\mathbb{P}} \right) \quad (\mathfrak{p} \in S), \tag{76}$$

where μ is a complex number, a *coupling* or weighting parameter. Similarly, for the indirect method, the equation is based on writing ϕ as a linear combination of a single and double layer potential

$$
\langle \boldsymbol{\uprho} \left( \boldsymbol{\uprho} \right) \rangle = \langle (L\_k + \mu M\_k) \boldsymbol{\uprho} \rangle (\boldsymbol{\uprho}) \ (\boldsymbol{\uprho} \in E). \tag{77}
$$

This returns the following boundary integral equations

$$\left\|\varphi(\mathbf{p}) = \left\|(L\_k + \mu(M\_k + \frac{1}{2}I))\eta\right\|\_{\mathcal{S}}(\mathbf{p}) \quad (\mathbf{p} \in S), \tag{78}$$

$$w(p) = \left| (M\_k^t - \frac{1}{2}I + \mu N\_k)\eta \right|\_S (p; \mathfrak{u}\_p) \quad (p \in S), \tag{79}$$

The issue of the determination of the values for the coupling parameter will be revisited in Section 6.2.

#### 4.3.5. Scattering

The boundary integral equations for the exterior acoustic problem are readily applicable to radiation problems. The equations can also be used for the scattering problem, but with extra terms involved that model the incident field. These same techniques can be applied in the interior problem. For example, for the indirect method, the exterior acoustic field (77) is presumed to be made up of the layer potentials, superposed with the (known) incident field:

$$\varphi(\mathcal{p}) = \varphi\_{\text{inc}}(\mathcal{p}) + \{ (L\_k + \mu M\_k)\eta \}(\mathcal{p}) \ (\mathcal{p} \in E)^2$$

This similarly adjusts Equations (78) and (79):

$$\begin{aligned} \varphi(p) &= \varphi\_{\text{inc}}(p) + \left\|(L\_k + \mu(M\_k - \frac{1}{2}I))\eta\right\|\_S(p) \quad (p \in S), \\\\ \upsilon(p) &= \upsilon\_{\text{inc}}(p) + \left\|(M\_k^t + \frac{1}{2}I + \mu N\_k)\eta\right\|\_S(p; \mu\_p) \quad (p \in S). \end{aligned}$$

#### **5. Extending the Boundary Element Method in Acoustics**

In this Section, extensions in the standard boundary element methods of the previous section are considered. These include the Rayleigh integral method for computing the acoustic field surrounding a vibrating plate set in an infinite reflecting baffle, as a case of the more general half-space methods. This Section also includes shell elements in which a revision in the boundary integral equation for the exterior problem returns a model for the acoustic field surrounding a thin screen. Through principles of continuity and superposition, hybrids of these models and the standard models also significantly extend the range of acoustic problems that come under the boundary element fold. Vibro-acoustic, aero-acoustic and inverse acoustic problems are also considered in this Section.

#### *5.1. Half-Space Methods*

In Section 3 it was stated that the boundary element solution of Laplace's equation was a useful entry to the BEM in acoustics. The Rayleigh integral method (RIM) is also a good a starting point, in that it required only one of the four Helmholtz operators, and, for Neumann problems, it is an integral, rather than an integral equation. In this Subsection, the Rayleigh integral method is defined and further developments are reviewed.

The Rayleigh integral method can be viewed as a solution to a half-space problem. If further boundaries are placed in the half-space, then the integral equation formulation, with a simple modification of the Green's function, forms the model with ϕ = 0 or ∂ϕ <sup>∂</sup>*<sup>n</sup>* = 0 on the plane. Further development of the Green's function have been researched in order that an impedance condition is modelled on the plane, a useful model in outdoor sound propagation.

On the other hand, if there is a cavity in the plane then the interior boundary element method can model the acoustic field within the cavity and this is coupled to the Rayleigh integral. The model is based on based on the interior formulation to model the cavity, applying a false flat boundary across the opening and coupling the interior formulation with that of the half-space. This method is a hybrid of the boundary element method and the Rayleigh integral method and is termed *BERIM*.

#### 5.1.1. The Rayleigh Integral Method

In the operator notation of this paper, the Rayleigh integral is as follows:

$$\varphi(\mathfrak{p}) = -2 \langle L\_k \upsilon \rangle\_{\varPi}(\mathfrak{p}) \quad \left(\mathfrak{p} \in E^+ \cup \varPi\right)\_\prime \tag{80}$$

where Π is the flat plate, radiating into the half-space *E*+. The solution of the Neumann problem, finding ϕ from *v*, is simply the evaluation of an integral. For the general boundary condition of the

form (7), the technique is again to solve the boundary integral Equation ((80) with *p* ∈ Π) to obtain *v* and then evaluating the integral to obtain ϕ at any point in *E*+. In the Rayleigh integral method [71], the plate is divided into elements, as discussed and applying collocation or, the equivalent for the integral, product integration. Substituting (80) into the boundary condition (7) gives

$$-2\alpha(\mathfrak{p})\langle L\_k v \rangle\_{\varPi}(\mathfrak{p}) + \beta(\mathfrak{p})v(\mathfrak{p}) = f(\mathfrak{p}).\tag{81}$$

Using the discrete notation of this paper, the equation following the application of collocation is as follows:

$$\left(-2D\_{\bowtie}L\_{\text{III}\upharpoonright} + D\_{\beta}\right)\underline{\upleftarrow}\_{\text{II}} = \underline{f}\_{\text{II}}.\tag{82}$$

There have been several reported implementations and developments of the Rayleigh integral method [244–254], including vibro-acoustics [197,255,256]. Applications of the method include sandwich panels [257–263], engine or machine noise [203], electrostatic speaker [264] and transducers [52,265,266].

#### 5.1.2. Developments on Half-Space Problems

Integral equations for scattering or radiating boundaries above an infinite plane can be developed through altering the Green's function [247] in order that it also satisfies the condition on the plane. For the perfectly reflecting plane the Green's function must satisfy the Neumann condition on the plane; <sup>∂</sup>*G*<sup>∗</sup> <sup>∂</sup>*<sup>n</sup>* = 0 and hence *G*∗(*p*, *q*) = *G*(*p*, *q*) + *G*(*p*, *q*∗), where *q*<sup>∗</sup> is the point that corresponds to *q*, when reflected through the plane. Similarly, if there is a homogeneous Dirichlet condition on the plane then the revised Green's function is *G*∗(*p*, *q*) = *G*(*p*, *q*) − *G*(*p*, *q*∗). More generally, the modified Green's function is *G*∗(*p*, *q*) = *G*(*p*, *q*) + *R G*(*p*, *q*∗), with *R* representing the reflection of the plane (−1 ≤ *R* ≤ 1). An impedance boundary condition is generally required in modelling environmental noise problems [220,267–270].

However, the methods in the previous paragraph are applicable when the acoustic scatterers or radiators are on or above the plane. Another set of methods apply if there is a cavity in the plane. Early examples of this type of problem have arisen in harbour modelling, in which the Helmholtz equation has been used to model the waves in a harbour (the cavity), which is open to the sea bounded by the coastline (the plane) [271,272]. The Boundary Element—Rayleigh Integral Method (BERIM) [169], is applicable to open cavity problems in acoustics. The interior boundary element method (Section 4.1) models the sealed interior and the Rayleigh integral method models the field exterior. The advantage in this model over the exterior model could be substantial; the mesh covers the interior of the cavity and the opening only, rather than the inside and outside. There is an important issue with the model, as it presumes that the cavity opens out on to an infinite reflecting baffle. However, this could be a small price to pay, and some problems—such as the loudspeaker problems considered by the author [169]—the mouth opens onto a front face of the cabinet. Motivated again by environmental noise problems, several methods have been developed for a cavity opening on to an impedance plane [273].

#### *5.2. Shell Elements*

The derivation of integral equations and methods for modelling thin shells (that is an open boundary modelling a discontinuity in the field) in the boundary element context takes us back to the works such as Krishnasamy [274], Gray [275], Terai [276] and Martin [277]. In these, and various other references, the shell is also termed a crack. For the Helmholtz equation, the derivation of the integral equations is set out in Warham [73] and Wu [72].

Shell problems may be solved using the standard boundary element methods already described in this paper. One method is to mesh the shell as a closed boundary, with a finite thickness. However, if the thickness of the shell is a fraction of the element size, then the boundary integral equations representing collocation points on either side of the shell are similar, and the equation approaches

degeneracy. Alternatively, many elements may be required, to ensure that the elements on the shell are not disproportionate in comparison with those along the edges.

A more productive method of using existing methods is to artificially extend the shell to form a complete boundary, with the interior and exterior BEM applied to the inner and outer domains, and continuity applied over the artificial boundary. However, this approach, in many cases, would be prohibitive, as the number of elements would be substantially greater. However, analytic test problems for shell problems are difficult to develop and using the more-established interior and exterior BEM in this way can provide comparative solutions. As with the Rayleigh integral method of the previous Subsection, the same Helmholtz operators re-occur with the shell model. Hence, the inclusion of shells in boundary element software significantly extends the functionality of the library, at little extra cost.

#### 5.2.1. Derivation of the BEM for Shells

In this model, ϕ is the solution of the exterior Helmholtz Equation (3) in the exterior domain, surrounding an open boundary Ω. The boundary Ω is presumed to be an infinitesimally thin discontinuity and so, at the points on the boundary, ϕ and its normal derivative *v* take two values, one at each side of the discontinuity. The two sides of the shell are labelled "+" and "-": it doesn't matter which way round this is, as long as it remains consistent. Hence, on the shell, four functions are modelled, ϕ+(*p*), ϕ−(*p*), *v*+(*p*) and *v*−(*p*) (*p* ∈ **Ω**), where the normal to the boundary, which orientates *v*<sup>+</sup> and *v*−, is taken to point outward from the '+' side of the shell. However, rather than working with these functions, it is more straightforward to work with the difference and average functions;

$$\begin{aligned} \delta(p) &= \wp\_+(p) - \wp\_-(p), \\ \phi(p) &= \frac{1}{2}(\wp\_+(p) + \wp\_-(p)), \\ \nu(p) &= \upsilon\_+(p) - \upsilon\_-(p), \\ V(p) &= \frac{1}{2}(\upsilon\_+(p) + \upsilon\_-(p)), \end{aligned}$$

for (*p* ∈ Ω) and where, for simplicity, it is presumed that the boundary is smooth.

The boundary condition is defined in a similar way as in Equation (7), but as there are double the number of unknown function, two boundary conditions are required (*p* ∈ Ω):

$$
\alpha(p)\delta(p) + \beta(p)\upsilon(p) = f(p), \tag{83}
$$

$$A(\mathfrak{p})\otimes(\mathfrak{p}) + B(\mathfrak{p})V(\mathfrak{p}) = F(\mathfrak{p}),\tag{84}$$

The integral equations that govern the field around the shell discontinuities derived from the exterior direct formulation (70) [73], by taking the limit as the boundary becomes thinner, and they are as follows:

$$\mathcal{O}(p) = \langle M\_k \delta \rangle\_\Gamma(p) - \langle L\_k \nu \rangle\_\Gamma(p) \ (p \epsilon \Gamma), \tag{85}$$

$$\mathcal{V}(p) = \langle \mathbb{N}\_k \delta \rangle\_{\Gamma}(p; \mathfrak{n}\_p) - \left\langle \mathbb{M}\_k^t \nu \right\rangle\_{\Gamma}(p; \mathfrak{n}\_p) \ (p \in \Gamma), \tag{86}$$

$$
\varphi(\mathfrak{p}) = \langle \mathbb{M}\_k \delta \rangle\_\Gamma(\mathfrak{p}) - \langle \mathbb{L}\_k \boldsymbol{\nu} \rangle\_\Gamma(\mathfrak{p}) \text{ ( $\mathfrak{p} \in \mathcal{E}$ )}.\tag{87}
$$

There are few tests of methods based on these equations in the literature. The only known issue is at the edge, as it is likely that the solution is singular there. Therefore, mesh grading, using smaller and smaller elements close to the edge may improve efficiency.

#### 5.2.2. Mixing Opem with Closed Boundaries

Again, using the superposition principle, discussed in Section 4.3.5, shell boundaries can be mixed with the traditional boundaries. For example, for the interior problem made up of a domain *D* with boundary *S* and with shell discontinuities Γ within the domain, the superposition of the direct Equations (48) and (61) with Equations (85)–(87) returns the following equations:

<sup>Φ</sup>(*p*) = {*Lkv*}*S*(*p*) <sup>−</sup> *Mk*ϕ *<sup>S</sup>*(*p*) + {*Mk*δ}Γ(*p*) − {*Lk*ν}Γ(*p*) (*p*Γ), *<sup>V</sup>*(*p*) <sup>=</sup> *Mt kv S* (*p*) <sup>−</sup> *Nk*ϕ *<sup>S</sup>*(*p*) + {*Nk*δ}<sup>Γ</sup> *p*; *np* − *Mt k*ν Γ *p*; *np* (*p*Γ), <sup>ϕ</sup>(*p*) = {*Lkv*}*S*(*p*) <sup>−</sup> *Mk*ϕ *<sup>S</sup>*(*p*) + {*Mk*δ}Γ(*p*) − {*Lk*ν}Γ(*p*) (*pD*), 1 2 <sup>ϕ</sup>(*p*) = {*Lkv*}*S*(*p*) <sup>−</sup> *Mk*ϕ *<sup>S</sup>*(*p*) + {*Mk*δ}Γ(*p*) − {*Lk*ν}Γ(*p*) (*pS*).

Methods based on this analysis and equations were developed and demonstrated by the author for the interior Laplace equation [85], the exterior acoustic/Helmholtz problem [278], and for the interior acoustic/Helmholtz problem [135,279].

#### *5.3. Vibro-Acoustics*

Problems that involve structural vibration, as well as an acoustic response, are termed vibro-acoustics. In this Subsection, the modelling by a domain method, such as the finite element method is outlined. The FEM can be applied to the structure and/or the acoustic/fluid domain, but, in the context of this Subsection, the FEM is used as the structural model. When the structure and the fluid significantly influence each other's dynamic response then they are modelled as coupled fluid-structure interaction.

#### 5.3.1. Discrete Structural or Acoustic (Finite Element) Model

The properties of the interior acoustic problem parallel the expected response an excited elastic structure, presuming no damping, or, more simply, simple harmonic motion. In this discussion, let us consider this further in order to bring context. With a mass M and a stiffness K, the equation of (unforced) simple harmonic motion is

$$M\,\bar{q} + K\,q = \underline{0}\,\,\prime$$

where *<sup>q</sup>* is the displacement and .. *q* is the acceleration. The same equation results from a system of masses or from the finite element method solution with *M* and *K* termed the mass and stiffness matrices and *<sup>q</sup>* and .. *q* are vectors of displacement and acceleration of the individual masses or nodes in the FEM. The phasor solution *q* = *Qej*ω*<sup>t</sup>* returns the following equation,

$$-\omega^2 M \underline{\underline{Q}} + K \underline{\underline{Q}} = 0 \ , \tag{88}$$

Hence, applying the finite element method to the structural or interior acoustic modal analysis problem returns a generalised algebraic eigenvalue problem (of the form of Equation (68)). The matrices are sparse and hence are amenable to more efficient methods of solution. Although the matrices are larger and the domain needs to be meshed, the BEM with its nonlinear eigenvalue problem struggles to compete with the FEM in acoustic modal analysis. Let ω<sup>2</sup> *<sup>j</sup>* and *Qj* for j = 1, 2, ... be the eigenvalues (natural frequencies) and corresponding eigenvectors (mode shapes) of Equation (88). It follows that *<sup>M</sup>*−<sup>1</sup>*KQj* = ω<sup>2</sup> *j Qj* .

Let us also now generalise the model in order to include and external excitation force *g*:

$$
\mathcal{M}\ddot{\mathfrak{q}} + \mathcal{K}\mathfrak{q} = \mathfrak{g}.\tag{89}
$$

Following the phasor substitution *g* = *Gej*ω*<sup>t</sup>* Equation (89) is modified as follows:

$$-\alpha^2 \underline{\mathcal{M}} \underline{\mathcal{Q}} + \underline{\mathcal{K}} \underline{\mathcal{Q}} = \underline{\mathcal{G}} \, , \text{ or}$$

$$-\alpha^2 \underline{Q} + M^{-1} K \underline{Q} = M^{-1} \underline{G} \dots$$

For the homogeneous equations (*G* = 0), the above may be cast as a generalised or standard algebraic eigenvalue problems, as discussed in Section 4.2.2. Let us write the response and the excitation in terms of the modal basis

$$
\underline{Q} = \sum\_{j} \gamma\_{j} Q\_{j} \text{ and } \mathcal{M}^{-1} \underline{\mathcal{G}} = \sum\_{j} a\_{j} Q\_{j}. \tag{90}
$$

Considering each eigen-solution in turn relates the coefficients, so that

$$\gamma\_j = \frac{a\_j}{\omega\_j^2 - \omega^2}.$$

In practice, a structure or an enclosed fluid experience damping, the simplest and usual model is to relate damping to the velocity:

$$M\,\,\underline{\ddot{q}} + \mathcal{C}\,\,\underline{\dot{q}} + \,\,\mathcal{K}\,\,\underline{q} = \underline{\underline{g}}$$

where *C* is termed the damping matrix. Following the phasor substitutions, Equation (89) is modified as follows:

$$-\omega^2 M \underline{Q} + j\omega \underline{C} \underline{Q} + \underline{K} \underline{Q} = \underline{G}.\tag{91}$$

The response of the system to an applied boundary condition across a frequency sweep is to have a smoothed peak at the resonance frequency and a more gradual phase change. In general, the response of the system can be modelled as a sum of weighted modes (90) with the coefficients that are relatable as follows:

$$
\gamma\_j = d\_j(\omega) a\_i. \tag{92}
$$

#### 5.3.2. Coupled Fluid-Structure Interaction

The finite element model in the previous Sub-subsection is directly applicable to a structure when there is no significant coupling between the structure and a fluid. Similarly, the acoustic analysis methods of Sections 4 and 5.1 and Section 5.2 are directly applicable when there is no significant coupling between the fluid and a structure. However, for many dynamical systems, it is appropriate to couple the structural model, of the FEM form outlined in the previous Sub-subsection, with the acoustic model of the fluid with which it is in contact. The finite element method can be applied in both domains, with appropriate properties. In the context of this work, the boundary element method provides the computational acoustic model. The models are coupled together, through continuity of the particle velocity at the interface, and the resultant forcing on the structure is affected by the sound pressure. The discrete coupling is applied at the elements that coincide with the boundary.

The fluid-structure interaction model with the boundary element method modelling the acoustic domain has been developed and applied over several decades. Expansions on the general method can be found in the following works and the references therein [183,267,280–283]. Applications include the interaction of plates with fluids [23,245,284], sandwich panels/lightweight structures [257,261–263,285], sound insulation [283], screens [286], the passenger compartment of an automobile [111], hydrophones [284,287] and in seismo-acoustics [288].

In Section 4.2, the acoustic modal analysis of an enclosed fluid was discussed. Similarly, in the previous Section, structural modal analysis via the finite element method was outlined. Of course, when coupling occurs, the eigenvalue analysis need to be applied to the coupled system. For example, a structure typically exhibits different resonant frequencies in vacuo than it does when immersed in a fluid. In the literature, these are often termed the *dry* and *wet* modes (at least for structures placed outside of and in water). Modal analysis via the finite element method returns a standard and sparse eigenvalue problem (88), modal analysis by the boundary element method yields a non-linear eigenvalue problem (69) and hence the hybrid coupled FEM-BEM system is also non-linear. There are several reported implementations of coupled fluid-structure modal analysis using the boundary element method [19,255–257,265,284,287,289–301]. The coupled matrices are generally much larger than for the boundary element method alone and this can be resolved by determining the response in the dry structural modal basis and coupling that with the boundary element model [299,301].

#### *5.4. Aero-Acoustics*

While the boundary element method has been applied to coupled fluid-structure interaction problems for almost as long as the method has been around, the BEM in aero-acoustics is a much more recent development. The standard exterior BEM in a vibro-acoustic settting has been applied in aircraft noise [210,211]. The wave Equation (1) does not include convective flow and is only an adequate model in aero-acoustics in the special case of insignificant flow and the formulations and methods outlined in this work are no longer directly applicable. The Navier-Stokes equations model fluid flow and their solution by numerical methods is termed *computational fluid dynamics* (CFD). There are reports of applying the standard acoustic BEM and CFD to aircraft [47,48,302–304] and, similarly, to underwater vehicles [305,306].

To model the noise from aircraft, the domains are typically large and significant resolution is required to capture the higher frequencies, and hence domain methods come with a high computational cost. Computational aero-acoustics [11] has arisen for developing and applying numerical methods in this particular area. The attraction of the BEM in computational aero-acoustics is the same as it is in standard acoustics, a significant reduction in meshing and hence the potential for faster computation.

Work on the adaption of the BEM to a wider scope of problems has been developed, for example by the dual and multiple reciprocity method, which has also been applied to variants of the Helmholtz equation [46–48,103,131,148]. A generalisation of the boundary element method in acoustics that includes convection, is applicable to problems with uniform flow, but this can also form a useful approximation method with the mean flow rate substituting the value for the uniform flow rate [307]. Similar to the approach in half-pace problems, discussed in Section 5.1.2, aero-acoustic problems are adapted for the BEM by revising the Green's function [308–310]. Recently the BEM in acoustics has been adapted to model viscous and thermal losses [311–313].

Methods in aero-acoustics that use the BEM generally involve setting a fictitious surface, enclosing the significant effects such as noise generation and turbulence, within which typical methods of CFD are used to model the Navier-Stokes equations; the sound generation and sound propagation are modelled separately. The boundary element method models the outer domain, but requiring the fictitious surface to be meshed, modelling the uniform flow and the radiation condition from the boundary and into the far-field. There are several reported implementations and aero-acoustic and related applications of these methods [307,314–326].

#### *5.5. Inverse Problems*

An inverse problem in acoustics, involves measuring properties of the acoustic field and processing that information in order to determine something about its origin. Over recent decades, the main text that guides inverse acoustic (and electromagnetic) (scattering) problems is that of Colton and Kress [16], now in its third edition. Colton and Kress provide a mathematical analysis of inverse problems, and, in acoustics, their focus is on determining the shape of the boundary of a scatterer from the (disturbed) sound pressure at points in the far-field. An example of an application of this is identifying bodies on the sea floor [327].

Inverse problems are ill-posed. This means a range of solutions can give rise to the same measurements, in contrast to the forward problems, studied so far, which are usually well-posed, with a unique solution. In practice, the linear system that is returned by the boundary element method (or any standard method) is significantly ill-conditioned. If a conventional method is used to solve the linear system then the solution will not be acceptable. In the author's previous work on the inverse diffusion problem (classically, the backward heat conduction problem) [328], also included in Visser [237]) it was discussed that there is effectively insufficient information in the data to determine the solution of an

inverse problem. The notion of an 'acceptable' solution, the bias of the observer, provides the final constraint that enables the ill-posed problem to be substituted by a nearby well-posed problem. This returns an acceptable solution, even though it is a less accurate solution of the discrete equations. In the literature, this technique is termed *regularisation*.

Acoustic holography—determining the sound field near the source from acoustic properties at a distance from the source—is one of the main application areas of the inverse boundary element method. In general, an array of pressure or velocity transducers provides the data and the inversion returns the acoustic properties on the surface. Acoustic holography can determine the sources of noise from a radiating structure, which can help guide a re-design. For example, starting with the discrete equivalent of Equation (71),

$$
\underline{\underline{\underline{\underline{\mu}}}}\_{E} = L\_{ES} \underline{\underline{\underline{\mu}}}\_{S} \, \underline{\underline{\underline{\mu}}}\_{S} \, \, \, \, \tag{93}
$$

the field data <sup>ϕ</sup>*<sup>E</sup>* (effectively the sound pressure, see Section 2.1) is related to the layer potential <sup>σ</sup>*S*. If we could find a discrete approximation to σ*S*, then the approximation to the surface potential (sound pressure) and velocity could be found by the matrix-vector multiplication of the discrete equivalent of Equations (53) and (54) and the surface intensity could then be found by the Equation (5). However, as discussed, the solution of the linear system will not yield acceptable results. Even if the number of data points massively exceeds the number of elements, the system is stll effectively underdetermined.

Perhaps the most popular method of resolving the ill-posedness is to use *Tikonov regularisation*. This involves minimising the residual in system, along with a penalty function that constrains the variability of the solution in some sense. For example, Equation (93) is replaced by

$$\min\_{\widetilde{\mathcal{Q}}, \widetilde{\mathcal{Q}}} \left\{ \left\| \underline{\mathcal{Q}}\_{E} - L\_{ES,k} \, \widetilde{\mathcal{Q}}\_{S} \right\| + \eta^{2} \left\| \left\| \underline{\mathcal{Q}}\_{S} \right\| \right\|^{2} \right\},\tag{94}$$

where η is a parameter that can be 'tuned' to achieve and acceptable solution and the norm is the 2-norm.

A related regularisation method is termed truncated singular value decomposition (TSVD). For example, the singular value decomposition (SVD) of the LES,k matrix factorises it as follows:

$$L\_{ES,k} = \mathcal{U}\Sigma V^H,\tag{95}$$

where *U* is a *nE* × *nS* matrix, *V* is *nS* × *nS* and Σ is a diagonal matrix containing *nS* singular values *si*, non-negative values in non-decreasing order. In Equation (95), the *H* denotes the complex conjugate transposed. Let *<sup>U</sup>* <sup>=</sup> # *u*1, *u*2, ... *un* \$ and *<sup>V</sup>* <sup>=</sup> # *v*1, *v*2, ... *vn* \$ , where the *ui* are the left singular vectors and the *vi* are the right singular vectors. For ill-posed problems, the final singular vectors are oscillatory, and the corresponding singular values are very small. Hence, on solution of a system like (93), the oscillations dominate. In TSVD, the offending singular values are simply removed or filtered. The solution may then be determined as

$$\widehat{\underline{\varrho}}\_{\mathcal{S}} = \sum\_{i=1}^{n^\*} \frac{\underline{u}\_i^H \cdot \underline{\varrho}\_E}{s\_i} \underline{v}\_i \,\,\,\underline{v}$$

where *n*<sup>∗</sup> < min (*nS*, *nE*).

There have been several reported implementation of acoustic holography using the methods discussed [237,329–334]. Similar methods have also been developed for finding the impedance of the surfaces in rooms from measured sound pressure data [335,336].

#### *5.6. Meshless Methods*

One of the main advantages of the boundary element method over domain methods, such as the finite element method is the reduced burden of meshing. With meshless methods, a mesh is not required at all. The simplest concept of a meshless method is the equivalent sources method (ESM) in which it is presumed that the acoustic field is equivalent to a field generated by a finite set of point sources:

$$\varphi(p) \approx \sum\_{j=1}^{m} \gamma\_j G\_k(p, q\_j)$$

where the *q<sup>j</sup>* are the positions of the sources that are outside of the acoustic domain and the γ*<sup>j</sup>* are the unknown source strengths and *Gk* is the appropriate Green's function. For the Dirichlet boundary condition, by matching the Dirichlet data on the boundary ϕ *pj* for *p<sup>j</sup>* ∈ *S*, gives a linear system of equations, provided there are at least as many source points as there are items of boundary data. Similarly, by differentiating the above equation with respect to the normal to the boundary

$$\frac{\partial \varphi(p)}{\partial n} \approx \sum\_{j=1}^{m} \gamma\_j \frac{\partial G\_k}{\partial n} (p\_\prime q\_j)\_\prime$$

then approximate solution can be found by matching the values with Neumann data, and using a combination. Usually, more internal points than items of boundary data are chosen and a least-squares solution is sought. The meshless methods avoid the problem of singular integration. However, the ESM is not based on a boundary integral equation formulation and hence it is an alternative to the boundary element method and is therefore beyond the scope of this paper. Lee [337] provides a recent review of the ESM in acoustics.

An alternative method for developing meshless methods has more in common with the standard boundary element method. The method can be applied to the standard interior or exterior problem. For the exterior problem, the methods relates to the Schenck or CHIEF method, of Section 4.3.3, but only the integral equations for the internal points are applied, and the integrals are approximated by using only the midpoint value

$$
\overline{M}\_{DS,k}\underline{\boldsymbol{\phi}}\_{\mathcal{S}} = \overline{L}\_{DS,k}\underline{\boldsymbol{\upmathbb{S}}}\_{\mathcal{S}}.
$$

With the number of interior points exceeding the number of boundary points the solution can be found in the least-squares sense. For the exterior problem, the Burton and Miller form is expected to have superior numerical properties:

$$(\overline{M}\_{DS,k} + \mu \overline{N}\_{DS,k})\underline{\boldsymbol{\varrho}}\_{\mathcal{S}} = (\overline{L}\_{DS,k} + \mu \overline{M}\_{DS,k}^{\mathbf{t}})\underline{\boldsymbol{\varrho}}\_{\mathcal{S}}.$$

These methods still avoid integration, particularly singular integration and can also be used on the modal problem. There are several reports of implementations of these methods [144,145,338–341].

As methods for solving acoustic problems, the meshless methods, are well behind the standard boundary element method in the sense of developing robust software. As with the Schenck or CHIEF method of Section 4.3.3, there is the added issue of determining the placement of the equivalent sources. The methods outlined in Section 3.4.2 could be used in determining whether source points are interior or exterior.

#### **6. Areas of Discussion**

Much of the development work on the boundary element method in acoustic has been on relatively simple shapes and relatively low frequencies. For practical problems, the existing methods must be applicable to more complicated domains and to high frequencies. In this Section, two significantly challenging areas in the acoustic BEM are surveyed. The first is that of efficiency, the relationship between the accuracy achieved and the computational effort. Error analysis of the acoustic/Helmholtz BEM has received significant attention [230,342,343]. Three categories of error arise in the BEM; the discretisation error due to the approximation of the boundary and boundary functions, the quadrature induced error resulting from the numerical integration method and the error in the solution of the linear system of equations. In general, efficiency is maximised by balancing the errors. Much of the

focus in the development of the BEM in acoustics is on improving its efficiency so that a full frequency sweep, particularly including the resolution required to model high frequencies, can be achieved with reasonable computer time and memory requirements. In Section 6.1 the efficiency of the acoustic BEM is discussed and methods for improving efficiency are reviewed.

The most valuable acoustic boundary element method—the standard exterior problem, outlined in Section 5.3—has also been found the most difficult to maintain. The combined integral equations of Section 4.3.4 have held the most promise. The earlier work on the acoustic BEM suggested that the coupling parameter was somewhat arbitrary. In terms of scaling up the method, however, the coupling seems to have become a significant issue and, in Section 6.2, the choice of coupling parameter is revisited.

#### *6.1. E*ffi*ciency*

It is often stated that the boundary element method had a significant efficiency advantage over domain methods, such as the finite element method. However, efficiency concerns remain critical, not least for the acoustic boundary element method, in which it has been a focus of research for many decades. Efficiency, relates the accuracy of the output to the computational resources required to achieve it. The computational resources include the computer processing time and/or the memory requirements. One approach for improving the execution time of numerical methods is to use parallel processing, so that instructions are executed simultaneously, rather than sequentially, and such techniques have been applied in the boundary element method [304,344].

In this Subsection, the efficiency of the acoustic BEM is analysed and techniques for improvement are reviewed. Although a range of classes of boundary element methods have been outlined in this paper, the analysis is fairly generic, and can be applied to the boundary-value problems that arise in acoustics. The modal analysis or eigenvalue problem is not considered, but much of the analysis carries over. For the extended problems considered in Section 5, the analysis in this Subsection is only relevant insofar as the BEM is implemented with the wider context.

#### 6.1.1. The Computational Cost of the Acoustic Boundary Element Method at One Frequency

In this Sub-subsection, the boundary element method in acoustics is developed in its most typical way for one frequency. The efficiency of the method is analysed and discussed. More disruptive methods for improving the efficiency are considered subsequently.

As discussed, the boundary element method is composed of two stages, the first stage involving determining the boundary functions and the domain solution in the second stage. In the first stage one or more *n*·*n* matrices are formed, where *n* is the number of boundary solution points (e.g., collocation points, or elements for simple constant elements). In the second stage, usually one *nP*·*n* matrix is formed, where *nP* is the number of domain points. Hence the storage requirement is *O n*<sup>2</sup> + *nPn* complex numbers.

The matrices are normally computed by numerical integration. In general, the number of quadrature points required on each panel would be varied with the size of the panel and the distance between the panel and the point [345]. For each quadrature point, the Green's function would be computed, together with the other geometric information, although the former is likely to incur the greater part of the computational cost. The discretisation could be carried out separately for each operator; however, for efficiency reasons, the computed values should be shared between the operators, as, for example, followed in the author's previous work [61,91]. The diagonal components of the square matrix (or matrices) in the initial stage may be the result of evaluating singular and hypersingular integrals, as discussed, and, although they require special treatment, it is usually more efficient if the quadrature points are similarly the same for all required discrete operators. Lumping together the matrices, for the *n* × *n* and *nP* × *n* components, let the average number of quadrature points be *NQ*

and let the average cost of the evaluation of the Green's function be *CG* and the average cost of the geometrical properties be *Cg*, then the total cost (or computer time) of computing the matrices is

$$N\_Q(n^2 + n\_P n)(C\_G + C\_{\mathcal{S}})\,.$$

Clearly *Cg* will be larger for the combined operators, used in the primary stage for exterior problems. Alternatively, for the Schenck method, the matrices are augmented by the discrete operators for the interior points.

Once, the matrix vector system is formed, the next step is to solve it. The most straightforward method of solution of a square system is the Gaussian elimination or LU factorisation method (re-writing a matrix as the product of a lower and upper triangular marices). The overall computer time for the acoustic boundary element method at one frequency can be summarised with the equation

$$N\_Q(n^2 + n\_P n) \left( C\_G + C\_{\mathcal{S}} \right) + \left( \frac{4}{3} \ n^3 + O(n^2) \right) C\_f = 0$$

where *Cf* represents the cost of a floating-point operation, such as multiplication or division. The non-square system that arises in the Schenck method can also easily be resolved as a square system through pre-multiplying both sides by the transpose of the matrix over which the solution is sought and this is equivalent to the least-squares solution.

In its earlier development, the cost of computing the matrices generally far outweighed the cost of solving the system. However, in the first stage of the boundary element method, the computational cost of setting up the linear system is *O n*2 , whereas the cost of solving it (using the methods stated) is *O n*3 . For this reason, it has also been understood for a long time that LU factorisation or related methods for solving the linear system was unsustainable, as progress is made towards the finer meshes that are particularly required for high frequency problems. However, it is important also to state the particular advantages that LU factorisation has. LU factorisation is a robust method and, in the case when there are a range of inputs to a problem with a fixed boundary and boundary condition, once the *O n*3 factorisation as been stored, it can be used repeatedly again with *O n*2 cost.

The fast multiple method (FMM) is a popular method of speeding up the computation of the matrices in the boundary element method and it has been applied to acoustic/Helmholtz problems [220,227,254,284,287,290,320,346]. In this method the Green's function is approximated by local polar expansions that can be translated through the domain, rather than re-computed. The FMM is often used with iterative methods.

In this Subsection alternatives to standard LU factorisation are reviewed. Faster solution methods include the use of hierarchical matrices [190] or panel clustering [347] and iterative methods are also considered. If the solution of the linear system is potentially dominant, in terms of computational cost, then the attention also re-focuses on the integral equation method. For it may be advantageous to choose an integral equation method that results in a linear system that has faster convergence properties, rather than one that is the easiest to apply or the most efficient (such as the collocation method that is highlighted in this work).

#### 6.1.2. The Frequency Factor, Multi-Frequency Methods and Wave Boundary Elements

Typically, in vibration and acoustics, the time-dependent signals are resolved into frequency components. For example, for air-acoustics, the audible range for human being is up to 20 kHz and typically his could be resolved into components with a 10Hz resolution; multiplying the computational cost, considered in the previous Sub-subsection, by around 2000. The prospect of solving acoustic problems over the full frequency range has led to another set of techniques focus on reducing the computational burden by solving a range of frequencies together and these are usually termed multi-frequency methods [297,348–355]. However, with the nature of an acoustic field, originating typically from structural vibration, is such that structural and acoustic resonances, and the driving

profile, dominate the total acoustic response. From a computational point of view, effort should therefore be concentrated on these 'peak' frequencies. Hence, this adaption returns us to the standard mono-frequency problem, or applying the multifrequency method, centred on each peak.

Moreover, in general the frequency of observation is matched in the acoustic and vibratory properties. In practice, as the frequency increases, a higher elemental resolution is required and hence *n* = *O*(*k*) for two-dimensional and axisymmetric problems and *n* = *O k*2 in three dimensions. Although this seems to imply a revision of the mesh at every frequency, it is more likely that one mesh that is suitable for the highest frequency is used throughout the range, or separate meshes are applied in significant sub-ranges of the frequency range, in order that the burden of meshing is proportionate within the overall project.

A review of the actual number of elements required to capture the solution is provided by Marburg [356–358]. In general, for the simple elements that are often used, it is considered that 6-10 elements per wavelength is a reasonable guideline. Obviously, higher-order boundary element methods [100,359–363] would often return the same accuracy with fewer elements. There is significant interest in isogeometric elements, in which the boundary and the boundary functions are modelled with the same basis functions (typically splines), so that the BEM can be more easily integrated with computer-aided design [364–371]. For air acoustics, at high frequencies reaching 20 kHz, the wavelength is less than 2 cm. For example, even for a 10cm cube, the number of elements required to model the high frequencies is in the tens of thousands and for a 1m cube, over a million elements would be required. A variation on the boundary element method has arisen for developing sinusoidal basis functions or wave boundary elements, in order to more accurately model the acoustic functions and the term PUBEM or partition of unity BEM is often used to identify such methods [267,371–374], which has similarity with the application of the Treffetz method [292–294].

#### 6.1.3. Iterative Methods and Preconditioning

The scalability of the standard BEM in acoustics has been in question for several decades. The computational estimate relates two areas of particular concern, the *O*(*n*3) nature of LU factorisation and the *O*(*n*2) cost of computing the matrices. Hence, much of the current research on the BEM in acoustics is focused on reducing its computational burden, so that the methods can be casually applied to more significant problems and high frequencies.

The conjugate gradient method was identified as a useful iterative solution method for the acoustic BEM [155,375,376]. The conjugate gradient method and related methods are generally termed Krylov subspace methods and these variants have been significantly tested on the linear systems arising from the boundary element method in acoustics [106,377–379]. In general, if the underlying operator is compact, corresponding to a clustering of the corresponding matrix eigenvalues, then iterative methods, such as the Krylov subspace methods, work well [380]. However, one of the operators, *Nk*, is not compact, the eigenvalues of its matrix equivalent are spread out, and hence the raw iterative methods are only applicable in particular cases, in which the hypersingular operator is not in play. Hence, in the spirit of generality that is sought in this work, the iterative methods are of limited value in solving the untreated equations. However, in the acoustic boundary element method, these iterative methods are usually applied following an intervention with the original linear system, termed preconditioning [381].

In the acoustic BEM, preconditioning has been advocated since its early days. The original concept can be derived from the integral equation formulations. By various substitutions with the boundary integral equations, operator identities can be derived. One of the most useful is the substitution of (74) and (75) into (70) (*p* ∈ *S*) (or substituting (63) or (64) into (48) (for *pS*)) giving

$$L\_k N\_k = \left(M\_k - \frac{1}{2}I\right)\left(M\_k + \frac{1}{2}I\right) = \left(M\_k + \frac{1}{2}I\right)\left(M\_k - \frac{1}{2}I\right) = M\_k^{-2} - \frac{1}{4}I\*$$

and hence the left-hand side of the Burton and Miller Equation (76) may be re-written as follows

$$\begin{split} \left\{ (\mathcal{M}\_{k} - \frac{1}{2}I + \mu L\_{k} N\_{k}) \varphi \right\}\_{\mathcal{S}} \left( p; n\_{p} \right) &= \left\{ (\mathcal{M}\_{k} - \frac{1}{2}I + \mu (\mathcal{M}\_{k} + \frac{1}{2}I)(\mathcal{M}\_{k} - \frac{1}{2}I)\_{k}) \varphi \right\}\_{\mathcal{S}} \left( p; n\_{p} \right) \\ &= \left\{ (L\_{k} + \mu L\_{k} (\mathcal{M}\_{k}^{t} + \frac{1}{2}I) \upsilon \right\}\_{\mathcal{S}} \left( p; n\_{p} \right). \end{split} \tag{96}$$

The equation is preconditioned, as the *Nk* operator has been eliminated and all the operators are compact. However, the motivation for this technique can be with the elimination of the hypersingular integration, with incidental preconditioning. Methods based on (96) along with *LkNk* = *M<sup>t</sup>* 2 *<sup>k</sup>* <sup>−</sup> <sup>1</sup> <sup>4</sup> , that is obtained through the substitution of (53) and (54) into (61) or the substitution of (71) and (72) into (73) (*p* ∈ *S*), are often termed the Calderon equations [382,383]. Methods that make use of this substitution have been developed [155,157,241]. Alternative approaches in developing well-conditioned integral formulations for iterative solution have also been the subject of research [384–388]. Langou [389] develops methods for the iterative solution of linear systems with multiple right-hand sides, echoing the stated advantage of the LU factorisation method, discussed in Section 6.1.1.

Matrix preconditioners are often based on constructing an approximate inverse and following a fixed point/ defect correction/ contraction method. Introducing a preconditioning matrix also introduces an *O*(*n*2) matrix multiplication at each step, but this can be reduced to *O*(*n*) if the approximate inverse matrix is sparse or banded [378,380]. Methods based on an incomplete LU factorisation [390] have been tested in the acoustic BEM [391]. A similar approach is the construction of a *DtN* (or *DN*) map (Dirichlet to Neumann) or *NtD* (or *ND*) map, also termed an on surface radiation condition (OSRC), as the approximate inverse [190,392–394].

#### *6.2. The Coupling Parameter*

The direct and indirect integral formulations of the exterior problem that that were free of the characteristic frequencies or non-uniqueness, were introduced in Section 4.3.4. The equations were defined with a coupling parameter μ, which weighs the contribution of derivative equation with the original equation in the direct method, or between the double layer potential and the single layer potential equations in the indirect method. Mathematically, the coupling parameter has a non-zero imaginary part that ensures that the equations are free of characteristic wavenumbers, and, in general, μ is presumed to be an imaginary number. In accepting that μ is imaginary, the starting point is −∞ < *Im*(μ) < ∞. If |μ| is very small or very large then one or the other of the original equations would dominate, and the issues with the dominant equation would be evident with the hybrid equation, and this narrows its range; 0 |μ| ∞.

The argument quickly shifted from introducing μ in order to avoid the potential catastrophic errors in the original equations of Section 4.3.1 to considering an optimal value. For several decades, based originally on the works of Kress [395–398] and the further works of Amini [399], there has been a strong recommendation, with supporting research, that μ = *<sup>i</sup> <sup>k</sup>* is reasonably close to the 'optimal' choice, as least for the simple boundaries, like spheres. Obviously *<sup>i</sup> <sup>k</sup>* is unsuitable for low wavenumbers, as the parameter would be large, violating the condition |μ| ∞, and the second equation would dominate and a cap on its value is recommended, for example,

$$\mathfrak{u} = \begin{cases} 1 \text{ if } k \le 1 \\ \frac{i}{k} \text{ if } k > 1 \end{cases}$$

The term 'optimal' in these paragraphs, refers to the condition of the combined operators over which the equation is being solved. The stronger the condition of the operator, the less able it is to magnify the solution, or magnify its error. The original equations, on their own, have 'spikes' of ill-conditioning around their respective eigenvalues. The rationale is that in combining the equations, these spikes are significantly reduced, and the condition of the combined operator steers an even keel

through the frequencies. It has been shown, for example for a sphere that *<sup>i</sup> <sup>k</sup>* provides a generally well-conditioned combined operator. A recent paper by Zheng [400] showed that each eigenvalue of the combined operator follows a loop-shaped locus as μ varies, joining the real axis when |μ| = 0, ∞, and supporting μ = *<sup>i</sup> <sup>k</sup>* , as this reasonably approximated the point on the locus that was furthest from the real axis. Zheng terms the inclusion of the derivative equation as 'adding damping', relating the analogy with, for example, Equation (91), wherein the damping term similarly shifts the eigenvalue off the real axis; the combined integral equations are not free from irregular frequencies, rather they are simply moved from the real axis into the complex plane.

There is a simple case for the parameter with an inverse proportionality with the wavenumber. Considering the three-dimensional case, the Green's function (52) is the kernel of the *Lk* operator, and its magnitude does not change with *k eikr* = <sup>1</sup> . Hovever its derivative,

$$\frac{\partial G\_k}{\partial r} = \frac{e^{ikr}}{4\pi r^2} (ikr - 1),\tag{97}$$

a factor in the kernel of the *Mk* and *M<sup>t</sup> <sup>k</sup>* operators, is *<sup>O</sup>*(*k*). Hence, combining the operators of *Lk* with *Mk* or *M<sup>t</sup> k* , as in Section 4.3.4, a parameter that is inversely proportional to *k* equalises the contribution from the two operators. Similarly, differentiating again, the kernel of the *Nk* is *O*(*k*2), and the same parameter has a similar function, when *Mk* or *Mt <sup>k</sup>* are combined with it.

A recent paper by Marburg [401] considered principally the sign of μ. Marburg noted that many papers had inadvertently used <sup>μ</sup> = <sup>−</sup> *<sup>i</sup> <sup>k</sup>* as the coupling parameter. Given that its value was not thought to be critical, it might be thought that this would work as well as μ = *<sup>i</sup> <sup>k</sup>* . However, Marburg demonstrated that *<sup>i</sup> <sup>k</sup>* usually returns significantly more accurate results, and much faster convergence with the iterative methods for solving the systems of equations. The parameter <sup>μ</sup> = <sup>−</sup> *<sup>i</sup> <sup>k</sup>* was not used out of choice in the papers reviewed therein, but rather through the confusion of the signs that underlie the basic definitions in the mathematical model. These results have been confirmed by Galkowski [116].

If *Lk* is coupled with its derivative, as it is in the Burton and Miller equation, then the following equation is obtained, for one side of the equation

$$\begin{split} \left\{ (L\_k + \mu (M\_k^t + \frac{1}{2}I)v )\_S \right\} (p; \mu\_p) &= \int\_S (G\_k + \mu \frac{\partial G\_k}{\partial r} \frac{\partial r}{\partial n\_p}) v(q) dS + \frac{1}{2} \mu v(p), \\ &= \int\_S \frac{e^{ikr}}{4\pi r^2} \Big( r + \mu ikr \frac{\partial r}{\partial n\_p} - \mu \frac{\partial r}{\partial n\_p} \Big) v(q) \, dS \, + \frac{1}{2} \mu v(p). \end{split}$$

With μ = *<sup>i</sup> <sup>k</sup>* there is significant cancellation in *<sup>r</sup>* <sup>+</sup> <sup>μ</sup>*ikr* <sup>∂</sup>*<sup>r</sup>* <sup>∂</sup>*np* when <sup>∂</sup>*<sup>r</sup>* <sup>∂</sup>*np* ≈ 1, possibly with a diagonalising effect on the operator, and this is also presumed for the other combinations of operators.

Another approach is to consider the condition of the matrices that arise in the BEM [402–404]. Earlier, it was discussed that the *Nk* operator was a pseudo-differential operator, rather than an integral operator. As a result of this, the *NSS*,*k* increases with the number of elements *ns*, whereas for the other integral operators, the norm of the matrix is independent of the number of elements. This is also supported by Equations (45) and (46); the diagonal components of the *NSS*,*<sup>k</sup>* matrix are *O* 1 *h* = *O*(*nS*) for two- dimensional and axisymmetric problems and *O* <sup>√</sup>*nS* for three-dimensional problems. It follows that *NSS*,*k* = *O*(*nS*) or *O* <sup>√</sup>*nS* . The scene is therefore set for two conflicting interests in choosing the coupling parameter μ = *<sup>i</sup> <sup>k</sup>* is near-optimal in conditioning the combined operator, as a result of mathematical analysis whereas μ = *<sup>i</sup> nS* or <sup>μ</sup> <sup>=</sup> *<sup>i</sup>* <sup>√</sup>*nS* (put simply), balancing the matrix norms, as a result of numerical analysis. Without the latter correction, the *NSS*,*<sup>k</sup>* matrix will apparently provide an increasing dominant potential for (quadrature-induced) error as the number of elements increases.

However, in Section 6.1, the modern form of the combined method was outlined, involving pre-conditioning and iterative methods for solving the linear systems. For example, with the preconditioner in (96), the *Nk* operator is avoided and, in such cases, the analysis of the previous two paragraphs is not applicable.

#### **7. Concluding Discussion**

The main purpose of this paper is to encapsulate the modern scope of the boundary element method in acoustics; how it can be adopted for general 2D, 3D and axisymmetric dimensions, interior, exterior and half-space problems, thin shells and modal analysis. The standard BEM can be directly applied to an enclosed domain with cavities or to multi-boundary problems in an exterior domain. Boundary element methods are founded on a boundary integral equations formulations and these can be fused to form additional methods: thin shells can be used to model discontinuities in an interior or exterior domain defined by closed boundaries, half-space problems can be modelled by the Rayleigh integral or modifying the Green's function, domain methods can be linked to the BEM so that vibro-acoustic and aero-acoustic problems can be tackled and—through regularization of the BEM equations—inverse problems can be addressed.

The boundary element method is a numerical method that only becomes a potentially useful tool for solving real-world problems in acoustic engineering when it has been implemented in software. Clearly, it is important that the executing code fits within the memory of the computer and completes the computation in reasonable time, and efficiency issues have been considered in Section 6.1. There are issues of generality and, in that regard, this work has focused on the generalised Robin boundary condition and general topologies. Reliability and maintenance are also very important issues in software development, and hence much of the focus has been on the standard exterior acoustic problem, which has had issues and solution approaches documented for more than half a century. The robustness of the BEM, particularly considering the validity of a defined elemental boundary, has had little attention in the field, but methods for assisting with this are summarised in Section 3.4.2.

There is much commonality across the various topologies to which the boundary element method can be applied in acoustic/Helmholtz problems; the equations for different dimensional problems are literally the same, with a change in the definition of the Green's function and line integrals become surface integrals when we move from 2D to 3D. All the core equations, whether they are interior, exterior, shell problems or the straightforward half-space problems, have the same essential operators within a chosen dimensional space. Hence there is significant scope for component-wise development of software, adopting a 'library' approach and unifying the method.

Significant issues have always surrounded the *Nk* operator. This operator was initially included in the combined integral equation formulations for the exterior problem, outlined in Section 4.3.4. Firstly, for surface (collocation) points its integral definition is hypersingular, which is difficult to interpret and evaluate, and perhaps for that reason alone, the alternative Schenck/CHIEF method is the preference of many. Once this significant issue is surpassed, the *Nk* is not compact, it is a pseudo-differential operator and this causes further issues for solving over by iterative methods and hence preconditioning has been introduced to circumvent this. Preconditioning can eliminate the hypersingular operator *Nk* for the standard exterior problem (96). However, *Nk* is still required for shell elements and hence it has to be included in a general library.

Much of the current research theme in the area, put simply, is to shunt the BEM in acoustics from its traditional comfort-zone of problems with <sup>∼</sup>10<sup>3</sup> elements to modern problems, to include high frequencies, with <sup>∼</sup>106 elements and from straightforward problems to complicated applications [190,238,405]. The computational bottle-neck in the traditional BEM is in the use of LU factorisation or similar methods for solving the linear system. Hence the LU factorisation is the first necessary casualty of this move and iterative methods are favoured, although this has been presaged for several decades. The shift to <sup>∼</sup> 106 elements could also have a significant demand on computer memory and hence the expectation of storing matrices must also be relaxed. With iterative methods, the effect of the approximation of the matrix components is more controllable and methods such as

the fast multipole method and panel clustering are more able to broker the issues of storage and computation time to this end.

Probably the majority of research on the BEM in acoustics is on the seemingly intractable exterior problem. Our expectation that the solution in its infinite domain can be thread through the boundary is ultimately a questionable one, as the method is scaled up. The combined methods, and most prominently, the Burton and Miller method, have held centre-stage in this endeavour, seen by most as numerically superior to the Schenck/CHIEF method [238]. However, the two erstwhile competing methods may have to be married in our best efforts to achieve <sup>∼</sup>106 elements and high frequencies; using the equations from interior points to provide more stability in the Burton and Miller equation, as they did originally with the elementary equations to form the CHIEF method. The system could be augmented further with directional derivates of Equation (70) for points in the interior *D*.

The combined methods throw up two issues, one is the inclusion of the *Nk* operator, as discussed, and the other is the coupling parameter. The systems of equations arising in the BEM using combined methods require that the *Nk* operator is preconditioned in order to achieve convergence with iterative methods, as discussed in Section 6.1. The choice of coupling parameter is discussed in Section 6.2. However, the approach to choosing the coupling parameter is already to potentially optimise the system and hence it is itself a pre-conditioner, as pointed out in Betcke et al. [190]. In Harris and Amini [406], the coupling parameter is generalised from a singular value to a function over the surface. Hence, another marriage is proposed. As the coupling parameter and the preconditioner have a similar purpose, then the Burton and Miller Equation (76), for example, could take the following form

$$\left\{ (\mathcal{M}\_k - \frac{1}{2}I + \mathcal{U}\_k \mathcal{N}\_k) \varphi \right\}\_S \Big( p; n\_p \Big) = \left\{ (L\_k + \mathcal{U}\_k (\mathcal{M}\_k^t + \frac{1}{2}I) \upsilon \right\}\_S \Big( p; n\_p \Big).$$

where *Uk* is the preconditioning operator, fusing the coupling parameter and the preconditioner into one entity. The objective then is to determine *Uk*, or an analogous preconditioning matrix, *USS*,*k*. For example, following the method in Equation (96), *Uk* = *<sup>i</sup> <sup>k</sup>Lk*; however, in general, the goal is to set *Uk* or *USS*,*<sup>k</sup>* to best-prepare the system for iterative solution.

**Acknowledgments:** The author would like to thank the anonymous reviewers for their enthusiasm for the paper. The author is also very grateful for the reviewers' many useful points and further references, that have helped to significantly improve on the original manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **The Role of Powered Surgical Instruments in Ear Surgery: An Acoustical Blessing or a Curse?**

#### **Tsukasa Ito \* , Toshinori Kubota, Takatoshi Furukawa, Hirooki Matsui, Kazunori Futai, Melinda Hull and Seiji Kakehata**

Department of Otolaryngology, Head and Neck Surgery, Yamagata University Faculty of Medicine; Yamagata 990-9585, Japan; t-kubota@med.id.yamagata-u.ac.jp (T.K.); t-furukawa@med.id.yamagata-u.ac.jp (T.F.); matsuihirooki@gmail.com (H.M.); kfutai@gmail.com (K.F.); melindahull@nifty.com (M.H.); seijik06@gmail.com (S.K.)

**\*** Correspondence: tuito@med.id.yamagata-u.ac.jp

Received: 23 January 2019; Accepted: 19 February 2019; Published: 21 February 2019

**Abstract:** Ear surgery in many ways lagged behind other surgical fields because of the delicate anatomical structures within the ear which leave surgeons with little room for error. Thus, while surgical instruments have long been available, their use in the ear would most often do more damage than good. This state of affairs remained the status quo well into the first half of the 20th century. However, the introduction of powered surgical instruments, specifically the electric drill used in conventional microscopic ear surgery (MES) and the ultrasonic aspirator, the Sonopet® Omni, in transcanal endoscopic ear surgery (TEES) marked major turning points. Yet, these breakthroughs have also raised concerns about whether the use of these powered surgical instruments within the confines of the ear generated so much noise and vibrations that patients could suffer sensorineural hearing loss as a result of the surgery itself. This paper reviews the intersection between the noise and vibrations generated during surgery; the history of surgical instruments, particularly powered surgical instruments, used in ear surgeries and the two main types of surgical procedures to determine whether these powered surgical instruments may pose a threat to postoperative hearing.

**Keywords:** noised-induced hearing loss; powered surgical instruments; ultrasonic aspirator; transcanal endoscopic ear surgery

#### **1. Introduction**

The internal anatomy of the ear is made up of extremely tiny, delicate, and interlocking anatomical structures that are surrounded by bone and muscle, with sound traveling through the external auditory canal as shown in Figure 1. In particular, the mastoid portion of the temporal bone lies behind the ear and serves as a solid, normally impenetrable, barrier protecting the internal ear. This bony barrier has made it particularly challenging to access the anatomical structures within the internal ear, which has primarily been accomplished by drilling straight through and removing the mastoid bone in a procedure called a mastoidectomy. This procedure has been the mainstay of ear surgery up until the turn of the 20th century. While most ear surgery procedures are performed today with the objective of either preserving or improving hearing, the potential exists for noise and vibration generating surgical instruments required in a mastoidectomy to damage hearing. Of special concern has been the adverse impact on hearing of the use of powered surgical instruments, particularly drills, and the more recently introduced ultrasonic devices which are used to remove bone and expose the internal anatomy of the ear. This paper discusses issues available in the literature which have reported on the effects on hearing of these powered surgical instruments.

**Figure 1.** The internal anatomy of the ear.

The majority of current ear surgery procedures can be broken down into two broad approaches: conventional microscopic ear surgery (MES) and the more recently developed transcanal endoscopic ear surgery (TEES). While powered surgical instruments are a standard part of MES, such instruments are only used in a subset of TEES procedures referred to as powered TEES. As was noted above, TEES has not totally eliminated the need for a mastoidectomy, because some surgeons are more comfortable with MES, and some procedures are not indicated for TEES, such as the mastoid air cells, inner ear, and skull base that are beyond its reach.

This paper presents an overview on types and causes of hearing loss in the context of ear surgery; a brief history of surgical instruments used to access the internal anatomy of the ear focusing on the powered surgical instruments, specifically electric drills and ultrasonic devices; an overview of the two main surgical procedures of MES and TEES used in ear surgery to access the internal anatomy of the ear; and a review of the literature on the potential for hearing loss caused by noise and vibrations produced by powered surgical instruments used in MES and TEES procedures, together with the presentation of previously unpublished data on the use of the ultrasonic aspirator, a powered surgical instrument, used in TEES.

#### **2. Types and Causes of Hearing Loss**

#### *2.1. Types of Hearing Loss*

Hearing loss falls into three broad categories: conductive hearing loss, sensorineural hearing loss and mixed hearing loss. Conductive hearing loss generally occurs when a physical impediment or barrier prevents the transmission of sound waves through the pathway from the outer ear through to the middle ear. Such impediments or barriers can range from a simple buildup of ear wax, accumulation of fluid within the ear due to an infection or an abnormal growth such as bony tissue or a tumor. The other type of hearing loss, sensorineural hearing loss, can be attributed to problems within the inner ear, primarily the cochlea and associated hair cells or the vestibulocochlear nerve (cranial nerve VIII). Sensorineural hearing loss can be caused by either intrinsic factors such as genetics resulting in congenital abnormalities, or extrinsic factors such as inner ear infections; ototoxic drugs such as aminoglycosides and cisplatin; or exposure to high noise levels both over an extended period of time such as in an industrial workplace, prolonged use of headphones or concert/entertainment venues or a single discrete event such as a blast of noise from equipment, gun shot, or bomb blast. The third type of hearing loss, mixed hearing loss, as the name implies, is a combination of the other two types of hearing loss [1].

#### *2.2. Hearing Loss and Powered Surgical Instruments*

Hearing loss as related to powered surgical instruments has primarily been studied from two perspectives: noise levels (air-conducted) and vibrations (bone-conducted). The vast majority of studies have focused on noised-induced hearing loss, which is a clear subcategory of sensorineural hearing loss, while vibrations, or more precisely, skull vibrations, have garnered much less attention until recently and deserve further study and consideration. Such hearing loss is measured based on the degree to which the hearing threshold sensitivity has risen and is classified as either a permanent threshold shift or a temporary threshold shift [2]. Most sensorineural hearing loss caused by powered surgical instruments fortunately falls into the temporary threshold shift category.

Separate from these two types of hearing loss caused by powered surgical instruments is the physical contact of an instrument with the ossicular chain. Such drill-to-bone contact results in vibrations which are transmitted via the ossicular chain to the cochlea. The subsequent damage to the cochlea generally results in permanent hearing loss. However, this hearing loss can be attributed to surgeon error rather than conventional use of the surgical instrument itself.

#### 2.2.1. Noised-Induced Hearing Loss

Noised-induced hearing loss can, as alluded to above, be caused by either chronic, accumulative, and gradual exposure, or an acute, one-time event. The chronic, accumulative, and gradual exposure is generally defined in terms of daily exposure over years. Specifically, the National Institute for Occupational Safety and Health (NIOSH) has set down the recommended exposure limit (REL) as 85 decibels, A-weighted (dBA) for an 8-h time-weighted average (TWA) [3] while the acute, one-time event is generally set at 140 dB or higher [4]. However, any noise-induced hearing loss caused by powered surgical instruments falls into an undefined category between these two defined categories. While the noise levels tend to fall within chronic category, surgery time is only measured in minutes up to a couple of hours on a single day, rather than accumulated hours over multiple days. Moreover, while the noise levels generated by powered surgical instruments fall below noise levels defined as dangerous for acute one-time events, the fact that these instruments are used directly within the anatomy of the ear needs to be factored in.

The connection between possible noise-induced hearing loss and powered surgical instruments used during ear surgery has long been a source of concern and study in the case of surgical drills [5–8], as well as a recent target of research in the case of ultrasonic devices [9].

#### 2.2.2. Vibration-Induced Hearing Loss

The second potential way that powered surgical instruments can cause sensorineural hearing loss is by skull vibrations generated by these instruments. This cause has not been totally ignored, but has not received the same amount of attention over the years as noise levels. Moreover, occupational standards are not widely codified, particularly in regard to hearing, as opposed to damage to the circulatory system through the use of hand-held heavy equipment which can result in what is known as hand-arm vibration syndrome (HAVS), a type of Raynaud's syndrome [10]. Vibrations have been posited to cause damage to hearing through inner ear damage and Seki et al. [11] and Miyasaka [12] have both posited that morphological changes, specifically permeability, occur in stria vascularis capillaries, when the auditory ossicles or mastoid are subjected to vibrations. Again, as was the case with noised-induced hearing loss, vibration-induced hearing loss has been studied with surgical drills [10,13,14], but less so with ultrasonic devices [15].

#### **3. History of Otological Surgical Instruments**

Progress in ear surgery has, like other surgical fields, been driven primarily by progress in the development of appropriate instruments. These instruments were often originally designed with another purpose in mind, often dentistry, but were eventually adapted by enterprising ear surgeons who co-opted them for their own purposes. Mudry has reported on the history of instruments in ear surgery and divided this history into three periods: trepans (Figure 2a); chisels and gouges (Figure 2b): and electrical drills (Figure 2c [16]. The introduction of powered TEES has added a fourth period of ultrasonic devices (Figure 2d) [17].

**Figure 2.** (**a**) Surgical instruments for trepanning. Engraving with etching by T. Jefferys. Credit: Wellcome Collection. CC BY (**b**) Tools for mastoidectomy [18]. (**c**) Visao® High-Speed Otologic Drill, Medtronic, Shinagawa, Tokyo, Japan). (**d**) Sonopet® Omni, (model UST-2001 Ultrasonic Surgical Aspirator, Stryker Corporation, Kalamazoo, Michigan, USA).

#### *3.1. Pre-Powered Surgical Instruments*

Though scattered references have been found throughout history that may be describing a mastoidectomy-like procedure such as by Galen of Pergamon (129 A.D−210 A.D.) [19]. Surgical instruments had only developed to the point where a mastoidectomy could seriously be contemplated in the 18th century. These developments were attributed to advances made in metallurgy and machine tooling that fueled the Industrial Revolution and underscore the critical importance of technological advances outside medicine, leading to breakthroughs in medical procedures.

#### *3.2. Powered Surgical Drills*

The powered surgical drill, as opposed to a hand drill type instrument such as a trepan, was first used in the 1880s. However, its use did not catch on, and like many new approaches, was a bit ahead of its time. Instead, the electric drill was reintroduced into otological practice in the late 1920s by Julius Lempert, who is commonly recognized as the father of the use of the electric drill in ear surgery [16]. The electrical drill has been the workhorse in ear surgery up until the present day.

#### *3.3. Ultrasonic Devices*

A new type of device technology began to make an appearance in the medical literature around the 1970s: ultrasonic devices. Broadly speaking, three types of ultrasonic devices used in surgery have appeared on the market: cavitation ultrasonic aspirators or CUSA devices; piezoelectric devices; and ultrasonic aspirators equipped with longitudinal and torsional oscillation.

A CUSA instrument selectively targets the fluid in soft tissue such as tumors, but leaves hard tissue such as bone untouched [20]. CUSAs are quite frequently used in neurosurgery [21] and renal surgery [22] and have enabled surgeons in these fields to remove tumor tissue located in areas with little room for error with less possibility of damaging anatomical structures in comparison to electric drills.

Piezoelectric devices are another type of ultrasonic device that can cut both bone and soft tissue depending on the frequency setting [23]. These devices are another example of borrowing a tool primarily used in dental surgery and expanding its reach into a wide variety of fields, including craniofacial surgery, to perform osteotomies [24] and ear surgery to perform

mastoidectomies [25]. However, some surgeons have found that piezoelectric devices compare unfavorably to electric drills in terms of bone cutting power and speed [9].

A third type of ultrasonic device offers the reverse of the CUSA by targeting hard tissue such as bone and leaving soft tissue relatively untouched. Both the CUSA and this new type of ultrasonic aspirator that removes only bone are often referred to in the literature as an ultrasonic bone curette (UBC) because they can both use the same handpiece to which a specialized tip is attached. These third types of ultrasonic aspirators remove bone using a longitudinal oscillation (L mode) or longitudinal and torsional oscillation (LT mode) to emulsify bone and these L mode/LT mode UBCs have been used in fields such as neurosurgery [26], paranasal sinus surgery [27], maxillofacial surgery [28], and spinal surgery [29].

Our original paper on powered TEES referred to the ultrasonic aspirator which we used as a Sonopet® UBC [17]. However, this ultrasonic aspirator that we used and still use today, emerged on the market in the early 2000s and has passed through a number of companies [30]. This ultrasonic aspirator is today more properly called the Sonopet® Omni and uses an H101 tip (model UST-2001 Ultrasonic Surgical Aspirator, Stryker Corporation, Kalamazoo, Michigan, USA) (Figure 3).

**Figure 3.** (**a**) Sonopet® Omni, (model UST-2001 Ultrasonic Surgical Aspirator, Stryker Corporation, Kalamazoo, Michigan, USA). (**b**) Bone removal view (**c**) Aspirator view.

It is the Sonopet® Omni together with its H101 tip that has allowed us to perform powered TEES. The ability of this Sonopet® Omni to target bone while generally avoiding soft tissue made it the perfect tool for working within the narrow confines of the external auditory canal and for removing bone from the ear canal. In addition, the Sonopet® Omni not only removes bone, but also has irrigation and aspiration functions that were critical in opening the door to powered TEES, as will be described below. However, little information has been reported on the safety of the Sonopet® Omni in terms of noise levels and vibrations [15], which will also be discussed below.

#### **4. History of Otological Surgical Procedures**

The primary otological surgical procedure has been the mastoidectomy. A mastoidectomy involves the removal of the bone behind the ear in order to access the internal anatomy of the ear. However, the tools used in performing a mastoidectomy and related objectives and considerations in terms of safety have changed dramatically over the years. Moreover, since around the turn of the 20th century, a new surgical approach, TEES, has emerged that circumvents the need for a mastoidectomy. It should be noted that TEES has not totally eliminated the need for a mastoidectomy; some surgeons have yet to adopt the procedure and some procedures are not indicated for TEES. Moreover, TEES has raised its own safety issues that need to be addressed.

#### *4.1. Pre-MES Procedures*

Once tools were available for removing bone, humans opened holes in skulls for various purposes. The objective in prehistoric and early historic times could have been to release evil spirits, while other

procedures may have been to relieve a buildup of pressure within the skull. This pressure buildup could be due to the presence of excess fluids caused by internal hemorrhaging or infection. Prior to the development of antibiotics and the appropriate tools, the main concern of otologists was the treatment of middle ear infections that resulted in pus accumulating within the ear, and having nowhere to go, were potentially life threatening. Such procedures often involved simple incisions of the tympanic membrane or abscesses located inside or outside of the ear. The French physician Ambroise Paré was the first to be credited in the 16th century with making a surgical incision to drain such an infection [19].

The first confirmed mastoidectomy is credited to the French surgeon Jean Louis Petit in 1736 [16] and thereafter the mastoidectomy became a standard, but not universally practiced part of ear surgery from around 1860 [19]. A mastoidectomy was sometimes performed to allow draining of the pus and cleaning out of the infected area. However, removing the mastoid bone using the hand-drill-like trepan or by chipping away at the bone using chisels and gouges did not offer the level of control needed. Some patients were reported to have developed post-operative meningitis and subsequently died, an outcome which would have been due to removing too much bone and tissue. The inadequacy of the pre-powered age surgical instruments led to the mastoidectomy falling out of favor.

It should be noted that since these mastoidectomies were often performed as a life saving measure with little to no regard to the potential damage to internal ear anatomy and hearing caused by the surgical instruments. It is not a stretch, however, to assume that those patients who survived probably did have hearing loss. Eventually, Gustave Bondy developed improvements to the mastoidectomy with the aim of preserving the middle ear anatomical structures and better results were achieved. Bondy's improvements and the introduction of the electric drill resulted in the mastoidectomy becoming a mainstay of ear surgery. This situation improved even further with the introduction of antibiotics in the 1940s, which dramatically cut down on the severity of middle ear infections [19].

#### *4.2. MES Procedures*

While the introduction of an electric drill was a major breakthrough in ear surgery, an equally important breakthrough was the introduction of the microscope. This breakthrough is generally attributed to Dr. William House in the mid-1950s [31]. The combination of the widespread use of antibiotics and the microscope led to the development of new surgical possibilities and techniques under the general category of microscopic ear surgery (MES). A majority of these procedures still involve the mastoidectomy as a first step to opening up the internal anatomy of the ear to the enhanced visualization afforded by the microscope. At the same time, higher standards emerged in terms of preserving hearing and preventing hearing loss. Concerns thus began to be raised, and are still raised today, about the potential for noised-induced hearing loss due to the use of these electric drills.

Figure 4a shows a close-up view of a mastoidectomy in progress while Figure 4b shows the typical set up of the operating room and positioning of surgeons when performing an MES procedure.

(**a**) (**b**) **Figure 4.** Microscopic ear surgery (MES) (**a**) A mastoidectomy in progress. (**b**) Surgeons performing MES.

#### *4.3. TEES Procedures*

While MES is still, by far, the common surgical approach employed in ear surgery, transcanal endoscopic ear surgery (TEES) has emerged as a viable alternative in the last twenty years. The microscope is not used in TEES, but instead an endoscope is employed to access and visualize the internal anatomy of the ear through the external auditory canal instead of the more invasive mastoidectomy approach. Even though the endoscope has long been a standard part of surgery in other surgical fields, ear surgeons have had to face a unique set of anatomical and technological challenges that delayed the use of the endoscope within the ear. Endoscopes equipped with cameras were first used together with the microscope when performing mastoidectomies in the 1990s as a way to get a better view of the internal anatomy of the ear [32]. The first surgeries performed completely via the external auditory canal with the endoscope alone were reported on by Dr. Muaaz Tarabichi in 1997 [33] and once again in 1999 [34]. TEES really took off after 2008 with the further development of the 3-charged-coupled device (CCD) camera connected to high definition (HD) monitors, which resulted in high-resolution images during surgery of the tiny structures of the ear [35,36]. Figure 5a illustrates how the endoscope and forceps can be simultaneously inserted into the external auditory canal to perform TEES, and Figure 5b shows the typical set up of the operating room and positioning of surgeons when performing a TEES procedure.

**Figure 5.** Transcanal endoscopic ear surgery (TEES) (**a**) Accessing the internal anatomy of the ear with an endoscope via the external auditory canal. (**b**) A surgeon performing TEES.

TEES offers a number of advantages over MES, including no need to perform an invasive mastoidectomy which requires bone loss, better visualization of the surgical field, the ability to see into deep recesses within the ear, no disfiguring retroauricular scarring, and a quicker recovery time. Moreover, most TEES procedures are performed entirely without powered surgical instruments, which eliminates the potential for sensorineural hearing loss resulting from noise levels or vibrations. However, the "conventional" or what we have unofficially dubbed "non-powered" TEES can be performed only so far into the middle ear, but a subset of TEES procedures use the Sonopet® Omni to remove bone within the middle ear and enable access to the antrum which we call powered TEES. Figure 6 illustrates the internal anatomy of the middle ear (Figure 6a); the scope of the indications for non-powered TEES which can only reach into the inferior portion of the attic (Figure 6b); and the scope of the indications for powered TEES that can reach into the antrum (Figure 6c). This figure underscores that the Sonopet® Omni is reaching deep into the middle ear and closer to the ossicular chain and cochlea, which raises the specter of sensorineural hearing loss resulting from noise levels or vibrations, and is addressed below.

**Figure 6.** Revised system for staging and classifying of middle ear cholesteatomas to be treated by non-TEES based on the PTAM System for Staging and Classification of Middle Ear Cholesteatomas as proposed by the Japan Otological Society [37] (**a**) The middle ear (P: protympanum; T: tympanic cavity; A: attic; An: antrum; M: mastoid cells). (**b**) Non-powered TEES (orange overlay). (**c**) Powered TEES (yellow overlay).

#### **5. Generation of Noise and Vibrations by Otological Surgical Instruments and Hearing Loss**

#### *5.1. Pre-MES Procedures*

One can reliably posit that the surgical instruments and procedures used in the pre-MES period were so crude that sensorineural hearing loss was common, but no data is available.

#### *5.2. MES Procedures*

#### 5.2.1. Drill Generated Noise Levels

Many researchers have conducted basic research targeting the problem of measuring the noise levels generated by surgical drills in a non-clinical setting. Among the earliest studies of drill generated noise levels and most frequently cited is that of Kylén and Arlinger [5], published in 1976. They measured vibrations generated by drills using isolated temporal bones and cadavers, whereafter the data was then converted into equivalent air-borne noise levels. They found that the isolated temporal bones produced lower noise levels than cadavers, and thus concluded that cadavers better simulated real surgical conditions. The maximum equivalent air-borne noise levels were found to be around 100 dB, which falls below the 130 to 140 dB threshold for causing permanent hearing loss from an acute one-time exposure. The same group reported additional results in 1977 [6], looking at different variables focusing on the size of the burr; the type of the burr; (sizes: 2-mm, 4-mm and 6-mm burrs; types: diamond versus cutting burrs); drill rotation speed; and location of drilling. They found that the smaller the burr, the lower the noise levels, as well as that lower noise levels were obtained with diamond burrs, as opposed to cutting burrs with the highest noise level of 108 dB discovered with a 6-mm cutting burr. They further found that drill rotation speed and location of drilling had little effect on noise levels.

Several researchers subsequently collected intraoperative in vivo data with the objective of getting a better picture of drill generated noise during real-life conditions. Holmquist et al. [38] reported higher in vivo noise levels in 1979 from the contralateral ear (the ear which was not being operated on) of six patients while undergoing a mastoidectomy on the ipsilateral ear (the ear which was being operated on). They reported noise levels in excess of the 108 dB recorded by the Kylén group at 116 dBA for 8-mm burrs; 109 dBA for 4-mm burrs; and surprisingly, a dangerously high 125 dBA for 2-mm burrs. Hickey and Fitzgerald O'Connor [39] conducted an in vivo study in 1991 in which they attempted to directly measure the drill-generated noise levels by monitoring intraoperative *electrocochleographic* (ECoG) noise levels and calculating these levels by using a masking technique. They were able to determine that drill-generated noise levels were present in excess of 90 dBHL at the level of the cochlea, but were unable to determine peak values.

Other researchers collected data on postoperative sensorineural hearing loss in the contralateral ear in order to attempt to eliminate surgical trauma as a possible cause of hearing loss from surgery, while at the same time acknowledging the added distance between the contralateral and ipsilateral ears. Man and Winerman [40] conducted a study on 62 patients and reported in 1985 that they found a minimal difference between the peak noise levels which did not exceed 84 dB in the ipsilateral ear and 82 dB in the contralateral ear. Moreover, they found no hearing loss in the contralateral ear, but did find loss in 16 out of the 62 patients in the ipsilateral ear. They suggested that these results indicated that drill-generated noise levels did not cause sensorineural hearing loss. daCruz et al. [41] reported on drill-induced hearing loss in the contralateral ear in 1997 based on determining outer hair cell (OHC) function using distortion-product otoacoustic emissions (DPOAE). They found a change in the amplitude of the intraoperative DPOAE in 2 out of 12 patients undergoing temporal bone surgery, indicating a transient OHC dysfunction that subsequently returned to normal. This transient but reversed dysfunction was attributed to drill-generated noise levels. Goyal et al. [42] reported on the effect of mastoid bone drilling on the contralateral ear in 2013 based on otoacoustic emissions (OAE). Their study looked at the results for 30 patients and they stated that 15 of these patients exhibited a reduction in postoperative OAE levels out of which only 10 out of 15 completely recovered. However, data was only collected for up to 72 h, which makes a definite conclusion that the hearing loss was permanent a bit premature. In contrast, in two similar studies, Baradaranfar et al. in 2015 [43] and Latheef et al. [44] in 2018 reported transient hearing losses that all had disappeared by 72 h, while Badarudeen et al. [45] in 2018 reported transient hearing losses that were still present on the 7th day after surgery.

The above findings thus paint a mixed picture on transient hearing loss in the contralateral ear caused by drilling.

#### 5.2.2. Drill Generated Vibrations

Recent studies which have examined skull vibrations have indicated that this factor should not be discounted. Suits et al. reported on a guinea pig model in 1993 which was used to measure both noise and vibration levels based on the auditory brainstem response (ABR). They did find that a temporary threshold shift occurred, but that the shift had disappeared by three weeks [13]. In 2001, Zou et al. also created a guinea pig model and compared the results for younger animals versus older animals when they were exposed to both noise and noise + vibrations. They reported that the older animals were more vulnerable to a threshold shift and that the sound-induced threshold shift was significantly less than the vibration + sound-induced threshold shift at three days after exposure [14]. The same group of researchers reported in 2007 that temporal bone vibration in a guinea pig model showed that vibrations at high frequencies caused more severe hearing loss than at lower frequencies, but that the threshold shift was generally temporary [10]. However, Hilmi et al. contended that high speed drills do not produce sufficient high levels of high frequency skull vibrations to result in damage to hearing [8].

#### 5.2.3. Hearing Loss after MES Procedures

A commonly accepted range for the incidence of sensorineural hearing loss in the ipsilateral ear after ear surgery is from 1.2% to 4.5% of patients. The higher figure of 4.5% is from a study of 1680 ear surgeries reported by Palva et al. [46] in 1973, and the lower figure of 1.2% is from a study of 2,303 ear surgeries performed from 1965 to 1980, reported by Tos et al. [47] in 1984. What is notable about these papers is that they are among the first large scale studies published in the literature and that they are from more than close to 35 years ago. Both authors attributed some of these hearing losses to the surgeon being too aggressive in the area of the ossicular chain, and Tos et al. [47] in particular, noted that the incidence of sensorineural hearing loss was lower in patients treated from 1974 to 1980. They attributed this drop to better technique and better drills. In 1989, Doménech et al. [48] reported what they characterize as an important sensorineural hearing loss after tympanoplasty in

a larger percentage of patients at 16.7% than previously reported in literature, but the hearing loss was restricted to over 8000 Hz, which is typically not measured. Urquhart et al. [7] reported in 1992 that they found no evidence of an even temporary threshold shift after ear surgery in a patient group of 40; however, they only tested up to 4000 Hz. A recent study by Kent et al. [49] in 2017 looked at factors of the experience of the surgeon and the use of a powered drill in hearing outcomes for patients undergoing a tympanoplasty which required drilling in the ear canal, a less invasive procedure than a mastoidectomy, and they found that neither factor exhibited a correlation with high-frequency hearing loss. In contrast, Al Anazy et al. [50] found in 2016 that the experience of the surgeon was a significant factor, but the use of a drill was not significant in the incidence of postoperative sensorineural hearing loss between 500 Hz to 4000 Hz after tympanoplasty which was 7% for residents, but only 1% for more experienced surgeons. However, they could not identify any obvious errors on the part of the residents.

Thus, the literature as it relates to the role of powered surgical instruments, specifically electric drills in postoperative sensorineural hearing loss, is inconclusive at best, and depends on the study design and definition of postoperative sensorineural hearing loss.

#### *5.3. TEES Procedures*

While a considerable amount of study has been done on noise-induced and vibration-induced hearing loss caused by surgical drills when performing MES, ultrasonic devices have not been studied in similar depth. CUSA devices are not used in ear surgery, and a search of the literature did not reveal any related studies from outside of ear surgery. Kramer et al. [9] look at the potential for noise trauma caused by piezoelectric devices in craniofacial osteotomies and concluded that piezoelectric devices offer no advantage over regular drills in acoustic properties. They ultimately recommended using a drill because of the slower speed of the piezoelectric device. Research on TEES and the Sonopet® Omni has, to our knowledge, only been conducted by our own group. The reason for this difference is that, as stated above, TEES procedures do not usually require the use of electric drills and only the small, but important subset of powered TEES, employ an ultrasonic aspirator together with a standard drill. The Sonopet® Omni is inserted via the external auditory canal and used to remove part of the canal wall to expose the antrum while the surgical drill is used to polish a bony shelf which is designed to protect the facial nerve.

We took it upon ourselves to collect data on noise levels and vibration levels produced by the Sonopet® Omni and compare it to data collected for standard surgical drills. The data on noise levels is presented herein for the first time while the data on skull vibrations was previously published in 2013 [15].

#### 5.3.1. Ultrasonic Aspirator Generated Noise Levels

Our study was designed to confirm that the noise levels generated by an ultrasonic aspirator during powered TEES fall within safe levels and should not induce sensorineural hearing loss. The study was conducted from September 2014 to February 2015 and data was collected during surgery from patients undergoing a powered TEES procedure to remove a cholesteatoma with a total of 14 patients (5 males/9 females) ranging in age from 15 to 76 and a median age of 50. All patients underwent a transcanal atticoantrotomy which was performed using a Sonopet® Omni ultrasonic aspirator (model UST-2001 Ultrasonic Surgical Aspirator, Stryker Corporation, Kalamazoo, Michigan, USA) in the LT (longitudinal-torsion) mode at 25 kHz (Figure 3a) and a high-speed drill with a curved burr at 80,000 rpm/1333 Hz (Visao® High-Speed Otologic Drill, Medtronic, Shinagawa, Tokyo, Japan) (Figure 7).

**Figure 7.** The high-speed drill with a curved burr at 80,000 rpm/1333 Hz (Visao® High-Speed Otologic Drill, Medtronic, Shinagawa, Tokyo, Japan) used for polishing bone during the transcanal atticoantrotomy.

The noise levels were measured at 10 cm (Figure 8a and 70 cm (Figure 8b) from the devices from 0.5 kHz to 16 kHz. These distances were selected because 10 cm is the closest that it was physically possible to measure noise levels generated within the external auditory canal and 70 cm represents the distance to the surgeon's ear. The noise level at 70 cm was measured in our original study because when new powered instruments are introduced, it is important to take into consideration the potential for hearing loss for anyone in the operating room due to long term noise exposure. Furthermore, we only measured up to 16 kHz despite the Sonopet® LT mode generating frequencies of up to 25 kHz, because standard noise measurement equipment can only measure up to 16 kHz, which is also close to the maximum auditory threshold of a normal adult. The results for both the Sonopet® Omni and the Visao® drill at 10 cm were at 85 dB or below with the noise levels lower for the Sonopet® Omni up to 2000 Hz and higher for the Sonopet® Omni from 4000 Hz to 16,000 Hz. The results for both the Sonopet® Omni and the Visao® drill at 70 cm were below 70 dB with the noise levels for the Sonopet® Omni lower or equal to the Visao® drill up to 8000 Hz and higher for the Sonopet® Omni from 8000 Hz to 16,000 Hz.

**Figure 8.** (**a**) Air-conducted noise levels measured at 10 cm from the tip of the tool; (**b**) Air-conducted noise levels measured at 70 cm from the tip of the tool.

5.3.2. Ultrasonic Aspirator Generated Vibration Levels

Our group conducted a study designed to determine the vibration levels generated by an ultrasonic aspirator and compare the ultrasonic aspirator vibration levels to those of two surgical drills: an Osteon® drill at 20,000 rpm/333 Hz (Zimmer Biomet, Warsaw, Indiana, USA) and a Visao® High-speed Otologic Drill at 40,000 rpm/667 Hz and 80,000 rpm/1333 Hz (Medtronic, Shinagawa, Tokyo, Japan) [15]. All measurements were taken during an MES mastoidectomy and the skull bone vibrations were measured with a polyvinylidene difluoride (PVDF) film taped to the forehead as shown in Figure 9. PVDF is a piezoelectric material and the charge builds up in the PVDF film in response to any applied mechanical stress.

**Figure 9.** The setup for measuring skull vibrations using PolyVinylidene DiFluoride film.

Figure 10 shows the mean values of the measured skull vibrations and the background noise level at four frequency bands. In the frequency bands of 500–2000 Hz and 2000–8000 Hz, the mean values of the Sonopet® Omni with an LT-vibration tip did not exceed the values for the Visao® revolution speed of 40k rpm or 80k rpm as well as the Osteon® drill. The peak values of skull vibrations by the Sonopet® Omni with an LT-vibration tip were lower than the vibrations of Visao® at 40 k rpm in the band of 500–2000 Hz; those of Visao® at 80k rpm in the bands of 500–2000 Hz and 2000–8000 Hz; and those of Osteon® in the bands of 500–2000 Hz and 2000–8000 Hz. No significant differences in the skull vibrations were observed among the three instruments below 500 Hz or above 8000 Hz.

**Figure 10.** (**a**) Mean values for vibration levels. (**b**) Peak values for vibration levels.

Figure 11 shows the power spectrum produced by the Sonopet® at 25 kHz in LT-mode versus background noise for the purpose of reference.

**Figure 11.** Comparison of the power spectrum produced by the Sonopet® at 25 kHz in LT mode versus background noise.

#### 5.3.3. Hearing loss and TEES procedures

Our group has performed powered TEES since approximately 2011, and we have yet to record any postoperative sensorineural hearing loss in any patient up to a frequency of 8000 Hz. Thus, powered TEES can be performed, when indicated, without the worry of postoperative sensorineural hearing loss up to the highest measurable frequency.

#### **6. Conclusions**

The early days of ear surgery often focused on saving the lives of patients with little regard to postoperative hearing loss. However, the end of the 19th century and into the 20th century saw the introduction and steady improvement of surgical approaches and tools, particularly powered surgical instruments in the form of electric drills for MES and, more recently, ultrasonic aspirators for powered TEES. Concerns were raised about the impact that these powered surgical instruments could have on sensorineural hearing loss, and research has been done looking at both noise-induced hearing loss and vibration-induced hearing loss caused by powered surgical instruments when performing MES. The results of this research are long and often contradictory, as related to conventional MES and postoperative sensorineural hearing loss, particularly as related to high frequency hearing loss. Thus, effort and research should continue to be expended toward improving both MES techniques and electric drill noise level specifications, because MES, as stated above, will continue to be a standard and essential part of ear surgery, particularly in the areas of the mastoid air cells, inner ear, and skull base.

The introduction of TEES in the late 20th century eliminated the issue of sensorineural hearing loss for a circumscribed set of middle ear procedures because non-powered TEES does not require the use of powered surgical instruments. However, the emergence of powered TEES in the 21th century, once again, required that the research be conducted anew related to the use of powered surgical instruments, specifically the ultrasonic aspirator and curved burr, directly within the confines of the EAC and postoperative sensorineural hearing loss. We presented our research herein on whether powered TEES could be shown to cause either noise-induced hearing loss or vibration-induced hearing loss, and we found that powered TEES is as safe as, if not safer, in regard to the potential for postoperative sensorineural hearing loss than MES based on the data which we collected on noise levels, vibration levels, and the occurrence of postoperative sensorineural hearing loss for powered TEES. Thus, while MES will continue to be an essential part of ear surgery, surgeons can now rest assured that, when indicated, powered TEES can be performed safely and does not present any more of a risk of postoperative sensorineural hearing loss than MES.

**Author Contributions:** Conceptualization, T.I.; Investigation, T.I., T.K., T.F., H.M. and K.F.; Project administration, S.K.; Writing—original draft, M.H.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18