^{1}, Oscar P. Bruno

^{2,a)}, Theresa Y. Cheung

^{3}and Robin O. Cleveland

^{4}

### Abstract

On the basis of recently developed Fourier continuation (FC) methods and associated efficient parallelization techniques, this text introduces numerical algorithms that, due to very low dispersive errors, can accurately and efficiently solve the types of nonlinear partial differential equation(PDE)models of nonlinear acoustics in hundred-wavelength domains as arise in the simulation of focused medical ultrasound. As demonstrated in the examples presented in this text, the FC approach can be used to produce solutions to nonlinear acousticsPDEsmodels with significantly reduced discretization requirements over those associated with finite-difference, finite-element and finite-volume methods, especially in cases involving waves that travel distances that are orders of magnitude longer than their respective wavelengths. In these examples, the FC methodology is shown to lead to improvements in computing times by factors of hundreds and even thousands over those required by the standard approaches. A variety of one-and two-dimensional examples presented in this text demonstrate the power and capabilities of the proposed methodology, including an example containing a number of scattering centers and nonlinear multiple-scattering events.

N.A. gratefully acknowledges the support of the NSF postdoctoral fellowship. N.A. and O.P.B. gratefully acknowledge support from NSF and AFOSR. T.Y.C. and R.O.C. gratefully acknowledge support from NSF through grant 0835795.

I. INTRODUCTION

II. FOURIER CONTINUATION FUNDAMENTALS: FC(GRAM) ALGORITHM

III. FC SOLVERS IN ONE-DIMENSIONAL SPACE

A. Block domain decomposition

B. Application of the FC derivative operator

IV. FC SOLVERS: ONE-DIMENSIONAL CASE STUDIES

A. Case study I: Linear advection equation

B. Case study II: Nonlinear viscous Burgers equation

V. FC SOLVERS IN MULTIDIMENSIONAL SPACE

A. Overlapping curvilinear patches

B. Parallelization techniques

VI. TWO-DIMENSIONAL APPLICATION: COMPARISON OF THE NAVIER-STOKES AND KZK MODELS

A. Model and equations

1. Ultrasound transducermodel

2. Nonlinear acousticsNavier-Stokes equations

3. KZK equation

B. KZK solver

C. Navier-Stokes solver

1. Reflection symmetry

2. Spatial and temporal windowing

D. Results of Navier-Stokes/KZK comparisons

VII. TWO-DIMENSIONAL APPLICATION: COMPLEX GEOMETRIES AND FULL-DOMAIN SIMULATIONS

VIII. DISCUSSION AND CONCLUSIONS

A. Three-dimensional simulation

B. Advanced tissue models

C. Heterogeneous media

D. Adaptive mesh refinement

E. Solutions containing shocks

### Key Topics

- Transducers
- 35.0
- Navier Stokes equations
- 18.0
- Acoustic waves
- 16.0
- Partial differential equations
- 16.0
- Acoustic modeling
- 14.0

##### A61B8/00

## Figures

An illustration of the main elements of the FC method: after producing a translated copy of the given values of *f* on the interval [*b,* 1 + *b*], the algorithm proceeds to (1) obtain a smooth function that blends the original function on the interval [0, 1] to zero in the interval [1, *b*], and, similarly, to obtain a smooth function that blends the original (translated) function on the interval [*b,* 1 + *b*] to zero in the interval [1, *b*]; and to (2) add these two blending functions to obtain a smooth continuation. The original function values are represented by solid circles and the additional continuation values are represented by open circles. The thick curves denote the polynomial approximations in the matching intervals. The thin gray lines show the blending functions that are summed to produce the continuation values. The Fourier continuation *f ^{c} * then results by a direct application of the FFT algorithm.

An illustration of the main elements of the FC method: after producing a translated copy of the given values of *f* on the interval [*b,* 1 + *b*], the algorithm proceeds to (1) obtain a smooth function that blends the original function on the interval [0, 1] to zero in the interval [1, *b*], and, similarly, to obtain a smooth function that blends the original (translated) function on the interval [*b,* 1 + *b*] to zero in the interval [1, *b*]; and to (2) add these two blending functions to obtain a smooth continuation. The original function values are represented by solid circles and the additional continuation values are represented by open circles. The thick curves denote the polynomial approximations in the matching intervals. The thin gray lines show the blending functions that are summed to produce the continuation values. The Fourier continuation *f ^{c} * then results by a direct application of the FFT algorithm.

A simplified illustration of the spatial discretization used by the FC-based solver for the advection Eq. (6). (a) *N* = 30 grid points are evenly divided into *M* = 3 subdomains. (b) The curve represents a function that is to be differentiated by means of the FC approach. The derivative values at the points in the right-most subdomain [black circles in (a)] are computed using the function values at these points as well as external fringe region points [gray circles in (a)]; the function values on this augmented set ofpoints is represented by the thicker portion of the curve in (b). This (non-periodic) function is differentiated by means of the FC differentiation algorithm. The derivative values in the remaining two subdomains are produced analogously.

A simplified illustration of the spatial discretization used by the FC-based solver for the advection Eq. (6). (a) *N* = 30 grid points are evenly divided into *M* = 3 subdomains. (b) The curve represents a function that is to be differentiated by means of the FC approach. The derivative values at the points in the right-most subdomain [black circles in (a)] are computed using the function values at these points as well as external fringe region points [gray circles in (a)]; the function values on this augmented set ofpoints is represented by the thicker portion of the curve in (b). This (non-periodic) function is differentiated by means of the FC differentiation algorithm. The derivative values in the remaining two subdomains are produced analogously.

(Top left) Initial conditions for the velocity *u* in the one-dimensional advection equation example presented in Sec. IV A. (Top right) Close-up view of the initial condition near *x* = 0. (Bottom left and right) Errors vs computing time and discretization size, respectively at time *t* = 200 for the advection equation example presented in Sec. IV A. The graphs indicate that the FC and CL(4) computing-time curves reach the 10^{−2}-accuracy level at about 2.7 s (using 1400 discretization points) and at 1300 s (using 116 000 cells), respectively.

(Top left) Initial conditions for the velocity *u* in the one-dimensional advection equation example presented in Sec. IV A. (Top right) Close-up view of the initial condition near *x* = 0. (Bottom left and right) Errors vs computing time and discretization size, respectively at time *t* = 200 for the advection equation example presented in Sec. IV A. The graphs indicate that the FC and CL(4) computing-time curves reach the 10^{−2}-accuracy level at about 2.7 s (using 1400 discretization points) and at 1300 s (using 116 000 cells), respectively.

(Top left) Initial conditions for the velocity *u* in the one-dimensional (non-linear) Burgers equation example presented in Sec. IV. (Top right) Comparison of the corresponding highly resolved, converged FC solution at time *t* = 30 (solid line, obtained numerically, as indicated in the text, with a relative accuracy of less than 10^{−5}) and the corresponding FC solution (shown in black dots), obtained using a 11 500-point discretization. The relative maximum error in the latter FC solution is 1.2 × 10^{−2}. (Bottom left and right) Errors vs computing time and discretization size respectively at time *t* = 30 for the Burgers equation example presented in Sec. IV. The FC and CL(4) computing-time curves reach the 10^{−2}-accuracy level at 106 s (using 12 100 discretization points) and at 2700 s (using 138 000 cells), respectively.

(Top left) Initial conditions for the velocity *u* in the one-dimensional (non-linear) Burgers equation example presented in Sec. IV. (Top right) Comparison of the corresponding highly resolved, converged FC solution at time *t* = 30 (solid line, obtained numerically, as indicated in the text, with a relative accuracy of less than 10^{−5}) and the corresponding FC solution (shown in black dots), obtained using a 11 500-point discretization. The relative maximum error in the latter FC solution is 1.2 × 10^{−2}. (Bottom left and right) Errors vs computing time and discretization size respectively at time *t* = 30 for the Burgers equation example presented in Sec. IV. The FC and CL(4) computing-time curves reach the 10^{−2}-accuracy level at 106 s (using 12 100 discretization points) and at 2700 s (using 138 000 cells), respectively.

(Color online) (a) An illustration of the domain decomposition used for a parallel FC solver. The dashed lines represent processor boundaries. In order to update the solution values in the central subdomain (represented by dark circles), the associated processor receives a thin layer of shared values (represented by light circles) from each of four neighboring processors. (b)The treatment of the mirror points across the line of symmetry *x* = 0. Values are copied from the mirror image values (on the same processor) above the axis according to the symmetry relations (20)–(22). Additionally, the vertical component of the velocity is set to 0 on the axis as required by (22). (c) A portion of a grid covered by overlapping curvilinear patches. (d) As indicated by solid filled regions, one of the patches subdivided into subpatches for distribution in a parallel computing environment.

(Color online) (a) An illustration of the domain decomposition used for a parallel FC solver. The dashed lines represent processor boundaries. In order to update the solution values in the central subdomain (represented by dark circles), the associated processor receives a thin layer of shared values (represented by light circles) from each of four neighboring processors. (b)The treatment of the mirror points across the line of symmetry *x* = 0. Values are copied from the mirror image values (on the same processor) above the axis according to the symmetry relations (20)–(22). Additionally, the vertical component of the velocity is set to 0 on the axis as required by (22). (c) A portion of a grid covered by overlapping curvilinear patches. (d) As indicated by solid filled regions, one of the patches subdivided into subpatches for distribution in a parallel computing environment.

A circular arc transducer with focal length *F,* aperture radius *a* and radius of curvature .

A circular arc transducer with focal length *F,* aperture radius *a* and radius of curvature .

(Color online) Stages of the FC focused transducer solver. In Stage I, the computational window is fixed in space next to the transducer (left image). The computation proceeds in Stage I until the signal is centered within the computational window (center image). At this point, the solver enters Stage II computation and the computational window begins moving away from the transducer (right image).

(Color online) Stages of the FC focused transducer solver. In Stage I, the computational window is fixed in space next to the transducer (left image). The computation proceeds in Stage I until the signal is centered within the computational window (center image). At this point, the solver enters Stage II computation and the computational window begins moving away from the transducer (right image).

Comparisons of the KZK and Navier-Stokes (DNS) solutions for the 5*λ* aperture radius transducer. (a) Peak positive on-axis pressure. (b) Peak negative on-axis pressure. Axial pressure waveform (c) 50*λ* and (d) 80*λ* from the transducer.

Comparisons of the KZK and Navier-Stokes (DNS) solutions for the 5*λ* aperture radius transducer. (a) Peak positive on-axis pressure. (b) Peak negative on-axis pressure. Axial pressure waveform (c) 50*λ* and (d) 80*λ* from the transducer.

(Color online) Comparisons of the amplitudes of the first five harmonics (top to bottom) for the 5*λ* aperture radius transducer. (Left) Results obtained from the Navier-Stokes FC simulation. (Right) Results obtained from the KZK model. Horizontal axes are measured in focal lengths. Vertical axes are measured in aperture radii.

(Color online) Comparisons of the amplitudes of the first five harmonics (top to bottom) for the 5*λ* aperture radius transducer. (Left) Results obtained from the Navier-Stokes FC simulation. (Right) Results obtained from the KZK model. Horizontal axes are measured in focal lengths. Vertical axes are measured in aperture radii.

Additional comparisons of the KZK and Navier-Stokes (DNS) solutions for transducers with various aperture radii. Peak positive on-axis pressure for the (a) 10*λ*, (b)20*λ*, and (c) 30*λ* aperture radius transducers. (d) Axial pressure waveform 50*λ* from the transducer for the 30*λ* transducer.

Additional comparisons of the KZK and Navier-Stokes (DNS) solutions for transducers with various aperture radii. Peak positive on-axis pressure for the (a) 10*λ*, (b)20*λ*, and (c) 30*λ* aperture radius transducers. (d) Axial pressure waveform 50*λ* from the transducer for the 30*λ* transducer.

(Color online) Comparisons of the amplitudes of the first five harmonics (top to bottom) for the 30*λ* aperture radius transducer. (Left) Results obtained from theNavier-Stokes FC simulation. (Right) Results obtained from the KZK model. Horizontal axes are measured in focal lengths. Vertical axes are measured in aperture radii.

(Color online) Comparisons of the amplitudes of the first five harmonics (top to bottom) for the 30*λ* aperture radius transducer. (Left) Results obtained from theNavier-Stokes FC simulation. (Right) Results obtained from the KZK model. Horizontal axes are measured in focal lengths. Vertical axes are measured in aperture radii.

(Color online) Pressure field for focusing transducer simulations. (a) Propagation and focusing in a homogeneous medium. (b) Propagation and focusing across an array of scatterers.

(Color online) Pressure field for focusing transducer simulations. (a) Propagation and focusing in a homogeneous medium. (b) Propagation and focusing across an array of scatterers.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content