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$.

  1. 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}$

  1. 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}$.

  1. 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 ?


1 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.

  • $\begingroup$ I still have a question on this. I use the "conservative variables". In the case of the diffusion, I would get $$ \frac{\partial}{\partial t}\left(\frac{3}{2}nT+\frac{1}{2}n v^2\right)=\frac{\partial}{\partial x}\left(D(T)\frac{\partial T}{\partial x} \right)$$. This requires a non-linear solver. This could be more convenient if rewritten as $$ n^* \frac{\partial T}{\partial t}=\frac{\partial}{\partial x}\left(D(T^*)\frac{\partial T}{\partial x}\right)$$ This would give a tridiagonal matrix, easier to solve. Is it possible to do it ? $\endgroup$ Sep 15, 2018 at 21:37
  • $\begingroup$ $n$ does not change in the diffusion step so thats fine. But using $D(T^*)$ would give you only first order accuracy in time. $\endgroup$
    – cfdlab
    Sep 16, 2018 at 4:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.