I gave up on this approach, when I could not find stable parameters for the explicit formulation of the scalar acoustic wave propagation in viscous material. Building the wave equation from basic operators lead to a rather complex condition for stability, as I will show below. I did not expect the stability analysis for the compounded wave operator to be more complex than the analysis for the standard implementation. Anyway, here is what I found until I abandoned the project:

The derivation of the acoustic wave equation is based on the
combination of three equations that have physical meaning by themselves.
The **continuity equation** describes the excess mass flow
due to the divergence of , the local density and
material velocity:

(46) |

(47) |

(48) |

The three equations do constitute the wave equation. The substitution of in the equation of state by the corresponding expression in the equation of continuity 48, yields

The partial time derivative of the expression is where I interchanged the order of differentiation assuming well-behaved solutions. Substitution of by an expression derived from the equation of motion 50 yields Finally, the partial time derivative of this expression and the divergence of the equation of motion 50 allow us to eliminate the velocity term and yield the viscous wave equation for the material's pressureInstead of implementing the wave equation directly, I decided to implement a system of three partial differential equations of which each equation has its own meaning and application. The elastic part of equation of state 49 and the continuity equation 48 combine to

(49) |

(50) |

The discrete version of the combined state and continuity equation 51 is

where are the components of the velocity vector in two dimensions, the superscript indicates the spatial sample, the subscript indicates the temporal sample, and and are the respective discrete increments between samples. In two dimensions, the equation of motion 50 has two componentsFigure 9

I programmed the three equations, and Figure 9 shows snapshots for a point source and progressing time. The later frames show the computation to be instable, a standard problem with explicit finite difference schemes. I first tried to guess parameters and for which the computation would be stable. I mostly increased the viscosity to counteract the exponential amplitude growth by exponential amplitude reduction. When I failed to find appropriate values, I decided to estimate the stability limits of the parameters by careful analysis of the Courant-Levy criteria.

The eigenmodes of the discrete equation system are

Substitution of the right hand side in the system of difference equations (now written as a matrix operator) yields where To find the stability limit, I need to find parameters and so that all the eigenvalues of the propagation matrix are smaller than 1. This is equivalent with ensuring that the values of for which the determinant of the matrix is zero are all smaller than 1. The limit has to hold for all wave numbers3/9/1999