Numerical instabilities in the viscous wave equation

I'm currently working out the wave equation with the introduction a viscosity term, s.t. $a(x)\ddot{u}(x, t) = \nabla (b(x)\nabla . u(x, t) + c(x) \nabla . \dot{u}(x, t))$ (where $u$ is the displacement field of a particle).

I use a spectral method for space discretization (such that the mass matrix is diagonal) and high order Runge-Kutta for time marching. I obtain very poor results, and some numerical instabilities appear depending on $a, b, c$. The wave amplitude blows up at some point.

I studied the behavior of the equation and, numerically, it shows that there's a strong bind between the value of the ratio $\frac{c}{b}(x)$ and the stability of the resolution. We end up with a upper bound for this ratio such that when overshooting it, it seems impossible to solve the equation ($\frac{c}{b} < D$ is required, $D \in \mathbb{R}$).

Obviously, this constant $D$ is dependent on the both space and time discretizations (scheme + time,space steps $\triangle t, \triangle x$). It's really likely that a CFL condition is involved, as $v(b, c)\frac{\triangle t}{\triangle x} < E$ where $v$ is the wave velocity depending on both coefficients $b, c$. Indeed, increasing $\triangle x$ makes the resolution feasible for larger values of $\frac{c}{b}$.

For some specific boundary conditions (periodic boundary conditions + dirichlet + free boundary), I can get an analytic solution and when $b, c$ are constants, it seems that the resolution works well and my numerical solution fits pretty well. Also, without the viscous term, the code works perfectly fine.

As an additional info, the phase velocity, tends to infinity as frequency increases. By the way, I asked myself if you choose a sufficiently large $\triangle t$, such that $1/\triangle t << f_{max}$ where $f_{max}$ is a chosen frequency, do you get, in your modelization, the "action" of frequencies beyond that $f_{max}$ ? (kind of Nyquist consideration)

I know that it's a very general consideration and you will probably ask for more details but I was wondering, if someone had an idea where, theoretically and numerically things could get wrong ?

Thanks

• $\nabla . \dot{u}(x, t)$ seems to be velocity gradient to me. Velocity should be $\dot{u}(x, t)$. – mike Jun 13 '14 at 8:38
• This is the divergence of the velocity field (velocity in time not in space). – Amzocks Jun 13 '14 at 8:45

• I quickly read the article but I don't see how, even though you have no causality (meaning that it depends on the history) for a solution to the viscous wave equation with constant parameters (section $6$), it can play any role in my case. – Amzocks Jun 13 '14 at 13:40