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
The problem is most likely that due to the added first order time derivative term (viscosity), Your numerical solution is not causal anymore. This is a common problem. Make sure the Kramers-Kronig relations are satisfied, and that will fix the discrepancies. Basically the Hilbert transform of the the 'attenuation' is related to the frequency and the phase velocity in a certain way. Make sure your coefficients obey this so your solution is physical (causal).