Finite-difference time domain (FDTD) techniques are widely used to model the propagation of viscoelastic waves through complex and heterogeneous structures. However, in the specific case of media mixing liquid and solid, attempts to model continuous media onto a Cartesian grid produces errors when the liquid-solid interface between different media do not align precisely with the Cartesian grid. The increase in spatial resolution required to eliminate this grid staircasing effect can be computationally prohibitive. Here, a modification to the Virieux staggered-grid FDTD scheme called the superposition method is presented. This method is intended to reduce this staircasing effect while keeping a manageable computational time. The method was validated by comparing low-spatial-resolution simulations against simulations with sufficiently high resolution to provide reasonably accurate results at any incident angle. The comparison of the root-mean-square of the stress amplitude maps showed that the amplitude of artifactual waves could be reduced by several orders of magnitude when compared to the Virieux staggered-grid FDTD method and that the superposition method helped to significantly reduce the staircasing effect in FDTD simulations.

Finite-difference time domain (FDTD) methods have been used extensively to model the propagation of ultrasound in heterogeneous media.1 The use of FDTD simulations has been validated for planning transcranial ultrasound treatments,2 which is an emerging modality for treatment of brain disorders including thalamotomy for treatment of essential tremor3 and chronic neuropathic pain,4 and thermal ablation of brain tumours.5 Applying transcranial focused ultrasound presents several challenges due to the presence of the skull, which is a problematic biological tissue to focus through due to its strong acoustic mismatch with water, conversion between compressive and transverse waves, scattering, and attenuative properties, which gives rise to phase aberrations, standing waves, and undesired tissue heating.6–9 

The Virieux-staggered grid FDTD scheme10 is widely used to model the viscoelastic wave propagation in heterogeneous media11,12 due to its computational efficiency, smaller memory requirement, and easy implementation,13 and has been validated for the case of shear mode conversion.14 Staggered grid methods provide greater accuracy and stability than traditional FDTD methods, and form the basis of reverse-time migration and full-waveform inversion.12,15

A source of error in FDTD solutions of the viscoelastic wave equation arises when the liquid-solid interface does not align with the Cartesian grid.16 This staircase effect was first described by Fornberg.11,17 Staircasing refers to the spatial approximation required to define continuous geometries onto a discrete Cartesian grid in two and three dimensions (2D and 3D), resulting in many vertices or edge positions that do not correspond to points on the grid.1 This error is exacerbated by end-point misregistration, which occurs when the angle of the interface does not form a Pythagorean triangle on the grid, leading to line-source endpoints that are not coincident with specific grid point positions. Both the staircase representation of a continuous line and the endpoint misregistration are considered aspects of the staircase effect.1 The discrete discontinuities in the representation of the interface produces circular wavefronts emanating from points where the discrete staircase approximation causes the interface to appear “jagged” or discontinuous on the rectangular grid, which act as a source of unphysical noise and artifactual waves.17 

Previous approaches to satisfy the boundary conditions at a material interface include the use of a hybrid grid-/finite-difference- (FD-) method, where the grid method18 is applied in the fluid for its flexibility in handling complex fluid-solid interfaces, while the FD method is applied in the solid due to its lower numerical dispersion.9 Other approaches may account for surface topography, such as the technique developed by Tessmer et al.,19 which is based on a pseudo-spectral method that maps a rectangular grid onto a curved surface, but requires significantly more time and memory and is unstable for large surface curvature. The FD method on a rotated staggered grid20 uses rotated finite-difference operators to align the grid to the material interface. Finite differences on deformed grids21–24 are used for their relative simplicity and preservation of stencil computations,25 but are incompatible with the use of the fully staggered Virieux scheme, and so require more computationally expensive approaches. Lisitsa et al.25 developed a combined Galerkin/FD method that exploits the polyhedral mesh of the discontinuous Galerkin method in the regions around the interface in order to more accurately handle complex topography, while employing the more efficient FD method in the rest of the computational domain. The use of finite element methods provides flexibility for their ability to accurately represent complex geometries, but construction of irregular polyhedral meshes is a problem of high complexity, and the use of global finite element methods is more computationally expensive than FDTD methods on deformed grids.25 

Previous examinations have demonstrated that this staircasing effect decreases with finer spatial sampling, and that reasonably accurate results may be obtained from the standard FDTD method by using grid point spacing of 30 points per minimum wavelength (PPW).26 For 2D simulations, the memory requirements increase proportional to Δx–2, where Δx is the spatial step in the FDTD implementation, while computational time requirements increase proportional to Δx–3, which imposes limits on the use of high-resolution FDTD algorithms in practical applications.

Due to these limitations of these numerical schemes, there is a demand for solutions that can reduce the staircasing effect with manageable increases to the memory and time requirements, while maintaining the simplicity, accuracy, and efficiency of FDTD methods. To address this problem, we present a modification to the standard Virieux staggered-grid FDTD scheme that increases the accuracy of simulations without the constraining resources required for increased spatial resolution. Conceptually, we applied the principle of superposition for linear systems27 to represent the liquid-solid interface as a superposition of different configurations. Our results demonstrate that this approach is especially advantageous to address the problem of staircase artifacts in FDTD modelling.

The proposed method is inspired by the principle of superposition, which states that for linear systems the overall state may be broken into a sum of separate pieces, where the equation can be solved in each separate piece, the pieces may be added back together to obtain the solution for the whole.27 The goal of this method is to address the ambiguity of the liquid-solid interface by having the locations of the step discontinuities distributed at different locations along the interface in different simulations by using different approximations of the interface based on the location of the interface relative to the grid positions, and then combine the results of those different simulations to obtain a solution that contains contributions each point along the interface.

The superposition method aims to reduce the artifactual wavefronts produced by the discrete staircase approximation of angled or curved interfaces by recombining the waves from multiple simulations, where each simulation exhibits differing approximations of the interface based on the location of the interface relative to the grid positions. The application of the superposition method will be considered analytically in 1D, and numerically in 2D.

The staircasing effect in FDTD simulations arises when the location of the material interface does not align with the discrete grid used in numerical simulations. In order to examine the application of the superposition method analytically in 1D, the transmission and reflection of waves is considered from an interface that can only be defined at discrete points that are offset from the location of the physical interface. To begin, supposed that a boundary between two media exists at x = 0. Incident wave ψi at frequency f0 is propagating toward increasing x, giving rise to the reflected wave ψr and transmitted wave ψt. This scenario is shown in Fig. 1. The wave speed c is given by

c(x)={c1,if x<0,c2,if x0.
(1)
FIG. 1.

Analytical scenario for the superposition method in 1D.

FIG. 1.

Analytical scenario for the superposition method in 1D.

Close modal

The incident, reflected, and transmitted waves will then be given, respectively, by

ψi=exp(ik1x),
(2)
ψr=rexp(ik1x),
(3)
ψt=texp(ik2x),
(4)

where the time-dependence has been suppressed. The wave numbers are given by

ki=2πλi=2πf0ci
(5)

and the reflection and transmission coefficients satisfy

ψi*ψi=ψr*ψr+ψt*ψt=|r|2+|t|2=1.
(6)

The staircasing effect in FDTD simulations arises when the location of the material interface does not align with the discrete grid used in simulations. To approximate the staircase problem analytically in 1D, the location of the interface is restricted to discrete points that are offset from x = 0. When these points do not align with the location of the interface, the interface can either be located at the point that immediately precedes the location of the interface, or the point immediately after. For the first case, the location of the material interface is considered to be at position x=(1β)Δx, while in the second case the location of the material interface is located at x=βΔx. The quantity β can range from 0 to 1, and represents the fraction of the distance between the two grid points that is occupied by the solid. The length Δx is analogous to the spatial resolution in numerical simulations. The transmitted wave in scenarios 1 and 2 are denoted ψt1 and ψt2, respectively, and are shown in Fig. 2.

FIG. 2.

Two difference cases where the boundary between two media is offset from zero, where the transmitted wave ψt will be reconstructed from the weighted superposition of the two transmitted waves ψt1 and ψt2.

FIG. 2.

Two difference cases where the boundary between two media is offset from zero, where the transmitted wave ψt will be reconstructed from the weighted superposition of the two transmitted waves ψt1 and ψt2.

Close modal

The goal of the superposition method is to use a linear combination of the basis functions ψt1 and ψt2 to create a transmitted wave ψt that is mathematically more similar to the exact transmitted wave ψt than either of the basis functions. For 1D, the resultant transmitted wave will be given by

ψt=βψt1+(1β)ψt2.
(7)

Compared to waves that are transmitted through the interface at the actual location ψt, the transmitted waves in the situations shown above experience a phase shift, such that

ψt1=exp(iΔk(1β)Δx)ψt
(8)

and

ψt2=exp(iΔkβΔx)ψt,
(9)

where

Δk=k1k2.
(10)

To quantify the correlation between the reconstructed transmitted wave ψt and the analytical transmitted wave ψt, consider an error function et that is given by the difference between the reconstructed wave and the analytical transmitted wave, such that

et=ψtψt.
(11)

The relative error for the transmitted wave, εt, can be given as the magnitude of the error function, relative to the magnitude of the transmitted wave, which can be found using

ϵt=et*etψt*ψt=(ψt*ψt*)(ψtψt)ψt*ψt,
(12)

which for the transmitted wave reconstructed using the superposition method yields

ϵtSPP=(1β)2+2β(1β)cos(ΔkΔx)+β22(1β)cos(βΔkΔx)2βcos[(1β)ΔkΔx]+1.
(13)

This can be compared with the error ϵtst calculated from the stepwise approximation to the transmitted wave, where the location of the interface is assigned to whichever grid point is closer to the location of the actual interface, and the transmitted wave is given by

ψt={ψt1,if 0β12,ψt2,if 12<β1.
(14)

Solving for the magnitude ϵtst yields

ϵtst={4sin2(βΔkΔx2),if 0β12,4sin2((1β)ΔkΔx2),if 12<β1.
(15)

Because the stepwise approximation preserves the amplitude of the reconstructed waves and differs from the exact solution only by a phase shift, all error for the waves in the stepwise approximation can be attributed purely to phase error and not amplitude error. The same analysis was applied to quantify the reconstructed reflected wave error ϵrSPP which yields

ϵrSPP=(1β)2+2β(1β)cos(2k1Δx)+β22(1β)cos(2βk1Δx)2βcos[2(1β)k1Δx]+1.
(16)

This can be compared against the error for the reconstructed reflected wave in the stepwise approximation ϵrst, which yields

ϵrst={4sin2(βk1Δx),if 0β12,4sin2((1β)k1Δx),if 12<β1.
(17)

While the error for the transmitted wave reconstructed using the superposition method does not depend on the direction of propagation, the reflected wave error depends on the wave number in the incident/reflected medium, which for a given interface will depend on the direction in which the wave is incident on the boundary. This same analysis may be applied for a wave that is reflected from the opposite side of the boundary due to an incident wave propagating towards –x to obtain the error, ϵr. For the reflected wave reconstructed using the superposition method, this yields

ϵrSPP=(1β)2+2β(1β)cos(2k2Δx)+β22(1β)cos(2βk2Δx)2βcos[2(1β)k2Δx]+1,
(18)

while for the wave reflected from the opposite side of the interface calculated using the stepwise approximation, this yields

ϵrst={4sin2(βk2Δx),if 0β12,4sin2((1β)k2Δx),if 12<β1.
(19)

To evaluate the magnitude of the error, wave speeds were selected based on the properties of water and acrylic, such that c1=1482m/s and c2=2848m/s. To be consistent with numerical comparisons, the magnitude of Δx is given by 1/6th of the smallest wavelength, and a frequency f0=1MHz is selected, giving a Δx value of 0.247 mm. The magnitude of the transmitted and reflected wave error for the superposition method and the stepwise approximation are shown in Fig. 3. The error is lower for the superposition method than for the stepwise approximation for all wave components and all values of β, and converge with the exact solutions as the value of β approaches 0 or 1.

FIG. 3.

(Color online) 1D error magnitude vs β.

FIG. 3.

(Color online) 1D error magnitude vs β.

Close modal

To evaluate the effects of amplitude and phase error separately, the waves reconstructed using the superposition method given in Eq. (7) allows for the reflected and transmitted waves to be written in terms of the corresponding analytical wave as

ψi=Aijψi,
(20)

where Aij is a complex amplitude (with i=t,r,r, and j=SPP,st). The amplitude and phase of A for the superposition and piecewise approximations are shown in Fig. 4. The piecewise approximations to the various wave components only differ from the analytical solutions by a phase difference, so that the absolute value of the amplitude coefficient is always equal to 1, and so are not included in the figures. While the superposition approximation shows reduction in the overall amplitude of the various waves, the overall reduction in total error by the superposition method compared to the stepwise approximation shown in Fig. 3 can be attributed to the significant reduction in phase error achieved by the superposition method.

FIG. 4.

(Color online) Relative amplitude of superposition amplitudes (left), phase of superposition amplitude (center), and phase of piecewise amplitude (right) compared to the analytical waves.

FIG. 4.

(Color online) Relative amplitude of superposition amplitudes (left), phase of superposition amplitude (center), and phase of piecewise amplitude (right) compared to the analytical waves.

Close modal

The Virieux scheme is a finite difference method of modelling viscoelastic wave propagation on a fully staggered grid with the P-SV formulation that defines the viscoelastic equations given as a first-order hyperbolic system according to10,14,28

vxt=1ρ(σxxx+σxzz),
(21)
vzt=1ρ(σzzz+σxzx),
(22)
σxxt=[(λ+2μ)(1+τl)2μ(1τs)](vxx+vzz)+2μ(1+τs)vxx+rxx,
(23)
σzzt=[(λ+2μ)(1+τl)2μ(1τs)](vxx+vzz)+2μ(1+τs)vzz+rzz,
(24)

and

σxzt=μ(1+τs)(vxz+vzx)+rxz,
(25)

where ρ is the density, and λ and μ are the Lamé coefficients. The rij values are memory variables used to include the standard linear solid model of viscoelasticity, which are given by

rzzt=1τσ{[(λ+2μ)τl2μτs](vxx+vzz)+2μτsvzz+rzz},
(26)
rxxt=1τσ{[(λ+2μ)τl2μτs](vxx+vzz)+2μτsvxx+rxx},
(27)

and

rxzt=1τσ{μτs(vzx+vxz)+rxz},
(28)

where τl and τs are the stress relaxation times for the longitudinal and shear stress, respectively, and τσ=1/2πf0 with f0 given by the central frequency. The full procedure for optimizing material parameters based on speed of sound and attenuation coefficient is given in Ref. 28.

The viscoelastic wave equation was solved using fourth order in space and second order in time29–31 finite-difference operators. Harmonic averaging of the bulk and shear moduli and arithmetic averaging of the density was applied in all simulations,32 which provides greater accuracy across fluid-elastic boundaries even when the interface is aligned with the grid.33 The full procedure for the implementation of the FDTD algorithm is given in Ref. 31.

The superposition method is a modification to the Virieux FDTD method. It consists of running multiple simulations of the physical space, which will be referred to as “configurations,” where each configuration contains a representation of the physical space constructed with different staircase approximations of the material interfaces. Measurements of the stress/velocity are based on the average values across all configurations.

For the 1D solution, the error in the location of the interface is related to the quantity β, and there are only two possible waves that can be recombined, allowing for the precise definition of the superposition weighting values used in Eq. (7). In contrast, the error in the location of the numerical interface in 2D numerical simulations can vary continuously for each cell that is intersected by the interface. The superposition weighting is approximated by averaging multiple simulations of the same domain with a different staircase approximation for each configuration. For a cell that is intersected by the interface, the fraction of configurations in which that cell is designated a solid is proportional to the amount of space in the cell occupied by the solid.

For FDTD simulations, each grid point must be defined as either a solid or a fluid, with the corresponding material properties. In 2D this is done by creating a map of material properties, designated χ(x,z), which is assigned a value of 0 to denote the fluid or 1 to denote the solid. The physical properties at each point in the grid are determined by the value of the material map at that point, such that

ρ(x,z)={ρfluid,if χ(x,z)=0,ρsolid,if χ(x,z)=1
(29)

with a similar relation for the Lamé parameters and relaxation coefficients, which determines the physical quantities (ρ,λ,μ,τ) at each point in space.

For a given cell position, the fraction of configurations in which a cell is designated as the material that lies beyond the interface is proportional to how much of the space within that cell is beyond the line defining the physical interface. Implementation is done using a weighting map β(x,z) defined at each point in space, where the value of β is given by the fraction of space within each cell that is beyond the line of the interface, and building material maps with a configuration-specific threshold that must be met for the cell to be designated as the material that lies beyond the interface. The grid coordinates are considered to be the location of the longitudinal stresses in the cell, such that the cell at coordinates (x0, z0) is considered to range from x0Δx/2x<x0+Δx/2 and z0z<z0+Δx. For regions where only a single material is present, the values are 0 or 1. At a material boundary that does not align with the Cartesian grid used for the simulation, the weighting values fall between 0 and 1. This is shown in Fig. 5(a).

FIG. 5.

(Color online) (a) Weighted map β(x,z) showing the relative area occupied by the solid within each cell near the interface for a 15° tilt angle. Tilted line (red in color) denotes the location of the physical interface. (b) The different material maps for N = 4 produced from the weighting map β(x,z).

FIG. 5.

(Color online) (a) Weighted map β(x,z) showing the relative area occupied by the solid within each cell near the interface for a 15° tilt angle. Tilted line (red in color) denotes the location of the physical interface. (b) The different material maps for N = 4 produced from the weighting map β(x,z).

Close modal

For FDTD simulations, only a single χ map is required. For the superposition method, a unique map is created for each separate configuration. For a solid/fluid interface where the solid lies beyond the line of the interface, each element in the nth(n=1,2,,N) material map χn will be designated as the solid if the weighting value (β) satisfies

χn(x,z)={1,if β(x,z)nN+1,0,if β(x,z)>nN+1,
(30)

where N is the total number of configurations. This threshold ensures that the number of configurations where a given cell is designated as a specific material is approximately equal to the relative area that that material occupies in that cell. For the example shown in Fig. 5(a) with N = 4 configurations, the corresponding χn maps are shown in Fig. 5(b). For N = 1, the superposition method reduces to the unaltered Virieux FDTD scheme, with a material map constructed using Eq. (30) with a threshold of 0.5.

The FDTD method is applied to each configuration independently. The resulting values of pressure/stress are given by the average value across all configurations, for example the stress components σij it yields

σij(x,z,t)=1Nn=1Nσij(x,y,t,n)
(31)

and root-mean-square (RMS) stress amplitude maps used for error quantification are calculated according to according to

σRMS(x,z)=1Tt=1Tσxx2(x,y,t)+σzz2(x,y,t).
(32)

To determine the accuracy of the superposition method, the model was tested under three different conditions, as shown in Fig. 6. The first scenario uses a highly simplified arrangement with a half space occupied by water and a lossless solid. The second scenario uses the same fluid/solid half space arrangement but introduces attenuation in the solid half-space. For the third and final test, the scenario consists of transmission through layered attenuating media, which provides a more practical application of the method.

FIG. 6.

Setup of the computational domain for lossless solid/fluid half-space and attenuating solid/fluid half-space (left), and layered media (right). Source plane is located at z = 0. The solid is located so that the center of the solid is 5 cm from the source. Tilt angle θ is the angle between the surface normal and the z-axis. Detectors are located at the boundaries of the domain. Lx=0.10m,Lz=0.11m.

FIG. 6.

Setup of the computational domain for lossless solid/fluid half-space and attenuating solid/fluid half-space (left), and layered media (right). Source plane is located at z = 0. The solid is located so that the center of the solid is 5 cm from the source. Tilt angle θ is the angle between the surface normal and the z-axis. Detectors are located at the boundaries of the domain. Lx=0.10m,Lz=0.11m.

Close modal

The simulations for various values of incident angles (θ), spatial resolutions (Δx), and number of configurations (N) were compared to high spatial-resolution simulations. Low spatial resolutions (6, 8, 10, 12, and 14 PPW) were compared against high spatial resolution simulations (50 PPW), which provide sufficient spatial resolution to produce reasonably accurate results for all incident angles.26 

The physical scenario is shown in Fig. 6. A solid/fluid interface is located 5 cm from the source plane. The solid has the longitudinal/shear speed of sound, longitudinal/shear attenuation (αl/αs) and density of acrylic as determined experimentally in Ref. 14. For the lossless half space simulation, the wave speeds and density of the solid are unchanged, and attenuation is set to zero. The surrounding medium is a lossless liquid with the density and longitudinal speed of sound of water. These quantities, as well as the corresponding relaxation values, are given in Table I. For trials with the layered attenuating media, the thickness of the solid layer was selected to be 8 mm, which is within the measured range of thickness of human skulls.34,35

TABLE I.

Viscoelastic and density parameters used for simulated solid and liquid.

LiquidSolid
ρ(kgm3) 1000 1190 
cl(ms1) 1482 2848 
cs(ms1) — 1026 
αl(Npm1) 34.5 
αs(Npm1) — 163 
τl(ms) 53.9 
τs(ms) — 93.5 
LiquidSolid
ρ(kgm3) 1000 1190 
cl(ms1) 1482 2848 
cs(ms1) — 1026 
αl(Npm1) 34.5 
αs(Npm1) — 163 
τl(ms) 53.9 
τs(ms) — 93.5 

The size of the domain was 10 cm along x and 11 cm along z. Spatial step Δx was set to minλs/PPW based on a frequency of 1 MHz and shear speed of sound given in Table I, while the time step Δt was calculated to incorporate the effects of attenuation on the stability criteria36 and provided a Courant-Friedrichs-Lewy number of 0.7. The domain is surrounded by a perfectly matched layer37 (PML) to absorb outgoing waves. Detectors were placed at the edges of the domain (just before the start of the PML) in order to detect the outgoing wave. To vary the incident angle θ, the material interface was tilted at various angle values with respect to the z-axis.

The source was located at z = 0 m and generated waves propagating towards the solid material. The source plane was aligned with the Cartesian grid to ensure that no artifacts arise from staircase approximation of the source plane, which is not addressed by this method. The impact of staircasing on acoustic line-sources has previously been examined in Refs. 43 and 44. An unfocused source with a Gaussian amplitude distribution along x was used. To minimize the spectral width, the source signal S(x, t) in time was a Gaussian wave packet modulated at a central frequency f0 of 1 MHz, given by

S(x,t)=exp((xa0)2)exp((t3ττ)2)sin(2πf0t),
(33)

where τ=4/f0 and a0=1cm.

To quantify the staircasing effect, the 2D RMS stress amplitude maps were compared against the RMS stress amplitude maps from high-resolution simulations. High resolution data were downsampled using a nearest-neighbour method to match the same spatial points as those from low-resolution simulations. Each RMS amplitude map σRMS(x,z) was compared with the downsampled high-resolution simulations σ̂RMS(x,z) using the structural similarity index (SSIM), which has been validated for spatial model assessment due to its inclusion of mean, variance, and correlation between two continuously valued maps.38 The structural similarity is calculated over the simulation domain using39 

SSIM(x,z)=(2μxμz+C1μx2+μz2+C1)(2ηxηz+C2ηx2ηz2+C2)(ηxz+C3ηxηz+C3),
(34)

where μx and μz are the average pixel values over x and z, ηx and ηz are pixel value standard deviations along x and z, and ηxz is the pixel value covariance along x and z. C1, C2, and C3 are constants introduced to stabilize instabilities when any denominator is close to zero. The mean SSIM index is calculated over all points in space to evaluate the similarity over the entire map. The error was then quantified using the structural dissimilarity index (DSSIM), which is related to the SSIM as40 

DSSIM=(1SSIM),
(35)

which has a range in the interval of [0,2], where two identical maps will have a DSSIM value of zero, while maps that are negatives of each other will have a DSSIM value of 2. The mean square error (MSE) was also calculated between the RMS amplitude maps and the downsampled high resolution maps using the formula

Ξ=1Kxy(σRMS(x,z)σ̂RMS(x,z))2,
(36)

where K is the total number of spatial points. The MSE error was calculated for different values of incident angle, spatial resolution, and number of configurations. While the DSSIM provides a relative metric for correlation between the two maps, the RMS MSE provides a direct quantitative metric of the absolute difference between the two RMS amplitude distributions.

The mean square error ΞS was also calculated for the detectors along the boundaries of the domain using

Ξs=1Mxt(σzz(x,t)σ̂zz(x,t))2,
(37)

where M is the total number of samples and x is the spatial coordinates. The high-resolution detector signals were downsampled to match the space and time coordinates of the low-resolution detector values.

Simulation code was implemented using the cuda language (V9.0.176, Nvidia, Santa Clara, CA, USA) and executed in a high-performance computing cluster counting with graphic processor units (GPUs) (GP100, Nvidia, Santa Clara, CA, USA) available at Compute Canada facilities. matlab software package (R2019a, MathWorks, Natick, Massachusetts, USA) was used to prepare and analyze data for simulations.

Examples of the reduction in scattering artifacts in the lossless solid/fluid half-space can clearly be seen in the RMS stress amplitude map as shown in Fig. 7. The Virieux method introduces multiple scattered wave artifacts that are significantly reduced when the superposition method is applied.

FIG. 7.

(Color online) Log-scale RMS stress amplitude map for the lossless half-space for 5° (top row) and 40° (bottom row) tilt angles with a low-resolution (left), low-resolution superposition method (center), and high-resolution (right) FDTD simulations. Low-resolution simulations used 6 PPW, and high-resolution used 50 PPW. Twenty domains were used for the superposition simulations. The tilt angle of 40° exceeds the longitudinal critical angle of 31°.

FIG. 7.

(Color online) Log-scale RMS stress amplitude map for the lossless half-space for 5° (top row) and 40° (bottom row) tilt angles with a low-resolution (left), low-resolution superposition method (center), and high-resolution (right) FDTD simulations. Low-resolution simulations used 6 PPW, and high-resolution used 50 PPW. Twenty domains were used for the superposition simulations. The tilt angle of 40° exceeds the longitudinal critical angle of 31°.

Close modal

The reduction of scattering artifacts is reflected in the reduced MSE error values as a function of N shown in Fig. 8, and the reduced DSSIM as a function of N shown in Fig. 9. These results show that the spatial distribution of the reflected and transmitted waves is improved with the application of the superposition method for all incident angles and spatial resolutions. The increase in error with increasing value of N shown for some incident angles may be attributed in part to the reduction of amplitude in the reflected wave with increasing value of N that is caused by destructive interference between waves in different configurations. As the value of N increases, more unique solutions are present which gives rise to more destructive interference between configurations that reduces the overall amplitude, which manifests as a magnitude error in the reflection and transmission coefficients. While this interference introduces error in the magnitude of the reflection and transmission coefficients, the 1D superposition method Sec. II A indicates it is also the source of the significant reduction in imaginary error that is the source of scattering artifacts. This is consistent with the analytical consideration of the error arising from the staircase approximation and effects the superposition are presented in  Appendixes A and  B, respectively.

FIG. 8.

(Color online) Logarithmic MSE for the lossless half-space 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW in the lossless half-space. Dotted lines show the error values for the Virieux method (N = 1).

FIG. 8.

(Color online) Logarithmic MSE for the lossless half-space 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW in the lossless half-space. Dotted lines show the error values for the Virieux method (N = 1).

Close modal
FIG. 9.

(Color online) Logarithmic DSSIM for lossless half-space 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW. Dotted lines show the error values for the Virieux method (N = 1).

FIG. 9.

(Color online) Logarithmic DSSIM for lossless half-space 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW. Dotted lines show the error values for the Virieux method (N = 1).

Close modal

The MSE calculated for the detectors around the boundaries of the domain show more limited reduction of error with increasing value of N, as shown in Fig. 10. This can be explained by the time delay that is introduced in low-resolution simulations due positional error in the detector positions and to the effects of numerical dispersion, to which this error metric is particularly sensitive,26 and which is not addressed by the superposition method.

FIG. 10.

(Color online) Logarithmic MSE for detectors around boundary of domain vs the number of configurations in the lossless half-space, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 points PPW in the lossless half-space. Dotted lines show the error values for the Virieux method (N = 1).

FIG. 10.

(Color online) Logarithmic MSE for detectors around boundary of domain vs the number of configurations in the lossless half-space, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 points PPW in the lossless half-space. Dotted lines show the error values for the Virieux method (N = 1).

Close modal

The RMS stress amplitude maps are shown in Fig. 11, where the effects of attenuation in the solid can be seen to cause rapid attenuation of the transmitted wave. The MSE as a function of N is shown in Fig. 12, which illustrates that the superposition method is less effective in reducing the error when attenuation is present in the solid. This can be explained with reference to the attenuative effects in the solid, which reduces the overall contribution of transmitted waves to the MSE and causes the MSE to be more heavily weighted towards the error in the reflected waves. When compared with the results of the lossless half space test, this lends further support to the conclusion drawn from the 1D analytical case, where both amplitude and phase errors from the interface are most significant in the waves that are reflected from the interface. While the RMS MSE shows some inconsistency in the ability of the superposition method to reduce scattering artifacts, the DSSIM metric shown in Fig. 13 for the attenuating half space illustrate how the application of the superposition method improves structural characteristics in the spatial distribution of the reflected and transmitted waves.

FIG. 11.

(Color online) Log-scale RMS stress amplitude map for the attenuating half-space for 5° (top row) and 40° (bottom row) tilt angles with a low-resolution (left), low-resolution superposition method (center), and high-resolution (right) FDTD simulations. Low-resolution simulations used 6 PPW, and high-resolution used 50 PPW. Twenty domains were used for the superposition simulations. The tilt angle of 40° exceeds the longitudinal critical angle of 31°.

FIG. 11.

(Color online) Log-scale RMS stress amplitude map for the attenuating half-space for 5° (top row) and 40° (bottom row) tilt angles with a low-resolution (left), low-resolution superposition method (center), and high-resolution (right) FDTD simulations. Low-resolution simulations used 6 PPW, and high-resolution used 50 PPW. Twenty domains were used for the superposition simulations. The tilt angle of 40° exceeds the longitudinal critical angle of 31°.

Close modal
FIG. 12.

(Color online) Logarithmic MSE for attenuating half-space 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW in the lossless half-space. Dotted lines show the error values for the Virieux method (N = 1).

FIG. 12.

(Color online) Logarithmic MSE for attenuating half-space 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW in the lossless half-space. Dotted lines show the error values for the Virieux method (N = 1).

Close modal
FIG. 13.

(Color online) Logarithmic DSSIM for attenuating half-space 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW. Dotted lines show the error values for the Virieux method (N = 1).

FIG. 13.

(Color online) Logarithmic DSSIM for attenuating half-space 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW. Dotted lines show the error values for the Virieux method (N = 1).

Close modal

As was the case with the lossless half space, the MSE for the detectors at the domain boundaries show limited reduction in error with the application of the superposition method, which may be attributed to the same explanation as the detector MSE in the lossless half-space.

Examples of the reduction in scattering artifacts in the lossless solid/fluid half-space can clearly be seen in the RMS of the stress amplitude map as shown in Figs. 14–16. The Virieux method introduces multiple scattered wave artifacts that are significantly reduced when the superposition method is applied, even in cases where multiple interfaces are present.

FIG. 14.

(Color online) Log-scale RMS stress amplitude map for layered attenuating media for 5° (top row) and 40° (bottom row) tilt angles with a low-resolution (left), low-resolution superposition method (center), and high-resolution (right) FDTD simulations. Low-resolution simulations used 6 PPW, and high-resolution used 50 PPW. Twenty domains were used for the superposition simulations. The tilt angle of 40° exceeds the longitudinal critical angle of 31°.

FIG. 14.

(Color online) Log-scale RMS stress amplitude map for layered attenuating media for 5° (top row) and 40° (bottom row) tilt angles with a low-resolution (left), low-resolution superposition method (center), and high-resolution (right) FDTD simulations. Low-resolution simulations used 6 PPW, and high-resolution used 50 PPW. Twenty domains were used for the superposition simulations. The tilt angle of 40° exceeds the longitudinal critical angle of 31°.

Close modal
FIG. 15.

(Color online) Logarithmic MSE for layered attenuating media 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW in the lossless half-space. Dotted lines show the error values for the Virieux method (N = 1).

FIG. 15.

(Color online) Logarithmic MSE for layered attenuating media 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW in the lossless half-space. Dotted lines show the error values for the Virieux method (N = 1).

Close modal
FIG. 16.

(Color online) Logarithmic DSSIM for the layered attenuating media 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW. Dotted lines show the error values for the Virieux method (N = 1).

FIG. 16.

(Color online) Logarithmic DSSIM for the layered attenuating media 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW. Dotted lines show the error values for the Virieux method (N = 1).

Close modal

The MSE was also calculated for the RMS of stress amplitude map exclusively in the regions beyond the solid plate, where only the transmitted wave is present. These results are shown in Fig. 17, which shows a reduction in the MSE with increasing value of N due to the reduction of artifactual transmitted waves.

FIG. 17.

(Color online) MSE for the transmission region of the layered attenuating media 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW. Dotted lines show the error values for the Virieux method (N = 1).

FIG. 17.

(Color online) MSE for the transmission region of the layered attenuating media 2D RMS stress amplitude maps vs the number of configurations, for incident angles ranging from 5° to 45° and spatial resolutions ranging from 6 to 14 PPW. Dotted lines show the error values for the Virieux method (N = 1).

Close modal

The application of the superposition method significantly reduces the effect of staircasing in viscoelastic FDTD simulations. The case of the fluid/solid half-space provides a simplified model for verification of the ability of the superposition method to reduce scattering artifacts in 2D that arise from the staircase approximation. The attenuating half space demonstrates that the superposition method is compatible with the inclusion of viscoelastic relaxation for modelling attenuative effects. The third case of transmission through layered media in 2D demonstrates the applicability of more practical and complex scenarios, like transcranial focusing, which may include multiple interfaces.

The RMS of the stress amplitude maps show that the profile of the reflected and transmitted waves are less affected by the diffractive effects from the staircase artifacts. The reduction in artifactual scattered waves is reflected in the reduced DSSIM for the RMS stress amplitude maps, which show a reduction in error with increasing number of configurations N. While the RMS MSE metric also shows a reduction in error with the application of the superposition method, the improvements are generally more limited, especially when the MSE is calculated over the whole simulation domain, which may be a product of the ability of SSIM to incorporate structural information from continuous valued maps, where such additional information is helpful in spatial model assessment.38 

The differences between the SSIM and the MSE reflect the improved distribution of energy from reflected and transmitted waves, at the expense of precision in the energy carried in those waves. The ability of the superposition method to improve the distribution of energy while reducing wave amplitude can be seen in the difference maps showing the absolute error, along with the corresponding high-resolution RMS map, shown in Fig. 18 for 5° and 40° interface angles.

FIG. 18.

(Color online) Absolute difference in root mean squared stress intensity maps for low resolution Virieux (left) and low resolution superposition (center), along with the high resolution RMS map for comparison (right), for a 5° (top row) and 40° tilt angles.

FIG. 18.

(Color online) Absolute difference in root mean squared stress intensity maps for low resolution Virieux (left) and low resolution superposition (center), along with the high resolution RMS map for comparison (right), for a 5° (top row) and 40° tilt angles.

Close modal

The 1D case presented in Sec. II A illustrates how the superposition method significantly reduces phase error, while introducing an error in the magnitude of the transmission and reflection coefficients, and that the error is highest in the reflected wave due to the magnitude of the wavenumber in the incident medium that produces a larger phase difference between configurations. As shown in Fig. 12, the superposition method produces RMS MSE error greater than uncorrected case for some interface angles when the transmitted wave is rapidly attenuated in the solid. It was explained in Sec. IV B that the attenuative effects in the solid reduces the overall contribution of transmitted waves to the MSE, causing the error to be more heavily weighted towards the error in reflected waves. The MSE metric is considered to be sensitive to changes in amplitude, compared to structural changes (such as the elimination of the diffraction artifacts), which was among the reasons why the DSSIM was included. The increase in error with increasing N shown in Fig. 12 can be attributed to the reduction of the reflection coefficient caused by the interference between different configurations, the attenuation of the transmitted wave that minimizes its contribution to the overall error, and the sensitivity to changes in amplitude that the MSE metric exhibits.

 Appendix A provides an analytical examination of how phase errors arise from staircase approximations.  Appendix B uses the same analytical framework to examine how the superposition method reduces the imaginary component of the error, while introducing an error in magnitude. This is consistent with the results shown in Fig. 18 where scattering artifacts caused by phase error are significantly reduced, but the amplitude of the reflected and transmitted waves are also reduced. The results of the analytical approach in  Appendix B provides a basis to conclude that as the value of N increases, the magnitude error caused by the application of the superposition method will converge to a finite value, which decreases with increased spatial resolution.

As mentioned in Sec. IV A, the detector MSE showed limited improvements with the application of the superposition method. This may be attributed to time delays that can arise due the effects of numerical dispersion during propagation as well as to positional errors in the detector positions, to which this error metric is particularly sensitive,26 and which is not addressed by the superposition method. This delay can be seen in the pressure vs time in Fig. 19, where the transmitted waves in the low-resolution simulations experiencing a slight delay. Additional methods are available for further reduction of numerical dispersion, but aside from the higher order FD approximation for differential operators, they have not been included here.41 While detector MSE metrics are presented here for completeness, dispersion analysis has previously been excluded from examinations of interface errors.41 Despite this limited improvement of the detector MSE, the RMS map error metrics are more effective at capturing the erroneous spatial distribution of waves that are produced by the staircase approximation, and which the superposition method is intended to address.

FIG. 19.

(Color online) Pressure vs time for Virieux (solid blue) and high resolution (dashed red) taken at (x,z) = (0, 150 mm) for a 0° tilt angle in layered media.

FIG. 19.

(Color online) Pressure vs time for Virieux (solid blue) and high resolution (dashed red) taken at (x,z) = (0, 150 mm) for a 0° tilt angle in layered media.

Close modal

The memory required to implement the superposition method scales proportionately to the number of configurations. Results indicate that improvements in accuracy are minimal beyond N10. An advantage of the superposition method is that each different simulation configurations can be run to completion, and the final signal can be progressively reconstructed as each configuration is simulated. This is in contrast to high-resolution simulations, where all stress, pressure, and memory variables at all points in space must be stored in memory at once. Additionally, while the time required for a standard FDTD 2D simulation scales with Δx3, the total time required to complete a simulation for the superposition method scales with NΔx3, with a larger value for Δx. For the 50 PPW simulations, each stress, velocity, and memory variable component must be calculated at 2.12×1012 different points in space and time, which requires 2.67×107 values be stored in memory for each variable array. For the lowest resolution (6 PPW) simulations with N = 8, accounting for all the different configurations, the values of the stress, velocity, and memory variable component must be calculated at 3.17×1010 different points in space and time, requiring only 4.13×105 values to be stored in memory for each variable array at a time.

Both the Virieux FDTD and superposition methods show a significant reduction in the transmitted amplitude of shear waves at high incident angles when compared to high-resolution simulations. This is evident by the decreased amplitude of the transmitted waves in 40° RMS figures shown in Fig. 14, as well as an increased amplitude in the waves reflected from the first material interface. This error is especially pronounced at incident angles above the longitudinal critical angle, indicating a reduced transmission coefficient for shear waves. Plane wave analysis has previously indicated that the numerical reflection and transmission coefficient error is purely imaginary, implying that the amplitude of the coefficients remain the same but introduce a phase shift.42 These results suggest that this analysis may break down for shear wave conversion and that an error in the magnitude of the reflection/transmission coefficients is introduced due to the staircase effect.

Some inconsistencies in the relationship between the error and the spatial resolutions, and specifically the increased error at higher spatial resolutions, as can be seen in Fig. 17, may be explained by diffractive effects that arise from the periodicity of the staircase approximation, the misalignment of the staircase representation of the liquid-solid interface1 and the compounded effect of modelling multiple interfaces simultaneously, which has been shown to produce higher errors than simulations which include only a single interface.1 

Lack of error reduction for angles approaching 45° is an expected characteristic of the superposition method. The number of unique material map representations is limited by the number of possible positions available for the step discontinuity at a given value of x, which is proportional to 1/tanθ for angles between 0° and 45°. As N increases, many of the new material maps are exact duplicates of already existing map in other configurations. This is apparent in the 6 PPW N = 20 simulation, where there are only three unique material maps configurations for the top interface between all 20 configurations, while for simulations employing low θ values, all 20 material map configurations are unique. Cases of identical configurations are only expected to occur when the interface angle is constant and close to 45°. In these cases, the limited ability of the superposition method to improve results may be mitigated by the natural error minima that occurs at 45° due to the interface aligning with grid points.26 In realistic simulations, it is expected that curved interfaces will result in a range of interface angles, making the duplication of exact configurations unlikely.

For cases in which identical material map configurations are produced, an increased value of N will lead to the proportion of configurations in which a given cells is designated a specific material approached the exact fraction of space within that cell that the material occupies. Increasing the value of N is therefore expected to provide improved accuracy even when there is a limited number of unique material map configurations, though the improvements in accuracy will be limited compared to increases in N that include unique material maps. Redundant calculations may be avoided without altering the basis of the superposition method by checking for configurations that are identical, ensuring that these configurations are only used in a single simulation, and using a weighted averaging when recombining the final field.

The superposition method is effective in reducing the staircasing artifacts, and does not require any modifications to the implementation of the Virieux FDTD algorithm in individual simulations. For this reason, the superposition method is compatible with the application of other other numerical modifications, such as the spatial averaging developed by Moczo et al.32 

While the Virieux method has been verified to provide reasonably accurate results at all incident angles with 30 PPW spatial resolution,26 some artifactual scattered waves are still present at low incident angles, as shown in the maps in Fig. 20, which illustrate the difference between the high resolution RMS maps and the various low-resolution RMS maps. These scattered waves correspond to the first order diffraction maxima that are produced by the diffraction grating-like structure of the interface arising due to the periodicity of the staircased approximation. While these artifacts remain in high-resolution simulations, they are not present in low-resolution simulations where the superposition method is applied. While other artifacts remain in the low-resolution superposition RMS amplitude map, this result may indicate that at small incident angles the superposition method is capable of reducing the diffractive effects that arise from the staircasing error effect more effectively than increasing the spatial resolution, which is consistent with the larger reduction in error that occurred with smaller tilt angles. For this reason, the spatial resolution for high-resolution simulations used for comparison in this paper was increased from 30 PPW to 50 PPW, which was the maximum resolution available with current resources, and where no artifactual waves were observed.

FIG. 20.

(Color online) Log-scale RMS stress amplitude difference maps for 1° tilt angle for 6 PPW low-resolution, 6 PPW low-resolution superposition method, 30 PPW resolution, and the corresponding 50 PPW high-resolution RMS stress amplitude map. Ten configurations were used for the superposition simulation. Artifactual waves remain in 30 PPW simulations but are not present in 50 PPW simulations.

FIG. 20.

(Color online) Log-scale RMS stress amplitude difference maps for 1° tilt angle for 6 PPW low-resolution, 6 PPW low-resolution superposition method, 30 PPW resolution, and the corresponding 50 PPW high-resolution RMS stress amplitude map. Ten configurations were used for the superposition simulation. Artifactual waves remain in 30 PPW simulations but are not present in 50 PPW simulations.

Close modal

Significant reduction in the staircasing artifacts in viscoelastic FDTD simulation can be achieved by the application of the superposition method. The examples provided demonstrate a significant reduction in the artifactual scattering of waves at the solid/fluid interface, without any alterations to the FDTD method applied in any of the individual simulations.

This research was enabled in part by support provided by Compute Ontario and Compute Canada. We acknowledge the support of the Discovery Program of the Natural Sciences and Engineering Research Council of Canada (NSERC).

In this appendix, the effect of staircasing on the sampling of a propagating plane wave is considered analytically. To begin, a sensor line shown as the dotted line in scenario 1 in Fig. 21.

FIG. 21.

Two staircased approximation (solid line) of the same angled sensor line (dashed line). Scenario 1 has the staircase approximation of the sensor above the actual sensor line, while scenario 2 has the staircase approximation below the actual sensor line.

FIG. 21.

Two staircased approximation (solid line) of the same angled sensor line (dashed line). Scenario 1 has the staircase approximation of the sensor above the actual sensor line, while scenario 2 has the staircase approximation below the actual sensor line.

Close modal

An incident wave with wave vector k is propagation along ẑ towards the sensor line given by z=xtanθ. The field at that line, S0(x), is given by

S0(x)=exp(ik·r)=exp(ikxtanθ),
(A1)

where k=2π/λ.

In scenario 1 shown in Fig. 1, the tilted sensor line is approximated by a discontinuous staircase approximation with step size 2a along ẑ and step width w along x̂, shown as a solid line. The effect of staircasing introduces a position-dependent phase error, due to the spatial misalignment of the actual sensor line and the staircase-approximated sensor. For the staircase approximation to the exact sensor line, this phase delay of the sensor staircasing can be described by

S(x)=S0(x)ϕ(x),
(A2)

where

|ϕ(x)|=1for w<x<w.
(A3)

Based on the diagram in Fig. 21, the periodicity of the staircasing error creates a periodic phase error that from wxw can be given by

ϕ(x)=exp(ikawxika),
(A4)

which satisfies the requirement based on Fig. 21 that ϕ(w)=1 and ϕ(w)=exp(2ika).

In the rotated coordinate system (x̂,ẑ) shown in Fig. 21 with x̂ parallel with the sensor line, such that x=xcosθ, Eqs. (A2) and (A4) can be written as

S0(x)=exp(iksinθx)
(A5)

and

ϕ(x)=exp(ika)exp(ikawx),for w<x<w,
(A6)

where w is the distance along the sensor line from x=w to x = w, with

w=w2+(2a)2.
(A7)

The saw-tooth-like phase error characteristic of ϕ(x) caused by the periodicity of the staircase error will produces multiple discrete k-vector components along the x axis, which can be found by writing ϕ(x) in terms of its Fourier series components, which yields

ϕ(x)=exp(ika)n=sinc(ka+nπ)exp(inπwx)
(A8)

with

S(x)=ϕ(x)S0(x),
(A9)
=exp(ika)n=sinc(ka+nπ)exp((iksinθinπw)x).
(A10)

Solving for the wave number along the x̂ axis yields

kx,n=ksinθn=ksinθi+nπw,
(A11)

which is equivalent to the grating equation with a periodicity equal to 2w.45 

This demonstrates that the staircased approximation of the sensor line gives rise to artifacts that can be expressed as a discrete set of k-vector components along x. In the case of the staircased approximation to the sensor, the exact solution corresponds to the n = 0 order, and it can be shown that the amplitude of the zeroth order is given by sinc(ka), with all other orders correspond to staircasing artifacts. This same behaviour has been observed in the diffraction-grating-like behaviour of reflected waves when the interface between two media is approximated by a staircase approximation. This is illustrated in Fig. 22, which show the RMS wave amplitude map from a 3° interface with the incident wave removed, and is consistent with the behaviour predicted by the analytical diffraction equation.

FIG. 22.

(Color online) Logarithmic RMS stress amplitude map showing diffraction pattern of reflected waves due to the periodicity of the staircase approximation of the angled interface. Lines denote analytical diffraction orders. Ray arising due to specular diffraction is shown in black, higher diffraction orders are shown in red (n=±1,2). Incident wave has been removed. Simulation used θi=3°, 8 points per wavelength.

FIG. 22.

(Color online) Logarithmic RMS stress amplitude map showing diffraction pattern of reflected waves due to the periodicity of the staircase approximation of the angled interface. Lines denote analytical diffraction orders. Ray arising due to specular diffraction is shown in black, higher diffraction orders are shown in red (n=±1,2). Incident wave has been removed. Simulation used θi=3°, 8 points per wavelength.

Close modal

Consider now the case where the field is sampled at two different staircased approximations of the sensor plane, separated by a distance 2a, as shown in Fig. 21. In scenario 1, the staircased approximation of the sensor is above the sensor line, while in scenario 2 the staircased approximation is below the sensor line. The field at these planes can be given by

S1(x)=S0(x)ϕ1(x)
(B1)

and

S2(x)=S0(x)ϕ2(x).
(B2)

The case for scenario 1 has already been examined in  Appendix A, which gives

ϕ1(x)=exp(ika)n=sinc(ka+nπ)exp(inπwx).
(B3)

ϕ1(x) and ϕ2(x) differ by a positional difference along z that produces a phase shift, which allows ϕ2(x) to be written as

ϕ2(x)=exp(ika)n=sinc(ka+nπ)exp(inπwx).
(B4)

To apply the superposition method, first a superposition weighting factor β(x) is introduced, which gives the relative distance from the first staircased sensor to the actual sensor line

β(x)=x2w+12.
(B5)

This weighting factor is used as coefficients to combine the two fields S1(x) and S2(x) into the superposition-weighted field, S(x), given by

S(x)=[1β(x)]S1(x)+β(x)S2(x)={[1β(x)]ϕ1(x)+β(x)ϕ(x)}S0(x)

or

S(x)=ϕ(x)S(x)
(B6)

with

ϕ(x)=[1β(x)]ϕ1(x)+β(x)ϕ(x).
(B7)

In contrast with the individual fields which satisfied |ϕ1|=|ϕ2|=1 for all values of x, in this case

|ϕ(x)|1.
(B8)

Using the values functions for ϕ1(x),ϕ2(x), and β(x), it can be shown that ϕ(x) can be represented in a Fourier series as

ϕ(x)=cos(ka)n={sinc(ka+nπ)tan(ka)ka+nπ×[cos(ka+nπ)sinc(ka+nπ)]}exp(inπwx).
(B9)

For the zeroth order term, we can determine that

a0=cos(ka)sinc(ka).
(B10)

The two staircase approximations giving in Fig. 21 represent the two basis functions for the superposition method. However, scenarios 1 consistently overestimate the position of the sensor position along z, and scenario 2 consistently underestimates the position along z. A more accurate staircase approximation would intersect the interface at x = 0, and would provide

ϕ(x)=n=sinc(ka+nπ)exp(inπwx)
(B11)

with the first term amplitude, a0, given by

a0,1=sinc(ka).
(B12)

The parameters in these equations can be expressed in terms of relevant FDTD variables by substitution a=Δx/2,Δx=λ/PPW, and k=2π/λ. For both the staircased approximation and the superposition approximation, the amplitude of the zeroth order is reduced. In the staircase approximation, the reduction of intensity in the central order is due purely to the energy scattered into higher diffraction order modes. With the application of the superposition method, the combination of the different basis functions leads to destructive interference that reduces the amplitude of all modes, but disproportionately affects the n > 0 modes. The absolute, real, and imaginary components of the amplitude for |n|10 modes are shown in Fig. 23(a), 23(b), and 23(c). The effect on the n = 0 mode amplitude as a function of spatial resolution is shown in Fig. 23(d).

FIG. 23.

(Color online) Absolute amplitude (a), real component of the amplitude (b), and imaginary component of the amplitude (c) for the staircase approximation (red) and the superposition approximation (blue) amplitude functions vs diffraction order with 8 PPW. Absolute amplitude of the zeroth order is shown as a function of wavelength in (d).

FIG. 23.

(Color online) Absolute amplitude (a), real component of the amplitude (b), and imaginary component of the amplitude (c) for the staircase approximation (red) and the superposition approximation (blue) amplitude functions vs diffraction order with 8 PPW. Absolute amplitude of the zeroth order is shown as a function of wavelength in (d).

Close modal
1.
J. L. B.
Robertson
,
B. T.
Cox
,
J.
Jaros
, and
B. E.
Treeby
, “
Accurate simulation of transcranial ultrasound propagation for ultrasonic neuromodulation and stimulation
,”
J. Acoust. Soc. Am.
141
(
3
),
1726
1738
(
2017
).
2.
F.
Marquet
,
M.
Pernot
,
J. F.
Aubry
,
G.
Montaldo
,
L.
Marsac
,
M.
Tanter
, and
M.
Fink
, “
Non-invasive transcranial ultrasound therapy based on a 3D CT scan: Protocol validation and in vitro results
,”
Phys. Med. Biol.
54
(
9
),
2597
2613
(
2009
).
3.
W. J.
Elias
,
N.
Lipsman
,
W. G.
Ondo
,
P.
Ghanouni
,
Y. G.
Kim
,
W.
Lee
,
M.
Schwartz
,
K.
Hynynen
,
A. M.
Lozano
,
B. B.
Shah
,
D.
Huss
,
R. F.
Dallapiazza
,
R.
Gwinn
,
J.
Witt
,
S.
Ro
,
H. M.
Eisenberg
,
P. S.
Fishman
,
D.
Gandhi
,
C. H.
Halpern
,
R.
Chuang
,
K.
Butts Pauly
,
T. S.
Tierney
,
M. T.
Hayes
,
G. R.
Cosgrove
,
T.
Yamaguchi
,
K.
Abe
,
T.
Taira
, and
J. W.
Chang
, “
A randomized trial of focused ultrasound thalamotomy for essential tremor
,”
New England J. Med.
375
(
8
),
730
739
(
2016
).
4.
D.
Jeanmonod
,
B.
Werner
,
A.
Morel
,
L.
Michels
,
E.
Zadicario
,
G.
Schiff
, and
E.
Martin
, “
Transcranial magnetic resonance imaging–guided focused ultrasound: Noninvasive central lateral thalamotomy for chronic neuropathic pain
,”
Neurosurg. Focus
32
(
1
),
E1
E11
(
2011
).
5.
N.
McDannold
,
G.
Clement
,
P.
Black
,
F.
Jolesz
, and
K.
Hynynen
, “
Transcranial MRI-guided focused ultrasound surgery of brain tumors: Initial findings in three patients
,”
Neurosugery
66
(
2
),
323
332
(
2010
).
6.
J.-F.
Aubry
,
M.
Tanter
,
M.
Pernot
,
J.-L.
Thomas
, and
M.
Fink
, “
Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans
,”
J. Acoust. Soc. Am.
113
(
1
),
84
93
(
2003
).
7.
G. T.
Clement
and
K.
Hynynen
, “
A non-invasive method for focusing ultrasound through the human skull
,”
Phys. Med. Biol.
47
(
8
),
1219
1236
(
2002
).
8.
G.
Pinton
,
M.
Pernot
,
E.
Bossy
,
J. F.
Aubry
,
M.
Muller
, and
M.
Tanter
, “
Mechanisms of attenuation and heating dissipation of ultrasound in the skull bone: Comparison between simulation models and experiments
,” in
Proceedings—IEEE Ultrasonics Symposium
(
2010
), pp.
225
228
.
9.
A.
Pulkkinen
,
B.
Werner
,
E.
Martin
, and
K.
Hynynen
, “
Numerical simulations of clinical focused ultrasound functional neurosurgery
,”
Phys. Med. Biol.
59
(
7
),
1679
1700
(
2014
).
10.
J.
Virieux
, “
P-SV wave propagation in heterogeneous media: Velocity–stress finite–difference method
,”
Geophysics
51
(
4
),
889
901
(
1986
).
11.
P.
Moczo
,
J. O.
Robertsson
, and
L.
Eisner
, “
The finite-difference time-domain method for modeling of seismic wave propagation
,”
Adv. Geophys.
48
,
421
516
(
2007
).
12.
W.
Liang
,
X.
Li
,
C.
Yang
,
Z.
Nashed
,
Y.
Wang
, and
G.
Liang
, “
Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method
,”
Geophysics
79
(
5
),
T277
T285
(
2014
).
13.
L.
Etemadsaeed
,
P.
Moczo
,
J.
Kristek
,
A.
Ansari
, and
M.
Kristekova
, “
A no-cost improved velocity-stress staggered-grid finite-difference scheme for modelling seismic wave propagation
,”
Geophys. J. Int.
207
(
1
),
481
511
(
2016
).
14.
S.
Pichardo
,
C.
Moreno-Hernández
,
R.
Andrew Drainville
,
V.
Sin
,
L.
Curiel
, and
K.
Hynynen
, “
A viscoelastic model for the prediction of transcranial ultrasound propagation: Application for the estimation of shear acoustic properties in the human skull
,”
Phys. Med. Biol.
62
(
17
),
6938
6962
(
2017
).
15.
Y.
Liu
, “
Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling
,”
Geophys. J. Int.
197
(
2
),
1033
1047
(
2014
).
16.
T. G.
Jurgens
,
A.
Taflove
,
K.
Umashankar
, and
T. G.
Moore
, “
Finite-difference time-domain modeling of curved surfaces
,”
IEEE Trans. Ant. Propag.
40
(
4
),
357
366
(
1992
).
17.
B.
Fornberg
, “
The pseudo-spectral method: Accurate representation in elastic wave calculations
,”
Geophysics
53
(
5
),
625
729
(
1988
).
18.
J.
Zhang
, “
Wave propagation across fluid-solid interfaces: A grid method approach
,”
Geophys. J. Int.
159
(
1
),
240
252
(
2004
).
19.
E.
Tessmer
,
D.
Kosloff
, and
A.
Behle
, “
Elastic wave propagation simulation in the presence of surface topography
,”
Geophys. J. Int.
108
(
2
),
621
632
(
1992
).
20.
E. H.
Saenger
and
T.
Bohlen
, “
Finite–difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid
,”
Geophysics
69
(
2
),
583
591
(
2004
).
21.
I.
Tarrass
,
L.
Giraud
, and
P.
Thore
, “
New curvilinear scheme for elastic wave propagation in presence of curved topography
,”
Geophys. Prospect.
59
(
5
),
889
906
(
2011
).
22.
B.
Lombard
,
J.
Piraux
,
C.
Gélis
, and
J.
Virieux
, “
Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves
,”
Geophys. J. Int.
172
(
1
),
252
261
(
2008
).
23.
S.
Hestholm
, “
Three-dimensional finite difference viscoelastic wave modelling including surface topography
,”
Geophys. J. Int.
139
(
3
),
852
878
(
1999
).
24.
W.
Zhang
,
Z.
Zhang
, and
X.
Chen
, “
Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids
,”
Geophys. J. Int.
190
(
1
),
358
378
(
2012
).
25.
V.
Lisitsa
,
V.
Tcheverda
, and
C.
Botter
, “
Combination of the discontinuous Galerkin method with finite differences for simulation of seismic wave propagation
,”
J. Comput. Phys.
311
,
142
157
(
2016
).
26.
R.
van Vossen
,
J. O. A.
Robertsson
, and
C. H.
Chapman
, “
Finite-difference modeling of wave propagation in a fluid–solid configuration
,”
Geophysics
67
(
2
),
618
(
2002
).
27.
R. P.
Feynman
,
R. B.
Leighton
, and
M.
Sands
, “
The Feynman lectures on physics: Mainly mechanics, radiation, and heat
,”
Phys. Today
17
(
8
),
45
(
1965
).
28.
J. O.
Blanch
,
J. O. A.
Robertsson
, and
W. W.
Symes
, “
Modeling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique
,”
Geophysics
60
(
1
),
176
184
(
2002
).
29.
A. R.
Levander
, “
Fourth–order finite–difference P-SV seismograms
,”
Geophysics
53
(
11
),
1425
1436
(
2002
).
30.
J. O. A.
Robertsson
,
J. O.
Blanch
, and
W. W.
Symes
, “
Viscoelastic finite–difference modeling
,”
Geophysics
59
(
9
),
1444
1456
(
2002
).
31.
T.
Bohlen
, “
Parallel 3-D viscoelastic finite-difference seismic modelling
,”
Comput. Geosci.
28
,
887
899
(
2002
).
32.
P.
Moczo
,
J.
Kristek
,
V.
Vavryčuk
,
R. J.
Archuleta
, and
L.
Halada
, “
3D heterogeneous staggered-grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities
,”
Bull. Seismol. Soc. Am.
92
(
8
),
3042
3066
(
2002
).
33.
O.
Podgornova
and
V.
Lisitsa
, “
Accuracy analysis of finite–difference staggered–grid numerical schemes for elastic–elastic and fluid–elastic interfaces
,” in
Society of Exploration Geophysicists
(
2010
), pp.
3087
3091
.
34.
O. P.
Murty
, “
Variability in thickness of human skull bones and sternum—An autopsy experience
,”
J. Forensic Med. Toxicol.
26
(
2
),
26
31
(
2009
).
35.
J.
Ruan
and
P.
Prasad
, “
The effects of skull thickness variations on human head dynamic impact responses
,” SAE Technical Paper Series (
2018
), Vol. 1(June).
36.
C.
Sun
,
Y.
Xiao
,
X.
Yin
, and
H.
Peng
, “
Stability condition of finite difference solution for viscoelastic wave equations
,”
Earthq. Sci.
22
(
5
),
479
485
(
2009
).
37.
F.
Collino
and
C.
Tsogka
, “
Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media
,”
Geophysics
66
(
1
),
294
(
2001
).
38.
C.
Robertson
,
J. A.
Long
,
F. S.
Nathoo
,
T. A.
Nelson
, and
C. C.
Plouffe
, “
Assessing quality of spatial models using the structural similarity index and posterior predictive checks
,”
Geograph. Anal.
46
(
1
),
53
74
(
2014
).
39.
G. P.
Renieblas
,
A. T.
Nogués
,
A. M.
González
,
N.
Gómez-Leon
, and
E. G.
del Castillo
, “
Structural similarity index family for image quality assessment in radiological images
,”
J. Med. Imag.
4
(
3
),
035501
(
2017
).
40.
R.
Karamichalis
,
L.
Kari
,
S.
Konstantinidis
, and
S.
Kopecki
, “
An investigation into inter- and intragenomic variations of graphic genomic signatures
,”
BMC Bioinf.
16
(
1
),
1
22
(
2015
).
41.
D.
Vishnevsky
,
V.
Lisitsa
,
V.
Tcheverda
, and
G.
Reshetova
, “
Numerical study of the interface errors of finite-difference simulations of seismic waves
,”
Geophysics
79
(
4
),
T219
T232
(
2014
).
42.
W. W.
Symes
and
T.
Vdovina
, “
Interface error analysis for numerical wave propagation
,”
Comput. Geosci.
13
(
3
),
363
371
(
2009
).
43.
Elliott S.
Wise
,
James L. B.
Robertson
,
Ben T.
Cox
, and
Bradley E.
Treeby
, “
Staircase-free acoustic sources for grid-based models of wave propagation
,” in
IEEE International Ultrasonics Symposium, IUS
(
2017
).
44.
Eleanor
Martin
,
Yan To
Ling
, and
Bradley E.
Treeby
, “
Simulating focused ultrasound transducers using discrete sources on regular Cartesian grids
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Contr.
63
(
10
),
1535
1542
(
2016
).
45.
Francisco J.
Duarte
,
Tunable Laser Optics
(
Elsevier
,
Amsterdam
,
2003
).