Time-domain simulation of ultrasound propagation in a tissue-like medium based on the resolution of the nonlinear acoustic constitutive relations

Time-domain simulation of ultrasound propagation in a tissue-like medium based on the resolution of the nonlinear acoustic constitutive relations

Noé Jiménez Instituto de Investigación para la Gestión Integrada de Zonas Costeras, Universitat Politècnica de València, Paranimf 1, 46730 Grao de Gandia, Spain    Francisco Camarena Instituto de Investigación para la Gestión Integrada de Zonas Costeras, Universitat Politècnica de València, Paranimf 1, 46730 Grao de Gandia, Spain    Javier Redondo Instituto de Investigación para la Gestión Integrada de Zonas Costeras, Universitat Politècnica de València, Paranimf 1, 46730 Grao de Gandia, Spain    Víctor Sánchez-Morcillo Instituto de Investigación para la Gestión Integrada de Zonas Costeras, Universitat Politècnica de València, Paranimf 1, 46730 Grao de Gandia, Spain    Yi Hou Department of Biomedical Engineering, Columbia University, 351 Engineering Terrace, mail code 8904, 1210 Amsterdam Avenue, New York, NY, USA    Elisa E. Konofagou Department of Biomedical Engineering, Columbia University, 351 Engineering Terrace, mail code 8904, 1210 Amsterdam Avenue, New York, NY, USA [
The present paper is a preprint, submited to J. Acoust. Soc. Am. in June 2015

A time-domain numerical code based on the constitutive relations of nonlinear acoustics for simulating ultrasound propagation is presented. To model frequency power law attenuation, such as observed in biological media, multiple relaxation processes are included and relaxation parameters are fitted to both exact frequency power law attenuation and empirically measured attenuation of a variety of tissues that does not fit an exact power law. A computational technique based on artificial relaxation is included to correct the non-negligible numerical dispersion of the numerical method and to improve stability when shock waves are present. This technique avoids the use of high order finite difference schemes, leading to fast calculations. The numerical code is especially suitable to study high intensity and focused axisymmetric acoustic beams in tissue-like medium, as it is based on the full constitutive relations that overcomes the limitations of the parabolic approximations, while some specific effects not contemplated by the Westervelt equation can be also studied. The accuracy of the method is discussed by comparing the proposed simulation solutions to one-dimensional analytical ones, to -space numerical solutions and also to experimental data from a focused beam propagating in a frequency power law attenuation media.

43.58.Ta, 43.80.Sh, 43.35.Wa, 43.35.Fj, 43.25.Ts

Also at: ]Department of Radiology, Columbia University, 351 Engineering Terrace, mail code 8904, 1210 Amsterdam Avenue, New York, NY, USA

I Introduction

The significant development of ultrasound technology in the medical field in recent years has led to the need for simulation tools increasingly realistic. Effects like absorption in biological media, nonlinear propagation, heterogeneities, strong focusing, streaming, resonances, multiple scattering or the presence of discontinuities due to tissue layers or rigid boundaries have to be taken into consideration. The most general approach for ultrasound simulation is to directly solve the constitutive relations of the nonlinear acoustics. It also allows the explicit calculation of the particle velocity, what can be used to compute important magnitudes as the vector components of the nonlinear acoustic intensity or the acoustic radiation force.

The numerical resolution of the nonlinear constitutive equations in tissue-like medium supposes a difficult problem due to the large size of the region of interest in relation to the size of the acoustic wavelength and the complexity of the model. Simplifying assumptions have been needed in the past for modeling beam patterns from ultrasound transducers, as one-way parabolic approximations, most based on the Khokhlov-Zabolotskaya-Kuznetsov equation (KZK) Aanonsen et al. (1984); Lee and Hamilton (1995); Cleveland, Hamilton, and Blackstock (1996); Yang and Cleveland (2005); Khokhlova et al. (2006); Jing and Cleveland (2007); Soneson and Myers (2007); Prieur and Holm (2011). To overcome the validity limitations of the parabolic approximations, i.e. for large aperture focused sound sources or modeling sound field near the acoustic source, many one-way numerical methods have been proposed, including phenomenological approaches Christopher and Parker (1991); Tavakkoli et al. (1998) with tissue attenuation Zemp, Tavakkoli, and Cobbold (2003), or based on the one-way formulation of the Westervelt equation Varslot and Taraldsen (2005); Yuldashev and Khokhlova (2011).

Tissue inhomogeneity can be modeled in these one-way models Jing and Cleveland (2007), like transmission though tissue layers with refraction, but they do not take into account backscattering and multiple reflections. More realistic models, e.g. those accounting for scattering from internal tissue structures, are based on the Westervelt-type full-wave equations Hallaj and Cleveland (1999); Hallaj, Cleveland, and Hynynen (2001); Pinton et al. (2009); Huijssen and Verweij (2010); Demi, Van Dongen, and Verweij (2011); Verweij, Demi, and van Dongen (2013); Jing, Wang, and Clement (2012). This full-wave equation has been validated for strongly focused sources Jing, Shen, and Clement (2011). However, due to the assumptions taken in the derivation of the Westervelt equation, the accuracy of this model is limited in practical situations as () the modeling of rigid boundaries where the thermo-viscous boundary layer effects are not-negligible, i.e. in general case where the particle velocity field becomes rotational, () situations where the second order Lagrangian density of acoustical energy not vanish, i.e. near the source or in situations where plane progressive waves does not exist and the acoustic field becomes complex due to multiple scattering, reverberation or resonances, () situations where the equilibrium-state particle velocity is not null, including the self generation of acoustic streaming. See Ref. Hamilton and Blackstock (1998) Chap. 3 for further discussion.

The recent development of computational capacity has made possible to consider the full constitutive relations (i.e. without the assumptions discussed above). Thus, for small-amplitude acoustic waves, the linearized pressure-velocity formulation of constitutive relations in inhomogeneous media was solved by means of Finite-Differences in Time-Domain (FDTD) methods with frequency independent losses by Manry et al. Manry and Broschat (1996), or using two-step MacCormack finite-differences scheme by Mast et al. Mast et al. (1997). Also, relaxation processes can be included in finite difference methods in an efficient way in order to model tissue attenuation and dispersion Yuan et al. (1999). Furthermore, -space numerical methods have also been applied to solve the linearized first order equations in lossless inhomogeneous media Mast et al. (2001). In order to account for soft tissue losses, the computational solution of the fractional Laplacian by -space spectral methods have demonstrated to be extremely efficient due to the spatial frequency domain representation of the acoustic field Treeby and Cox (2010).

In the case of nonlinear constitutive relations models, the evolution of the acoustic magnitudes have been simulated in time-domain by means of finite differences schemes such as Dispersion Relation Preserving method (DRP) in ideal fluids and axisymmetric domains Ginter et al. (2002). Thermo-viscous losses in finite-differences methods have been widely used, see Sparrow et al. Sparrow and Raspet (1991). In order to introduce tissue attenuation in the governing equations time-dependent fractional derivatives can be included by convolutional operators. Thus, in Ref. Liebler et al. (2004) an efficient method has been presented, but although the memory requirements can be strongly reduced compared to direct convolutions the algorithm employs up to ten auxiliary fields and a memory buffer of three time steps. Furthermore, construction of specific causal memory functions that models soft tissue attenuation and dispersion in Navier-Stokes equations is also possible Lobanova, Lobanov, and Khokhlova (2014), but certain time history must be stored in memory and in this case the computational domain was restricted to one dimensional propagation.

In order to overcome those numerical limitations, recently -space and pseudo-spectral numerical methods have been applied to constitutive relations in nonlinear regime to solve fractional Laplacian operators efficiently Treeby et al. (2012). Furthermore, in the case of domains of hundreds of wavelengths, when the cumulative phase error due to numerical dispersion of standard finite-difference schemes can not be neglected, those spectral numerical methods have reported an improvement in accuracy of the numerical solution. This two factors, i.e. the negligible numerical dispersion and the efficient resolution of fractional Laplacian operators, have led the spectral methods to be widely used in practical applications. However, their main limitation is that the implementation of natural space discontinuities due to tissue layers or rigid boundary conditions leads to errors in the reconstruction of the spectral information due to the poor convergence of Fourier series at jumps, i.e. the well-known Gibbs oscillations. Preventing this kind of errors is typically achieved by filtering the spatial spectrum Jing, Wang, and Clement (2012), so the theoretical spatial minimum sampling of two point per wavelength becomes larger. In addition, these errors propagate globally and affect to the accuracy all over the domain, in contrast with locally propagating errors in finite differences methods. On the other hand, taking into account the spatial discontinuity due to symmetry boundary condition, axisymmetric domains becomes not feasible by standard -space methods, and full 3D domains must be employed even for axisymmetric configurations. Those errors can be prevented by means of the recently developed Fourier Continuation (FC) method Albin et al. (2012). However, the discontinuities formed due to shock propagation are still not solved by FC methods and other additional numerical treatment must be applied for correctly describe shock formation, e.g. intensive computations by high order accurate weighted essentially non-oscillatory schemes (FC/WENO) Shahbazi et al. (2011). Unlikely, the computational times increases by using those intensive computational techniques and the multiresolution analysis to detect discontinuities in the domain.

The aim of the present work is to present a generalization of the constitutive relations of nonlinear acoustics including multiple relaxation processes in a non-convolutional formulation that allows the time-domain numerical solution by an explicit finite differences numerical scheme. Frequency power law attenuation based in relaxation have been applied in the same way than it has been applied to generalized Burgers equation Cleveland, Hamilton, and Blackstock (1996), Khokhlov-Zabolotskaya-Kuznetsov (KZK) model Cleveland, Hamilton, and Blackstock (1996); Yang and Cleveland (2005) and Westervelt equation Pinton et al. (2009). The relaxation parameters have been fitted to both exact frequency power law attenuation and empirically measured attenuation of a variety of tissues that does not fit an exact power law. Two processes have been enough to model tissue attenuation with acceptable accuracy over a frequency range covering about 4 octaves, as it was demonstrated by Yang et al. Yang and Cleveland (2005). A numerical technique based on artificial relaxation is included to control the non-negligible numerical dispersion of the FDTD method and improve stability when shock waves are present in the solution. The method includes backscattering and arbitrary propagation direction of finite amplitude beams, and can be specially suitable in axisymmetric configurations where the computational resources for full 3D -space methods are prohibitive.

The paper is organized as follows: in Sec. II the model equations that describes the problem are exposed, Sec. III describes the computational method presented in this work and in Sec. IV the method is validated comparing the numerical results with analytic solutions for linear, smooth and discontinuous nonlinear waves. In Sec. V solutions for frequency power law media are presented and compared with analytic and numerical solutions obtained by -space methods as benchmark case. Furthermore, an experimental test was presented where a focused beam propagating in castor oil has been modeled. Finally, a high intensity focused source in the high nonlinear regime focusing in soft-tissue is modeled.

Ii Generalized nonlinear acoustics model for multiple relaxation media

ii.1 Full-wave modeling

The principles of mass and momentum conservation lead to the main constitutive relations for nonlinear acoustic waves, which for a fluid can be expressed as Naugolnykh and Ostrovsky (1998)




where is the total density field, v is the particle velocity vector, is the pressure, and and are the coefficients of shear and the bulk viscosity respectively. The acoustic waves described by this model exhibit viscous losses with quadratic power law dependence on frequency. In order to include a power law frequency dependence on the attenuation, a multiple relaxation model will be added into the time domain equations.

The basic mechanism for energy loss in relaxing media is the appearance of a phase shift between the pressure and density fields. This behavior is commonly modeled as a time dependent relation at the fluid state equation, that for a fluid retaining the material nonlinear effects up to second order an be expressed as Naugolnykh and Ostrovsky (1998); Rudenko and Soluian (1977):


where is the density perturbation over the stationary density , is the nonlinear parameter, is the small amplitude sound speed, and is the kernel associated with the relaxation mechanism. The first two terms of the right hand side of Eq.  (3) describe the instantaneous response of the medium, where the convolutional third term accounts for the “memory time” of the relaxing media. Thus, by choosing an adequate time function for the kernel the model can present an attenuation and dispersion response that fits the experimental data of the heterogeneous media. However, the direct resolution of the constitutive relations Eqs. (1-3) in this integral form is a complex numerical task due to the convolutional operator. Thus, instead of describe with a specific time domain waveform, the response of the heterogeneous medium can be alternatively described by a sum of relaxation processes with exponential time dependence as:


with the -th order relaxation kernel expressed as


where is the Heaviside piecewise function , , is the characteristic relaxation time and the relaxation parameter for the -th order process. This last dimensionless parameter controls the amount of attenuation and dispersion for each process as , where is the sound speed in the high frequency limit associated to -th order relaxation process, also known as the speed of sound in the “frozen” state Pierce (1989). In order to describe relaxation without the need of including a convolutional operator, we shall define a state variable for each process as


Thus, using the convolutional property , the time derivative of the relaxation state variable obeys the following relation for the -th order process:


where is the Dirac delta function. Using the Eq. (6) this relation becomes a simple ordinary differential equation for each process as:


Using again convolutional properties, we can substitute Eq. (8) into Eq. (4), and the relaxing nonlinear state Eq. (3) becomes:


Moreover, if “frozen” sound speed for mechanisms is defined as , Eq. (9) leads to:


Due to the smallness of the relaxation parameter, , i.e. when weak dispersion is modeled, the sound speed in the high frequency limit reduces to Naugolnykh and Ostrovsky (1998):


Note Eq. (10) for a mono-relaxing media is equivalent to that can be found in literature Rudenko and Soluian (1977)


Thus, the constitutive relations to solve by means of the numerical method in the nonlinear regime are the continuity Eq. (1), the motion Eq. (2) and the second order fluid state relaxing Eq. (10), where the state variable obeys the relation Eq. (8) for the -th order relaxation process. Although the aim of this work is to model biological media, the generalized formulation presented here can be used to describe the attenuation and hence the dispersion observed in other relaxing media, as the relaxation processes of oxygen and nitrogen molecules in air or the relaxation associated with boric acid and magnesium sulfate in seawater Pierce (1989).

ii.2 Small amplitude modeling

If small amplitude perturbations are considered, an equivalent derivation of this model can be expressed for multiple relaxation media Yuan et al. (1999). Thus, for an homogeneous inviscid relaxing fluid, the linearized continuity and motion Eq. (1-2) reduces to




and linearizing the fluid state Eq. (10) we obtain:


These equations can be solved directly in this form, however, if expressed in pressure-velocity formulation the density field is no longer necessary and computational effort can be reduced. Thereby, assuming a linear “instantaneous” compressibility , and substituting Eq. (15) into Eq. (13) yields


Then, taking the time derivative of the state variable Eq. (8) we get


Finally, substituting again the linearized state Eq. (15) and arranging terms the linearized continuity equation leads to


On the other hand, the state evolution equation can be expressed as a function of the acoustic pressure as


Thus, the linearized governing Eq. (14, 18) for a relaxing media are expressed in a pressure-velocity formulation and can be solved together with the coupled state evolution Eq. (19) by means of standard finite differences numerical techniques Yuan et al. (1999). In this way, lossless linear acoustics equations can be obtained by setting or in the limit when the relaxation times . The relaxation behavior described by this linearized model is achieved too by the formulation described by (Yuan et al., 1999), where the relaxation coefficients and the relaxation variable are defined in a different, but analogous way.

Iii Numerical solution by finite-difference time-domain

In this section the numerical techniques for solving the complete set of equations (continuity Eq. (1), momentum Eq. (2), state Eq. (10) and the relaxation Eq. (8)) are presented. The numerical method is based on a second order FDTD method where multiple relaxation processes are included in order to: first, modeling physical attenuation and dispersion at the frequencies of interests and second, correct the numerical dispersion and include artificial attenuation to guarantee convergence in nonlinear regime. Moreover the inclusion of relaxation processes in the presented formulation require only one extra field per relaxation process and no memory buffer is needed.

iii.1 Discretization

Cylindrical axisymmetric coordinate system is considered in this work, however, the method can be derived in other coordinate systems. As in the standard acoustic FDTD method Botteldooren (1996), the particle velocity fields are discretized staggered in time and space respect to the density and pressure fields. As shown in Fig. 1 uniform grid is considered, where , , , with and as the radial and axial spatial steps, and is the temporal step.

Figure 1: Spatial staggered discretization. The pressure and the -th order relaxation process state fields are evaluated at same discrete location as the density . Particle velocity fields are discretized staggered in both space and time respect to the density, pressure and the -th order state fields.

Centered finite differences operators are applied over the partial derivatives of the governing equations. Thus, spatial interpolation is needed over the off-center grid variables in order to fulfill the conservation principles over each discrete cell of the domain LeVeque (1992). The component of Eq. (2) is expressed in a cylindrical axisymmetric system as


Each term of the above expression is approximated by centered finite differences evaluated at , , . This equation can be solved obtaining an update equation for . In the same way, the component of the motion Eq. (2) is expressed as


This equation is approximated by centered finite differences and evaluated at , , . An update equation is obtained solving this equation for . Equation (1) in cylindrical axisymmetric coordinate system is expressed as


Following the same procedure, each term of the above expression is approximated by centered finite differences and evaluated at , , , and the update equation is obtained solving this expression for . A leap-frog time marching is applied to solve Eq. (20-22) for each time step until the desired simulation time is reached. Finally, Eq. (8) is locally solved for by and explicit fourth-order Runge-Kutta method and then Eq. (10) is used for update the pressure field.

iii.2 Boundary conditions

Figure 2: Reflection coefficient of the perfectly matched layer (PML) versus layer thickness for different wave amplitudes.

The staggered grid is terminated on velocity nodes, so the boundary conditions are applied on these external nodes, allowing to prevent the singularity of the cylindrical coordinate system: due to the staggered grid, the only variable discretized at is , and axisymmetric condition is applied there. Furthermore, to solve spatial differential operators at boundaries some “ghost” nodes must be created with the conditions: , , and .

Perfectly matched layers (PML) Liu and Tao (1997) were placed in the limits of the domain ( and ) to avoid spurious reflections from the limits of the integration domain. Inside the PML domains linearized acoustic equations were solved using the complex coordinate screeching formulation Liu (1999). For a layer of 30 elements and a broadband incident wave with 1 MHz central frequency and non-normal incidence angle, these absorbent boundary conditions have reported a reflection coefficient of  dB. However, the performance of the PML is amplitude dependent as long as the nonlinear terms are uncoupled to the PML domains. The amplitude dependence of the reflection coefficient is shown in Fig. 2, where a PML of 25 layers have reported reflection coefficients  dB for waves in linear regime and highly nonlinear waves including shocks.

iii.3 Minimizing numerical dispersion

The stability for the lossless linear FDTD algorithm follows the Courant-Friedrich-Levy (CFL) condition, that for uniform grid the maximum duration of the time step is limited by where is the number of dimensions (i.e. in cylindrical axisymmetric coordinate system). That condition essentially states that for a single time step information can not propagate in the numerical grid a distance longer that one cell. However, if relaxation is included numerical instabilities have been observed when . Due to this empirical relation, the maximum values for relaxation frequencies are limited too by the chosen spatial discretization by the simple relation where is the maximum relaxation frequency for all processes, is the number of spatial samples per wavelength and the frequency of the propagating wave.

On the other hand, nonlinear effects induce the progressively growing of harmonics of the fundamental frequency of the initial wave. The diffusive viscous terms in Eq. (2), attenuates the small-amplitude high-spatial frequencies, damping the “node to node” numerical oscillations and ensuring numerical stability in weakly nonlinear regime. Thus, for a smooth solution the numerical algorithm shows consistency when , so if stability is achieved by the CFL condition, the convergence is guaranteed. However, in strongly nonlinear regime, i.e. when sharp waveforms or even shocks are present in the solution, extra numerical techniques must be employed to guarantee convergence. Artificial viscosity can be added when shock waves are present in the solution where a common implementation follows a fourth order spatial filtering Sparrow and Raspet (1991); Ginter et al. (2002). Thus, the artificial attenuation retrieved by this spatial operator is fourth power of frequency: the low frequency components of the solution remains quasi-undamped, while the higher spatial frequencies are strongly attenuated. In this way, the solution is smoothed and shock thickness depends on the artificial viscosity coefficient.

Figure 3: Normalized dispersion relation for of a FDTD lattice of Courant number and anomalous dispersion relation for a mono relaxing acoustic media. The straight line represents the reference nondispersive case.

However, the main drawback for finite difference methods is numerical dispersion, where the analytic dispersion relation can be expressed in 1D as , with the Courant number . In this way, numerical dispersion reduces phase speed for high frequency components so traveling sharp solutions develop tail oscillations: the low wavenumbers travels fast and left behind high spatial frequencies. In nonlinear regime, is well-known that the combined effects of nonlinearity and strong dispersion can lead rich phenomena, e.g. beatings in the generated harmonics, pulsations on the vertex of a sawtooth wave or soliton formation in strong dispersive media Rudenko and Soluian (1977). In this way, the numerical dispersion by discreteness of the FDTD methods couple to the physical nonlinearity can lead to a great variety of non-physical or even unstable solutions.

In order to overcome those two limitations, i.e. the generation of harmonics over the discrete limit and the numerical dispersion, we propose the use of artificial relaxation. As Fig. 3 shows, physical relaxation processes introduces anomalous dispersion, i.e. the phase speed increases in the high frequency regime, opposite to the numerical (lattice) dispersion of the finite differences scheme. Thus, by introducing a collection of relaxation processes and choosing its adequate relaxation parameters the high frequency numerical dispersion can be compensated. As a consequence, introduction of those relaxation processes in the high frequency lead to the inevitable inclusion of artificial attenuation. However, this numerical attenuation is then exploited to limit the growing of higher harmonics in a similar way than artificial viscosity Sparrow and Raspet (1991). It is worth noting here that, due to the attenuation using artificial relaxation is, at maximum, only second power of frequency, the low frequency range of the solution is therefore also attenuated. Thus, the proposed method is restricted to lossy media.

The adequate relaxation parameters that corrects the numerical dispersion have been found by multi-objective optimization techniques, where two cost function are proposed: one for dispersion and other for attenuation. In the first case, the error between the goal (ideal) dispersion relation and the retrieved numerical dispersion corrected by multiple artificial relaxation and is evaluated in the high frequency regime. Finally, the second cost function is the error between the desired (ideal) attenuation and the artificial attenuation evaluated in the low frequency regime.

Iv Validation

iv.1 Single relaxation process

Figure 4: (Color online) Pareto front retrieved by the multi-objective genetic optimization. Square marker (A) is the solution those relaxation parameters minimizes the numerical dispersion. The best fit in the attenuation are the parameters that provided the solution marked by the triangle marker (B). A compromise between attenuation and dispersion is achieved at the individuals around the center of the Pareto front, as shows the sample marked by the circle (C). The inset shows the normalized dispersion relation retrieved by the individuals (A) dotted line, (B) dashed-dotted line, (C) continuous line. Dashed black line shows the numerical dispersion relation of the FDTD method for a Courant number of 0.94 without artificial relaxation.
Figure 5: (Color online) Retrieved attenuation (top) and dispersion (bottom) by the inclusion of artificial relaxation for the individuals A (dotted), B (dotted-dashed) and C (dashed) of the individuals marked in Fig. 4. Continuous lines represent the frequency range included by the optimization, while dashed lines shows the not optimized frequency range. As a consequence of correcting dispersion, attenuation increases in the high frequency range, as shown in the top subplot. This extra attenuation is used as artificial attenuation for numerical-nonlinear stability.

A canonical case of a physical single relaxation process is presented. In order to correct numerical dispersion the parameters of two extra artificial relaxation processes have been found using the multi-objective genetic algorithms provided by the optimization toolbox in MATLAB R2014a v8.03. Linear propagation was considered and simulation parameters were set to typical values for water: m/s, kg/m, , Pas. A single physical relaxation process was included, with a characteristic relaxation time of and MHz, and relaxation modulus of that leads to a frozen sound speed of m/s. In this case, the numerical parameters were set to m and s. A plane wave traveling in direction was considered.

Thus, the theoretical attenuation for the relaxation processes and including viscosity can be expressed as


and the theoretical phase velocity can be predicted as Pierce (1989)


In order to compute the attenuation and dispersion of the numerical method, simulated pressure was recorded at two locations and , and attenuation and phase velocity were estimated from the spectral components over the bandwidth of the input signal. The numerical attenuation was calculated as


where is the Fourier transform of the measured pressure waveforms at points and . On the other hand, the phase velocity was computed as


where correct phase unwrapping is needed in the function.

In this way, Fig. 4 shows the retrieved Pareto front of the optimization, where 3 different areas can be distinguished. The first area, marked as (A) in Fig. 4, represents individuals whose numerical dispersion is minimal but attenuation is not optimal. On the other hand, the individuals around area (B) represent a set of relaxation parameters that provides the best agreement between numerical and physical attenuation. A good compromise between both situations can be obtained in the central area of the Pareto front (C), where retrieved attenuation and dispersion in the low frequency band shows good agreement with the physical, and the numerical dispersion has been corrected over a wide frequency range. However, as can be seen in the inset of Fig. 4, the dispersion relation retrieved by all the cases corrects the FDTD lossless numerical dispersion relation. The phase speed of those tree individuals is shown in Fig. 5 b), where it can be seen that in the frequency range selected for the optimization the numerical phase speed is corrected for all the individuals, where the best fit is obtained for individuals in the Pareto front area A. On the other hand, the inclusion of artificial relaxation leads to an increasing of the attenuation in the high frequency range, as is shown in Fig. 5 a). In this way, as the phase speed error is reduced the effect of artificial attenuation increases. Although this increasing can be seen as a non-desired counterpart, the appearance of this attenuation is useful in order to control the harmonic growing in nonlinear regime in the same way as artificial viscosity spatial operators Sparrow and Raspet (1991).

iv.2 Nonlinear steady solution for single relaxation process

Figure 6: (Color online) (a) Analytical (thick gray line), numerical using artificial relaxation with corrected dispersion (continuous line) and using artificial viscosity Sparrow and Raspet (1991) (dotted line) for the nonlinear steady state solution for . (b) Nonlinear steady state solution for . Inset shows detailed shock numerical solution for artificial relaxation ( markers) and artificial viscosity ( markers). (c) Nonlinear steady state solution for =0.01.

In order to validate the method in the nonlinear regime a full-wave simulation was developed in a mono-relaxing media using above parameters. Thus, the analytical (inverted) solution for the steady solution with for , for and for , for the retarded time reads Hamilton and Blackstock (1998)


where measures the ratio of relaxation effects to nonlinear effects. For , where no shock is present, the solution retrieved by FDTD algorithm shows good agreement with analytical and no artificial attenuation is needed. However, for a discontinuity is present in the solution and convergence is only possible with the inclusion of extra numerical techniques.

Thus, Fig. 6 (a-c) shows the analytical and numerical solutions including artificial relaxation and artificial viscosity, where excellent agreement is achieved in all cases. In the case of artificial viscosity, higher harmonics are strongly attenuated and by reducing grid step convergence can be achieved. Due to artificial viscosity operator is essentially a low pass spatial filter, a smoothed version of the shock is achieved. However, the phase speed of the higher spatial frequencies present in the shock is modified due to numerical dispersion, and for (Fig. 6 (b, c)) oscillations appears in the tail of the discontinuity, leading to the appearance of non-physical solutions.

On the other hand, the proposed method of artificial attenuation by relaxation also limits the harmonic growing so a smooth version of the shock appears. Moreover, artificial relaxation also corrects phase velocity so all the spatial frequencies travels at same speed and no oscillatory tail appears. The case of is shown in Fig. 6 (c), where nonlinear effects strongly dominates over attenuation. In this case, tail oscillations provided by artificial viscosity increases in amplitude. In contrast, by including artificial relaxation a smoothed version of the shock is captured and accuracy is maintained.

V Results

v.1 Frequency power law attenuation

Figure 7: Attenuation retrieved by the numerical algorithm (markers) and target frequency power law attenuation (gray lines). By using the optimization algorithm the relaxation times and modulus were optimized for minimize the relative error between the target power laws of and the attenuation retrieved.

Using methodology described above, the optimal relaxation parameters were obtained in order to fit the multiple-relaxation numerical attenuation to frequency power law attenuation in the form


where is the power law exponent and the power law coefficient in Np (rad/s)m. Moreover, the numerical dispersion was corrected by means of artificial relaxation in order to fit the corresponding frequency power law dispersion, where its analytical form satisfying causality can be expressed as Waters, Mobley, and Miller (2005)


This expression is valid for with , and an alternate equation can be found Waters, Mobley, and Miller (2005) in the limit for . Here, simulation parameters were m/s, kg/m, , Pas, MHz, m, s; that leads to 26 elements per wavelength and a CFL number of 0.9. Only two independent relaxation processes were employed in this section to obtain the target frequency power laws.

Figure 8: Phase speed retrieved by the numerical algorithm (markers) and target frequency power law attenuation (gray lines) for .

Following the above procedure, the relaxation times and relaxation modulus were optimized for different frequency power laws covering the range of that observed in tissues . The attenuation coefficient was chosen for each power law to present an attenuation dB/cm/MHz. The fitting was developed over the typical frequency range for medical ultrasound applications, i. e. 1 to 20 MHz for both attenuation and phase speed. The results for the attenuation and phase speed curves are plotted in Fig. 7 and Fig. 8, where the theoretical and the numerical predictions agree over the frequency range used for the fitting.

v.2 Fitting attenuation for tissue experimental data

Figure 9: Experimental attenuation data for some tissues adapted from Ref. Hill, Bamber, and ter Haar (2004) (lines), and obtained by the numerical method (markers) by fitting the parameters of 2 relaxation processes. The numbers above the curves show the exponent of the frequency power law for each frequency region (i.e. the slope of the curve).

Although a frequency power law dependence can describe the ultrasound attenuation over a finite frequency range, the attenuation data of some particular examples shows variation of the exponent over the entire frequency range Hill, Bamber, and ter Haar (2004). Thus, Fig. 9 shows experimental attenuation data curves for some tissues where the local slope of the power law changes over the measured frequency range. This behavior can be modeled by a sum of relaxation processes by optimizing the relaxation parameters as described above. Thus, the results show that most tissues with locally variable can be fitted by only a pair of relaxation processes, as the same way that for constant-slope frequency power law attenuationCleveland, Hamilton, and Blackstock (1996).

In this way, Table 1 shows the error of the numerical attenuation relative to the experimental data. The percent relative error was computed as , where is the experimental attenuation data, and define the frequency range of the measurement.

Tissue Power law
(local slope)
Skin 6.67 0.167 0.136 0.120
Liver 7.62 0.517 0.404 0.165
Blood , 8.34 0.349 0.330 0.310
Breast , 5.20 0.216 0.209 0.205
Skull bone , , 10.60 10.54 8.628 5.189
Table 1: Error of the optimized attenuation response relative to the experimental data for total relaxation processes.

As expected, the goodness of fit grows as the number of relaxation processes included increases. However, only two processes are enough to obtain relative errors below 1% for tissues with . In the case of tissues where a local value of has been observed, the fitting procedure fails, like in the skull bone in the 2 MHz range Hill, Bamber, and ter Haar (2004). The maximum slope achieved by single relaxation and thermo-viscous losses is for any frequency, so a tissue showing that slope cannot be accurately modeled in this frequency region with the method proposed in this work. From another point of view, Eq. (29) states that frequency power law medium with presents standard dispersion Waters, Mobley, and Miller (2005), opposite to anomalous dispersion for media falling in the range . Therefore, the dispersion relation of media with cannot be modeled by a sum of relaxation processes as long relaxation includes only anomalous dispersion.

Tissue Numerical Analytical
Skull bone 80.737 70.720
Skin 10.148 2.460
Breast 2.323 2.455
Liver 3.118 2.339
Blood 0.865 0.907
Table 2: Variation of sound speed observed numerically for the modeled tissues by means of two relaxation processes and analytical using the Kramers-Kronig relations.

Using Kramers-Kronig relations O’Donnell, Jaynes, and Miller (1981), the variations of sound speed can be predicted by the frequency dependent attenuation. Table 2 shows the variation of sound speed observed in the numerical solution over the fitted frequency range. The magnitude of these variations are of the order of magnitude of those measured experimentally in this frequency range, and the frequency dependence observed for the variation is roughly linear as observed in real tissueHill, Bamber, and ter Haar (2004). As expected from the relations between dispersion and absorption O’Donnell, Jaynes, and Miller (1981), the magnitude of the variation in sound speed increases as the total variation of the absorption increases for a given frequency range.

v.3 Nonlinear one-dimensional propagation in tissue-like media

v.3.1 Non-dispersive media

Figure 10: (Top) Waveforms at and for thermo-viscous attenuation (). Mendousse analytical solution (black line), -space (gray line) and FDTD numerical solution (markers). (Bottom) Spatial distribution of the first ten harmonics for Mendousse analytical solution (black line), -space (gray line) and FDTD numerical solution (markers).
Figure 11: (Top) Waveforms at and in a tissue-like media with frequency power law . -space (gray line) and FDTD numerical solution (markers). (Bottom) Spatial distribution of the first ten harmonics for -space (gray line) and FDTD numerical solution (markers).

In order to study the convergence of the numerical calculations to an analytical solution of the model in the nonlinear regime, a medium with frequency squared dependence attenuation is implemented using the adequate relaxation times and relaxation modulus as explained above. The numerical solution is compared with the analytical solution for a plane wave traveling through a thermo-viscous fluid proposed by Mendousse Pierce (1989):


where is the Gol’dberg number, defined for power law media as , with absorption length , shock formation distance and normalized distance ; with the parameter of nonlinearity and the acoustic Mach number , with is the source particle velocity and the wavenumber.

Figure 10 (top) presents the simulated waveforms at and . The wave steepening due to nonlinear processes in the absence of dispersion are well resolved by the numerical method presented here. In order to study the accuracy of the algorithm, the amplitude of the first ten harmonics has been extracted for numerical and analytic solutions and plotted versus in Fig. 10 (bottom). The observed relative error of the computational method decreases due to grid coarsening by a square law (i.e. the numerical scheme is second order accuracy). The magnitude of the error mainly depends on the number of elements per wavelength but, due to not ideal dispersion, also on the traveled propagated distance. Including the correction of dispersion by artificial relaxation, for a path length of 100 , a grid of 26 elements per wavelength was needed to obtain a relative error below 1 % for the third harmonic. Obviously, the relative error of the first and second harmonics will be always lower, i.e. the fundamental component error was 0.072 %.

In addition, the solution was compared also to the obtained by a -space method applied to the constitutive relations, i.e. the -wave algorithm Treeby et al. (2012). This method was selected due to the low numerical dispersion and the possibility of including frequency power law attenuation. The result of both computational methods and Eq. (30) agree over all the spectral components analyzed, showing convergence to the analytic solution.

v.3.2 Dispersive media

In the case of frequency power law attenuation media with no general analytic solution exist in nonlinear regime for monochromatic progressive waves. Thus, in order to study convergence in this regime, the proposed FDTD solution was compared with the solution obtained by -space methods Treeby et al. (2012). By using same physical and grid parameters in both methods, the solutions agrees for different power laws. Thus, Fig. 11 (top) shows the good agreement for the waveforms measured at and obtained for , where the characteristic asymmetry effect of media with anomalous dispersion (e.g. relaxing, boundary layer effects)Hamilton and Blackstock (1998) is observed: the shock front after the rarefaction phase is followed by a rounded positive compression profile. The spatial distribution for each harmonic is shown in Fig. 11 (bottom), where it is observed that the proposed FDTD solution with optimized attenuation and dispersion converges to the obtained by pseudo-spectral methods up to ten harmonics. As in the case of frequency squared media, grid refinement numerical tests have reported a second order accuracy of the FDTD method in nonlinear regime.

v.4 Nonlinear propagation in tissue-like media including diffraction

v.4.1 Experimental validation

An experiment was designed to test the validity of the algorithm for intense beams in frequency power law attenuation media. The source was formed by a plane single element piezoceramic crystal (PZ 26, Ferroperm Piezoceramics, Denmark) mounted in a custom designed steel housing and a polymethyl methacrylate (PMMA) focusing lens with aperture mm and radius of curvature mm. The source was driven with a sinusoidal pulse burst of frequency MHz and cycles using a function generator (14 bits, 100 MS/s, model PXI5412, National Instruments) and a linear RF amplifier (ENI 1040L, 400W, 55dB, ENI, Rochester, NY). The pressure waveforms were acquired with a HNR 500 m needle PVDF hydrophone (Onda Corp, CA), and a digitizer (64 MS/s, model PXI5620, National Instruments) was used. A three-axis micropositioning system (OWIS GmbH, Germany) was used to move the hydrophone in three orthogonal directions with an accuracy of 10 m. The amplitude frequency response of the hydrophone was compensated in post-processing but not in phase due to the absence of phase calibration for this equipment.

The source was completely immersed in a castor oil tank ( mm). We select this frequency power law attenuation media due to the low variability of its acoustic properties along existent literature Liebler et al. (2004); Treeby et al. (2009). Using a sound speed inside the bulk of the lens m/s, and a sound speed of the castor oil of m/s (at 26º C room temperature), the effective lens geometrical focal is estimated as mm, leading to a linear lossless gain of .

On the other hand, a nonlinear simulation including diffraction and frequency power law attenuation with same parameters was carried out in a workstation (20 cores Intel Xeon E5-2680 CPU, 2.8GHz with 256 GB RAM). The boundary conditions were implemented for a spherical focused ultrasound source. The castor oil parameters at 26º C room temperature Treeby et al. (2009), were m/s, kg/m, dB/cm/MHz, , . The grid parameters were m and ns, leading to a CFL number and elements per wavelength at fundamental frequency, note that this grid leads to for third harmonic.

Figure 12: Spatial distribution of the fundamental (), second () and third () obtained by the numerical (continuous line) and experimental methods (markers) for a focused transducer immersed in castor oil.

The balance between nonlinear effects and power law attenuation can be estimated by using the Gol’dberg ratio, where is the shock formation distance and the media attenuation characteristic length. Thus, the amplitude of the source were selected to obtain a Gol’dberg ratio of in order to let the frequency power law attenuation effects slightly dominate over nonlinear effects. On the other hand, the ratio between diffraction effects and nonlinear effects can be described by the so called Khokhlov number Naugolnykh and Ostrovsky (1998) as , where is the diffraction length and the source radius. For the proposed test a Khokhlov number of was selected to let the nonlinear effect slightly dominate over diffraction effects. The selected excitation pressure amplitude was kPa.

The results are summarized in Fig. 12, where axial pressure distribution for the fundamental, second and third harmonic are presented. A good agreement is found between simulations and the experimental test. Only far to the focal point the amplitude there exist differences between computations and experiments, that can be caused by nonuniform vibration of the source Canney et al. (2008), boundary effects of the PMMA lens, or miss-alignment of the source axis and micro-positioning system orthogonal directions along the 100 mm axial measurement.

The maximum amplitudes of the first harmonic were MPa for the experiment and for the numeric. The second harmonic peak pressure was kPa and kPa, and the third harmonic peak pressure and kPa for the experimental and numeric respectively. The relative errors between numerical and experimental results are % for the fundamental frequency, % for the second harmonic and % for the third. No error estimation was done for the peak pressure due to the absence of a phase calibration of the hydrophone.

v.4.2 Highly focused beam

Figure 13: Axial spatial distribution of the peak compression () and minimum rarefaction () pressure for a focused beam propagating through a liver tissue layer. (Top) weakly nonlinear propagation () and (bottom) strong nonlinear effects (). The insets show the waveforms recorded at the geometrical focal.
Figure 14: (Color online) Spatial distribution of the peak compression () and minimum rarefaction () pressure for a focused beam propagating through a liver tissue layer. (boundary marked with dashed line) for . Colorbars are in units.

In order to test the algorithm in the very high nonlinear regime with realistic tissue parameters a focused bowl of geometrical focal mm and aperture mm, driven at MHz was numerically tested. These parameters leads to a source gain and a -number, showing that source is beyond the paraxial limit. The media consist in two layers. The first layer, where the source was located, was water at 20º C with parameters m/s, kg/m, , dB/cm MHz, . At a distance mm, a layer of human liver tissue was placed, therefore the focal spot is located inside tissue at a depth of 15 mm. Liver tissue parameters Hill, Bamber, and ter Haar (2004) were , kg/m, , dB/cm MHz, . In this case, the numerical grid parameters were m and ns, leading to a CFL number and elements per wavelength at fundamental frequency ( for third harmonic).

Figure 13 shows the spatial peak pressure profiles for different excitation amplitudes. In Fig. 13 (top) the pressure of the source was MPa, while in the bottom figure was increased to MPa. Thus, for the selected parameters, Fig. 13 (top) present results for and , so the attenuation effects dominates over nonlinear effects and nonlinearity slightly dominates over diffraction effects. In this way, low asymmetry is observed between the positive compression peak, , and the minimum rarefaction pressure distribution, . The calculated waveform at , shown in the inset of Fig. 13 (top), is weakly distorted. However, there exist differences between its normalized peak amplitude and , and the source characteristic linear gain, . They are caused, in one hand by the attenuation effects, where the value of the lossy linear gain observed was measured at , i.e. the amplitude at the focal was reduced to 87.2% of the lossless amplitude. On the other hand, the differences due to the asymmetry between compression and rarefaction cycles are caused by the combined effect of nonlinearity and focusing.

If source amplitude is increased to MPa, as shown in Fig. 13 (bottom), the ratio between attenuation and nonlinear effects is increased to . In this regime, nonlinear effects are almost of the same order of attenuation effects. On the other hand, increasing source amplitude while keeping same transducer parameters implies also the Khokhlov number changes to , so the nonlinearity clearly dominates over diffraction effects. In this regime, highly asymmetric pressure distribution is observed, where the values at focal point are MPa and MPa.

Other typical nonlinear phenomena characteristic of high intensity focused sources can also be observed: formation of sharp shock front and its corresponding harmonic generation, or, as Fig. 14 shows, narrowing of the beam for and broadening for pressure distributions. In addition, nonlinear focal shift, i.e. displacement of the peak pressure relative to the position of the linear peak pressure can be also predicted for tissue propagation. In the case of it was observed a nonlinear focal shift mm and mm for the and pressure distribution respectively.

Vi Summary

A general model based on the full constitutive relations of nonlinear acoustics in relaxing media have been presented in a time-domain formulation which does not require convolutional operators. A numerical solution by means of finite-differences in time-domain have been obtained, showing that the theoretical attenuation and dispersion due to relaxation processes can be achieved by the numerical method with accuracy. These results can be also used to model typical relaxation process of other relaxing media (e. g. the processes observed in air, associated to the molecules of oxygen and nitrogen, or in seawater, associated to the relaxation of boric acid and magnesium sulfate).

Moreover, a method for modeling frequency power law attenuation by means of multiple relaxation has been implemented in the constitutive relations. The proposed method can describe local variations of the exponent of the frequency power law, so an arbitrary attenuation curve in the range can be modeled by means of the proper optimization of the relaxation coefficients. This feature of the presented method is an advantage when compared with most fractional derivatives methods, where the attenuation follows an exact but unique frequency power law over the entire frequency range. A broad range of human tissues have been modeled and the goodness of the fit using from two to four relaxation processes has been discussed.

Furthermore, a computational technique that exploits the anomalous dispersion of the relaxation processes is employed to mitigate the numerical dispersion of the finite-differences scheme. Thus, while phase speed is corrected by including artificial relaxation processes, its corresponding artificial attenuation is used to improve stability in the nonlinear regime. In this way, smooth and stable versions of shock waves have been obtained and compared with its analytic solution. Furthermore, the validity of the algorithm including diffraction have been tested with experimental measurements of a focused beam in castor oil.

Due to the model is developed from the constitutive relations for nonlinear acoustics, most wave phenomena is captured. As a difference from the one-way models the proposed model implicitly includes multiple wave direction, and, due to the Lagrangian density of acoustic energy is implicitly included in the computation, multiple scattering and strong resonance effects can be accurately described. Moreover, unlike KZK and other parabolic approximations, the proposed model captures the diffraction exactly, so for simulation of acoustic beams the field is not approximated only to the beam axis, but also in the near field and far to the beam axis and thus high focused devices can be simulated.

The code has shown to be particularly appropriate if the problem to simulate presents axisymmetry, because the constitutive relations for nonlinear acoustics are solved in a computational 2D domain, while standard -space methods need to employ full 3D domains due to the poor convergence of the Fourier series at discontinuities (). This is the case, for example, of the focused ultrasound transducer simulated in Sec. V.4.2, where a full 3D solution will require huge computational resources and calculation times. Finally, due to the particle velocity vector and the acoustic density fields are solved implicitly by the code, this information can be used to estimate other relevant magnitudes as the full nonlinear intensity vector, the nonlinear acoustic radiation forces in these absorbing media or the acoustic streaming generated in frequency power law attenuation fluids.


The authors acknowledge financial support from the FPI program of the Universitat Politècnica de València.


  • Aanonsen et al. (1984) S. I. Aanonsen, T. Barkve, J. N. Tjotta,  and S. Tjotta, “Distortion and harmonic generation in the nearfield of a finite amplitude sound beam,” J. Acoust. Soc. Am. 75, 749–768 (1984).
  • Lee and Hamilton (1995) Y.-S. Lee and M. F. Hamilton, “Time-domain modeling of pulsed finite-amplitude sound beams,” J. Acoust. Soc. Am. 97, 906–917 (1995).
  • Cleveland, Hamilton, and Blackstock (1996) R. O. Cleveland, M. F. Hamilton,  and D. T. Blackstock, “Time-domain modeling of finite-amplitude sound in relaxing fluids,” J. Acoust. Soc. Am. 99, 3312–3318 (1996).
  • Yang and Cleveland (2005) X. Yang and R. O. Cleveland, “Time domain simulation of nonlinear acoustic beams generated by rectangular pistons with application to harmonic imaging,” J. Acoust. Soc. Am. 117, 113–123 (2005).
  • Khokhlova et al. (2006) V. Khokhlova, A. Ponomarev, M. Averkiou,  and L. Crum, “Nonlinear pulsed ultrasound beams radiated by rectangular focused diagnostic transducers,” Acoustical Physics 52, 481–489 (2006).
  • Jing and Cleveland (2007) Y. Jing and R. O. Cleveland, “Modeling the propagation of nonlinear three-dimensional acoustic beams in inhomogeneous media,” J. Acoust. Soc. Am. 122, 1352–1364 (2007).
  • Soneson and Myers (2007) J. E. Soneson and M. R. Myers, “Gaussian representation of high-intensity focused ultrasound beams,” The Journal of the Acoustical Society of America 122, 2526–2531 (2007).
  • Prieur and Holm (2011) F. Prieur and S. Holm, “Nonlinear acoustic wave equations with fractional loss operators,” J. Acoust. Soc. Am. 130, 1125–1132 (2011).
  • Christopher and Parker (1991) P. T. Christopher and K. J. Parker, “New approaches to nonlinear diffractive field propagation,” J. Acoust. Soc. Am. 90, 488–499 (1991).
  • Tavakkoli et al. (1998) J. Tavakkoli, D. Cathignol, R. Souchon,  and O. A. Sapozhnikov, “Modeling of pulsed finite-amplitude focused sound beams in time domain,” The Journal of the Acoustical Society of America 104, 2061–2072 (1998).
  • Zemp, Tavakkoli, and Cobbold (2003) R. J. Zemp, J. Tavakkoli,  and R. S. Cobbold, “Modeling of nonlinear ultrasound propagation in tissue from array transducers,” J. Acoust. Soc. Am. 113, 139–152 (2003).
  • Varslot and Taraldsen (2005) T. Varslot and G. Taraldsen, “Computer simulation of forward wave propagation in soft tissue,” Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on 52, 1473–1482 (2005).
  • Yuldashev and Khokhlova (2011) P. Yuldashev and V. Khokhlova, “Simulation of three-dimensional nonlinear fields of ultrasound therapeutic arrays,” Acoustical physics 57, 334–343 (2011).
  • Hallaj and Cleveland (1999) I. M. Hallaj and R. O. Cleveland, ‘‘Fdtd simulation of finite-amplitude pressure and temperature fields for biomedical ultrasound,” J. Acoust. Soc. Am. 105, L7–L12 (1999).
  • Hallaj, Cleveland, and Hynynen (2001) I. M. Hallaj, R. O. Cleveland,  and K. Hynynen, “Simulations of the thermo-acoustic lens effect during focused ultrasound surgery,” The Journal of the Acoustical Society of America 109, 2245–2253 (2001).
  • Pinton et al. (2009) G. F. Pinton, J. Dahl, S. Rosenzweig,  and G. E. Trahey, “A heterogeneous nonlinear attenuating full-wave model of ultrasound,” Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on 56, 474–488 (2009).
  • Huijssen and Verweij (2010) J. Huijssen and M. D. Verweij, “An iterative method for the computation of nonlinear, wide-angle, pulsed acoustic fields of medical diagnostic transducers,” J. Acoust. Soc. Am. 127, 33–44 (2010).
  • Demi, Van Dongen, and Verweij (2011) L. Demi, K. Van Dongen,  and M. Verweij, “A contrast source method for nonlinear acoustic wave fields in media with spatially inhomogeneous attenuation,” J. Acoust. Soc. Am. 129, 1221–1230 (2011).
  • Verweij, Demi, and van Dongen (2013) M. D. Verweij, L. Demi,  and K. W. van Dongen, “Computation of nonlinear ultrasound fields using a linearized contrast source method,” J. Acoust. Soc. Am. 134, 1442–1453 (2013).
  • Jing, Wang, and Clement (2012) Y. Jing, T. Wang,  and G. T. Clement, “A k-space method for moderately nonlinear wave propagation,” Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on 59, 1664–1673 (2012).
  • Jing, Shen, and Clement (2011) Y. Jing, D. Shen,  and G. T. Clement, “Verification of the westervelt equation for focused transducers,” Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on 58, 1097–1101 (2011).
  • Hamilton and Blackstock (1998) M. Hamilton and D. Blackstock, Nonlinear Acoustics (Academic Press, San Diego, CA 92101, USA, 1998) p. 455.
  • Manry and Broschat (1996) C. W. Manry and S. L. Broschat, “Fdtd simulations for ultrasound propagation in a 2-d breast model,” Ultrasonic imaging 18, 25–34 (1996).
  • Mast et al. (1997) T. D. Mast, L. M. Hinkelman, M. J. Orr, V. W. Sparrow,  and R. C. Waag, “Simulation of ultrasonic pulse propagation through the abdominal wall,” J. Acoust. Soc. Am. 102, 1177–1190 (1997).
  • Yuan et al. (1999) X. Yuan, D. Borup, J. Wiskin, M. Berggren,  and S. A. Johnson, ‘‘Simulation of acoustic wave propagation in dispersive media with relaxation losses by using fdtd method with pml absorbing boundary condition,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 46, 14–23 (1999).
  • Mast et al. (2001) T. D. Mast, L. P. Souriau, D.-L. Liu, M. Tabei, A. I. Nachman,  and R. C. Waag, “A k-space method for large-scale models of wave propagation in tissue,” Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on 48, 341–354 (2001).
  • Treeby and Cox (2010) B. E. Treeby and B. Cox, “Modeling power law absorption and dispersion for acoustic propagation using the fractional laplacian,” J. Acoust. Soc. Am. 127, 2741–2748 (2010).
  • Ginter et al. (2002) S. Ginter, M. Liebler, E. Steiger, T. Dreyer,  and R. E. Riedlinger, “Full-wave modeling of therapeutic ultrasound: nonlinear ultrasound propagation in ideal fluids.” J. Acoust. Soc. Am. 111, 2049–2059 (2002).
  • Sparrow and Raspet (1991) V. W. Sparrow and R. Raspet, ‘‘A numerical method for general finite amplitude wave propagation in two dimensions and its application to spark pulses,” J. Acoust. Soc. Am. 90, 2683–2691 (1991).
  • Liebler et al. (2004) M. Liebler, S. Ginter, T. Dreyer,  and R. E. Riedlinger, “Full wave modeling of therapeutic ultrasound: efficient time-domain implementation of the frequency power-law attenuation.” J. Acoust. Soc. Am. 116, 2742–2750 (2004).
  • Lobanova, Lobanov, and Khokhlova (2014) E. G. Lobanova, S. V. Lobanov,  and V. A. Khokhlova, ‘‘Counterpropagation of waves with shock fronts in a nonlinear tissue-like medium,” Acoustical Physics 60, 387–397 (2014).
  • Treeby et al. (2012) B. E. Treeby, J. Jaros, A. P. Rendell,  and B. Cox, “Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method,” J. Acoust. Soc. Am. 131, 4324–4336 (2012).
  • Albin et al. (2012) N. Albin, O. P. Bruno, T. Y. Cheung,  and R. O. Cleveland, “Fourier continuation methods for high-fidelity simulation of nonlinear acoustic beams,” J. Acoust. Soc. Am. 132, 2371–2387 (2012).
  • Shahbazi et al. (2011) K. Shahbazi, N. Albin, O. P. Bruno,  and J. S. Hesthaven, “Multi-domain fourier-continuation/weno hybrid solver for conservation laws,” J. Comput. Phys. 230, 8779 – 8796 (2011).
  • Naugolnykh and Ostrovsky (1998) K. Naugolnykh and L. Ostrovsky, Nonlinear Wave Processes in Acoustics, Cambridge Texts in Applied Mathematics (Cambridge University Press, 40 West 20th Street, New York, NY 10011-4211, USA, 1998) p. 298.
  • Rudenko and Soluian (1977) O. Rudenko and S. Soluian, Theoretical foundations of nonlinear acoustics, Studies in Soviet science (Consultants Bureau, New York, NY 10011, USA, 1977) p. 274.
  • Pierce (1989) A. Pierce, Acoustics: An Introduction to Its Physical Principles and Applications (Acoustical Society of America, Melville, New York, NY 11747-4502, USA, 1989) p. 678.
  • Botteldooren (1996) D. Botteldooren, ‘‘Numerical model for moderately nonlinear sound propagation in three-dimensional structures,” J. Acoust. Soc. Am. 100, 1357–1367 (1996).
  • LeVeque (1992) R. LeVeque, Numerical Methods for Conservation Laws, edited by L. in Mathematics. ETH Zurich (Birkhauser Verlag, 133 CH-4010 Basel, Switzerland, 1992) p. 214.
  • Liu and Tao (1997) Q.-H. Liu and J. Tao, “The perfectly matched layer for acoustic waves in absorptive media,” J. Acoust. Soc. Am. 102, 2072–2082 (1997).
  • Liu (1999) Q. H. Liu, “Perfectly matched layers for elastic waves in cylindrical and spherical coordinates,” J. Acoust. Soc. Am. 105, 2075–2084 (1999).
  • Waters, Mobley, and Miller (2005) K. R. Waters, J. Mobley,  and J. G. Miller, “Causality-imposed (kramers-kronig) relationships between attenuation and dispersion,” Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on 52, 822–823 (2005).
  • Hill, Bamber, and ter Haar (2004) C. Hill, J. C. Bamber,  and G. ter Haar, Physical Principles of Medical Ultrasonics, 2nd ed. (John Wiley & Sons Ltd, Chichester, West Sussex PO19 8SQ, England, 2004) p. 528.
  • O’Donnell, Jaynes, and Miller (1981) M. O’Donnell, E. T. Jaynes,  and J. G. Miller, “Kramers–kronig relationship between ultrasonic attenuation and phase velocity,” J. Acoust. Soc. Am. 69, 696–701 (1981).
  • Treeby et al. (2009) B. E. Treeby, B. T. Cox, E. Z. Zhang, S. K. Patch,  and P. C. Beard, “Measurement of broadband temperature-dependent ultrasonic attenuation and dispersion using photoacoustics,” Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on 56, 1666–1676 (2009).
  • Canney et al. (2008) M. S. Canney, M. R. Bailey, L. A. Crum, V. A. Khokhlova,  and O. A. Sapozhnikov, “Acoustic characterization of high intensity focused ultrasound fields: A combined measurement and modeling approach,” The Journal of the Acoustical Society of America 124, 2406–2420 (2008).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description