Sound field prediction of ultrasonic lithotripsy in water with spheroidal beam equations*
Zhang Lue, Wang Xiang-Da, Liu Xiao-Zhou, Gong Xiu-Fen
Key Laboratory of Modern Acoustics, Institute of Acoustics, Nanjing University, Nanjing 210093, China

Corresponding author. E-mail: xzliu@nju.edu.cn

Project supported by the National Basic Research Program of China (Grant Nos. 2012CB921504 and 2011CB707902), the National Natural Science Foundation of China (Grant No. 11274166), the State Key Laboratory of Acoustics, Chinese Academy of Sciences (Grant No. SKLA201401), and the China Postdoctoral Science Foundation (Grant No. 2013M531313).

Abstract

With converged shock wave, extracorporeal shock wave lithotripsy (ESWL) has become a preferable way to crush human calculi because of its advantages of efficiency and non-intrusion. Nonlinear spheroidal beam equations (SBE) are employed to illustrate the acoustic wave propagation for transducers with a wide aperture angle. To predict the acoustic field distribution precisely, boundary conditions are obtained for the SBE model of the monochromatic wave when the source is located on the focus of an ESWL transducer. Numerical results of the monochromatic wave propagation in water are analyzed and the influences of half-angle, fundamental frequency, and initial pressure are investigated. According to our results, with optimization of these factors, the pressure focal gain of ESWL can be enhanced and the effectiveness of treatment can be improved.

Keyword: 43.35.+d; 43.25.+y; spheroidal beam equation; extracorporeal shock wave lithotripsy; transducer with wide aperture angle
1. Introduction

As an effective and non-intrusive medical technique proposed in the last century, extracorporeal shock wave lithotripsy (ESWL) is applied worldwide for the treatment of calculi in the human body, such as kidney stones.[1] This technique has been used to pursue the complete damage of calculi with less side effects. For ESWL, the shock wave acts on calculi and causes extrusion, fatigue, and breaking directly while the cavitation effect caused by negative pressure helps facilitate the breakage of calculi.[2] Since the original Dornier HM-3 was introduced as a lithotripter in clinical applications in the 1980s, continuous development has been taken on the systems to produce a lithotripter field with optimized peak pressure and beam size.[1] According to Jian et al., in order to comminute a kidney stone in clinic treatments, the peak pressure at focus was set from 40 MPa to 120 MPa.[3] Nevertheless, the sizes, depths, and tensile strengths of kidney stones vary in different cases, so a comprehensive understanding of the working mechanisms of ESWL is still necessary. In order to avoid harm to the normal biological tissue, selective action on target is required and the lesions within the focal area need to be controlled precisely by predicting the acoustic field distribution of ESWL transducers.

Although Khokhlov– Zabolotaskaya– Kuznetsov (KZK) equation has been employed to predict the acoustic fields of transducers with numerical approaches, the half aperture angle of the transducer has to be limited within 16.6° to avoid large errors due to parabola approximation.[4, 5] The transducers with wide aperture angles have performed with remarkable superiority in acoustic focusing effect, focusing area, and transmissivity, therefore it is better to use transducers with wide aperture angles in ESWL treatment.

As far as nonlinearity and its influence on waveform distortion are concerned for high-frequency focal systems with wide aperture angles, Tjø tta et al. proposed an equation that can hold for 2nd Mach number to modeling acoustic propagation in fluids including viscous and thermal losses.[4] The spheroidal beam equation (SBE) is parabolic, which can precisely predict the propagation of a high-frequency focused beam emitted from a circular aperture.[6] Considering the nonlinearity, diffraction, and absorption effects, SBE model can effectively predict the distribution of the focusing acoustic field for a transducer with a wide aperture angle.[7, 8] In this paper, the acoustic field distribution of ESWL transducer is analyzed based on the SBE model.

2. Nonlinear SBE model

In this section, nonlinear SBE model is used in the ellipsoidal coordinate to illustrate the acoustic beam that propagates in water, the boundary conditions are obtained by equivalent transformation when the source is located at the focus of an ESWL transducer, and the frequency domain algorithm is employed to solve the SBE model on the basis of Fourier expansion.

The ellipsoidal coordinate system (σ , η , φ ) is shown in Fig. 1. In the xz plane (or the yz plane), the definition domain of σ is (− ∞ , + ∞ ), the contour curves of its absolute values are ellipsoids; the definition domain of η is [0, 1], its contour curves are hyperbolas with common focuses f1 and f2 of interfocal length 2b0; φ denotes the angle that these contour curves revolve around the z axis with a range of [0, 2π ]. For − ∞ < σ < ∞ , 0 ≤ η ≤ 1, 0 ≤ ϕ < 2π ; , the transformation between ellipsoidal coordinate (σ , η , φ ) and rectangular coordinate (x, y, z) is

Fig. 1. The ellipsoidal coordinate.

For an elliptical transducer, half-angle α 0 is determined by tan α 0 = af/df, where af and df denote the circular aperture radius and the distance from transducer incision surface to its right focus f2 respectively. As shown in Fig. 2, the propagation zone can be divided into two parts: the area near the reflecting surface with σ < σ 0 < 0 and the area near the focal point with σ 0σ . The wave equation for each area is derived by defining an observing point P whose distance to the origin O is.

Fig. 2. A focused ESWL transducer in σ θ (η ) coordinate plane.

In order to obtain the nonlinear SBE model, the following assumptions are made for simplification: 1) the particles in fluid near source vibrate along the radial direction from the origin; and 2) the acoustic impedance is related to acoustic pressure and particle velocity. Thereafter, the simplified 2nd-order nonlinear wave equation is then obtained, known as the Westervelt equation[9]

where p, c0, and t denote sound pressure, sound velocity, and time respectively; ρ 0, δ , and β denote density, nonlinear coefficient, and diffusion coefficient of thermal viscous medium respectively.

To ensure the accuracy of calculation of Eq. (2), the aperture angles should be limited to no more than 30° , otherwise, the influence of the nonlinearity on harmonic excitation would become significant and the above simplification of Eq. (2) cannot hold any longer.[8]

For the area of σ < σ 0 < 0 where the spherical wave is dominant, by defining delay time ts = t + R/c0, the wave propagation is expressed by[8]

where ɛ = 1/2kb0, p = p/p0, and τ s = ω ts when k, ω , and p0 denote the wave number, angular frequency, and the pressure amplitude of the source at the left focus f1; the angular variable θ is determined by η = cos θ ; is acoustic attenuation coefficient and is the discontinuous distance of plane wave. Since the sound pressure of transducers discussed in this paper is symmetric, the coordinate φ does not appear in Eq. (3).

For the area of σ σ 0 and σ 0< 0 where the planar wave is dominant, the delay time is defined as tp = tb0σ η /c0. With similar simplification of Eq. (3), in paraxial region, the wave propagation is expressed by[8]

where τ p = ω tp.

The coupling condition for Eqs. (3) and (4) at σ = σ 0 is

Equations (3)– (5) are the expressions of nonlinear SBE model, whose solution can be obtained in frequency domain by Fourier series expansion.

For a focused ESWL transducer in σ θ (η ) coordinate plane, the boundary conditions need to be determined to solve the SBE equations. In Fig. 2, with a, b, and c denoting the semi-major length, semi-minor length, and semi focal of the ellipse, the distance from transducer focus f2 to its vertex is d= a+ c. As shown in Fig. 3, for σ max = d/b0, the source of acoustic wave locates at f1, the vertex of reflecting surface locates at σ = − σ max and the focal point f2 falls at the origin. The surface of σ = − σ max is not identical to the real surface of transducer reflection.

Fig. 3. Equivalent transformation of the source.

In order to focus at f2, the monochromatic wave radiated from f1 is equivalent to be radiated from the surface σ = − σ max. For a spherical wave sourced from f1, its phase when reaching the surface of the ellipsoid transducer can be expressed as

where is the distance from f1 to a point Ps on the surface of the transducer.

For the propagation of a spherical wave, the normalized acoustic field distribution of the transducer surface can be obtained as

where dop = ac is the distance from transducer vertex to the left focal point.

The sound beams emitting from f1 focus at f2 after the reflection of the transducer. Since the transducer surface is rigid, the reflection will cause a phase change of constant π , which does not affect the solution and will be neglected in the following discussion of this paper.

Based on the phase compensation, the phase of Pσ max, the equivalent point of Ps on the surface σ = − σ max, can be expressed as

where is the distance from right focal point f2 to Pσ max, is the distance from right focal point f2 to Ps. The points Ps, Pσ max, and O are collinear.

While equivalently transformed from the transducer surface to the surface of σ = − σ max, the distribution of acoustic pressure amplitude remains the same; the difference lies in the phase compensating term of kΔ R0 in Eq. (8).

With the above analysis, after the equivalent transformation between the transducer surface and the surface of σ = − σ max, the normalized distribution of acoustic field can be expressed as

Since e− j2ka is constant, equation (9) can be further simplified to

The boundary condition for SBE model of an ellipsoid transducer is

where θ 0 is the angle in ellipsoidal coordinate after the equivalent transformation of half-angle α 0 to the surface of σ = − σ max.

For a transducer on ellipsoid (z+ c)2/a2+ x2/b2 = 1 in the xz plane of Fig. 1, according to Eq. (1), the following holds for a point (− df, af) on the transducer edge

Therefore, the angle θ 0 can be obtained as

On the other hand, for a point Ps (z, x) on the surface of transducer in the xz plane, as shown in Fig. 3, with coordinate transformation, the following can be obtained from Eq. (12):

with η = cos θ , σ , and η can be obtained by Eq. (14) for a given θ . Consequently, is determined by

The 2nd-order nonlinear term in the model leads to distortions, thus the distorted waveform can be expanded in Fourier series as[8]

where Fourier expansion coefficients , , , and are the functions of σ and θ , which represent the amplitude components of the fundamental and harmonic waves.

Combining Eqs. (3), (4), and (16), a pair of coupled nonlinear partial differential equations can be obtained. To solve the nonlinear SBE model accurately with the finite difference method, the expansion order should be high enough and the step sizes of σ and θ should be involved. The integration step for θ is taken as 0.2° , the step size for σ is proportionate to 1+ σ 2, and σ 0 = − 1.

The boundary conditions for fundamental wave of both linear and nonlinear SBE models are identical, while the harmonics are zero at the boundary of σ = − σ max. The relations for fundamental and harmonics at σ = σ 0 are both determined by

where . Pressure amplitudes for fundamental and harmonics are obtained by taking the square root of the square summation of corresponding Fourier expansion coefficients.

To derive nonlinear SBE model with finite difference method, the following expressions can be obtained by disregarding the superscripts of Fourier expansion coefficients:

Submitting Eq. (18) into Eqs. (3) and (4), the expansions for nonlinear spheroidal beam equations are

where x3 = σ (1 + σ 2), , E = (σ 2 + cos2θ )/(1 + σ 2), and α represents the absorption. The absorption coefficient for an acoustic wave propagating in a medium is

where α * is the absorption coefficient with the reference frequency f* .

Combining the differential equations of Eq. (20) with the boundary conditions for fundamental and harmonics, the numerical solution of nonlinear SBE model can be obtained by the finite difference method.

3. Acoustic field distribution of ESWL transducer in water

In this section, analyses are carried out for an ESWL transducer of Dornier HM-3 type on its acoustic-field distribution. The length of ellipsoid semi-major axis is a = 14 cm and the length of semi-minor axis is b = 8 cm. Necessary parameters for water are: density ρ 0 = 1000 kg/m3, acoustic velocity c0 = 1500 m/s, nonlinear coefficients β = 3.5, μ = 2, acoustic attenuation coefficient of 1 MHz is α * = 0.0253 Np/m.

3.1. Comparison of linear and nonlinear SBE models

With initial pressure p0 = 0.5 MPa, source exciting frequency f0 = 0.5 MHz, transducer half-angle α 0 = 20° and corresponding σ max = 2.0, since σ max = d/b0, θ 0 = 17° can be obtained according to Eqs. (12) and (13), θ 0 is slightly smaller than α 0.

Since the linear SBE model is a particular case of the nonlinear model by setting δ and β as zero, [4] the axial pressure distributions for the fundamental wave of linear and nonlinear models are shown in Fig. 4, where ξ = z/d is the normalized axial distance, pf = p = p/p0 is the normalized axial pressure, i.e., the pressure gain. The trends of fundamental pressure amplitude for both linear and nonlinear fundamental propagations are alike. However, while approaching the focus, the dashed curve that considers nonlinearity and absorption exhibits a significant deduction with respect to the linear case. At the focal point, the focusing gain of linear and nonlinear propagations are 15.0 and 13.2 respectively. This deduction is caused by not only the energy loss of acoustic absorption but also the energy transduction of nonlinearity, which transfers the energy from fundamental to harmonics. Meanwhile, for both linear and nonlinear models, the oscillations before the focal point are severe while the oscillations after the focal point are relatively steady.

Fig. 4. Axial pressure distributions for the fundamental wave of linear and nonlinear models.

Figure 5 shows the axial pressure distributions of the fundamental and the harmonic waves. It can be noticed that the focal gain of the 2nd and 3rd harmonics at the focal point are about 50% and 34% of the fundamental respectively. In addition to the satisfying focal gain, the focal areas of the harmonics are smaller and there are no oscillations beside the focus, therefore, the resolution of lithotripsy positioning can be improved and the biological tissue around the calculi will not be harmed when the harmonics are used for treatment. In practice, the harmonics can be extracted for optimized lithotripsy with the fundamental wave filtered out. Although the energy of fundamental is lost, the effect of treatment and the safety of patients can be ensured. Besides, figure 5 also indicates that the focal point is moving to the right with regard to the increase of the harmonic order.

Fig. 5. Axial pressure distributions of the 1st, 2nd, and 3rd harmonics.

Figure 6 is the pressure waveform at the focal point f2 for the linear and nonlinear models, which implies that the time domain waveform of linear propagation remains sinusoidal while the waveform of nonlinear propagation distorts severely, since the principles of lithotripsy are the fatigue crack of calculi caused by forces and the cavitation effect caused by negative pressure during nonlinear propagation. In Fig. 6, during nonlinear propagation, the pressure at focal point increases quickly to a maximal value, which can lead to a stronger pressure impact on the calculi more than the slow action in linear propagation; the pressure falls to negative rapidly after the maximal value, this instantaneous drop of pressure causes remarkable damage to the calculi. Meanwhile, the negative pressure zone of nonlinear propagation is wider and smoother than that of the linear propagation, which indicates that the nonlinear cavitation effect acts on the target longer and cracking the stone can be promoted. In fact, the effect of acoustic nonlinearity on the focused beam is notable in lithotripsy applications. The existence of shock fronts induced by nonlinearity has been experimentally confirmed in real media.[10] Therefore, the involvement of nonlinearity is essential for the precise prediction of the ESWL acoustic field.

Fig. 6. Pressure waveform at focal point f2 for linear and nonlinear models.

3.2. Influence of half-angle α 0

With initial pressure p0= 0.5 MPa, source exciting frequency f0= 0.5 MHz, transducer half-angle α 0 are taken 10° , 20° , and 30° . Corresponding σ max and θ 0 are listed in Table 1.

Table 1. Values of σ max and θ 0 with respect to α 0.

The axial pressure distributions for the fundamental wave with different α 0 are shown in Fig. 7. According to the results, the focal point of the fundamental is moving to the right while α 0 is increasing. Meanwhile, the reflecting surface is enlarged with the increase of α 0, which leads to the growth of reflected energy and the focal gain at focus.

Fig. 7. Axial pressure distributions for the fundamental wave with different α 0.

In order to reveal the influence of half-angle α 0 on the energy transfer of harmonics, the axial pressure distributions of 1st, 2nd, and 3rd harmonics with α 0 = 10° , 20° , and 30° are shown in Fig. 8. The maximum pressure gains of 2nd and 3rd harmonics are 45.6% and 29% of the fundamental when α 0 = 10° . For α 0 = 20° and 30° , the ratios are relatively stable at about 50% and 33.4% respectively. In addition, with the growth of α 0, the focusing areas are narrowed down, which is essential to avoid the damage of biological tissue nearby. Hence, the half-angle of ESWL transducer should be chosen relatively large in practice.

Fig. 8. Axial pressure distributions of 1st, 2nd, and 3rd harmonics with α 0 = 10° (a), 20° (b), and 30° (c).

With different values of α 0, the pressure waveform at the focal point f2 are shown in Fig. 9. Both the impact caused by positive pressure and cavitation effect caused by negative pressure are enhanced with the growth of α 0, which can be used to facilitate the treatment of lithotripsy. Since tan α 0 = af/df, for α 0 = 10° , 20° , and 30° , corresponding aperture radius is 2.4 cm, 5.0 cm, and 7.9 cm respectively when df remains the same. Therefore, the transducer geometric parameter can be optimized in order to obtain satisfactory treatment.

Fig. 9. Pressure waveform at the focal point f2 with different α 0.

3.3. Influence of fundamental frequency f0

With initial pressure p0 = 0.5 MPa, transducer half-angle α 0 = 20° , and corresponding σ max and θ 0 as listed in Table 1, the fundamental frequencies f0 are taken as 0.25 MHz, 0.5 MHz, and 1 MHz.

The axial pressure distributions of fundamentals with different f0 are shown in Fig. 10. It can be noticed that while the frequency f0 increases, the focal point of fundamental is moving to the right and its focal gain grows.

Fig. 10. Axial pressure distributions of fundamentals with different f0.

The axial pressure distributions of 1st, 2nd, and 3rd harmonics with f0 = 0.25 MHz, 0.5 MHz, and 1 MHz are presented in Fig. 11. The ratios of maximum pressure gains between 2nd, 3rd harmonics and the fundamental are 26.6% and 10.41% when f0 = 0.25 MHz. The ratios become 51% and 34.4% for f0 = 0.5 MHz, and 72.7% and 69.5% for f0 = 1 MHz. Therefore, with the increase of fundamental frequency f0, the ratio of maximum pressure gain between the harmonics and the fundamental grows, which illustrates that a higher fundamental frequency can improve the efficiency of harmonics. Besides, the focal area is narrowed down with the increase of fundamental frequency.

Fig. 11. Axial pressure distributions of 1st, 2nd, and 3rd harmonics with f0 = 0.25 MHz (a), 0.5 MHz (b), and 1 MHz (c).

In Fig. 12, the time domain distributions of pressure at the focal point with f0= 0.25 MHz and 0.5 MHz are presented. Obviously, both the impacts caused by positive pressure and cavitation effect caused by negative pressure are enhanced with the growth of f0, which can be used to facilitate the treatment of lithotripsy.

Fig. 12. Pressure waveform at the focal point f2 with different fundamental frequencies.

3.4. Influence of initial pressure p0

With fundamental frequency f0= 0.5 MPa, transducer half-angle α 0= 30° , and corresponding σ max, and θ 0 listed in Table 1, the initial pressures p0 are taken as 0.25 MPa, and 0.5 MPa, corresponding axial pressure distributions of fundamentals, 2nd and 3rd harmonics are presented in Fig. 13. With the increase of p0, although the focal gain of the fundamental decreases, the focal gains of harmonics grow and the growth rate of the 3rd harmonic is larger, which indicates that the increase of initial pressure p0 helps the energy to transfer to harmonics. This trend also agrees with the experimental results that increasing p0 leads to the amplification of the nonlinearity effect with creation of the higher harmonics.[10] Moreover, the previous experiment has approved that, with increasing the pressure amplitude, the effective energy of lithotripsy can be enhanced and the degree of disintegration can be significantly improved.[11] In real treatment, the initial pressure p0 is taken as a relatively large value so that the efficiency of harmonics is improved while the damage of sidelobes to the biological tissue around the calculi is reduced. Even though the decline of focal gain is induced, with the increase of p0, sufficient energy will concentrate at the calculi for lithotripsy.

Fig. 13. Axial pressure distributions of the (a) fundamental, (b) 2nd harmonic, and (c) 3rd harmonic with different p0.

4. Conclusion

In summary, to solve the nonlinear SBE model of an ESWL transducer by Fourier expansion, equivalent transformation is taken for a source located at the transducer focal point to obtain the boundary conditions for monochromatic waves. With nonlinear SBE model, the acoustic field distributions of monochromatic waves in water are investigated for both fundamental and harmonics, the influences of half-angle α 0, fundamental frequency f0, and initial pressure p0 are analyzed subsequently. According to the numerical results, by increasing α 0, f0, and p0, the squeezing and cavitation effect acting upon targets (calculi) are stronger while the harmonic pressure focal gain of ESWL is enhanced in general. Therefore, with the optimization of these factors, the effectiveness of treatment can be improved.

Reference
1 Qin J, Simmons W N, Sankin G and Zhong P 2010 J. Acoust. Soc. Am. 127 2635 DOI:10.1121/1.3308409 [Cited within:2] [JCR: 1.646]
2 Skolarikos A, Alivizatos G and de la Rosette J 2006 Eur. Urol. 50 981 DOI:10.1016/j.eururo.2006.01.045 [Cited within:1] [JCR: 10.476]
3 Jian X Q, Li W L, Liu D S, Shi Q D and Tan Z 2008 Progress in Modern Biomedicine 8 122(in Chinese) [Cited within:1] [CJCR: 0.3122]
4 Tjøtta J N and Tjøtta S 1993 Acta. Acustica. 1 69 [Cited within:3] [JCR: 0.714]
5 Fan T B, Zhang D, Zhang Z, Ma Y and Gong X F 2008 Chin. Phys. B 17 3372 DOI:10.1088/1674-1056/17/9/038 [Cited within:1] [JCR: 1.148] [CJCR: 1.2429]
6 Kamakura T 2004 Jpn. J. Appl. Phys. 43 2808 DOI:10.1143/JJAP.43.2808 [Cited within:1] [JCR: 1.067]
7 Fan T B, Liu Z B, Zhang Z, Zhang D and Gong X F 2009 Chin. Phys. Lett. 26 084302 DOI:10.1088/0256-307X/26/8/084302 [Cited within:1] [JCR: 0.811] [CJCR: 0.4541]
8 Kamakura T, Ishiwata T and Matsuda K 2000 J. Acoust. Soc. Am. 107 3035 DOI:10.1121/1.429332 [Cited within:5] [JCR: 1.646]
9 Norton G V and Purrington R D 2009 Journal of Sound and Vibration 327 163 DOI:10.1016/j.jsv.2009.05.031 [Cited within:1] [JCR: 1.613]
10 Tavakkoli J, Cathignol D, Souchon R and Sapozhnikov O A 1998 J. Acoust. Soc. Am. 104 2061 DOI:10.1121/1.423720 [Cited within:2] [JCR: 1.646]
11 Forssmann B 2006 838 291 DOI:10.1063/1.2210364 [Cited within:1]