Multi-steps method for Navier-stokes equations with strongly nonlinear diffusion

Question

I am trying to solve a particular form of the Euler / Navier-Stokes equations in 1D, with very strong and non-linear diffusion coefficients.

My system of equations is

\begin{cases} \frac{\partial n}{\partial t} + \frac{\partial}{\partial x}\left(nv \right ) = S_n\\ \frac{\partial \left( n v\right)}{\partial t} + \frac{\partial}{\partial x}\left(nv^2 + nT \right ) = S_M + \frac{\partial}{\partial x}\left(\nu(T)\frac{\partial v}{\partial x} \right)\\ \frac{\partial}{\partial t}\left(\frac{3}{2}nT + \frac{1}{2}n v^2 \right ) + \frac{\partial}{\partial x}\left(\frac{5}{2}nTv + \frac{1}{2}n v^3 \right ) - \frac{\partial}{\partial x}\left(D(T)\frac{\partial T}{\partial x} \right) - \frac{\partial}{\partial x}\left(v\nu(T)\frac{\partial v}{\partial x} \right) = S_E \end{cases} where $\nu(T)$ and $D(T)$ are strongly nonlinear coefficients. For instance, in some versions of my model, $\nu(T), D(T) \propto T^3$. This system is associated with a set of boundary conditions and initial conditions that I think are not essential to describe my issue. The temperature diffusion term ($D$) is by far (at least one order of magnitude) the fastest timescale in my system.

I've tried solving this using a purely time-explicit scheme, but this leads to very small time-steps (if the diffusion terms are removed from the system, I can use higher time-steps, so I guess they are the culprits). I also tried a fully implicit method (Crank-Nicholson or Backward-Euler, implemented with Petsc), but the non-linear Newton solve-step is quite slow. I would like to improve this.

I wondered if a scheme like the following could be more efficient:

1. Solve the advection part of the equation

using $Q = \begin{pmatrix} n \\ nv \\ \frac{3}{2}nT + \frac{1}{2}nv^2 \end{pmatrix}$ and $F(Q) = \begin{pmatrix} nv \\ nv^2+nT \\ \frac{5}{2}nTv + \frac{1}{2}nv^3 \end{pmatrix}$,

we solve $\frac{\partial Q}{\partial t} + \frac{\partial F}{\partial x} = S$ (where $S$ is a vector of the source terms). This is a quite "classical" problem, that can be solved explicitly (I've tried a FD-WENO3 that I was happy with). This gives us $Q^\star$, from which we can compute $n^\star$, $v^\star$ and $T^\star$.

2. We "push" $v$ implicitly, assuming other variables are frozen. We write

$n^\star \frac{\partial v}{\partial t} = \frac{\partial}{\partial x}\left(\nu(T^\star)\frac{\partial v}{\partial x} \right)$.

This is linear in $v$, so can be solved very effectively, and get $v^{\star\star}$

3. We "push" $T$ implicitly, assuming other variables are frozen. $n^\star \frac{\partial T}{\partial t} = \frac{\partial}{\partial x}\left(D(T^\star)\frac{\partial T}{\partial x} \right) + \frac{\partial}{\partial x}\left(v^{\star\star}\nu(T^\star)\frac{\partial v^{\star\star}}{\partial x} \right) $.

which again should not be complicated to solve. We get $T^{\star\star}$.

4. Solve again the explicit part using $n^\star, v^{\star\star}, T^{\star\star}$ as initial solutions.

This is then the end of the iteration, and we have as a solution at $t=t+\Delta t$, using for steps 1, 2, 3, 4 the time-steps $\Delta t /2 $, $\Delta t$,$\Delta t$,$\Delta t /2$, respectively.

As you can guess, I am not an expert in numerical methods, and I have a hard time determining if this idea makes any sense. Is it consistent? I have the feeling that it is plain wrong. If it makes sense, is the order of the steps correct ? Otherwise, would you have any ideas of methods I could explore to solve this problem ?

Answer 1

If your diffusion is strong, then it is a strongly elliptic problem and I would have thought this should be easier to solve.

You can try operator splitting approach which is basically what you are trying to describe. To get second order accuracy, you can use Strang splitting. If you search for "strang splitting" on this forum you will find many posts. The idea is simple.

$$ u_t = A(u) + B(u) $$

is solved in three steps. Let $E(A,\Delta t)$ and $E(B,\Delta t)$ denote the second order accurate numerical schemes to advance the split equations $u_t=A(u)$ and $u_t = B(u)$ by time step $\Delta t$.

Let $u^n$ be current solution.

Solve $u_t = A(u)$ for half time step

$$ u_1 = E(A, \Delta t/2) u^n $$

Solve $u_t = B(u)$ for full time step

$$ u_2 = E(B, \Delta t) u_1 $$

Solve $u_t = A(u)$ for half time step

$$ u^{n+1} = E(A, \Delta t/2) u_2 $$

where you now have the solution at next time step.