Nonlinear acoustics plays an important role in both diagnostic and therapeutic applications of biomedical ultrasound and a number of research and commercial software packages are available. In this manuscript, predictions of two solvers available in a commercial software package, pzflex, one using the finite-element-method (FEM) and the other a pseudo-spectral method, spectralflex, are compared with measurements and the Khokhlov-Zabolotskaya-Kuznetsov (KZK) Texas code (a finite-difference time-domain algorithm). The pzflex methods solve the continuity equation, momentum equation and equation of state where they account for nonlinearity to second order whereas the KZK code solves a nonlinear wave equation with a paraxial approximation for diffraction. Measurements of the field from a single element 3.3 MHz focused transducer were compared with the simulations and there was good agreement for the fundamental frequency and the harmonics; however the FEM pzflex solver incurred a high computational cost to achieve equivalent accuracy. In addition, pzflex results exhibited non-physical oscillations in the spatial distribution of harmonics when the amplitudes were relatively low. It was found that spectralflex was able to accurately capture the nonlinear fields at reasonable computational cost. These results emphasize the need to benchmark nonlinear simulations before using codes as predictive tools.

## I. INTRODUCTION

In medical ultrasound, nonlinear acoustic propagation effects are present in both diagnostic and therapeutic applications.^{1,2} In diagnostic ultrasound, harmonics generated during propagation can be used for high resolution imaging but the presence of harmonics means derating water measurements based on a single absorption value is not robust.^{3} In therapeutic applications, due to the higher attenuation of acoustic energy at higher frequencies, the generation of harmonics can also result in additional acoustic absorption in tissue, which can alter both the rate and the spatial distribution and of heat deposition, and have the potential to enhance heating but also impact patient safety.^{2}

The fundamental equations describing sound propagation in a fluid are conservation of mass, momentum and the equation of state. For small, but finite-amplitude fluctuations, a second order wave equation can be derived which accounts for arbitrary wave fields, nonlinearity and diffusivity.^{4} In cases where the Lagrangian density is negligible, the equation can be simplified to the Westervelt equation.^{5} For cases of sound beams propagating predominantly in one direction, the Westervelt equation can be reduced to the Khokhlov-Zabolotskaya-Kuznetzov (KZK) equation,^{6,7} which has been shown to agree well with hydrophone measurements of highly nonlinear ultrasound fields.^{8} The KZK equation can be written in the following form:^{7}

where *p* is pressure, z is the main direction of sound propagation, *τ = t − z/c _{0}* is retarded time,

*c*is small-signal sound speed, $\u2207\u22a52=\u22022/\u2202x2+\u22022/\u2202y2$ is the Laplacian operating in the plane perpendicular to the z axis,

_{0}*β*is the coefficient of nonlinearity,

*ρ*is density, and

_{0}*δ*is diffusivity of sound.

Several algorithms have been developed to solve the KZK equation. Aanonsen and Tjøtta^{9} employed a Fourier series expansion and solved the equation in the frequency domain. Mixed time frequency domains approaches have also been employed including capturing discontinuities.^{10} Lee and Hamilton^{11} developed a time-domain method for the propagation of pulsed finite-amplitude waves. The form given in Eq. (1) is for homogeneous fluids with frequency squared attenuation, however, other attenuation laws can be employed, for relaxation processes,^{12} and terms accounting for inhomogeneous medium added.^{13} One potential weakness of the KZK equation is that the parabolic approximation is considered appropriate for cases where the propagation is within about 20°,^{14} a condition that is violated by many therapeutic sources. Nevertheless, many researchers have used the KZK equation to model transducers with lower f-numbers and achieved good agreement with hydrophone measurements within the on-axis and focal region.

In order to capture full diffraction a number of approaches have been reported. One approach is to relax the parabolic approximation either in the frequency domain^{15} or the time domain.^{16} A second is to solve the Westervelt equation directly.^{17} A third approach is to solve the underlying nonlinear hydrodynamic equations, for example, with high order finite-difference time domain methods^{18} or a k-space approach.^{19} Full diffraction models typically have very high computational costs but modelling some nonlinear ultrasound problems can be tractable with the help of parallel and cluster computing.^{20,21}

Finite element methods (FEM) have also been adopted to simulate nonlinear wave propagation and they are well suited to handling complex geometries such as occur in the body.^{22,23} The commercial FEM software pzflex (Weidlinger Associates, Inc., Mountain View, CA) is one such full solver that is widely used in industrial and biomedical applications of ultrasound, because it provides a convenient platform with a graphical user interface and optimised computing kernel. The developers demonstrated that the solution for a focused high intensity focused ultrasound (HIFU) source matches the growth of harmonics to within 10%,^{24} and that an alternative pseudo-spectral method can model nonlinear pulses from imaging transducers.^{25} We are unaware of any independent assessments of the accuracy of the simulation results for nonlinear sound propagation.

In this paper, we compare measurements of the nonlinear wave propagation from a single-element HIFU transducer to simulations using a time-domain solver of the KZK equation and the two solvers available in the pzflex software. The main motivation is to independently validate spectralflex and pzflex for modelling therapeutic ultrasound.

## II. METHOD

The KZK equation solver used in the present work is the time-domain axisymmetric “Texas” code developed by Lee,^{11} which uses an operator splitting scheme to solve each part of the equation—attenuation, diffraction, and nonlinearity—in separate steps using the finite difference method. The KZK equation solver marches in the z direction and only needs to update *p*(*r, t*) at each step.^{11} For the simulations shown here 270 temporal-points per cycle and a temporal duration of 20 cycles were employed. In the radial direction 2500 points were used across the transducer and the outer boundary was placed at 4 times the transducer radius. These values lead to a grid of 5400 × 10 000= 54 000 000 points. The code was initialised by specifying *p*(*r, t*) in the source plane, z = 0, with focusing accounted for by time-shifting the waveform with radial distance. The marching step in the z-direction was fixed at 5 × 10^{−4} relative to the focal length, that is, there were 2000 steps from the source to the focus.

Rather than solving a wave equation both pzflex and spectralflex solve the underlying conservation laws and constitutive equations. The form of the momentum equations employed is^{26}

where u is the particle displacement. The equation of state employs a second order expansion of the Taylor series relationship between pressure and density,^{25}

where K is the bulk modulus and B/A the parameter of nonlinearity. Although this approximation will result in errors that accrue through the solution,^{25} it has been suggested that propagation on the order of hundreds of wavelengths is practical for pzflex and that spectralflex should be much less sensitive to dispersion errors.^{23,24,26}

The propagation of a 1D plane wave was simulated in pzflex first to find the optimal spatial discretisation for the simulation. The analytical solution for finite amplitude plane wave propagation was used as the benchmark. In the 1D plane wave model, the velocity source with amplitude 1 m/s (equivalent to 1.5 MPa) at 5 MHz was set as the boundary condition. The corresponding shock formation distance in water is *x*_{shock} = 20.5 mm and the computational domain extended to *σ = x/x*_{shock} = 3. The developers of pzflex suggest that at least 15 elements per wavelength are necessary for “negligible” numerical dispersion^{23} and here the value was varied from 60 (which should give sufficient accuracy at the third harmonic) to 240 points per wavelength.

In the pzflex FEM model, the focused transducer was modelled as a curved surface with a normal velocity condition. Hexahedral elements with eight nodes were used for regularly meshing the model. Because pzflex uses an explicit algorithm which is conditionally unstable, the time step must be smaller than the CFL limit. In this paper, the value of Δ*t* was chosen to be 80% of the critical value based on the CFL limit for stability. The models were built in axisymmetric coordinates and the total number of grid points depended on the transducer being modelled and on its driving frequency. The calculations were carried out using a Linux workstation (CPU: Intel E5-2650@2.6 GHz, RAM: 189.1 GB).

The other pzflex package, spectralflex, was also used to model nonlinear wave propagation for comparison.^{25,27} This package utilizes a pseudospectral algorithm which applies Fourier transforms over the spatial domain to simplify the calculation of the spatial derivatives of the functions. For the rest of this manuscript, “spectralflex” refers to the pseudo-spectral solver, whilst pzflex will be used to refer to the FEM solver only. The input to spectralflex was defined as pressure condition along a radial line in the source plane and rotated to give an axisymmetric coordinates with the same dimensions as the source used in the KZK solver. The focusing effect was introduced by applying a phase shift to the input waveform of each grid node. The pseudospectral method results in less numerical dispersion error^{28} and here 30 elements per wavelength at the fundamental frequency was used. A perfect matching layer boundary condition was applied around the computational domain to reduce the wraparound error caused by periodicity of the fast Fourier transform.

The HIFU transducer investigated had a focal length of 62.7 mm, an outer diameter of 64 mm (f-number = 0.98) and a central opening of diameter of 20 mm and here was driven at 3.3 MHz (Model H102B-35, Sonic Concepts Ltd., Woodinville, WA). For the simulations the nominal diameter was changed to 62 mm for the outer diameter as it resulted in a better match to the experimental measurements for both the gain and the position of the nulls. The pressure field was measured with a Fabry-Perot Optical Hydrophone (FPOH, Precision Acoustics, Ltd., Dorchester, UK) using the setup shown in Fig. 1. The therapeutic transducer was driven by a 20 cycle sinusoidal waves from a signal generator (33220A, Agilent Technologies, Santa Clara, CA) which was amplified by 55 dB (A300, E&I, Inc., Rochester, NY) and applied to the transducer by means of a matching circuit provided by the manufacturer. The hydrophone was attached to a three-dimensional positioning system and pressure waveforms measured at various locations in the field.

To obtain the correct input value for the simulations, the pressure fields of the transducers driven at low input voltage were first scanned by the hydrophone. The input values used by both codes were adjusted so that the simulations in the linear regime matched the experimental results. Measurements were then taken at a higher source voltage, where measurable harmonics began to appear in the spectrum of the waveform, and the input values for the numerical solvers were scaled by the same factor. The medium in all simulations was taken to be water with small-signal sound speed *c* = 1500 m/s, density *ρ* = 1000 kg/m, frequency square absorption $\alpha 1(f/1\u2009MHz)2$, where the attenuation at 1 MHz was $\alpha 1=0.0253\u2009Np/m$ and the coefficient of nonlinearity *β* = 3.5.

## III. RESULTS

The results of simulations of finite amplitude plane wave were normalized to the initial pressure in each case and depicted in Fig. 2. The pzflex results fit the analytical solution well until the propagation distance exceeds the shock formation distance. Beyond the shock formation distance there is generation of large numbers of harmonics and although at 240 points per wavelength pzflex should be able to model up to 16 harmonics this is not sufficient and errors induced in the higher harmonics affect the lower frequency components. These data suggest that pzflex is not well suited of modelling highly shocked waveforms. However, for the HIFU fields investigated in this study shock formation was not induced and therefore we define the error only for distances *σ < 1.* The root mean square error (RMSE), normalised to the analytical solution, is shown in Fig. 2(b) as a function of the number of grid points. It can be seen that once the number of grid points is greater than 120 points per wavelength, the RMSE of the first three harmonics are all less than 0.2%. Because of the limitation of the available memory the maximum number of points per wavelength for modelling of the HIFU transducer was 160 but based on the results in Fig. 2 this should be sufficient to capture waveforms up to the shock formation distance.

Figure 3 shows the axial pressure distribution of the transducer at low pressure as measured by the hydrophone and predicted by the lossless analytical O'Neil solution^{16,29} and pzflex. There is good agreement in the location of the nulls for the transducer. The central hole in the transducer was suitably captured in the linear analytical solution by subtracting the field of a transducer of the same size as the hole from the solution of an intact transducer and in the simulations by setting the source excitation to zero within the hole. For the KZK solver there were discrepancies with the exact solution in the near field (data not shown), which is expected due to the paraxial approximation, but good agreement with the linear Fresnel solution which employs the same paraxial approximation.

The transmit settings for the HIFU transducer were increased such that the pressure amplitude was close to 7 MPa and harmonic signals were apparent. Figure 4 shows the measurements and simulations of the fundamental, second harmonic and third harmonic. All the simulation methods captured the amplitude and profile of the fundamental frequency well. The KZK and spectralflex results showed good agreement with the measurements of the second and third harmonics as well. In comparison pzflex required 120 or 160 points per wavelength to capture the second harmonic correctly but even at 160 points per wavelength discretisation the third harmonic was more than 10% less than predicted by the other simulations. We note that using pzflex with 60 points per wavelength, which should be sufficient for the spatial discretisation according to the manual, the amplitude of the second harmonic is 62% of the other simulations and third harmonic only 33%. Although increasing discretisation improved the results of pzflex it came at a large cost, with a grid density of 160 points per wavelength requiring 155 GB memory and a computation time of almost 10 days, while the KZK solver required a few gigabytes and less than 10 h to finish the simulation. Finally, it can be seen that the output of pzflex has ripples in the focal region of the harmonics which are not expected.

To further investigate the performance of pzflex, the same geometry transducer was driven at 1, 1.5, 2, and 2.5 MHz was then modelled with pzflex at different grid densities and compared to the KZK code. A comparison of the axial distribution of the first three harmonics for each source frequency is shown in Fig. 5. The pzflex solutions show good agreement with the KZK solutions at the fundamental frequency, while the amplitude of the higher harmonics calculated by pzflex is lower than those obtained from KZK code. However, the non-physical oscillations alluded to above are more pronounced, increase at a lower frequency, and appear to be independent of grid spacing. For example, for the 1 MHz simulation in Fig. 5(a), the curves are almost identical suggesting that increasing the points per wavelength does not alleviate the oscillations. The effect of changing CFL was considered by varying CFL from 0.2 to 0.8 for the 1 MHz case. No change in the oscillation was observed (data not shown) suggesting that the effect is independent of time stepping. As the frequency is increased the oscillations are less pronounced, however, the amplitude of the harmonics decreases. Figure 6 shows the harmonic amplitude as a function of frequency and at higher frequencies the results diverge from the KZK solutions even for 160 points per wavelength.

The two-dimensional (2D) pressure fields calculated by KZK and pzflex are shown in Fig. 7 for the 1 MHz case. The field patterns of the fundamental component agree well for both numerical methods both on and off axis, with some discrepancies are in the near field where the paraxial approximation of the KZK equation is expected to break down. For the second harmonic, the general shape of the focus is the same for both methods, however, ripples appear in the pzflex results at 1 MHz which is consistent with the oscillations shown in Fig. 5. The 2D distribution of the third harmonic in the 1 MHz case calculated by pzflex shows the most pronounced ripple artefact and its distribution mirrors the shape of the fundamental component with lower amplitude, including the side lobe structure. This suggests that the oscillations are an artefact of signal from the fundamental bleeding through in to the harmonics.

## IV. DISCUSSION

The KZK equation has been widely used to model nonlinear acoustic beams including diagnostic and therapeutic transducers and the results here are consistent with prior work suggesting that the KZK equation can capture the acoustic field of the focal area for some tightly focused transducers. For example, the transducer investigated here had a linear gain of 102.3 and the agreement was good for the three harmonics considered. The results of the spectralflex were also in good agreement with the measurements up to the third harmonic, even with only 30 grid points per wavelength, which is consistent with the validation study done by Wojcik.^{23} The pzflex FEM solver was also able to capture the distributions of harmonics of the HIFU pressure field fairly well when a dense spatial discretization with 160 points per wavelength was applied, but at a relatively high computational cost: the axisymmetric simulation took 155 GB memory and 10 days of run time. Employing pzflex for a full 3D problem is likely prohibitively computationally expensive for the foreseeable future.

Spatial discretization had a larger effect on the calculated harmonics in pzflex than for the KZK solver. For pzflex using a discretisation of 60 points per wavelength resulted in more than a 65% reduction in the amplitude of the third harmonic where as an equivalent simulation with the KZK solver resulted in a difference of less than 10% (data not shown). Changing the number of radial points in the KZK simulation also had a minimal effect on the results.

Another issue with pzflex was that it predicted non-physical oscillations in the distributions of the harmonics and this phenomenon also occurred for spectralflex (data not shown). The oscillations were more pronounced in higher harmonics when the source pressure and frequency were low, conditions of low focusing and weak nonlinear generation of harmonics. However, as nonlinearity becomes more significant, e.g., by increasing the source level or driving frequency, the relative amplitude of the oscillation was reduced. This suggests that the oscillation originates from numerical errors which dominate when the real nonlinear signal is relatively low. Further evidence of this phenomenon was that the 2D distribution of the third harmonics of the 1 MHz simulation shows basically the same shape of the fundamental component with ripple patterns. This similarity implies that a small amount of energy from the fundamental is leaking into the higher harmonics in error and then causing the oscillations. The oscillations were not sensitive to changes in the CFL or the grid discretisation. Provided the same embedded function was used to derive the spectra from the time series for both spectralflex and pzflex, then we speculate that the numerical error resulting in the oscillations could have resulted from the embedded frequency analysis.

Capturing the levels of the harmonics is important in biomedical ultrasound for both imaging and therapy. In diagnostic ultrasound the second harmonic is used to form images and in order to understand and optimise the second harmonic its accurate prediction is therefore necessary. In therapeutic ultrasound the higher harmonics are absorbed more readily and result in enhanced temperature rise and spatial shift compared to the results of linear models.^{8,33} In order to use simulation tools for HIFU treatment planning, the amplitude of the harmonics need to be accurate in order to determine temperature in the tissue. For example, underestimating the heating rate could result in the used of higher source levels which might bring other unexpected phenomena, such as cavitation, boiling or undesirable pre-focal heating.

A key strength of pzflex lies in modelling the piezoelectric response and mechanical behaviour of different types of materials and ultrasonic transducers. Several groups have used pzflex for design and characterisation of ultrasonic devices^{21,30,31} and corresponding comparisons between pzflex and experiments were done to validate the performance of the software on piezoelectric modelling.^{30,32} The agreement we see with the fundamental frequency is consistent with these prior results. However, there is little validation of pzflex for modelling nonlinear wave propagation. Wojcik *et al*.^{25} used the FEM in pzflex for the model of the transducer dynamics, but then used spectralflex for the subsequent propagation, which showed good agreement with the experiments. The results in the paper suggest that the FEM solver of pzflex can simulate the nonlinear propagation well with a dense spatial discretization in the absence of shocks, but in the term of the computational efficiency, spectralflex outperforms the FEM of pzflex on modelling nonlinear propagation, given the accuracy provided by the pseudospectral method with fewer grid points.^{28}

In conclusion, results from three simulations packages have been considered for HIFU applications. Although theoretically, the KZK equation is only valid for focused transducers with a half-aperture angle of less than 16°, many papers use the KZK solver to model transducers with lower f-numbers and achieve good agreement with hydrophone measurements.^{33} Our results with the KZK solver are consistent with the literature, confirming that it captures the nonlinear pressure field from therapeutic transducers well enough to make it suitable for HIFU treatment planning. The KZK Texas code required fewer computational resources to run than either spectralflex or pzflex, however, the KZK equation cannot capture complex media where multiple reflections occur for example, near bones or gas pockets in the body. In comparison, pzflex has no restrictions on geometry and it was found here that it can be used to do linear acoustic simulations with good accuracy but did not correctly model the nonlinearity of ultrasound propagation for the cases we considered with an acceptable computational efficiency. The spectralflex algorithm was able to capture linear and nonlinear fields well, with a caution in interpreting the harmonics, and our results suggest it is well suited to modelling nonlinear propagation in the body where significant reflections, such as with bone or gas bodies are present.

## ACKNOWLEDGMENTS

The authors would like to thank Dr. Christian Coviello for his help with pzflex software. This work was supported by the UltraSpine award from the UK Engineering and Physical Sciences Research Council (Grant No. EP/K020757/1).