Navier-Stokes solver: How to adjust the time step based on non-linear terms?

Question

My code solves the incompressible Navier-Stokes equation in a conducting fluid, together with the induction equation:

$ \partial_t u + u \nabla u + 2\Omega \times u = -\nabla p + \nu \Delta u + (\nabla \times b) \times b \\ \partial_t b = \nabla \times (u \times b) + \eta \Delta b$

where $u$ and $b$ are the velocity and magnetic fields, $\eta$ and $\nu$ are diffusion constants, $p$ is the pressure, $\Omega$ the rotation rate of the reference frame. I solve this PDE system in spherical geometry using finite differences in the radial direction and spherical harmonic expansion.

Let's focus on the Navier-Stokes equation which contains the Lorentz force $(\nabla \times b) \times b$. I use a Crank-Nicholson semi-implicit scheme for the diffusive terms, while the non-linear terms $u\nabla u$ and $(\nabla \times b)\times b$ are treated explicitely:

$Au_{n+1} = Bu_{n} + NL(u_n,b_n)$

where $u_n$ and $b_n$ are respectively my velocity and magnetic field at time-step $n$, and $A$ and $B$ are matrices containing the time-step. (Note that $NL(u_n,b_n)$ actually also includes $2\Omega\times u$ for convenience, and that $\nabla p$ is eliminated by taking the curl)

I would like to adjust dynamically my time-step. For Navier-Stokes, the usual CFL criterion can be used, but what about the Lorentz force ? I can proably come up with some specific criterion, but here is my question:

Is there a way to adjust the time-step $\Delta t$ knowing only the non-linear term $NL(u_n,b_n)$ and ignoring its analytical form or physical meaning ? Some sort of generalised CFL criterion ?

I've thought about something like $NL \Delta t < c\, U$ (where $c$ is a constant). When $NL = u\nabla u \sim \frac{U^2}{\Delta x}$, it reduces to the CFL $\frac{U\Delta t}{\Delta x} < c$. Doest it makes sense ? Is there a better way to think about this ? Are there references on the subject ?

Answer 1

It would be helpful if you would write down the full system of PDEs you are considering, since (to me, at least) there is still some ambiguity. It would also be helpful if you would describe the discretization more precisely. But I will provide a rough answer based on the information given.

First we should distinguish between two concepts: (1) the strict CFL condition that is a necessary condition for convergence; and (2) the restriction on the CFL number that is both necessary and sufficient for stability. I think you are interested in the latter, which is the one with more practical significance.

Your thinking is correct -- for a typical explicit discretization of any kind of transport PDE, the restriction on the CFL number that ensures stability will have the form $$\Delta t \le c \frac{\Delta x}{U_{max}}$$ where $U_{max}$ is the largest characteristic speed and $c$ is a constant depending on the numerical discretization. The reason the stability condition scales this way is that the eigenvalues of the spatial discretization scale this way. For the convective term you mention, $U_{max}=2\max(U)$.

In general, determining $U_{max}$ can be complicated if the solution is not smooth or the flux is non-convex. But for smooth solutions and convex flux, $U_{max}$ is simply the largest eigenvalue of the jacobian of the flux function. To learn more, get a good book on hyperbolic PDEs. I would recommend either the second or the third book on this page. For a very simple introduction, look at Chapter 10 of the first book on that page.