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:
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 pressure
Instead 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
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 components
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
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 numbers kx and ky. The determinant is a third order polynomial in and the estimation of the upper limit requires very careful accounting of the sign of all the exponential expressions. When I first investigated the problem, I was stumped. But now I realize that I can find the boundaries by computing the zeroes of the corresponding derivatives and by checking the particular boundary values. The computation is tiresome but feasible. I, however, had already moved on and implemented the simpler acoustic wave equation.