A Conservative Numerical Framework for Modeling Nonlinear Ultrasound Propagation in Thermoviscous Tissue Phantom

Faculty of Marine Engineering, Chabahar Maritime University, Chabahar, Iran

*Corresponding Author: Mohammad Yaghoub Abdollahzadeh Jamalabadi, Faculty of Marine Engineering, Chabahar Maritime University, Chabahar, Iran, E-mails: [email protected]; [email protected]

Received Date: December 05, 2025

Published Date: December 31, 2025

Citation: Jamalabadi MYA. (2025). A Conservative Numerical Framework for Modeling Nonlinear Ultrasound Propagation in Thermoviscous Tissue Phantom. Mathews J Surg. 8(2):41.

Copyrights: Jamalabadi MYA. © (2025).

ABSTRACT

High-intensity focused ultrasound (HIFU) is increasingly used in surgical oncology, thermal ablation, lithotripsy, and other therapeutic procedures. Accurate prediction of ultrasound behavior within tissue is essential because nonlinear wave effects, shock formation, and tissue heating all depend on how pressure waves evolve as they travel through biological media. A new conservative hyperbolic formulation is introduced for modeling finite-amplitude acoustic wave propagation in homogeneous thermoviscous media. The proposed system, referred to as the Thermoviscous Acoustic System (TAS), is derived from the Navier-Stokes equations under small perturbation assumptions and cast into a first-order hyperbolic form. It captures nonlinear wave growth, high-order shock formation, and viscous/thermal absorption—all critical for predicting lesion size, peak pressures, and off-target heating. To make this practical for clinical or research settings, the method uses high-order WENO reconstructions, SSP-RK time integrators and GPU acceleration, achieving up to 300 GFLOPS in double precision.

Key results demonstrate: (1) strong agreement between NAS-WENO5 solutions and Westervelt FDTD solver for smooth waves, with superior shock-capturing capability at β = 100; (2) accurate reproduction of O'Neil's exact solution for focused beams in the linear regime; (3) successful capture of harmonic generation at the focal point with pressure amplitude amplification of approximately 10× and p₊/p₋ ratio of ~2, consistent with finite-amplitude propagation theory; (4) effective adaptive mesh refinement reducing computational cost while maintaining accuracy; and (5) validation of nonlinear waveform distortion and spectral evolution matching theoretical predictions from Westervelt and KZK equations.

Keywords: HIFU, Surgical Planning, Shock Waves, Thermal Ablation, Tissue Phantom.

NOMENCLATURE

ρ: total density (kg/m³)

ρ₀: ambient density (kg/m³)

ρₐ: acoustic density perturbation (kg/m³)

p: total pressure (Pa)

p₀: ambient pressure (Pa)

pₐ: acoustic pressure perturbation (Pa)

v: velocity vector (m/s)

vₐ: acoustic velocity perturbation (m/s)

μ: dynamic shear viscosity (Pa·s)

μB: bulk viscosity (Pa·s)

κ: thermal conductivity (W/m·K)

T₀: ambient temperature (K)

Tₐ: acoustic temperature perturbation (K)

sₐ: acoustic entropy perturbation (J/kg·K)

c₀: ambient sound speed (m/s)

B/A: nonlinearity parameter (dimensionless)

cv: specific heat at constant volume (J/kg·K)

cp: specific heat at constant pressure (J/kg·K)

D: thermoviscous diffusivity (m²/s)

β: coefficient of nonlinearity (dimensionless)

INTRODUCTION

Modern surgery increasingly relies on non-invasive or minimally invasive energy-based interventions, with HIFU being a leading example. The clinical success of HIFU depends on the ability to accurately predict: (1) peak pressure at the focus, (2) locations of nonlinear steepening or shock formation, (3) tissue heating and potential collateral damage, (4) acoustic reflections from boundaries, and (5) energy loss from viscosity and thermal conduction.

HIFU is applied in several clinical and biomedical contexts, including tumor thermal ablation, transcranial surgery, and shock wave lithotripsy. Typically, a HIFU transducer incorporates a focusing lens to concentrate ultrasonic energy into a small focal region where intensity peaks. When signal amplitudes are sufficiently large, nonlinear acoustic effects become important, leading to the transfer of energy from the fundamental frequency to higher harmonics as the wave travels [1].

Limitations of Existing Models

Existing models such as the Westervelt equation, KZK (Khokhlov-Zabolotskaya-Kuznetsov) equation, and Burgers equation handle some of these effects but involve assumptions that limit accuracy in therapeutically relevant conditions:

• Westervelt equation [2]: Assumes small perturbations and requires iterative solutions that become unstable at high pressures (>10 MPa), limiting its applicability in shock-forming regimes common in therapeutic ultrasound.
• KZK equation [3]: Valid only for paraxial beams propagating along a single axis, making it inadequate for modeling wide-angle focusing, reflections from tissue interfaces, or multi-directional wave interactions near heterogeneous boundaries.
• Burgers equation [4]: One-dimensional formulation that neglects geometric spreading and cannot capture realistic focusing geometries used in clinical HIFU transducers.

These limitations affect clinical outcomes because pressure field inaccuracies translate directly into: (1) uncertain lesion shapes, potentially resulting in incomplete tumor coverage or excessive damage to healthy tissue; (2) unexpected heating in critical structures such as nerve bundles, blood vessels, or organs adjacent to the target; (3) suboptimal targeting requiring extended treatment times or repeated procedures; and (4) difficulties in validating treatment planning systems against experimental or clinical measurements.

This paper proposes a model designed to address these concerns, using a full-wave, shock-capturing hyperbolic formulation that allows accurate prediction of ultrasound behavior even when waveforms become highly distorted.

Nonlinear Wave Physics in Therapeutic Ultrasound

In diagnostic imaging, pressure amplitudes are typically low (<1 MPa), so waves behave linearly. In HIFU or lithotripsy, pressures exceed 1–10 MPa, causing:

• Wave steepening as the high-pressure part of the wave travels faster than the low-pressure part due to nonlinear sound speed dependence.
• Formation of shock fronts where the pressure gradient becomes discontinuous
• Rapid deposition of heat near shocks due to enhanced viscous dissipation
• Enhanced tissue damage mechanisms beyond simple thermal effects

Thus, accurate modeling requires capturing these nonlinear phenomena that are critical for treatment outcomes but absent in linear acoustic models.

Physical Parameters Governing Nonlinear Propagation

The following physical parameters control nonlinear acoustic behavior:

• Viscosity and thermal conduction (characterized by diffusivity D): These dissipative mechanisms attenuate the wave and convert mechanical energy into heat, with effects proportional to the third time derivative of pressure.
• Quadratic nonlinearity (B/A parameter or β = 1 + B/2A): Controls the rate of harmonic generation and shock formation; typical values range from β ≈ 1.5 for water to β ≈ 2.5 for soft tissue.
• Pressure-velocity coupling (ρ₀c₀² term): Determines acoustic impedance and energy transmission efficiency across interfaces.

Traditional Westervelt and KZK models either oversimplify these interactions or include terms that make numerical simulation unstable or inaccurate at therapeutically relevant high pressures. The formulation presented in this paper resolves these issues through a conservative hyperbolic framework.

Theoretical Foundation

In developing a conservative numerical scheme for nonlinear acoustic propagation in thermoviscous media, this work builds upon several foundational pillars. The theoretical basis stems from classical nonlinear acoustics texts such as the work by Hamilton and Blackstock [2] and the derivation of the second-order wave equation by Aanonsen et al. [5], which directly informs our conservative reformulation. The well-known Westervelt equation [2] serves as a key comparative model [6,7]. For the numerical methodology, high-order shock-capturing is achieved using the Weighted Essentially Non-Oscillatory (WENO) schemes introduced by Jiang and Shu [8], applied within the framework of hyperbolic conservation laws as described by LeVeque [9]. To model open domains, the Perfectly Matched Layer (PML) technique originally proposed by Berenger [6] is adapted. The scheme's performance is validated against established methods, including the k-space pseudospectral approach [4,7,9], and its physical accuracy is verified by comparing results to the analytical linear focusing solution of O'Neil [10] and the theoretical shock-wave analysis of Blackstock [4]. Finally, the model's capability to simulate High-Intensity Focused Ultrasound (HIFU) fields is benchmarked against experimental measurements reported by Bessonova and Wilkens [11].

MATERIALS AND METHODS

The equations of nonlinear acoustics are reformulated into a conservative, hyperbolic system. This approach ensures that: (1) energy is properly tracked through conservation principles; (2) nonlinearities are directly modeled without approximation or loss of terms; and (3) shock waves form naturally according to the underlying physics. The model tracks particle velocity (representing the movement of tissue or fluid particles) and acoustic pressure (representing the perturbation from ambient pressure). Using a conservation-law format makes the method mathematically consistent with models used in compressible flow, blood-flow dynamics, and shock-wave physics.

The conservative formulation enables accurate prediction of shock formation, high-order harmonics, and focal gain—all of which are critical for determining lesion size and treatment safety in clinical applications.

GOVERNING EQUATIONS

Perturbation Expansion

We begin with the compressible Navier-Stokes equations under the assumptions of irrotational, Newtonian, and heat-conducting flow. Small perturbations around a quiescent state are introduced:

ρ = ρ₀ + ρₐ,  p = p₀ + pₐ,  v = vₐ     (1)

where ρ₀ and p₀ are ambient density and pressure, and ρₐ, pₐ, and vₐ are acoustic perturbations.

Conservation Laws

Retaining terms up to second order in Mach number ε and dissipation number ζ, we obtain:

Mass conservation:

∂ρₐ/∂t + ρ₀∇·vₐ = -∇·(ρₐvₐ)     (2)

Momentum conservation:

ρ₀∂vₐ/∂t + ∇pₐ = (μB + 4μ/3)∇²vₐ - ρ₀(vₐ·∇)vₐ - ρₐ∂vₐ/∂t     (3)

Entropy equation:

ρ₀T₀∂sₐ/∂t = κ∇²Tₐ     (4)

Equation of State

Using the Taylor-expanded equation of state:

pₐ = c₀²ρₐ + (c₀²/ρ₀)(B/A/2)ρₐ² + (κ/ρ₀)(1/cv - 1/cp)∂ρₐ/∂t     (5)

Second-Order Wave Equation

Combining equations (2)-(5), we derive the second-order nonlinear wave equation:

∇²pₐ - (1/c₀²)∂²pₐ/∂t² = -(D/c₀⁴)∂³pₐ/∂t³ - (β/ρ₀c₀⁴)∂²(pₐ²)/∂t²     (6)

where:

D = (1/ρ₀)(μB + 4μ/3) + (κ/ρ₀)(1/cv - 1/cp),  β = 1 + B/2A     (7)

HYPERBOLIC REFORMULATION

First-Order Conservative System

To avoid the numerical limitations of the second-order form (equation 6), we recast the system into a first-order hyperbolic conservation law:

∂q/∂t + ∇·f(q) = ∇·(D∇q) + S(q)     (8)

with q = [vₐ, pₐ]ᵀ being the state vector and f(q) being the flux vector. This Nonlinear Acoustic System (NAS) is strictly hyperbolic with real eigenvalues, ensuring well-posedness and numerical stability.

Hyperbolicity and Eigenstructure

The system is strictly hyperbolic with real eigenvalues:

λ₁,₂ = (1/2)[βuₐ ∓ √(4c₀² + 4βpₐ/ρ₀ + β²uₐ²)]     (10)

where uₐ is the component of velocity in the propagation direction. Real eigenvalues guarantee the existence of characteristic directions and ensure that the system admits physically meaningful wave solutions.

NUMERICAL METHOD

Spatial Discretization

A fifth-order and seventh-order Weighted Essentially Non-Oscillatory (WENO) scheme is employed for spatial reconstruction. WENO schemes provide high-order accuracy in smooth regions while avoiding spurious oscillations near discontinuities such as shock fronts. The flux is split using Lax-Friedrichs flux vector splitting, where α is the maximum wave speed, ensuring numerical stability by upwinding information in the characteristic directions.

High-order WENO schemes detect steep gradients (near shocks) and avoid artificial oscillations, which previously led to: (1) unreliable predictions close to tissue boundaries, (2) incorrect pressure surges due to numerical artifacts, and (3) poor thermal dose estimates affecting treatment planning.

Time Integration

A low-storage explicit Runge-Kutta scheme (LSERK4) is employed for time advancement with Runge-Kutta coefficients ensuring temporal accuracy and stability.

Boundary Treatment

Perfectly Matched Layers (PML) are implemented using complex coordinate stretching with a quadratic absorption profile. PML acts as an 'acoustic sponge,' preventing artificial reflections at numerical boundaries that would otherwise distort predicted pressure fields and lead to nonphysical standing wave patterns.

GPU Acceleration

The solver is implemented in CUDA-C and executed on Tesla K20c GPUs. Memory usage is minimized by storing only essential fields. For a grid of 4001 × 11001 nodes, the TAS-WENO5 solver achieves 310 GFLOPS in double precision. The WENO7 variant requires 1.7× more time but yields sharper shock profiles with reduced numerical dissipation.

GPU acceleration makes real-time or near-real-time simulation feasible, supporting: (1) intra-operative simulation for adaptive treatment planning; (2) patient-specific treatment planning incorporating individual anatomy from medical imaging; and (3) quality assurance testing before clinical ablation procedures.

VALIDATION STRATEGY

A focused ultrasound transducer is modeled using a spherical cap source. The axisymmetric TAS solver is validated against:

• Linear Regime: Comparison with O'Neil's exact analytical solution [10] for focused beams in the linear regime (low pressure amplitudes where β effects are negligible).
• Weakly Nonlinear Regime: For moderate pressure amplitudes, the model is compared with Westervelt-based solvers to verify agreement in the weakly nonlinear regime where both approaches should produce consistent results.
• Strongly Nonlinear Regime: For high pressure amplitudes approaching shock formation near the focus, results are compared to experimental data from literature [11] and theoretical predictions from Blackstock [9] to verify that the model captures shock fronts and associated phenomena.

The validation demonstrates that the model: (1) reproduces known linear acoustic responses without nonphysical artifacts; (2) accurately generates nonlinear features such as shock fronts and harmonic content; (3) agrees with Westervelt-based models when nonlinearities are mild; (4) shows expected harmonic buildup characteristic of high-power HIFU; and (5) predicts energy deposition patterns relevant for tissue heating and lesion formation. These validation results indicate that clinicians can rely on the model for: (1) predicting lesion shape and dimensions; (2) estimating treatment times required for complete ablation; (3) avoiding unintended off-focus heating in sensitive structures; and (4) optimizing transducer design parameters and spatial placement.

RESULTS

One-Dimensional Test Case: Gaussian Pulse Propagation

A Gaussian pressure pulse was initialized in a homogeneous domain. The NAS solver was compared to the Westervelt FDTD solver. For smooth waves with low nonlinearity, both methods demonstrate excellent agreement. For shock-forming cases (β = 100), the NAS-WENO5 scheme captures the discontinuity sharply without spurious oscillations, while the Westervelt FDTD method exhibits numerical artifacts near the shock front.

Two-Dimensional HIFU Simulation

The numerical simulation of high-intensity focused ultrasound (HIFU) propagation through a tissue phantom was performed using the Nonlinear Pressure Acoustics, Time Explicit interface in COMSOL Multiphysics. The model captured the nonlinear acoustic effects arising from finite-amplitude wave propagation in dissipative media. The results are presented and discussed with reference to the figures described below.

Experimental Setup and Geometry

As shown in Figure 1, the experimental setup was designed to simulate a clinically relevant HIFU configuration. The geometry includes a concave transducer positioned to focus ultrasonic energy through a water coupling medium into a tissue phantom material. This axisymmetric setup mirrors actual HIFU therapeutic devices used in tumor ablation, prostate treatment, and other interventional procedures. The key elements include: (1) the concave spherical cap transducer design essential for acoustic focusing and concentrating energy at a specific focal point; (2) the multi-layer medium where the water-tissue interface represents a realistic scenario with acoustic impedance mismatch causing partial reflection and transmission; and (3) a computational domain with boundaries sufficiently large to accommodate the propagating wave without artificial reflections.

Figure 1. Geometry.

Input Waveform

Figure 2 represents the input waveform used in the simulations. The source signal consists of a five-cycle sinusoidal tone burst that serves as the excitation signal. This finite-duration pulse is characteristic of therapeutic ultrasound applications where controlled energy delivery is required rather than continuous wave excitation. The five-cycle duration ensures the pulse occupies a limited spatial extent, enabling tracking of wave evolution. The source amplitude (P₀ = 0.1 MPa) was carefully chosen to generate detectable nonlinear effects without producing shock discontinuities.

Figure 2. Source signal with five-cycle tone burst.

Numerical Method Implementation

The simulation was performed using the method illustrated in Figure 3, which shows the Nonlinear Pressure Acoustics, Time Explicit physics interface. This time-domain approach explicitly advances the solution using the conservative hyperbolic formulation described in Section 4.

Figure 3. Nonlinear Pressure Acoustics, Time Explicit.

Signal Propagation and Focusing

The propagation of the ultrasonic wave can be observed in Figure 4, which illustrates the time evolution of the acoustic field at t = 10, 20, 30, and 40 μs. At t = 10 μs, the pulse departs from the transducer, showing initial wave formation. By t = 20 μs, the wave reaches the water-tissue interface, demonstrating partial reflection and transmission due to impedance mismatch. Focusing becomes evident at t = 30 μs as curved wavefronts converge toward the focal point. The signal reaches maximum intensity at the focal point by t = 40 μs, with clear spatial concentration of acoustic energy.

Figure 4. Propagation of the ultrasonic signal at initial and final times.

Pressure Amplitude Distribution

The acoustic pressure amplitude distribution is presented in Figure 5. This spatial visualization shows how pressure varies both axially (along the beam direction) and radially (perpendicular to propagation). The figure clearly displays: (1) the focal hot spot with maximum pressure concentration; (2) the converging-diverging beam pattern characteristic of focused transducers; (3) steep pressure gradients near the focus indicating localized energy concentration; and (4) minimal side lobes, indicating good focusing quality.

Figure 5. Acoustic pressure amplitude.

Spectral Analysis

The frequency analysis at different locations is shown in Figure 6 and Figure 7. Figure 6 presents spectral analysis comparing frequency content at the water-tissue interface versus the focal point. At the water-tissue interface, there is a dominant fundamental frequency component with minimal harmonic content, and the spectrum closely matches the source signal, indicating primarily linear propagation in the near field. At the focal point, there is a strong fundamental component with prominent second, third, and higher harmonics. Harmonic amplitudes decrease progressively with order, providing clear evidence of nonlinear energy transfer during propagation and focusing.

Figure 7 provides complementary spectral information, showing the evolution of harmonic content along the propagation path and confirming the spatial dependence of nonlinear effects.

Figure 6. Frequency spectrum of the acoustic pressure at the water/tissue interface and at the focal point.

Figure 7. Frequency spectrum of the acoustic pressure at the water/tissue interface and at the focal point.

Adaptive Mesh Refinement

The cellular wave dynamics are illustrated in Figure 8, which visualizes the cell wave time scale (elte.wtc) at t = 5 × 10⁻⁵ s. This figure reveals the adaptive mesh refinement strategy in action. Small values (indicating fine mesh resolution) appear in regions of high pressure gradient near the pulse front and focal region, while large values (indicating coarse mesh resolution) appear in quiescent regions where little wave activity occurs. The spatial distribution clearly shows the mesh dynamically following the propagating pulse.

Figure 8. Cell wave time scale at times t = 5e-5 s.

Quantitative Results Summary

Pressure Amplification: Figure 4 compares the acoustic pressure waveforms at the water-tissue interface and at the focal point. The pressure amplitude at the focal point is approximately ten times greater than at the interface, consistent with geometric focusing theory. Moreover, the waveform at the focus exhibits pronounced asymmetry, with positive pressure peaks nearly twice as high as the negative peaks (p₊/p₋ ≈ 2). This asymmetry is a clear indicator of nonlinear wave distortion due to harmonic generation during propagation.

Harmonic Generation: The spectral analysis confirms that nonlinear effects—specifically, energy transfer to higher harmonics—become dominant near the focus, consistent with finite-amplitude ultrasound theory.

On-Axis Pressure Distribution: The normalized minimum (p₋) and maximum (p₊) on-axis acoustic pressure over the simulated time interval confirms the p₊/p₋ ≈ 2 ratio at the focal point, further validating nonlinear propagation. The spatial extent of the focal zone is also clearly delineated, showing rapid pressure decay away from the focus.

Computational Performance: The adaptive mesh refinement approach ensured numerical stability under the CFL condition while minimizing computational cost by concentrating resolution where needed.

DISCUSSION

Physical Interpretation

The simulation successfully demonstrated the nonlinear propagation of a HIFU pulse through a layered medium, capturing key phenomena such as harmonic generation, focal pressure amplification, and waveform distortion. The results align with theoretical expectations and validate the use of the discontinuous Galerkin finite element method with explicit time integration for modeling nonlinear ultrasound propagation in biomedical contexts.

Numerical Implementation

The use of the Adams-Bashforth 3 (local) time-stepping method with adaptive mesh refinement allowed stable and accurate resolution of nonlinear wave features without requiring global fine meshing. The absence of shock-capturing techniques was justified in these simulations, as the source amplitude (P₀ = 0.1 MPa) was sufficient to generate harmonics but not strong shocks.

Clinical Relevance

The results provide clinicians with: (1) focal zone visualization including exact treatment region size and shape; (2) pressure predictions with expected peak pressures for lesion estimation; (3) safety assessment with identification of potential off-target heating; and (4) time estimates with wave propagation timing for pulse sequencing.

Engineers can use these results for: (1) transducer design optimizing geometry for desired focal patterns; (2) frequency selection balancing penetration and absorption; (3) power settings determining safe and effective amplitude ranges; and (4) quality control verifying device performance predictions.

Scientists benefit from: (1) mechanism understanding with clear visualization of nonlinear effects; (2) model validation as benchmark for new computational methods; (3) parameter studies as framework for systematic investigation; and (4) theory testing as comparison platform for analytical predictions.

Limitations and Future Directions

While comprehensive, the current results also reveal areas for future work: (1) extension to heterogeneous media with spatially varying tissue properties; (2) implementation of full 3D simulations beyond axisymmetric geometry; (3) systematic parameter variation studies across therapeutic operating ranges; and (4) direct comparison with experimental measurements from calibrated hydrophone data.

CONCLUSION

We have introduced a novel conservative hyperbolic formulation for nonlinear acoustics in thermoviscous media. The model:

1. Is derived from first principles and retains full nonlinearity and dissipation without ad hoc assumptions.

2. Allows use of high-order shock-capturing schemes (WENO5/7) that provide accurate resolution of steep gradients without spurious oscillations.

3. Is validated across linear, weakly nonlinear, and strongly nonlinear regimes through comparison with analytical solutions, established numerical methods, and theoretical predictions.

4. Is efficiently implemented on GPUs for large-scale simulations, achieving 310 GFLOPS and enabling practical clinical application.

Future work will extend the model to heterogeneous media with spatially varying material properties and include local nonlinear effects relevant to cavitation and tissue damage mechanisms.

Clinical Impact

For surgical applications, the major outcomes are:

1. More accurate HIFU treatment planning: The model captures real nonlinear wave physics, allowing better prediction of peak focal pressures, thermal dose distribution, and shock-induced heating.

2. Safer surgical use of high-power ultrasound: Better modeling means improved ability to avoid nerve bundles, blood vessels, and sensitive organs near the target.

3.  Support for patient-specific simulation: Fast GPU computation makes personalized planning feasible within clinical timeframes.

4. Foundational model for next-generation HIFU systems: The conservative formulation is robust, generalizable, and compatible with 3-D and heterogeneous tissue extensions.

ACKNOWLEDGEMENTS

None.

SOURCE OF FUNDING

This work was not supported by any funding.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

REFERENCES

1. Jamalabadi MYA. (2025). A Comprehensive Review of Computational Modeling in Tumor and Brain Disorder with a Focus on Ablation Therapies. J Biomed Res Environ Sci. 6(10):1481-1491.
2. Westervelt PJ. (1963). Parametric acoustic array. The Journal of the Acoustical Society of America. 35(4):535-537.
3. Hamilton MF, Blackstock DT. (1998). Nonlinear Acoustics. Academic Press, San Diego.
4. Blackstock DT. (1966). Connection between the Fay and Fubini solutions for plane sound waves of finite amplitude. The Journal of the Acoustical Society of America. 39(6):1019-1026.
5. Aanonsen SI, Barkve T, Tjøtta JN, Sigve T. (1984). Distortion and harmonic generation in the nearfield of a finite amplitude sound beam. The Journal of the Acoustical Society of America. 75(3):749-768.
6. Berenger JP. (1994). A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics. 114(2):185-200.
7. Treeby BE, Jaros J, Rendell AP, Cox BT. (2012). Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method. J Acoust Soc Am. 131(6):4324-4336.
8. Jiang GS, Shu CW. (1996). Efficient implementation of weighted ENO schemes. Journal of Computational Physics. 126(1):202-228.
9. LeVeque RJ. (2002). Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, UK.
10. O'Neil H. (1949). Theory of focusing radiators. The Journal of the Acoustical Society of America. 21(5):516-526.
11. Bessonova OV, Wilkens V. (2013). Membrane hydrophone measurement and numerical simulation of HIFU fields up to developed shock regimes. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control. 60(2):290-300.

APPENDIX: KEY FORMULAS SUMMARY

Equation of State:

pₐ = c₀²ρₐ + (c₀²/ρ₀)(B/A/2)ρₐ² + (κ/ρ₀)(1/cv - 1/cp)∂ρₐ/∂t     (16)

Conservative System:

∂q/∂t + ∇·f(q) = ∇·(D∇q)     (17)

Eigenvalues:

λ₁,₂ = (1/2)[βuₐ ∓ √(4c₀² + 4βpₐ/ρ₀ + β²uₐ²)]     (18)

PML Damping:

σ(x) = log(1/R)(3c₀/2ΔPML)(x/ΔPML)²     (19)