Fast Pressure Field Simulations of Large Ultrasound Phased Arrays

Results

Introduction

Spherical ultrasound phased arrays comprised of circular transducers are being evaluated in simulation studies for thermal therapy. These arrays, which have a diameter of 16cm, are harmonically excited at 1MHz. The apertures of these arrays are sampled every 1.5 wavelength to avoid grating lobes. The pressure fields are sampled twice per wavelength to avoid spatial aliasing during temperature calculations. These result in a compact array model with more than ten thousand sources and over ten million observation points in a three dimensional field. The Fast Near-field Method (FNM) is combined with the Angular Spectrum Approach (ASA) to quickly and accurately compute the acoustic pressure fields generated by these arrays. The fields are propagated through the computational grids using a 2D FFT, where the source fields are computed from FNM. This approach achieves a significant reduction in the computation time. The numerical errors in simulations using ASA are determined by the sampling and extent of the spatial grids and by the propagation distance. The errors are then characterized in terms of the corresponding changes in the computed temperature response encountered in bio-heat transfer simulations.

Array Model

A spherical phased array of 16cm diameter and 15cm curvature is populated with 11,845 circular sources with diameter λ/2. The center-center spacing between adjacent elements is λ. The excitation frequency is 1MHz.


Fig 1. Top view of the spherical phased array	Fig 2. Side view of the spherical phased array

Modeling methods

Fast near field method (FNM)

The formula to computes the pressure field generated by a harmonically excited circular transducer with radius a is

The 1D integral over ψ is evaluated by mid-point Riemann sum.

Angular Spectrum Approach (ASA)

This approach takes the known two-dimensional field generated by the array at plane A, evaluates the Fourier transform, multiplies by a propagation term H, and obtains the spectrum in plane B parallel to plane A, assuming the wave propagation satisfies the Helmholtz equation. An inverse Fourier transform of the result reconstructs the field at plane B. H also contains the linear attenuation effect through the trigonometric function in the exponent.

where Δz=z-z0, P(kx, ky, z) and P0(kx, ky, z0) are the Fourier transform of pressure at arbitrary transverse plane p(x, y, z) and at source plane p0(x, y, z0). H(kx, ky, Δz) is the spectral propagator.

Figure 3. Source plane in ASA calculation. The red circle illustrates the array. The dark blue grid indicates the square plane A where the pressure field is known. The light blue region is the zero-padding area, discretized at the same rate as in plane A.

The pressure is evaluated at MXM points at source plane A. The dimension L of the source plane determines how much spatial information is included. Zero padding prevents the wrap-round error caused by circular convolution and increases the spatial samples to N in each direction. Small sampling δ (L/M) is required to avoid spatial aliasing.

The discretized angular frequencies are

Results

The ASA error increases as the field propagates away from the source plane. By including multiple source planes the error is reduced. Two methods are used to evaluate the error at each plane parallel to z axis. The results are plotted in Figs 7 and 8. One is the peak normalized error, which is the absolute value of the difference between the ASA result and the reference normalized by the peak value in 3D field. The other is the root mean square error, which is the root mean square error in each plane normalized by the peak value in that plane.

Fig 4. Pressure field generated by the spherical Fig 5. Temperature field generated by the spherical

phased array evaluated in focal plane (z=15). phased array evaluated in focal plane (z=15).

Fig 6. Pressure field generated by the spherical Fig 7. Temperature field generated by the spherical

phased array evaluated in depth plane (y=0). phased array evaluated in depth plane (y=0).

Fig 8. The peak normalized error of the pressure Fig 9. The peak normalized error of temperature

in the depth plane (y=0) using five source planes. in the depth plane (y=0) using five source planes.

The error is restricted to 4%. The error is restricted to 2%.

Fig 10. ASA error in pressure field computation. Left: rms error, right: peak normalized error. The error increases in ASA as the distance from the source plane increases. 5 source planes are placed at z = 7, 10.45, 13.1, 14.6 and 19.7cm, and the error is evaluated in each plane normal to the z direction. The rms error and the normalized peak error are both restricted to 4%.

Fig 11. ASA error in temperature field computation. Left: rms error, right: peak normalized error. The 3D pressure field computed with the ASA using multiple source planes is used as the input to the BHTE. The error is evaluated in each plane normal to the z direction. The rms error and the normalized peak error are both under 2%.

In a 16cmX16cmX17cm computational volume with 217X217X227 field points, computing five source planes using the fast nearfield method takes 1.6 hours. Reconstructing the 3D fields using angular spectrum approach takes 2 seconds. The error is 4% compared with the reference computed with the FNM. The errors of temperature response calculated by the bio-heat transfer equation with different thermal parameters are consistently smaller than 2%.

Conclusion

Simulations of the three-dimensional pressure field and the temperature response generated by large ultrasound phased arrays are very time-consuming. By combining the fast near-field method and the angular spectrum approach, the simulation time is dramatically reduced. In the modeling of a spherical phased array of 11845 circular sources with 10 million observation points, the combined approach is 56 times faster than the direct integral approach, while the pressure result and the temperature response only differ by 4% and 2% compared with the reference.

Reference

[1] R. J. McGough, J. F. Kelly and T. V. Samulski, An efficient grid sectoring method for calculations of the nearfield pressure generated by a circular piston. J. Acoust.Soc.Am 115(5), 1942-1954 (2004).

[2] G. T. Clement, K. Hynynen, Field characterization of the therapeutic ultrasound phased arrays through forward and backward planar projection. J. Acoust.Soc.Am 108(1), 441-446 (2000).

[3] P. T. Christopher and K .J. Parker, New approaches to the linear propagation of acoustic fields. J. Acoust.Soc.Am 90(1), 507-521 (1991).

[4] P. R. Stepanishen and K. C. Benjamin, Forward and backward projection of acoustic fields using FFT method. J. Acoust.Soc.Am 71(4), 803-812 (1982).

[5] J. W. Goodman, Introduction to Fourier optics, 2nd edition. McGraw-Hill, 1988.

[6] P. Wu, R. Kazys and T. Stepinski, Analysis of the numerically implemented angular spectrum approach based on the evaluation of two-dimensionalacoustic fields. Part I. Errord due to the discrete Fourier transform and discretization, J. Acoust.Soc.Am 99(3), 1339-1348 (1996).

[7] P. Wu et al, -, Part II Characteristics as a function of angular range. J. Acoust.Soc.Am 99(3), 1349-1359 (1996).