- Email: [email protected]
Relationship between Acoustic Power and Acoustic Radiation Force on Absorbing and Reflecting Targets for Spherically Focusing Radiators
Ultrasound in Med. & Biol., Vol. 41, No. 3, pp. 832–844, 2015 Crown Copyright Ó 2015 Published by Elsevier Inc. on behalf of World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/$ - see front matter
http://dx.doi.org/10.1016/j.ultrasmedbio.2014.09.021
d
Original Contribution RELATIONSHIP BETWEEN ACOUSTIC POWER AND ACOUSTIC RADIATION FORCE ON ABSORBING AND REFLECTING TARGETS FOR SPHERICALLY FOCUSING RADIATORS PIERRE GELAT*y and ADAM SHAW* * National Physical Laboratory, Teddington, United Kingdom; and y Department of Mechanical Engineering, University College London, London, United Kingdom (Received 20 June 2014; revised 11 September 2014; in final form 16 September 2014)
Abstract—Total acoustic output power is an important parameter required by standards for most ultrasonic medical equipment including high-intensity focused ultrasound (HIFU) systems. Radiation force balances are routinely used; however, radiation force is not strictly dependent on the ultrasound power but, rather, on the wave momentum resolved in one direction. Consequently, measurements based on radiation force become progressively less accurate as the ultrasound wave deviates further from a true plane wave. HIFU transducers can be very strongly focused with F-numbers less than one: under these conditions, the uncertainty associated with use of the radiation force method becomes very significant. International Standards IEC 61161 and IEC 62555 suggest plane-wave correction factors for unfocused transducers radiating onto an ideal absorbing target and focusing corrections for focused transducers radiating onto ideal absorbing targets and onto conical reflecting targets (IEC 61161). Previous models have relied on calculations based on the Rayleigh integral, which is not strictly correct for curved sources. In the work described here, an approach combining finite element methods with a discretization of the Helmholtz equation was developed, making it possible to model the boundary condition at the structure/fluid interface more correctly. This has been used to calculate the relationship between radiation force and total power for both absorbing and conical reflecting targets for transducers ranging from planar to an F-number of 0.5 (hemispherical) and to compare with the recommendations of IEC 61161 and IEC 62555. (E-mail: [email protected]) Crown Copyright Ó 2015 Published by Elsevier Inc. on behalf of World Federation for Ultrasound in Medicine & Biology. Key Words: Acoustic radiation force, Acoustic power measurement, Ultrasound metrology, Focused ultrasonic fields, Finite element modeling, Boundary element modeling.
localized, non-invasive tissue ablation and offers the ability to treat deep-seated tumors locally with potentially fewer side effects than more invasive treatment modalities such as surgical resection, chemotherapy and ionizing radiation. Its efficacy has been reported in the treatment of a range of cancers, including those of the prostate (Ahmed et al. 2012; Sanghvi et al. 1996), liver (Leslie et al. 2012), kidney (Ritchie et al. 2010) and breast (ter Haar and Coussios 2007). It has also received U.S. Food and Drug Administration approval for the treatment of uterine fibroids, and its safety and efficacy continue to be established (Voogt et al. 2012). Two IEC standards relevant to HIFU were published in 2013 (IEC 606012-62 2013a; IEC 60601-2-62 2013a) and one Technical Specification (IEC 62556 2014). In these documents, as well as in longer established standards for ultrasonic devices such as imaging and physiotherapy systems, determination of acoustic power in water is required (IEC 60601-2-5 [IEC 2009]; IEC 61161 [IEC 2013b]).
INTRODUCTION Medical applications of ultrasound employ both focused and unfocused transducers. Ultrasound physiotherapy typically uses unfocused transducers, whereas focused transducers are widely used in ultrasonic imaging. In imaging, focusing is generally achieved by electronic beamforming with F-numbers (i.e., ratio of focal distance to aperture width) in the range 2–10. The emergence of high-intensity focused ultrasound (HIFU) over the past decade has led to an increased use of very strongly focused sources (F-number in the range 0.7–2), many of which are focused simply by virtue of having a concave radiating surface. HIFU enables highly
Address correspondence to: Pierre Gelat, Department of Mechanical Engineering, University College London, Torrington Place, London WC1E 6BT, UK. E-mail: [email protected] 832
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
The current primary standard method for measuring acoustic power radiated by an ultrasonic transducer is through measurement of the time-averaged acoustic radiation force using a gravimetric balance (IEC 61161 [IEC 2013b] and IEC 62555 [IEC 2013c]). A target is inserted in the field, and the time-averaged axial component of the radiation force vector is measured by means of the balance. The preferred target is of the absorbing type and is generally positioned orthogonal to the transducer axis. Another common target type has a right-angled convex conical reflecting surface aligned along the same axis as the transducer. Although the conical geometry is not recommended for focused fields, this type of target has the advantage that, by reflecting the ultrasonic energy, it is not susceptible to thermal damage, which can occur in absorbing targets. Acoustic radiation force is the result of a change in the momentum flux of traveling waves when absorbed or reflected by the target. It is given by the surface integral of the Langevin radiation pressure over the area of the target. If the target is perfectly absorbing, and its area is large enough, the radiation force is a characteristic of the ultrasonic field, like the acoustic power. The radiation force may be related to the acoustic power radiated by the transducer for a number of field types. For a plane progressive wave, the radiation force on a perfectly absorbing target is equal to P/c, where P is the acoustic power, and c is the speed of sound in the propagating medium. This statement is also true if the target is a right-angled convex reflector. Although this is often a good approximation to fields radiated by planar transducers (Beissner 1984a), there is evidence that large uncertainties are introduced for some transducer geometries when calculating the acoustic power from the axial component of the radiation force vector. This is particularly true in the case of strongly focused fields, as the acoustic beam may deviate much further from a plane wave than for a planar device, leading to potentially large measurement uncertainties when using a radiation force balance (Shaw 2006). An experimental method based on measuring the change in buoyancy generated by the thermal expansion of castor oil located inside a target suspended in a water bath was proposed to overcome this (IEC 62555 2013c; Shaw 2008). This technique shows promise in circumventing the limitations associated with radiation force balances and reducing measurement uncertainties for focused transducers (Rajagopal and Shaw 2012). However, because radiation force balances are still widely in use, there remains a requirement to produce satisfactory mathematical models of the measurement process and to estimate ensuing uncertainties for a range of radiators and targets. For arbitrary 3-D fields, the relationship between power and radiation force is complex and relies on a
833
tensorial treatment of the acoustic radiation stresses (Brillouin 1925). For circular rigidly vibrating plane piston sources, radiation force calculations have been the topic of much literature (Beissner 1984a, 1984b, 1987), and the relationship between power and radiation force may be characterized analytically (Beissner 1984a). For spherically curved bowl-type radiators, various attempts were made to relate radiation force to power. Corrections to the plane wave expression based on high-frequency approximations (Beissner 1987) and the geometric ray tracing approach (Shou et al. 1998) were proposed for curved radiators for absorbing targets, as well as for reflecting conical targets (IEC 62555 2013c; IEC 61161 2013b). This approach is valid only in the limit of large ka values, where k is the acoustic wavenumber and a is the source radius, and, as such, treats diffraction in an approximate way and does not address effects of multiple scattering. The radiation force from focused sources in the presence of absorbing targets was investigated as a function of target size for various ka values and focus half-angle values using a Rayleigh integral approach (Beissner 2010). The size of targets required to capture 97%, 98% and 99% of the radiation force that would be experienced by a target of infinite diameter was determined for a range of ka values, target distances and transducer apertures. The Rayleigh integral is a mathematical formulation of the Huygens–Fresnel principle, where all points on the surface of the radiator are local sources of secondary waves. The acoustic pressure generated by the radiator at any point within the field is then a superposition of point sources with amplitudes proportional to the component of the particle velocity normal to the vibrating surface. For a planar radiator mounted in an infinite baffle, this formulation is exact. For curved transducer surfaces, hemispherical elementary waves are initially diffracted by the radiator itself, and their structure is changed. Although various formulations have been proposed in the literature (Hasegawa et al. 1987; Lucas and Muir 1982; Wu 1995), a Rayleigh integral approach is not fully valid when investigating concave radiators with low F-numbers. The latter is defined as the ratio of the focal length of the source to its diameter. Various approaches have been suggested to overcome this problem (Cathignol and Sapozhnikov 1999; Sapozhnikov and Sinilo 2002). Nevertheless, these approaches are not generalizable to arbitrary sources or targets, and it is likely that a numerical solution to the problem will be more versatile in addressing a wider range of configurations. Although they can be computationally intensive, finite element/boundary element (FE/BE) methods are suitable for modeling a wide range of problems in
834
Ultrasound in Medicine and Biology
acoustics and vibration. Unlike ray tracing and Rayleigh integral approaches, FE/BE methods are exact to within discretization error. FE/BE methods overcome the limitations of the Rayleigh integral for curved radiators in that they account for both the pressure and velocity terms in the Kirchhoff–Helmholtz integral (Pierce 1981). The Rayleigh integral accounts only for the velocity term, which strictly is valid only for planar radiators. An FE/BE approach was proposed to model the field radiated by a plane rigidly vibrating circular piston in the presence of a polyurethane aperture, and radiation force calculations were carried out for a range of configurations (Gelat et al. 2005). As part of this work, the above approach was extended to feature a range of curved spherical concave radiators rigidly vibrating in an infinite baffle in the presence of both perfectly absorbing and conical reflecting targets. The dynamics of the conical target were accounted for. This approach takes into consideration full coupling between the radiator, the propagating medium and the target, hence implementing the correct boundary condition at the interfaces between the fluid and the vibrating structures. As such, effects of scattering, diffraction and multiple reflections may be accurately modeled. Acoustic radiation force calculations were carried out from first-order acoustic pressure and particle velocities for each configuration. The acoustic power radiated by each source was also calculated and the ratio Fc/P determined so as to estimate the uncertainties introduced by plane wave and approximate focusing corrections given in IEC 61161 (IEC 2013b) and, more recently, IEC 62555 (IEC 2013c). Currently, the most established method of measuring acoustic power radiated by an ultrasonic source is by means of the time-averaged force acting on a sufficiently large target positioned in the acoustic field. Assuming plane wave propagation in the direction normal to the target, the radiation force on a perfectly absorbing target is related to the acoustic power as follows. F5
P c
(1)
For an unfocused source represented by a rigidly vibrating plane circular piston mounted in an infinite baffle, it may be shown that if the target is perfectly absorbing (Beissner 1984a), F5
P 12J20 ðkaÞ 2J21 ðkaÞ : c 12J1 ð2kaÞ=ðkaÞ
(2)
The second term on the right-hand side is often called the ‘‘plane wave correction.’’ For values of ka $35, eqn (2) deviates from the plane wave expression
Volume 41, Number 3, 2015
in eqn (1) by ,2%. As such, many ultrasonic transducers can, in practice, be characterized using eqn (1) (Beissner 1984a). The force arising from fields generated by non-planar ultrasonic transducers will in general differ from the plane wave expression. For a spherical bowl transducer radiating onto a perfectly absorbing target, a relationship based on a high-frequency asymptotic approximation may be derived (IEC 61161 [IEC 2013b]; Shou et al. 1998). F5 where g is case of a right-angled cone tip is its centre of given by
P ð11cos gÞ c 2
(3)
the transducer focus half-angle. In the spherical transducer radiating onto a convex conical reflecting target where the positioned between the transducer and curvature, the corresponding expression is
F5
P 12cos 2g12g2sin 2g c 4ð12cos gÞ
(4)
where g must be expressed in radians. The expressions in eqns (3) and (4) are of practical significance when characterizing a subset of HIFU transducers, but nevertheless rely on high-frequency approximations to describe the radiation of sound from the transducer and its effect on the target: they also assume that the behavior of the real transducer is a good approximation to a spherical piston and that the wave reflected from the target is not incident on the transducer. Hence, they are of limited value for transducers with small or moderate ka values or of uncertain construction. Furthermore, in transducers with low F-numbers, effects caused by multiple scattering off the surface of the device become increasingly important and will not be accounted for by eqn (4). Finally, effects of multiple scattering of sound between the surfaces of the transducer and the conical target are not considered in the theory underpinning the plane wave or simple focusing corrections: of course, multiple reflections should be avoided for accurate measurement, but it is relevant to examine the effect theoretically to understand properly its importance. In this article, an approach based on the discretization of the Helmholtz equation in its integral form is used, which, although more computationally involved than high-frequency approximations, is more suitable for describing the above phenomena and is likely to help assess the range of validity of eqns (3) and (4) as a function of ka and g By combining this approach with finite element methods, it has been possible to model the boundary condition at the structure/fluid interface more correctly.
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
METHOD Radiation force theory Radiation force theory is now briefly reviewed. Consider a homogeneous, non-absorbing, non-dispersive, sound propagating fluid. Assume that the primary field quantities are acoustic pressure p and particle velocity ! u and are both functions of space and time. The medium is characterized by an equilibrium density (r) and sound speed (c). Let the position vector in Eulerian coordinates be characterized by ! r ðx; y; zÞ in a Cartesian coordinate system. The linear acoustics result for the Brillouin stress, which is the acoustic radiation stress tensor (S) in Eulerian coordinates, is given by (Beissner 1998) (5) Sij 5 T1 1T2 1T3 2V dij 2rui uj ; where i, j 5 1, 2 and 3; 1, 2 and 3 correspond to the global (x, y and z) directions in a Cartesian axis set, respectively; the bars denote temporal averages; Ti 5
rjui j 4
2
835
the radiation force is independent of the target distance because of the momentum conservation law (Beissner 1998). Acoustic radiation force theory for a perfectly absorbing circular target It is assumed that the main axis of propagation is along the positive Cartesian x-axis. For a circular, planar, perfectly absorbing target whose receiving area is parallel to the y-z plane, the radiation force vector is expressed as 0 1 0 1 0 1 Sxx Sxy Sxz 21 F ! @ xA F 5 Fy 5 %A @ Syx Syy Syz A$@ 0 AdA; (7) Szx Szy Szz 0 Fz where A is the surface area of the circular target. For an axisymmetric field, where x is the axialpcoordinate and ffiffiffiffiffiffiffiffiffiffiffiffi R is the radial coordinate such that R 5 y2 1z2 , it may be shown that (Beissner 1984a) ðb Fx 5 2p V1T x 2T R R dR; (8) 0
is the time averaged kinetic acoustic energy density;
where b is the target radius.
2
V5
jpj 4rc2
is the time-averaged potential acoustic energy density; dij is the Kronecker delta (dij 5 1 if i 5 j) and (dij5 0 if i s j); r is the medium equilibrium density; ui is the instantaneous (complex) particle velocity along the i direction; and r is the instantaneous (complex) acoustic pressure. 1 0 Sxx Fx ! @ A F 5 Fy 5 %A1 @ Syx Szx Fz 0
Sxy Syy Szy
Acoustic radiation force theory for a reflecting conical target The main axis of propagation is once again assumed to be the positive Cartesian x-axis. This is also assumed to be the axis of the conical target. If the axis of the cone makes an angle a to the generatrix, the radiation force vector is given by
1 0 1 0 Sxx Sxz 2sin a Syz A$@ cos a cos q AdA1 1%A2 @ Syx Szz Szx cos a sin q
The radiation stress tensor (S) expresses the momentum flux density (Brillouin 1925), and the above expression has been found to be valid to second order when r and ! u are first-order quantities. The radiation force vector ! ( F ) is given by integrating the acoustic radiation stress tensor over a surface A enclosing the target (Beissner 1984a) 0 1 Sxx Sxy Sxz ! n dA; (6) F 5 %A @ Syx Syy Syz A$! Szx Szy Szz where ! n is an outward pointing unit vector normal to the surface A. In a lossless fluid, for an infinitely large target,
Sxy Syy Szy
1 0 1 Sxz 1 Syz A$@ 0 AdA2 ; Szz 0
(9)
where tan q 5 z/y, A1 corresponds to the area defined by the surface of the conical region and A2 is the circular area defined by the base of the cone. A radiation force balance responds to the axial component of the radiation force vector Fx. For an axisymmetric field, where x is the axial coordinate pffiffiffiffiffiffiffiffiffiffiffiffiand R is the radial coordinate such that R 5 y2 1z2 , it may be shown that by substituting the expressions for the components of the stress tensor into eqn (9) and by applying the appropriate coordinate transformation, ðh Fx 5 22p tan2 a 0
T R 2T x 2V x dx
836
Ultrasound in Medicine and Biology
ðh 22p tan a
ð R1 rux uR x dx22p
0
Volume 41, Number 3, 2015
V1T x 2T R R dR;
0
(10) where h is the height of the cone and R1 is the radius of its base and ux uR 5
1 ! ! ux ð r ÞuR ð r Þ1uR ð! r Þux ð! rÞ ; 4
(11)
r Þ is the spatial component of the axial where ux ð! component of the complex particle velocity vector, r Þ is the spatial component of the radial component uR ð ! of the complex particle velocity vector, and the asterisk denotes the complex conjugate. Numerical modeling: Perfectly absorbing target The methodology employed here is similar to that used by Gelat et al. (2005). Circular symmetry about the Cartesian x-axis was assumed for both reflecting and absorbing target configurations. The spherical bowl source was modeled as a vibrating, rigid, curved piston mounted in an infinite baffle and radiating into a semi-infinite non-absorbing fluid. Three-noded quadratic thin axisymmetric shells of revolution were used, and all motions other than translational movement along the x-axis were restrained. A unit radial velocity magnitude and zero phase were assigned to one of the piston nodes, and this was repeated throughout all the nodes on the surface of the radiator. A combination of acoustic finite element and boundary element patches were used to model the fluid. A region of pressure-based, quadratic fluid finite elements in front of the curved piston was meshed to a distance of a quarter of a wavelength from the rim using a combination of six-noded triangular and eight-noded quadrilateral elements. A boundary element was then applied at the outer edge of fluid finite element region for the perfectly absorbing target calculations, as illustrated in Figure 1. The targets were positioned two wavelengths in the fluid away from the rim of the transducer. At positions along the radial direction from the rim of the transducer and outside its radius, a zero acoustic pressure normal derivative boundary condition was applied in the plane of the transducer, extending to infinity. This boundary condition effectively simulates a rigid baffle. A sound speed of 1487 m s21 and a density of 998 kg m23 were assumed for the acoustic medium. The frequency of excitation for all configurations was 1 MHz. Thirteen source radii were investigated: 1, 2, 3, 4, 6, 8, 10, 15, 20, 25, 30, 40 and 50 mm. For each of these configurations, 11 linearly spaced values for sin g were investigated ranging from 0 (plane piston case) to 1
Fig. 1. Representative curved radiator with absorbing target configuration revealing the finite element mesh and the boundary element. Curved radiator radius 5 10 mm. sin g 5 0.6.
(hemispherical radiator) in steps of 0.1. The F-numbers corresponding to the values of sin g are displayed in Table 1. This leads to a total of 143 absorbing target configurations. An analysis of how the target radius affects the estimate of the radiation force and measured acoustic power was not carried out, as this has been already been addressed by Beissner (2010). Instead, a large enough target radius b was considered to ensure that the value of Fx was beyond 0.98 of the limit of Fx as b (Beissner 2010). This was achieved by selecting b=a 5 10 for all configurations except for the sources of 1-mm radius, where b=a 5 100 was used, as using b=a 5 10 for the smallest, most highly focused hemispherical configurations resulted in underestimation of the acoustic power on the target by approximately 30%.
Table 1. sin g values and corresponding F-numbers sin g 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 F-number N 5 2.5 1.67 1.25 1 0.83 0.71 0.625 0.56 0.5
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
In boundary element methods, the discretization of the Helmholtz integral equation leads to a fully populated asymmetric system of linear equations. It is well known that the conventional Helmholtz boundary integral equation (BIE) formulation for exterior wave problems yields non-unique solutions at the eigenfrequencies associated with the (interior) Dirichlet problem (Burton and Miller 1971). When the coefficient matrix resulting from discretization of the Helmholtz boundary integral equation is inverted at frequencies approaching these (fictitious) eigenfrequencies, the problem becomes ill-conditioned. A way in which this may be overcome is to use the combined Helmholtz integral equation formulation (CHIEF) method (Schenk 1968). This is based on the fact that at any frequency, there is a unique surface pressure function that satisfies both the interior and surface Helmholtz integral equations. Additional collocation points were specified within regions surrounding the BE. These supplement the degrees of freedom on the BE surface and result in an overdetermined set of linear equations. Five collocation points were positioned in the region interior to the boundary element to overcome problems associated with non-uniqueness. The PAFEC (Programme for Automatic Finite Element Calculations) software (PACSYS, Nottingham, UK) was used to carry out the calculations. PAFEC has the capability of dealing with the coupling between the elastic structure and the acoustic wave in a fluid by using a combined finite element/boundary element method. For quadratic elements, meshing at three elements per wavelength is normally sufficient, but six elements per wavelength were nonetheless used in all meshes. PAFEC enables the calculation of acoustic pressures and particle velocities at user-specified locations within the boundary element exterior domain. These first-order quantities were used to determine Fx for both the reflecting and absorbing targets as expressed by eqns (8) and (10), in which numerical quadrature was used to carry out the integrations. The acoustic power radiated into the fluid by the source was computed in PAFEC as the integral of the normal component of the time-averaged acoustic intensity over the surface defined by the fluid–structure interface as
837
For a plane circular piston rigidly vibrating sinusoidally in an infinite baffle, the acoustic power radiated may be evaluated analytically as (Kinsler et al. 1982)
1 2J1 ð2kaÞ ; (13) P 5 U02 pa2 rc 12 2 2ka where U0 is the maximum piston velocity. There is no equivalent exact analytical solution in the literature for curved pistons. Numerical modeling: Reflecting conical target The reflecting conical targets were modeled using a three-noded quadratic thin shell of revolution elements, but no restraints other than those required for axisymmetric motion were applied. The reflecting target was assumed to be made of a steel outer shell 0.2 mm thick enclosing an air volume (modeled as a vacuum). This is similar to radiation force targets used previously on the National Physical Laboratory Primary Standard Radiation Force Balance. The mechanical properties of the steel were: Young’s modulus 5 215 GPa, density 5 7800 kg m23 and Poisson’s ratio 5 0.293. The boundary element was extended to encompass the region surrounding the cone, as illustrated in Figure 2. The region of air normally present inside conical reflecting targets was modeled as a vacuum. Continuity of normal particle velocity was assumed to couple the structure and the fluid. The cone angle was 90 , which is what is commonly used experimentally. The radius h of the base of the cone was assumed to be twice that of the source radius a which, for the conical target calculations, was fixed at 20 mm. This value for the source radius was chosen as it was toward the middle of the range of radii investigated in the absorbing target case. Investigation of other source radii was beyond the scope of the work described in this article. To increase the degrees of freedom on the BE surface, 419 additional collocation points were used inside the cone. Eleven linearly spaced sin g values ranging from 0 to 1 were investigated. The
1 P 5 %Asource ½pð! r Þ! u ð! r Þ1p ð! r Þ! u ð! r Þ:! n dAsource ; 4 (12) wherepð! r Þ is the spatial component of the complex acoustic pressure, ! u ð! r Þ is the spatial component of the complex particle velocity vector; Asource is the transducer radiating area; ! n is the unit vector normal to the area Asource ; and the asterisk denotes the complex conjugate.
Fig. 2. Representative curved radiator with reflecting target configuration revealing the finite element mesh and the boundary element. Curved radiator radius 5 20 mm. sin g 5 0.7.
838
Ultrasound in Medicine and Biology
Volume 41, Number 3, 2015
Fig. 3. Acoustic power radiated by a 1 MHz planar transducer. Percentage difference between analytical and finite element/boundary element (FE/BE) solutions. Inset: Comparison of FE/BE approach (solid line) with the analytical solution, eqn (13) (dotted line).
apex of the conical target was positioned two wavelengths into the fluid away from the rim of the transducer. Note that in the case of sin g 5 1, the cone apex lies beyond the focus so the target does not intercept the whole field; this most extreme condition may be of limited practical relevance, but is included for completeness. RESULTS Comparison against analytical solutions To confirm that the FE/BE model is reliable, acoustic power and the ratio Fc=P can be compared with theory. Figure 3 illustrates the acoustic power calculated using the FE/BE approach against the analytical solution in eqn (13) for a 1 m s21 piston normal velocity: agreement is within 0.1% for all cases. For convenience, the normal component of the radiation force vector will be referred to as F. Figure 4 illustrates the ratio Fc=P for the FE/BE calculations on the plane circular piston case compared against the analytical solution in eqn (2). The agreement between both sets of results is within 1.1%. Because only 11 ka values were used in these calculations, small oscillations caused by the Bessel function in eqn (13) are not resolved: nevertheless, Figures 3 and 4 taken together serve to validate the FE/BE model to better than 1%. Perfectly absorbing target calculations The ratio Fc=P was computed for the 143 source and target configurations described in the previous section.
The results are summarized in Figure 5(a, b). The focusing correction factor taken from IEC 61161 (IEC 2013b) expressed in eqn (3) is overlaid for each case as horizontal dotted lines for each value of g. Figure 5(a, b) illustrates that, for larger values of the source radius a, the IEC 61161 solution for the ratio Fc=P is very similar to the FE/BE results: the percentage difference is displayed as a function of a in Figure 6. The case of sin g 5 1 has been omitted from Figure 6, as the percentage difference with the IEC 61161 expression is large (between 241% and 131%). Conical reflecting target calculations The ratio Fc=P was computed for the 11 transducer aperture angles described in the previous section. The results are illustrated in Figure 7(a) and compared against the recommended correction factor provided by IEC 61161 (IEC 2013b) for reflecting conical targets described by eqn (4). Figure 7(b) illustrates the same data for values of sin g between 0 and 0.7. Figure 8 illustrates the acoustic pressure field between the transducer and the target for the case that gives an enhanced ratio of Fc=P. In this case, the boundary condition at the surface of the source was chosen as a radial displacement of 1028 m. The corresponding conical target displacements are illustrated in Figure 9, where the source radial displacement has been scaled by a factor of 105 relative to that in Figure 8 so as to better visualize the deflections of the target. The maximum target displacement amplitude is in fact 10.7 3 1028 m.
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
839
Fig. 4. 1 MHz planar transducer radiating on a perfectly absorbing target. Comparison of relationship between radiation force axial component and acoustic power for FE/BE approach (solid line) with that for analytical solution, eqn (2) (dotted line). FE/BE 5 finite element/boundary element.
Figure 8 illustrates that multiple reflections are occurring between the radiator and the target for this focusing angle, leading to a standing wave pattern. Figure 9 illustrates that the target displays modal behavior. DISCUSSION The IEC 61161 (2013b) standard typically requires that individual sources of uncertainty contribute less
than 1% to the uncertainty of the overall measured power. In IEC 62555 (IEC 2013c), which is specific to high-intensity therapeutic ultrasound sources, this criterion is relaxed to 2%. In the case of a perfectly absorbing target, analysis of the data in Figure 4 reveals that the agreement of the FE/BE calculations is within 0.1% of the analytical solution. It can therefore be concluded that, for the range of planar radiators investigated, the radiation force
Fig. 5. 1 MHz transducer radiating on a perfectly absorbing target. Comparison of relationship between radiation force axial component and acoustic power for FE/BE approach (solid line) with that for IEC 61161 (IEC 2013b, 2013c) solution, eqn (3) (dotted line). (a) 0 # sin g # 0.5. (b) 0.6 # sin g # 1. FE/BE 5 finite element/boundary element.
840
Ultrasound in Medicine and Biology
Volume 41, Number 3, 2015
Fig. 6. 1 MHz transducer radiating on a perfectly absorbing target. Percentage difference between radiation force axial component and acoustic power using eqn (3) (IEC 61161), relative to finite element/boundary element approach and as a function of a.
normal to the target predicted by the FE/BE method is within 1% of its value for an infinite target. The discrepancies are likely due to a combination of insufficient sampling of points on the target and the use of a target of finite size. The relationship between force and power given by eqn (3) for an absorbing target is based on a Rayleigh integral calculation for a spherical section radiator with uniform radial motion. The methodology employed in this article, also applied to a spherical section radiator with uniform radial motion, overcomes the limitations of the Rayleigh integral approach by including both the monopole and dipole terms in the Kirchhoff–Helmholtz integral. At 1 MHz, the ratio Fc=P provided by the IEC 61161 standard and expressed in eqn (3), is within approximately 1% of the improved calculation for transducers with F-number $0.83 and of radius $10 mm. It is consistent to within 2% for sin g # 0.9 as long as the source radius is larger than approximately 4 mm, that is, ka $ 17. Many HIFU devices fall within these ranges. These results provide additional confidence in the use of the IEC 61161 correction factor for spherically focusing sources. The IEC 61161 results diverge from the FE/BE solution as the source becomes increasingly focused and as its radius decreases. A similar set of calculations (not included in this article) for a spherical section radiator excited in a direction parallel to the x-axis gave comparable results for this range of source parameters. In Figure 5(a), the dashed lines
represent the focusing correction factor taken from IEC 61161 expressed in eqn (3). For the unfocused case, g 5 0, the dashed line is drawn at Fc=P 5 1, and is strictly only correct in the limit of large transducer radii. For plane piston transducers of finite size, diffraction occurs, leading to deviation form plane-wave theory, and the ratio Fc=P should be less than 1. Indeed, IEC 61161 provides a formula to calculate a plane-wave correction (which is approximately 2.4% for a 5-mm-radius source at 1 MHz, decreasing to 1.1% for 10-mm-radius and 0.2% for 50-mm-radius sources). It is possible that a weak degree of focusing may compensate for diffractive spreading resulting in a higher ratio of Fc=P than in the plane piston case: there is evidence here to support this possibility. Note in Figure 5(a) that the Fc=P ratios are comparable in the plane piston case and in the case of the curved radiator featuring sin g 5 0.1. This is true for the range of piston radii considered in this article. However, as sin g increases, the effects of focusing clearly overcome those of beam diffraction. Finally, it is worth noting that the IEC 61161 and IEC 62555 standards [IEC 2013b, 2013c] do not address the question of how to combine corrections for a finite-size transducer with corrections for focusing. Although the data in Figure 5(a, b) may serve as a starting point in establishing such a correction factor, only a subset of radiation force correction factors would be represented, as limited source radii have been investigated at only one excitation frequency (1 MHz).
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
841
a
b
Fig. 7. (a) 1 MHz 20-mm-radius curved transducer radiating on a conical reflecting target of 90 aperture angle and 40-mm base radius. Comparison of relationship between radiation force axial component and acoustic power for FE/BE approach (solid line) with that for IEC 61161 (dotted line). (b) 1 MHz, 20-mm-radius curved transducer radiating on a conical reflecting target of 90 aperture angle and 40 mm base radius. Comparison of relationship between radiation force axial component and acoustic power for FE/BE approach (solid line) with that for IEC 61161 (dotted line) for 0 # sin g # 0.7. FE/BE 5 finite element/boundary element.
More extensive numerical investigations on curved radiators excited over a range of frequencies could be carried for a wider range of source radii. This may establish further sets of radiation force correction factors,
but this is beyond the scope of the work described in this article. For the conical reflecting target configurations investigated as part of this work, it can be seen in
842
Ultrasound in Medicine and Biology
Fig. 8. Acoustic pressure field in megapascals from FE/BE calculations for a 1 MHz, 20-mm-radius curved transducer radiating on conical reflecting target of 90 aperture angle and 40-base radius. sin g 5 0.8. Source radial displacement 5 1028 m. FE/BE 5 finite element/boundary element.
Figure 7b that the IEC 61161 (2013b) correction factor agrees with the FE/BE model results within 2% for values of sin g between 0 and 0.6 (which is an F-number of 0.83): so, again, the IEC formula is sufficiently accurate for most practical purposes. However, with a reflecting target, there is the possibility of interaction between the reflected beam and the transducer, which can lead to a strong modulation of the radiation force. For sin g 5 0.8 and the target distance chosen (3 mm, which corresponds to two wavelengths in water at 1 MHz), the numerical modeling calculation of Fc=P is a factor of 5 higher than the result predicted by
Fig. 9. Conical target displacement profile resulting from FE/BE calculations for a 1 MHz, 20-mm-radius curved transducer radiating on conical reflecting target of 90 aperture angle and 40-mm base radius. sin g 5 0.8. Source radial displacement 5 1028 m, resulting in a maximum displacement of the target of 10.7 3 1028 m. For visualization purposes, displacements have been scaled by a factor of 105. FE/BE 5 finite element/boundary element.
Volume 41, Number 3, 2015
the IEC 61161 expression (see Fig. 7a). Equation (4) considerably underestimates the radiation force-toacoustic power ratio for this configuration. Although a 3mm separation distance might seem rather close, it is generally recommended that output power measurements be made close to the transducer surface, and conical targets do not show such enhanced behavior when used very close to planar transducers. The value of g, at which this interaction is strongest, will depend on the distance of the target; this type of enhancement has been observed experimentally (Shaw and Hodnett 2008) where an approximate doubling in the radiation force on a conical target was reported at distances between 3–15 mm from a focused transducer with F-number 1.5. The results in Figure 7(a, b) further highlight the limitations of conical targets and why their use should be avoided, particularly when dealing with focused fields (Shaw 2008). These simulations suggest that, in addition to the focusing errors introduced, the dynamic behavior of the target together with multiple reflections between the target and the radiator has the potential to cause large uncertainties in the power measurement. CONCLUSIONS A combination of finite element and boundary element modeling was used to calculate the relationship between radiation force and total power for both absorbing and conical reflecting targets for 1 MHz transducers ranging from planar to an F-number of 0.5 (hemispherical). This method overcomes the limitations of the Rayleigh integral for curved radiators in that it accounts for both the pressure and velocity terms in the Kirchhoff–Helmholtz integral (Pierce 1981). The Rayleigh integral accounts only for the velocity term, which is strictly valid for only planar radiators. The findings were compared with the recommendations of IEC 61161 and IEC 62555 (IEC 2013b, 2013c). The modeling described here, in which the radiation force was calculated by first determining the acoustic stress tensor, is more versatile than a Rayleigh integral or analytical approach because it allows acoustic interaction between the transducer and target to occur. Furthermore, the approach can be extended to more complex transducer models, including piezo-electric effects, and to more realistic targets of non-zero reflectivity. The modeling reported here is axisymmetric and limited to a single frequency of 1 MHz: long run times on some configurations will make it difficult to generalize to three dimensions on a desktop PC, but future work could address large 3-D problems using cluster computing. Within the restricted assumption of ideal spherical section transducers operating in an infinite baffle, it
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
was found that the correction factor from IEC 61161 and IEC 62555 (IEC 2013b, 2013c) for planar absorbing targets is consistent with the results presented in this article to within 1% for transducers with F-number $0.83 (i.e., sin g 5 0.6) and of radius $10 mm; it is consistent to within 2% for sin g , 0.9 as long as the source radius is greater than about 4 mm (ka $ 17); for smaller radii, a larger correction factor is required. The same correction factor can be applied irrespective of whether the surface displacement is radial or parallel to the x-axis. Finally, it is worth noting that the IEC 61161 and IEC 62555 standards do not presently address combining radiation force correction factors for focusing transducers of finite dimensions, which could be significant for small transducers or at low frequencies. The FE/BE method employed in the work described here to calculate the radiation force of focused radiators on perfectly absorbing targets may go some way toward addressing this, although simulations across a much greater range of devices would be required. With the same restrictions, this time considering the case of a 20-mm-radius transducer (ka 5 84) and a rightangled conical reflecting target, it was found that the focusing correction from IEC 61161 is again in good agreement with the FE/BE results for sin g , 0.6. However, care should be taken not to move the transducer too closely to the target, as reflections may greatly increase the radiation force and therefore overestimate the acoustic power. To the authors’ knowledge, this is the first time a full acoustic model of the radiation force on a reflecting target has been carried out. It is planned that this work be extended to more realistic transducer models in the future. Finally it is worth noting that as these calculations were carried out at an excitation frequency of 1 MHz, care should be taken in extrapolating from these results to other frequencies. Further calculations over a wider range of frequencies relevant to HIFU devices are desirable to generalize these conclusions. Acknowledgments—The authors gratefully acknowledge funding for this work provided by the United Kingdom Department for Business, Innovation and Skills, through the Acoustics and Ionising Radiation Metrology Programme.
REFERENCES Ahmed HU, Cathcart P, Chalasani V, Williams A, McCartan N, Freeman A, Kirkham A, Allen C, Chin J, Emberton M. Whole-gland salvage high-intensity focused ultrasound therapy for localized prostate cancer recurrence after external beam radiation therapy. Cancer 2012;118:3071–3078. Beissner K. Acoustic radiation pressure in the near field. J Sound Vib 1984a;93:537–548. Beissner K. Minimum target size in radiation force measurements. J Acoust Soc Am 1984b;76:1505–1510. Beissner K. Radiation force calculations. Acustica 1987;62:255–263.
843
Beissner K. The acoustic radiation force in lossless fluids in Eulerian and Lagrangian coordinates. J Acoust Soc Am 1998;103:2321–2332. Beissner K. Minimum radiation force target size for power measurements in focused ultrasonic fields with circular symmetry. J Acoust Soc Am 2010;128:3355–3362. Brillouin L. Sur les tensions de radiation [On the radiation stresses]. Ann Phys (Paris) 1925;4:528–586. Burton AJ, Miller GF. The application of integral equation methods to the numerical solution of some exterior boundary-value problems. Proc R Soc London Ser A 1971;323:201–210. Cathignol D, Sapozhnikov OA. On the application of the Rayleigh integral to the calculation of the field of a concave focusing radiator. Acoust Phys 1999;45:735–742. Gelat PN, Zeqiri B, Hodnett M. A finite element model of the aperture method for determining the effective radiating area of a physiotherapy treatment head. Ultrasonics 2005;43:321–330. Hasegawa T, Inoue N, Matsuzawa K. A new theory for the radiation from a concave piston source. J Acoust Soc Am 1987;82:706–708. International Electrotechnical Commission (IEC). IEC 60601-2-5: Medical electrical equipment: Particular requirements for the basic safety and essential performance of ultrasonic physiotherapy equipment. 3rd ed. Geneva: 2009. International Electrotechnical Commission (IEC). IEC 60601-2-62: Medical electrical equipment: Particular requirements for the basic safety and essential performance of high intensity therapeutic ultrasound (HITU) equipment. 1st ed. Geneva: 2013a. International Electrotechnical Commission (IEC). IEC 61161: Ultrasonics—Power measurements—Radiation force balances and performance requirements. 3rd ed. Geneva: 2013b. International Electrotechnical Commission (IEC). IEC 62555: Ultrasonics—Power measurement—Output power measurement for high intensity therapeutic ultrasound (HITU) transducers and systems. 1st ed. Geneva: 2013c. International Electrotechnical Commission (IEC). IEC 62556: Ultrasonics—Surgical systems—Specification and measurement of field parameters for high intensity therapeutic ultrasound (HITU) transducers and systems. Technical Specification. Geneva, 2014. Kinsler LE, Frey AR, Coppens AB, Sanders JV. Fundamentals of acoustics. 3rd ed. New York: Wiley; 1982. p. 192. Leslie T, Ritchie R, Illing R, ter Haar G, Phillips R, Middleton M, Wu F, Cranston D. High-intensity focused ultrasound treatment of liver tumours: post-treatment MRI correlates well with intraoperative estimates of treatment volume. Br J Radiol 2012;85: 1363–1370. Lucas BG, Muir TG. The field of a focusing source. J Acoust Soc Am 1982;72:1289–1296. Pierce AD. Acoustics: An introduction to its physical principles and applications. New York: McGraw–Hill; 1981. p. 214. Rajagopal S, Shaw A. The buoyancy method—A potential new primary ultrasound power standard. Metrologia 2012;49:327–339. Ritchie RW, Leslie T, Phillips R, Wu F, Illing R, ter Haar G, Protheroe A, Cranston D. Extracorporeal high intensity focused ultrasound for renal tumours: A 3-y follow-up. BJU Int 2010;106: 1004–1009. Sanghvi NT, Fry FJ, Bihrle R, Foster RS, Syrus J, Zaitsev AV, Hennige CW. Noninvasive surgery of prostate tissue by highintensity focused ultrasound. IEEE Trans Ultrason Ferroelectr Freq Control 1996;43:1099–1110. Sapozhnikov OA, Sinilo TV. Acoustic field produced by a concave radiating surface with allowance for the diffraction. Acoust Phys 2002;48:720–727. Schenck HA. Improved integral formulation for acoustic radiation problems. J Acoust Soc Am 1968;44:41–58. Shaw A. How to measure HIFU output power properly. In: Proceedings, 5th International Symposium on Therapeutic Ultrasound (2006). New York: IEEE; 2006. p. 628–632. Shaw A. A buoyancy method for the measurement of total ultrasound power generated by HIFU transducers. Ultrasound Med Biol 2008; 34:1327–1342. Shaw A, Hodnett M. Calibration and measurement issues for therapeutic ultrasound. Ultrasonics 2008;48:234–252.
844
Ultrasound in Medicine and Biology
Shou W, Wang Y, Qian D. Calculation for radiation force of focused ultrasound and experiment measuring HIFU power. Tech Acoust 1998;17:145. Ter Haar G, Coussios C. High-intensity focused ultrasound: Past, present and future. Int J Hyperthermia 2007;23:85–87. Voogt MJ, Trillaud H, Kim YS, Mali WPTh M, Barkhausen J, Bartels LW, Deckers R, Frulio N, Rhim H, Lim HK, Eckey T, Nieminen HJ,
Volume 41, Number 3, 2015 Mougenot C, Keserci B, Soini J, Vaara T, K€ohler MO, Sokka S, Maurice AJ, van den Bosch A. Volumetric feedback ablation of uterine fibroids using magnetic resonance-guided high intensity focused ultrasound therapy. Eur Radiol 2012;22:411–417. Wu J. Calculation of acoustic radiation force generated by focused beams using the ray acoustics approach. J Acoust Soc Am 1995; 97:2747–2750.
http://dx.doi.org/10.1016/j.ultrasmedbio.2014.09.021
d
Original Contribution RELATIONSHIP BETWEEN ACOUSTIC POWER AND ACOUSTIC RADIATION FORCE ON ABSORBING AND REFLECTING TARGETS FOR SPHERICALLY FOCUSING RADIATORS PIERRE GELAT*y and ADAM SHAW* * National Physical Laboratory, Teddington, United Kingdom; and y Department of Mechanical Engineering, University College London, London, United Kingdom (Received 20 June 2014; revised 11 September 2014; in final form 16 September 2014)
Abstract—Total acoustic output power is an important parameter required by standards for most ultrasonic medical equipment including high-intensity focused ultrasound (HIFU) systems. Radiation force balances are routinely used; however, radiation force is not strictly dependent on the ultrasound power but, rather, on the wave momentum resolved in one direction. Consequently, measurements based on radiation force become progressively less accurate as the ultrasound wave deviates further from a true plane wave. HIFU transducers can be very strongly focused with F-numbers less than one: under these conditions, the uncertainty associated with use of the radiation force method becomes very significant. International Standards IEC 61161 and IEC 62555 suggest plane-wave correction factors for unfocused transducers radiating onto an ideal absorbing target and focusing corrections for focused transducers radiating onto ideal absorbing targets and onto conical reflecting targets (IEC 61161). Previous models have relied on calculations based on the Rayleigh integral, which is not strictly correct for curved sources. In the work described here, an approach combining finite element methods with a discretization of the Helmholtz equation was developed, making it possible to model the boundary condition at the structure/fluid interface more correctly. This has been used to calculate the relationship between radiation force and total power for both absorbing and conical reflecting targets for transducers ranging from planar to an F-number of 0.5 (hemispherical) and to compare with the recommendations of IEC 61161 and IEC 62555. (E-mail: [email protected]) Crown Copyright Ó 2015 Published by Elsevier Inc. on behalf of World Federation for Ultrasound in Medicine & Biology. Key Words: Acoustic radiation force, Acoustic power measurement, Ultrasound metrology, Focused ultrasonic fields, Finite element modeling, Boundary element modeling.
localized, non-invasive tissue ablation and offers the ability to treat deep-seated tumors locally with potentially fewer side effects than more invasive treatment modalities such as surgical resection, chemotherapy and ionizing radiation. Its efficacy has been reported in the treatment of a range of cancers, including those of the prostate (Ahmed et al. 2012; Sanghvi et al. 1996), liver (Leslie et al. 2012), kidney (Ritchie et al. 2010) and breast (ter Haar and Coussios 2007). It has also received U.S. Food and Drug Administration approval for the treatment of uterine fibroids, and its safety and efficacy continue to be established (Voogt et al. 2012). Two IEC standards relevant to HIFU were published in 2013 (IEC 606012-62 2013a; IEC 60601-2-62 2013a) and one Technical Specification (IEC 62556 2014). In these documents, as well as in longer established standards for ultrasonic devices such as imaging and physiotherapy systems, determination of acoustic power in water is required (IEC 60601-2-5 [IEC 2009]; IEC 61161 [IEC 2013b]).
INTRODUCTION Medical applications of ultrasound employ both focused and unfocused transducers. Ultrasound physiotherapy typically uses unfocused transducers, whereas focused transducers are widely used in ultrasonic imaging. In imaging, focusing is generally achieved by electronic beamforming with F-numbers (i.e., ratio of focal distance to aperture width) in the range 2–10. The emergence of high-intensity focused ultrasound (HIFU) over the past decade has led to an increased use of very strongly focused sources (F-number in the range 0.7–2), many of which are focused simply by virtue of having a concave radiating surface. HIFU enables highly
Address correspondence to: Pierre Gelat, Department of Mechanical Engineering, University College London, Torrington Place, London WC1E 6BT, UK. E-mail: [email protected] 832
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
The current primary standard method for measuring acoustic power radiated by an ultrasonic transducer is through measurement of the time-averaged acoustic radiation force using a gravimetric balance (IEC 61161 [IEC 2013b] and IEC 62555 [IEC 2013c]). A target is inserted in the field, and the time-averaged axial component of the radiation force vector is measured by means of the balance. The preferred target is of the absorbing type and is generally positioned orthogonal to the transducer axis. Another common target type has a right-angled convex conical reflecting surface aligned along the same axis as the transducer. Although the conical geometry is not recommended for focused fields, this type of target has the advantage that, by reflecting the ultrasonic energy, it is not susceptible to thermal damage, which can occur in absorbing targets. Acoustic radiation force is the result of a change in the momentum flux of traveling waves when absorbed or reflected by the target. It is given by the surface integral of the Langevin radiation pressure over the area of the target. If the target is perfectly absorbing, and its area is large enough, the radiation force is a characteristic of the ultrasonic field, like the acoustic power. The radiation force may be related to the acoustic power radiated by the transducer for a number of field types. For a plane progressive wave, the radiation force on a perfectly absorbing target is equal to P/c, where P is the acoustic power, and c is the speed of sound in the propagating medium. This statement is also true if the target is a right-angled convex reflector. Although this is often a good approximation to fields radiated by planar transducers (Beissner 1984a), there is evidence that large uncertainties are introduced for some transducer geometries when calculating the acoustic power from the axial component of the radiation force vector. This is particularly true in the case of strongly focused fields, as the acoustic beam may deviate much further from a plane wave than for a planar device, leading to potentially large measurement uncertainties when using a radiation force balance (Shaw 2006). An experimental method based on measuring the change in buoyancy generated by the thermal expansion of castor oil located inside a target suspended in a water bath was proposed to overcome this (IEC 62555 2013c; Shaw 2008). This technique shows promise in circumventing the limitations associated with radiation force balances and reducing measurement uncertainties for focused transducers (Rajagopal and Shaw 2012). However, because radiation force balances are still widely in use, there remains a requirement to produce satisfactory mathematical models of the measurement process and to estimate ensuing uncertainties for a range of radiators and targets. For arbitrary 3-D fields, the relationship between power and radiation force is complex and relies on a
833
tensorial treatment of the acoustic radiation stresses (Brillouin 1925). For circular rigidly vibrating plane piston sources, radiation force calculations have been the topic of much literature (Beissner 1984a, 1984b, 1987), and the relationship between power and radiation force may be characterized analytically (Beissner 1984a). For spherically curved bowl-type radiators, various attempts were made to relate radiation force to power. Corrections to the plane wave expression based on high-frequency approximations (Beissner 1987) and the geometric ray tracing approach (Shou et al. 1998) were proposed for curved radiators for absorbing targets, as well as for reflecting conical targets (IEC 62555 2013c; IEC 61161 2013b). This approach is valid only in the limit of large ka values, where k is the acoustic wavenumber and a is the source radius, and, as such, treats diffraction in an approximate way and does not address effects of multiple scattering. The radiation force from focused sources in the presence of absorbing targets was investigated as a function of target size for various ka values and focus half-angle values using a Rayleigh integral approach (Beissner 2010). The size of targets required to capture 97%, 98% and 99% of the radiation force that would be experienced by a target of infinite diameter was determined for a range of ka values, target distances and transducer apertures. The Rayleigh integral is a mathematical formulation of the Huygens–Fresnel principle, where all points on the surface of the radiator are local sources of secondary waves. The acoustic pressure generated by the radiator at any point within the field is then a superposition of point sources with amplitudes proportional to the component of the particle velocity normal to the vibrating surface. For a planar radiator mounted in an infinite baffle, this formulation is exact. For curved transducer surfaces, hemispherical elementary waves are initially diffracted by the radiator itself, and their structure is changed. Although various formulations have been proposed in the literature (Hasegawa et al. 1987; Lucas and Muir 1982; Wu 1995), a Rayleigh integral approach is not fully valid when investigating concave radiators with low F-numbers. The latter is defined as the ratio of the focal length of the source to its diameter. Various approaches have been suggested to overcome this problem (Cathignol and Sapozhnikov 1999; Sapozhnikov and Sinilo 2002). Nevertheless, these approaches are not generalizable to arbitrary sources or targets, and it is likely that a numerical solution to the problem will be more versatile in addressing a wider range of configurations. Although they can be computationally intensive, finite element/boundary element (FE/BE) methods are suitable for modeling a wide range of problems in
834
Ultrasound in Medicine and Biology
acoustics and vibration. Unlike ray tracing and Rayleigh integral approaches, FE/BE methods are exact to within discretization error. FE/BE methods overcome the limitations of the Rayleigh integral for curved radiators in that they account for both the pressure and velocity terms in the Kirchhoff–Helmholtz integral (Pierce 1981). The Rayleigh integral accounts only for the velocity term, which strictly is valid only for planar radiators. An FE/BE approach was proposed to model the field radiated by a plane rigidly vibrating circular piston in the presence of a polyurethane aperture, and radiation force calculations were carried out for a range of configurations (Gelat et al. 2005). As part of this work, the above approach was extended to feature a range of curved spherical concave radiators rigidly vibrating in an infinite baffle in the presence of both perfectly absorbing and conical reflecting targets. The dynamics of the conical target were accounted for. This approach takes into consideration full coupling between the radiator, the propagating medium and the target, hence implementing the correct boundary condition at the interfaces between the fluid and the vibrating structures. As such, effects of scattering, diffraction and multiple reflections may be accurately modeled. Acoustic radiation force calculations were carried out from first-order acoustic pressure and particle velocities for each configuration. The acoustic power radiated by each source was also calculated and the ratio Fc/P determined so as to estimate the uncertainties introduced by plane wave and approximate focusing corrections given in IEC 61161 (IEC 2013b) and, more recently, IEC 62555 (IEC 2013c). Currently, the most established method of measuring acoustic power radiated by an ultrasonic source is by means of the time-averaged force acting on a sufficiently large target positioned in the acoustic field. Assuming plane wave propagation in the direction normal to the target, the radiation force on a perfectly absorbing target is related to the acoustic power as follows. F5
P c
(1)
For an unfocused source represented by a rigidly vibrating plane circular piston mounted in an infinite baffle, it may be shown that if the target is perfectly absorbing (Beissner 1984a), F5
P 12J20 ðkaÞ 2J21 ðkaÞ : c 12J1 ð2kaÞ=ðkaÞ
(2)
The second term on the right-hand side is often called the ‘‘plane wave correction.’’ For values of ka $35, eqn (2) deviates from the plane wave expression
Volume 41, Number 3, 2015
in eqn (1) by ,2%. As such, many ultrasonic transducers can, in practice, be characterized using eqn (1) (Beissner 1984a). The force arising from fields generated by non-planar ultrasonic transducers will in general differ from the plane wave expression. For a spherical bowl transducer radiating onto a perfectly absorbing target, a relationship based on a high-frequency asymptotic approximation may be derived (IEC 61161 [IEC 2013b]; Shou et al. 1998). F5 where g is case of a right-angled cone tip is its centre of given by
P ð11cos gÞ c 2
(3)
the transducer focus half-angle. In the spherical transducer radiating onto a convex conical reflecting target where the positioned between the transducer and curvature, the corresponding expression is
F5
P 12cos 2g12g2sin 2g c 4ð12cos gÞ
(4)
where g must be expressed in radians. The expressions in eqns (3) and (4) are of practical significance when characterizing a subset of HIFU transducers, but nevertheless rely on high-frequency approximations to describe the radiation of sound from the transducer and its effect on the target: they also assume that the behavior of the real transducer is a good approximation to a spherical piston and that the wave reflected from the target is not incident on the transducer. Hence, they are of limited value for transducers with small or moderate ka values or of uncertain construction. Furthermore, in transducers with low F-numbers, effects caused by multiple scattering off the surface of the device become increasingly important and will not be accounted for by eqn (4). Finally, effects of multiple scattering of sound between the surfaces of the transducer and the conical target are not considered in the theory underpinning the plane wave or simple focusing corrections: of course, multiple reflections should be avoided for accurate measurement, but it is relevant to examine the effect theoretically to understand properly its importance. In this article, an approach based on the discretization of the Helmholtz equation in its integral form is used, which, although more computationally involved than high-frequency approximations, is more suitable for describing the above phenomena and is likely to help assess the range of validity of eqns (3) and (4) as a function of ka and g By combining this approach with finite element methods, it has been possible to model the boundary condition at the structure/fluid interface more correctly.
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
METHOD Radiation force theory Radiation force theory is now briefly reviewed. Consider a homogeneous, non-absorbing, non-dispersive, sound propagating fluid. Assume that the primary field quantities are acoustic pressure p and particle velocity ! u and are both functions of space and time. The medium is characterized by an equilibrium density (r) and sound speed (c). Let the position vector in Eulerian coordinates be characterized by ! r ðx; y; zÞ in a Cartesian coordinate system. The linear acoustics result for the Brillouin stress, which is the acoustic radiation stress tensor (S) in Eulerian coordinates, is given by (Beissner 1998) (5) Sij 5 T1 1T2 1T3 2V dij 2rui uj ; where i, j 5 1, 2 and 3; 1, 2 and 3 correspond to the global (x, y and z) directions in a Cartesian axis set, respectively; the bars denote temporal averages; Ti 5
rjui j 4
2
835
the radiation force is independent of the target distance because of the momentum conservation law (Beissner 1998). Acoustic radiation force theory for a perfectly absorbing circular target It is assumed that the main axis of propagation is along the positive Cartesian x-axis. For a circular, planar, perfectly absorbing target whose receiving area is parallel to the y-z plane, the radiation force vector is expressed as 0 1 0 1 0 1 Sxx Sxy Sxz 21 F ! @ xA F 5 Fy 5 %A @ Syx Syy Syz A$@ 0 AdA; (7) Szx Szy Szz 0 Fz where A is the surface area of the circular target. For an axisymmetric field, where x is the axialpcoordinate and ffiffiffiffiffiffiffiffiffiffiffiffi R is the radial coordinate such that R 5 y2 1z2 , it may be shown that (Beissner 1984a) ðb Fx 5 2p V1T x 2T R R dR; (8) 0
is the time averaged kinetic acoustic energy density;
where b is the target radius.
2
V5
jpj 4rc2
is the time-averaged potential acoustic energy density; dij is the Kronecker delta (dij 5 1 if i 5 j) and (dij5 0 if i s j); r is the medium equilibrium density; ui is the instantaneous (complex) particle velocity along the i direction; and r is the instantaneous (complex) acoustic pressure. 1 0 Sxx Fx ! @ A F 5 Fy 5 %A1 @ Syx Szx Fz 0
Sxy Syy Szy
Acoustic radiation force theory for a reflecting conical target The main axis of propagation is once again assumed to be the positive Cartesian x-axis. This is also assumed to be the axis of the conical target. If the axis of the cone makes an angle a to the generatrix, the radiation force vector is given by
1 0 1 0 Sxx Sxz 2sin a Syz A$@ cos a cos q AdA1 1%A2 @ Syx Szz Szx cos a sin q
The radiation stress tensor (S) expresses the momentum flux density (Brillouin 1925), and the above expression has been found to be valid to second order when r and ! u are first-order quantities. The radiation force vector ! ( F ) is given by integrating the acoustic radiation stress tensor over a surface A enclosing the target (Beissner 1984a) 0 1 Sxx Sxy Sxz ! n dA; (6) F 5 %A @ Syx Syy Syz A$! Szx Szy Szz where ! n is an outward pointing unit vector normal to the surface A. In a lossless fluid, for an infinitely large target,
Sxy Syy Szy
1 0 1 Sxz 1 Syz A$@ 0 AdA2 ; Szz 0
(9)
where tan q 5 z/y, A1 corresponds to the area defined by the surface of the conical region and A2 is the circular area defined by the base of the cone. A radiation force balance responds to the axial component of the radiation force vector Fx. For an axisymmetric field, where x is the axial coordinate pffiffiffiffiffiffiffiffiffiffiffiffiand R is the radial coordinate such that R 5 y2 1z2 , it may be shown that by substituting the expressions for the components of the stress tensor into eqn (9) and by applying the appropriate coordinate transformation, ðh Fx 5 22p tan2 a 0
T R 2T x 2V x dx
836
Ultrasound in Medicine and Biology
ðh 22p tan a
ð R1 rux uR x dx22p
0
Volume 41, Number 3, 2015
V1T x 2T R R dR;
0
(10) where h is the height of the cone and R1 is the radius of its base and ux uR 5
1 ! ! ux ð r ÞuR ð r Þ1uR ð! r Þux ð! rÞ ; 4
(11)
r Þ is the spatial component of the axial where ux ð! component of the complex particle velocity vector, r Þ is the spatial component of the radial component uR ð ! of the complex particle velocity vector, and the asterisk denotes the complex conjugate. Numerical modeling: Perfectly absorbing target The methodology employed here is similar to that used by Gelat et al. (2005). Circular symmetry about the Cartesian x-axis was assumed for both reflecting and absorbing target configurations. The spherical bowl source was modeled as a vibrating, rigid, curved piston mounted in an infinite baffle and radiating into a semi-infinite non-absorbing fluid. Three-noded quadratic thin axisymmetric shells of revolution were used, and all motions other than translational movement along the x-axis were restrained. A unit radial velocity magnitude and zero phase were assigned to one of the piston nodes, and this was repeated throughout all the nodes on the surface of the radiator. A combination of acoustic finite element and boundary element patches were used to model the fluid. A region of pressure-based, quadratic fluid finite elements in front of the curved piston was meshed to a distance of a quarter of a wavelength from the rim using a combination of six-noded triangular and eight-noded quadrilateral elements. A boundary element was then applied at the outer edge of fluid finite element region for the perfectly absorbing target calculations, as illustrated in Figure 1. The targets were positioned two wavelengths in the fluid away from the rim of the transducer. At positions along the radial direction from the rim of the transducer and outside its radius, a zero acoustic pressure normal derivative boundary condition was applied in the plane of the transducer, extending to infinity. This boundary condition effectively simulates a rigid baffle. A sound speed of 1487 m s21 and a density of 998 kg m23 were assumed for the acoustic medium. The frequency of excitation for all configurations was 1 MHz. Thirteen source radii were investigated: 1, 2, 3, 4, 6, 8, 10, 15, 20, 25, 30, 40 and 50 mm. For each of these configurations, 11 linearly spaced values for sin g were investigated ranging from 0 (plane piston case) to 1
Fig. 1. Representative curved radiator with absorbing target configuration revealing the finite element mesh and the boundary element. Curved radiator radius 5 10 mm. sin g 5 0.6.
(hemispherical radiator) in steps of 0.1. The F-numbers corresponding to the values of sin g are displayed in Table 1. This leads to a total of 143 absorbing target configurations. An analysis of how the target radius affects the estimate of the radiation force and measured acoustic power was not carried out, as this has been already been addressed by Beissner (2010). Instead, a large enough target radius b was considered to ensure that the value of Fx was beyond 0.98 of the limit of Fx as b (Beissner 2010). This was achieved by selecting b=a 5 10 for all configurations except for the sources of 1-mm radius, where b=a 5 100 was used, as using b=a 5 10 for the smallest, most highly focused hemispherical configurations resulted in underestimation of the acoustic power on the target by approximately 30%.
Table 1. sin g values and corresponding F-numbers sin g 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 F-number N 5 2.5 1.67 1.25 1 0.83 0.71 0.625 0.56 0.5
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
In boundary element methods, the discretization of the Helmholtz integral equation leads to a fully populated asymmetric system of linear equations. It is well known that the conventional Helmholtz boundary integral equation (BIE) formulation for exterior wave problems yields non-unique solutions at the eigenfrequencies associated with the (interior) Dirichlet problem (Burton and Miller 1971). When the coefficient matrix resulting from discretization of the Helmholtz boundary integral equation is inverted at frequencies approaching these (fictitious) eigenfrequencies, the problem becomes ill-conditioned. A way in which this may be overcome is to use the combined Helmholtz integral equation formulation (CHIEF) method (Schenk 1968). This is based on the fact that at any frequency, there is a unique surface pressure function that satisfies both the interior and surface Helmholtz integral equations. Additional collocation points were specified within regions surrounding the BE. These supplement the degrees of freedom on the BE surface and result in an overdetermined set of linear equations. Five collocation points were positioned in the region interior to the boundary element to overcome problems associated with non-uniqueness. The PAFEC (Programme for Automatic Finite Element Calculations) software (PACSYS, Nottingham, UK) was used to carry out the calculations. PAFEC has the capability of dealing with the coupling between the elastic structure and the acoustic wave in a fluid by using a combined finite element/boundary element method. For quadratic elements, meshing at three elements per wavelength is normally sufficient, but six elements per wavelength were nonetheless used in all meshes. PAFEC enables the calculation of acoustic pressures and particle velocities at user-specified locations within the boundary element exterior domain. These first-order quantities were used to determine Fx for both the reflecting and absorbing targets as expressed by eqns (8) and (10), in which numerical quadrature was used to carry out the integrations. The acoustic power radiated into the fluid by the source was computed in PAFEC as the integral of the normal component of the time-averaged acoustic intensity over the surface defined by the fluid–structure interface as
837
For a plane circular piston rigidly vibrating sinusoidally in an infinite baffle, the acoustic power radiated may be evaluated analytically as (Kinsler et al. 1982)
1 2J1 ð2kaÞ ; (13) P 5 U02 pa2 rc 12 2 2ka where U0 is the maximum piston velocity. There is no equivalent exact analytical solution in the literature for curved pistons. Numerical modeling: Reflecting conical target The reflecting conical targets were modeled using a three-noded quadratic thin shell of revolution elements, but no restraints other than those required for axisymmetric motion were applied. The reflecting target was assumed to be made of a steel outer shell 0.2 mm thick enclosing an air volume (modeled as a vacuum). This is similar to radiation force targets used previously on the National Physical Laboratory Primary Standard Radiation Force Balance. The mechanical properties of the steel were: Young’s modulus 5 215 GPa, density 5 7800 kg m23 and Poisson’s ratio 5 0.293. The boundary element was extended to encompass the region surrounding the cone, as illustrated in Figure 2. The region of air normally present inside conical reflecting targets was modeled as a vacuum. Continuity of normal particle velocity was assumed to couple the structure and the fluid. The cone angle was 90 , which is what is commonly used experimentally. The radius h of the base of the cone was assumed to be twice that of the source radius a which, for the conical target calculations, was fixed at 20 mm. This value for the source radius was chosen as it was toward the middle of the range of radii investigated in the absorbing target case. Investigation of other source radii was beyond the scope of the work described in this article. To increase the degrees of freedom on the BE surface, 419 additional collocation points were used inside the cone. Eleven linearly spaced sin g values ranging from 0 to 1 were investigated. The
1 P 5 %Asource ½pð! r Þ! u ð! r Þ1p ð! r Þ! u ð! r Þ:! n dAsource ; 4 (12) wherepð! r Þ is the spatial component of the complex acoustic pressure, ! u ð! r Þ is the spatial component of the complex particle velocity vector; Asource is the transducer radiating area; ! n is the unit vector normal to the area Asource ; and the asterisk denotes the complex conjugate.
Fig. 2. Representative curved radiator with reflecting target configuration revealing the finite element mesh and the boundary element. Curved radiator radius 5 20 mm. sin g 5 0.7.
838
Ultrasound in Medicine and Biology
Volume 41, Number 3, 2015
Fig. 3. Acoustic power radiated by a 1 MHz planar transducer. Percentage difference between analytical and finite element/boundary element (FE/BE) solutions. Inset: Comparison of FE/BE approach (solid line) with the analytical solution, eqn (13) (dotted line).
apex of the conical target was positioned two wavelengths into the fluid away from the rim of the transducer. Note that in the case of sin g 5 1, the cone apex lies beyond the focus so the target does not intercept the whole field; this most extreme condition may be of limited practical relevance, but is included for completeness. RESULTS Comparison against analytical solutions To confirm that the FE/BE model is reliable, acoustic power and the ratio Fc=P can be compared with theory. Figure 3 illustrates the acoustic power calculated using the FE/BE approach against the analytical solution in eqn (13) for a 1 m s21 piston normal velocity: agreement is within 0.1% for all cases. For convenience, the normal component of the radiation force vector will be referred to as F. Figure 4 illustrates the ratio Fc=P for the FE/BE calculations on the plane circular piston case compared against the analytical solution in eqn (2). The agreement between both sets of results is within 1.1%. Because only 11 ka values were used in these calculations, small oscillations caused by the Bessel function in eqn (13) are not resolved: nevertheless, Figures 3 and 4 taken together serve to validate the FE/BE model to better than 1%. Perfectly absorbing target calculations The ratio Fc=P was computed for the 143 source and target configurations described in the previous section.
The results are summarized in Figure 5(a, b). The focusing correction factor taken from IEC 61161 (IEC 2013b) expressed in eqn (3) is overlaid for each case as horizontal dotted lines for each value of g. Figure 5(a, b) illustrates that, for larger values of the source radius a, the IEC 61161 solution for the ratio Fc=P is very similar to the FE/BE results: the percentage difference is displayed as a function of a in Figure 6. The case of sin g 5 1 has been omitted from Figure 6, as the percentage difference with the IEC 61161 expression is large (between 241% and 131%). Conical reflecting target calculations The ratio Fc=P was computed for the 11 transducer aperture angles described in the previous section. The results are illustrated in Figure 7(a) and compared against the recommended correction factor provided by IEC 61161 (IEC 2013b) for reflecting conical targets described by eqn (4). Figure 7(b) illustrates the same data for values of sin g between 0 and 0.7. Figure 8 illustrates the acoustic pressure field between the transducer and the target for the case that gives an enhanced ratio of Fc=P. In this case, the boundary condition at the surface of the source was chosen as a radial displacement of 1028 m. The corresponding conical target displacements are illustrated in Figure 9, where the source radial displacement has been scaled by a factor of 105 relative to that in Figure 8 so as to better visualize the deflections of the target. The maximum target displacement amplitude is in fact 10.7 3 1028 m.
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
839
Fig. 4. 1 MHz planar transducer radiating on a perfectly absorbing target. Comparison of relationship between radiation force axial component and acoustic power for FE/BE approach (solid line) with that for analytical solution, eqn (2) (dotted line). FE/BE 5 finite element/boundary element.
Figure 8 illustrates that multiple reflections are occurring between the radiator and the target for this focusing angle, leading to a standing wave pattern. Figure 9 illustrates that the target displays modal behavior. DISCUSSION The IEC 61161 (2013b) standard typically requires that individual sources of uncertainty contribute less
than 1% to the uncertainty of the overall measured power. In IEC 62555 (IEC 2013c), which is specific to high-intensity therapeutic ultrasound sources, this criterion is relaxed to 2%. In the case of a perfectly absorbing target, analysis of the data in Figure 4 reveals that the agreement of the FE/BE calculations is within 0.1% of the analytical solution. It can therefore be concluded that, for the range of planar radiators investigated, the radiation force
Fig. 5. 1 MHz transducer radiating on a perfectly absorbing target. Comparison of relationship between radiation force axial component and acoustic power for FE/BE approach (solid line) with that for IEC 61161 (IEC 2013b, 2013c) solution, eqn (3) (dotted line). (a) 0 # sin g # 0.5. (b) 0.6 # sin g # 1. FE/BE 5 finite element/boundary element.
840
Ultrasound in Medicine and Biology
Volume 41, Number 3, 2015
Fig. 6. 1 MHz transducer radiating on a perfectly absorbing target. Percentage difference between radiation force axial component and acoustic power using eqn (3) (IEC 61161), relative to finite element/boundary element approach and as a function of a.
normal to the target predicted by the FE/BE method is within 1% of its value for an infinite target. The discrepancies are likely due to a combination of insufficient sampling of points on the target and the use of a target of finite size. The relationship between force and power given by eqn (3) for an absorbing target is based on a Rayleigh integral calculation for a spherical section radiator with uniform radial motion. The methodology employed in this article, also applied to a spherical section radiator with uniform radial motion, overcomes the limitations of the Rayleigh integral approach by including both the monopole and dipole terms in the Kirchhoff–Helmholtz integral. At 1 MHz, the ratio Fc=P provided by the IEC 61161 standard and expressed in eqn (3), is within approximately 1% of the improved calculation for transducers with F-number $0.83 and of radius $10 mm. It is consistent to within 2% for sin g # 0.9 as long as the source radius is larger than approximately 4 mm, that is, ka $ 17. Many HIFU devices fall within these ranges. These results provide additional confidence in the use of the IEC 61161 correction factor for spherically focusing sources. The IEC 61161 results diverge from the FE/BE solution as the source becomes increasingly focused and as its radius decreases. A similar set of calculations (not included in this article) for a spherical section radiator excited in a direction parallel to the x-axis gave comparable results for this range of source parameters. In Figure 5(a), the dashed lines
represent the focusing correction factor taken from IEC 61161 expressed in eqn (3). For the unfocused case, g 5 0, the dashed line is drawn at Fc=P 5 1, and is strictly only correct in the limit of large transducer radii. For plane piston transducers of finite size, diffraction occurs, leading to deviation form plane-wave theory, and the ratio Fc=P should be less than 1. Indeed, IEC 61161 provides a formula to calculate a plane-wave correction (which is approximately 2.4% for a 5-mm-radius source at 1 MHz, decreasing to 1.1% for 10-mm-radius and 0.2% for 50-mm-radius sources). It is possible that a weak degree of focusing may compensate for diffractive spreading resulting in a higher ratio of Fc=P than in the plane piston case: there is evidence here to support this possibility. Note in Figure 5(a) that the Fc=P ratios are comparable in the plane piston case and in the case of the curved radiator featuring sin g 5 0.1. This is true for the range of piston radii considered in this article. However, as sin g increases, the effects of focusing clearly overcome those of beam diffraction. Finally, it is worth noting that the IEC 61161 and IEC 62555 standards [IEC 2013b, 2013c] do not address the question of how to combine corrections for a finite-size transducer with corrections for focusing. Although the data in Figure 5(a, b) may serve as a starting point in establishing such a correction factor, only a subset of radiation force correction factors would be represented, as limited source radii have been investigated at only one excitation frequency (1 MHz).
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
841
a
b
Fig. 7. (a) 1 MHz 20-mm-radius curved transducer radiating on a conical reflecting target of 90 aperture angle and 40-mm base radius. Comparison of relationship between radiation force axial component and acoustic power for FE/BE approach (solid line) with that for IEC 61161 (dotted line). (b) 1 MHz, 20-mm-radius curved transducer radiating on a conical reflecting target of 90 aperture angle and 40 mm base radius. Comparison of relationship between radiation force axial component and acoustic power for FE/BE approach (solid line) with that for IEC 61161 (dotted line) for 0 # sin g # 0.7. FE/BE 5 finite element/boundary element.
More extensive numerical investigations on curved radiators excited over a range of frequencies could be carried for a wider range of source radii. This may establish further sets of radiation force correction factors,
but this is beyond the scope of the work described in this article. For the conical reflecting target configurations investigated as part of this work, it can be seen in
842
Ultrasound in Medicine and Biology
Fig. 8. Acoustic pressure field in megapascals from FE/BE calculations for a 1 MHz, 20-mm-radius curved transducer radiating on conical reflecting target of 90 aperture angle and 40-base radius. sin g 5 0.8. Source radial displacement 5 1028 m. FE/BE 5 finite element/boundary element.
Figure 7b that the IEC 61161 (2013b) correction factor agrees with the FE/BE model results within 2% for values of sin g between 0 and 0.6 (which is an F-number of 0.83): so, again, the IEC formula is sufficiently accurate for most practical purposes. However, with a reflecting target, there is the possibility of interaction between the reflected beam and the transducer, which can lead to a strong modulation of the radiation force. For sin g 5 0.8 and the target distance chosen (3 mm, which corresponds to two wavelengths in water at 1 MHz), the numerical modeling calculation of Fc=P is a factor of 5 higher than the result predicted by
Fig. 9. Conical target displacement profile resulting from FE/BE calculations for a 1 MHz, 20-mm-radius curved transducer radiating on conical reflecting target of 90 aperture angle and 40-mm base radius. sin g 5 0.8. Source radial displacement 5 1028 m, resulting in a maximum displacement of the target of 10.7 3 1028 m. For visualization purposes, displacements have been scaled by a factor of 105. FE/BE 5 finite element/boundary element.
Volume 41, Number 3, 2015
the IEC 61161 expression (see Fig. 7a). Equation (4) considerably underestimates the radiation force-toacoustic power ratio for this configuration. Although a 3mm separation distance might seem rather close, it is generally recommended that output power measurements be made close to the transducer surface, and conical targets do not show such enhanced behavior when used very close to planar transducers. The value of g, at which this interaction is strongest, will depend on the distance of the target; this type of enhancement has been observed experimentally (Shaw and Hodnett 2008) where an approximate doubling in the radiation force on a conical target was reported at distances between 3–15 mm from a focused transducer with F-number 1.5. The results in Figure 7(a, b) further highlight the limitations of conical targets and why their use should be avoided, particularly when dealing with focused fields (Shaw 2008). These simulations suggest that, in addition to the focusing errors introduced, the dynamic behavior of the target together with multiple reflections between the target and the radiator has the potential to cause large uncertainties in the power measurement. CONCLUSIONS A combination of finite element and boundary element modeling was used to calculate the relationship between radiation force and total power for both absorbing and conical reflecting targets for 1 MHz transducers ranging from planar to an F-number of 0.5 (hemispherical). This method overcomes the limitations of the Rayleigh integral for curved radiators in that it accounts for both the pressure and velocity terms in the Kirchhoff–Helmholtz integral (Pierce 1981). The Rayleigh integral accounts only for the velocity term, which is strictly valid for only planar radiators. The findings were compared with the recommendations of IEC 61161 and IEC 62555 (IEC 2013b, 2013c). The modeling described here, in which the radiation force was calculated by first determining the acoustic stress tensor, is more versatile than a Rayleigh integral or analytical approach because it allows acoustic interaction between the transducer and target to occur. Furthermore, the approach can be extended to more complex transducer models, including piezo-electric effects, and to more realistic targets of non-zero reflectivity. The modeling reported here is axisymmetric and limited to a single frequency of 1 MHz: long run times on some configurations will make it difficult to generalize to three dimensions on a desktop PC, but future work could address large 3-D problems using cluster computing. Within the restricted assumption of ideal spherical section transducers operating in an infinite baffle, it
Acoustic power and acoustic radiation force d P. GELAT and A. SHAW
was found that the correction factor from IEC 61161 and IEC 62555 (IEC 2013b, 2013c) for planar absorbing targets is consistent with the results presented in this article to within 1% for transducers with F-number $0.83 (i.e., sin g 5 0.6) and of radius $10 mm; it is consistent to within 2% for sin g , 0.9 as long as the source radius is greater than about 4 mm (ka $ 17); for smaller radii, a larger correction factor is required. The same correction factor can be applied irrespective of whether the surface displacement is radial or parallel to the x-axis. Finally, it is worth noting that the IEC 61161 and IEC 62555 standards do not presently address combining radiation force correction factors for focusing transducers of finite dimensions, which could be significant for small transducers or at low frequencies. The FE/BE method employed in the work described here to calculate the radiation force of focused radiators on perfectly absorbing targets may go some way toward addressing this, although simulations across a much greater range of devices would be required. With the same restrictions, this time considering the case of a 20-mm-radius transducer (ka 5 84) and a rightangled conical reflecting target, it was found that the focusing correction from IEC 61161 is again in good agreement with the FE/BE results for sin g , 0.6. However, care should be taken not to move the transducer too closely to the target, as reflections may greatly increase the radiation force and therefore overestimate the acoustic power. To the authors’ knowledge, this is the first time a full acoustic model of the radiation force on a reflecting target has been carried out. It is planned that this work be extended to more realistic transducer models in the future. Finally it is worth noting that as these calculations were carried out at an excitation frequency of 1 MHz, care should be taken in extrapolating from these results to other frequencies. Further calculations over a wider range of frequencies relevant to HIFU devices are desirable to generalize these conclusions. Acknowledgments—The authors gratefully acknowledge funding for this work provided by the United Kingdom Department for Business, Innovation and Skills, through the Acoustics and Ionising Radiation Metrology Programme.
REFERENCES Ahmed HU, Cathcart P, Chalasani V, Williams A, McCartan N, Freeman A, Kirkham A, Allen C, Chin J, Emberton M. Whole-gland salvage high-intensity focused ultrasound therapy for localized prostate cancer recurrence after external beam radiation therapy. Cancer 2012;118:3071–3078. Beissner K. Acoustic radiation pressure in the near field. J Sound Vib 1984a;93:537–548. Beissner K. Minimum target size in radiation force measurements. J Acoust Soc Am 1984b;76:1505–1510. Beissner K. Radiation force calculations. Acustica 1987;62:255–263.
843
Beissner K. The acoustic radiation force in lossless fluids in Eulerian and Lagrangian coordinates. J Acoust Soc Am 1998;103:2321–2332. Beissner K. Minimum radiation force target size for power measurements in focused ultrasonic fields with circular symmetry. J Acoust Soc Am 2010;128:3355–3362. Brillouin L. Sur les tensions de radiation [On the radiation stresses]. Ann Phys (Paris) 1925;4:528–586. Burton AJ, Miller GF. The application of integral equation methods to the numerical solution of some exterior boundary-value problems. Proc R Soc London Ser A 1971;323:201–210. Cathignol D, Sapozhnikov OA. On the application of the Rayleigh integral to the calculation of the field of a concave focusing radiator. Acoust Phys 1999;45:735–742. Gelat PN, Zeqiri B, Hodnett M. A finite element model of the aperture method for determining the effective radiating area of a physiotherapy treatment head. Ultrasonics 2005;43:321–330. Hasegawa T, Inoue N, Matsuzawa K. A new theory for the radiation from a concave piston source. J Acoust Soc Am 1987;82:706–708. International Electrotechnical Commission (IEC). IEC 60601-2-5: Medical electrical equipment: Particular requirements for the basic safety and essential performance of ultrasonic physiotherapy equipment. 3rd ed. Geneva: 2009. International Electrotechnical Commission (IEC). IEC 60601-2-62: Medical electrical equipment: Particular requirements for the basic safety and essential performance of high intensity therapeutic ultrasound (HITU) equipment. 1st ed. Geneva: 2013a. International Electrotechnical Commission (IEC). IEC 61161: Ultrasonics—Power measurements—Radiation force balances and performance requirements. 3rd ed. Geneva: 2013b. International Electrotechnical Commission (IEC). IEC 62555: Ultrasonics—Power measurement—Output power measurement for high intensity therapeutic ultrasound (HITU) transducers and systems. 1st ed. Geneva: 2013c. International Electrotechnical Commission (IEC). IEC 62556: Ultrasonics—Surgical systems—Specification and measurement of field parameters for high intensity therapeutic ultrasound (HITU) transducers and systems. Technical Specification. Geneva, 2014. Kinsler LE, Frey AR, Coppens AB, Sanders JV. Fundamentals of acoustics. 3rd ed. New York: Wiley; 1982. p. 192. Leslie T, Ritchie R, Illing R, ter Haar G, Phillips R, Middleton M, Wu F, Cranston D. High-intensity focused ultrasound treatment of liver tumours: post-treatment MRI correlates well with intraoperative estimates of treatment volume. Br J Radiol 2012;85: 1363–1370. Lucas BG, Muir TG. The field of a focusing source. J Acoust Soc Am 1982;72:1289–1296. Pierce AD. Acoustics: An introduction to its physical principles and applications. New York: McGraw–Hill; 1981. p. 214. Rajagopal S, Shaw A. The buoyancy method—A potential new primary ultrasound power standard. Metrologia 2012;49:327–339. Ritchie RW, Leslie T, Phillips R, Wu F, Illing R, ter Haar G, Protheroe A, Cranston D. Extracorporeal high intensity focused ultrasound for renal tumours: A 3-y follow-up. BJU Int 2010;106: 1004–1009. Sanghvi NT, Fry FJ, Bihrle R, Foster RS, Syrus J, Zaitsev AV, Hennige CW. Noninvasive surgery of prostate tissue by highintensity focused ultrasound. IEEE Trans Ultrason Ferroelectr Freq Control 1996;43:1099–1110. Sapozhnikov OA, Sinilo TV. Acoustic field produced by a concave radiating surface with allowance for the diffraction. Acoust Phys 2002;48:720–727. Schenck HA. Improved integral formulation for acoustic radiation problems. J Acoust Soc Am 1968;44:41–58. Shaw A. How to measure HIFU output power properly. In: Proceedings, 5th International Symposium on Therapeutic Ultrasound (2006). New York: IEEE; 2006. p. 628–632. Shaw A. A buoyancy method for the measurement of total ultrasound power generated by HIFU transducers. Ultrasound Med Biol 2008; 34:1327–1342. Shaw A, Hodnett M. Calibration and measurement issues for therapeutic ultrasound. Ultrasonics 2008;48:234–252.
844
Ultrasound in Medicine and Biology
Shou W, Wang Y, Qian D. Calculation for radiation force of focused ultrasound and experiment measuring HIFU power. Tech Acoust 1998;17:145. Ter Haar G, Coussios C. High-intensity focused ultrasound: Past, present and future. Int J Hyperthermia 2007;23:85–87. Voogt MJ, Trillaud H, Kim YS, Mali WPTh M, Barkhausen J, Bartels LW, Deckers R, Frulio N, Rhim H, Lim HK, Eckey T, Nieminen HJ,
Volume 41, Number 3, 2015 Mougenot C, Keserci B, Soini J, Vaara T, K€ohler MO, Sokka S, Maurice AJ, van den Bosch A. Volumetric feedback ablation of uterine fibroids using magnetic resonance-guided high intensity focused ultrasound therapy. Eur Radiol 2012;22:411–417. Wu J. Calculation of acoustic radiation force generated by focused beams using the ray acoustics approach. J Acoust Soc Am 1995; 97:2747–2750.