Numerical methods for nonlinear wave propagation in ultrasound

The intensities associated with the propagation of diagnostic and therapeutic ultrasound pulses are large enough to require a nonlinear description. As a nonlinear wave propagates it distorts, creating harmonics and eventually acoustic shocks. Harmonics can be used to generate images with improved spatial resolution and less clutter. The energy from nonlinear waves is deposited in a different way than in the linear case which modifies predictions for in situ acoustic exposure. Tissue heating and radiation force depend on this intensity. High intensity shock waves are essential for stone communition with lithotripsy because it depends on the shear gradients caused by the pressure differentials and on the peak negative pressures for cavitation. The work presented in this dissertation investigates numerical simulations that solve nonlinear ultrasonic wave propagation in both the strongly nonlinear regime, where shocks develop, and the weakly nonlinear regime, where the acoustic attenuation prevents the formation of pressure discontinuities. The Rankine-Hugoniot relation for shock wave propagation describes the shock speed of a nonlinear wave. This dissertation investigates time domain numerical methods that solve the nonlinear parabolic wave equation, or the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation, and the conditions they require to satisfy the Rankine-Hugoniot relation. Two numerical methods commonly used in hyperbolic conservation laws are adapted to solve the KZK equation: Godunov's method and the monotonic upwind scheme for conservation laws (MUSCL). It is shown that they satisfy the Rankine-Hugoniot relation regardless of attenuation. These two methods are compared with the current implicit solution based method. When the attenuation is small, such as in water, the current method requires a degree of grid refinement that is computationally impractical. All three numerical methods are compared in simulations for lithotripters and high intensity focused ultrasound (HIFU) where the attenuation is small compared to the nonlinearity because much of the propagation occurs in water. The simulations are performed on grid sizes that are consistent with present-day computational resources but are not sufficiently refined for the current method to satisfy the Rankine-Hugoniot condition. It is shown that satisfying the Rankine-Hugoniot conditions has a significant impact on metrics relevant to lithotripsy (such as peak pressures), and HIFU (intensity). Because the Godunov and MUSCL schemes satisfy the Rankine-Hugoniot conditions on coarse grids they are particularly advantageous for three dimensional simulations. The propagation of focused and intense ultrasound beams is determined by nonlinearity, diffraction, and absorption. Most descriptions of nonlinear wave propagation in ultrasound, such as the KZK equation, rely on quadratic nonlinearity. At diagnostic and some therapeutic amplitudes the quadratic, or B/A, term dominates the nonlinear term. However, when the amplitudes are sufficiently large, such as in shock wave lithotripsy, the cubic, or C/A, term becomes significant. Conventionally the parabolic wave equation has only included the quadratic terms. This dissertation establishes a time domain numerical method that solves the parabolic wave equation with cubic nonlinearity in an attenuating medium. The differences between solutions of the quadratic and cubic equations for a focused lithotripter in a water bath are investigated. A study of numerical solutions to the linear full-wave equation and the KZK or parabolic wave equation is presented. Finite difference time domain methods are used to calculate the acoustic field emitted from a diagnostic ultrasound transducer. Results are compared to Field II, a simulation package that has been used extensively to linearly model transducers in ultrasound. The simulation of the parabolic equation can accurately predict the lateral beamplot for large F-numbers but exhibits errors for small F-numbers. It also overestimates the depth at which the focus occurs. It is shown that the finite difference solution of the full-wave equation is accurate for small and large F-numbers. The lateral beamplots and axial intensities are in excellent agreement with the Field II simulations. For these reasons the KZK equation is abandoned in favor of the full-wave equation to describe nonlinear propagation for ultrasound imaging. A full-wave equation that describes nonlinear propagation in a heterogeneous attenuating medium is solved numerically with finite differences in the time domain (FDTD). Three dimensional solutions of the equation are verified with water tank measurements of a commercial diagnostic ultrasound transducer and are shown to be in excellent agreement in terms of the fundamental and harmonic acoustic fields, and the power spectrum at the focus. The linear and nonlinear components of the algorithm are also verified independently. In the linear non-attenuating regime solutions match simulations from Field II to within 0.3 dB. Nonlinear plane wave propagation is shown to closely match results from the Galerkin method up to four times the fundamental frequency. In addition to thermoviscous attenuation we present a numerical solution of the relaxation attenuation laws that allows modeling of arbitrary frequency dependent attenuation, such as that observed in tissue. A perfectly matched layer (PML) is implemented at the boundaries with a novel numerical implementation that allows the PML to be used with high order discretizations. A -78 dB reduction in the reflected amplitude is demonstrated. The numerical algorithm is used to simulate a diagnostic ultrasound pulse propagating through a histologically measured representation of human abdominal wall with spatial variation in the speed of sound, attenuation, nonlinearity, and density. An ultrasound image is created in silico using the same physical and algorithmic process used in an ultrasound scanner: a series of pulses are transmitted through heterogeneous scattering tissue and the received echoes are used in a delay-and-sum beamforming algorithm to generate images. The resulting harmonic image exhibits characteristic improvement in lesion boundary definition and contrast when compared to the fundamental image. We demonstrate a mechanism of harmonic image quality improvement by showing that the harmonic point spread function is less sensitive to reverberation clutter. Numerical solutions of the nonlinear full-wave equation in a heterogeneous attenuating medium are used to simulate the propagation of diagnostic ultrasound pulses through a measured representation of the human abdomen with heterogeneities in speed of sound, attenuation, density, and nonlinearity. Conventional delay-and-sum beamforming is used to generate point spread functions (PSF) from a point target located at the focus. These PSFs reveal that, for the particular imaging system considered, the primary source of degradation in fundamental imaging is due to reverberation from near-field structures. Compared to the harmonic PSF the mean magnitude of the reverberation clutter in the fundamental PSF is 26 dB higher. An artificial medium with uniform velocity but unchanged impedance characteristics is used to show that for the fundamental PSF the primary source of degradation is phase aberration. Ultrasound images are created in silico and these beamformed images are compared to images obtained from convolution of the PSF with a scatterer field to demonstrate that a very large portion of the PSF must be used to accurately represent the clutter observed in conventional imaging. Conventional delay-and-sum beamforming is used to generate images of an anechoic lesion located beneath the abdominal layer for various transducer configurations. Point spread functions (PSF) and estimates of the contrast to noise ratio (CNR) are used to quantify and determine the sources of improvement between harmonic and fundamental imaging. Simulations indicate that reducing the pressure amplitude at the transducer surface has no discernible effect on image quality. It is shown that when the aperture is reduced there is an increase in the image degradation due to reverberation clutter in the fundamental and an increase in the effects of reverberation and phase aberration in the harmonic. A doubling of the transmit frequency shows that the harmonic lesion CNR becomes worse than the fundamental CNR due to increases in pulse lengthening and phase aberration. Acoustic Radiation Force Impulse (ARFI) imaging uses brief, high intensity, focused ultrasound pulses to generate a radiation force that displaces tissue. Nonlinear propagation of acoustic pulses transfers energy to higher frequencies where it is preferentially absorbed by tissue. The radiation force is proportional to the absorbed energy. This dissertation examines the effects of nonlinearity on the displacements induced by radiation force with various ultrasound transducer configurations. A three dimensional numerical method that simulates nonlinear acoustic propagation is used to calculate the intensity and absorption losses for typical ARFI pulses. It is demonstrated that nonlinearity has a relatively small effect on the intensity but increases estimates of the loss by up to a factor of 20. The intensity fields obtained from the acoustic simulations are used as an input to a finite element method (FEM) model of the mechanical tissue response to a radiation force excitation. These simulations show that including nonlinearity in the acoustic intensity significantly reduces predictions of the displacement without having a significant impact on the lateral and elevation resolution.Dissertatio