Stability analysid with pde, analogous fourrier number

Question

I am working with a team to discretize a PDE. The PDE has the following form:

$$ \frac{\partial ^2 u}{ \partial R^2} + \left( \frac{2}{R} + \frac{1}{k} \frac{\partial k}{\partial R} \right) \frac{\partial u}{\partial R} = \frac{1}{k} \frac{\partial u}{\partial t} $$

We have run into stability issues. I'm trying to figure out an expression
for something like the Fourier number, i.e. I need the equivelant expression for: $$ Fo = \alpha \frac{\Delta t}{\Delta x^2} $$ for the equation: $$ \frac{\partial u}{\partial t} = \alpha \frac{\partial ^2 u}{\partial x^2} $$

I am not looking for the solution per se, but documentation or advice on the matter would help greatly.

Answer 1

See Von Neumann Stability Analysis. The derivation is shown in Wikipedia and is

$$\qquad r=\frac{\alpha \Delta t}{\left( \Delta x \right)^2} \leq \frac{1}{2} $$

It is also covered in Chapters 9 and 10 of Finite Difference Methods for Ordinary and Partial Differential Equations by Randall LeVeque.

Von Neumann Stability Analysis is a standard technique that you (or someone on your team) should invest time in learning.

Answer 2

You can sort of pattern match the condition you gave for this equation. For the numerical diffusion term to be stable, we must have $$\max{k}\cdot\frac{\Delta t}{\Delta R^2}<\frac{1}{2}.$$ There is a similar condition for numerical advection, i.e. the first derivative term. It is sometimes called the CFL condition. For this problem, it is $$\max{\left(\frac{2}{R}+\frac{1}{k}\frac{\partial k}{\partial R}\right)}\cdot\frac{\Delta t}{\Delta R}<1.$$

If both these conditions are satisfied for $R$ in the computational domain, you should be fine. However, notice that in the second condition, the factor multiplying the deltas becomes infinitely large if $R$ approaches $0$ or if $k$ approaches $0$. This will cause instability and may be what's happening with your problem.

Some of this might be avoided if you could reformulate your problem into conservation form, i.e. not expanding the spatial derivatives and keeping it as something like $\frac{1}{R}(Rf(R)u_R)_R$

Note that these results only apply for explicit schemes as I have written them