A good finite difference for the continuity equation

Question

What would be a good finite difference discretization for the following equation:

$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho u\right)=0$?

We can take the 1D case:

$\frac{\partial \rho}{\partial t} + \frac{d}{dx}\left(\rho u\right)=0$

For some reason all schemes I can find is for the formulation in Lagrangian coordinates. I came up with this scheme for the time being (disregard the j index):

$ \frac{\rho^{n+1}_{i,j}-\rho^{n}_{i,j}}{\tau} + \frac{1}{h_x}\left(\frac{\rho^{n+1}_{i+1,j}+\rho^{n+1}_{i,j}}{2}u^{n}_{_xi+1/2,j}- \frac{\rho^{n+1}_{i,j}+\rho^{n+1}_{i-1,j}}{2}u^{n}_{_xi-1/2}\right)=0$

But seems to be really unstable or have some horrible stability condition. Is that so?

The velocity is actually calculated through the darcy law $u=-\frac{k}{\mu}\nabla p$. Plus we have the equation of state. The full system consists also of an energy equation and the equation of state for the ideal gas. The velocities can turn negative.

Answer 1

You are looking at the mass conservation equation:

$\dfrac{dm}{dt}=0$

When considering mass evolution per unit volume, this boils down to the density advection equation in flux form:

$\dfrac{\partial \rho}{\partial t} = -\nabla \cdot (\rho u)$

Good thing about this is that it is just the advection equation of an arbitrary scalar field (in our case, this happens to be density $\rho$) and it is (relatively) easy to solve, provided adequate time and space differencing schemes, and initial and boundary conditions.

When designing a finite differencing scheme, we worry about convergence, stability and accuracy. A scheme is converging if $\dfrac{\Delta A}{\Delta t} \rightarrow \dfrac{\partial A}{\partial t}$ when $\Delta t \rightarrow 0$. Stability of the schemes ensures that the quantity $A$ remains finite when $t \rightarrow \infty$. Formal accuracy of the scheme tells where the truncation error in the Taylor expansion series of the partial derivative lies. Look into a CFD textbook for more details on these fundamental properties of a differencing scheme.

Now, the simplest approach is to go straight to 1st order upstream differencing. This scheme is positive-definite, conservative and computationally efficient. The first two properties are especially important when we model the evolution of a quantity which is always positive (i.e. mass or density).

For simplicity, let's look at 1-D case:

$\dfrac{\partial \rho}{\partial t} = -\dfrac{\partial(\rho u)}{\partial x}$

It is convenient now to define the flux $\Phi = \rho u$, so that:

$\dfrac{\partial(\rho u)}{\partial x} = \dfrac{\partial \Phi}{\partial x} \approx \dfrac{\Delta \Phi}{\Delta x} \approx \dfrac{\Phi_{i+1/2}-\Phi_{i-1/2}}{\Delta x}$

Here's a schematic of what we are simulating:

            u           u
|          -->         -->          |
|    rho    |    rho    |    rho    |
x-----o-----x-----o-----x-----o-----x
     i-1  i-1/2   i   i+1/2  i+1

We are evaluating the evolution of $\rho$ at cell $i$. The net gain or loss comes from the difference of what comes in, $\Phi_{i-1/2}$ and what goes out, $\Phi_{i+1/2}$. This is where we start to diverge from Paul's answer. In true conservative upstream differencing, the quantity at the cell center is being carried by velocity at its cell edge, in the direction of its motion. In other words, if you imagine you are the advected quantity and you are sitting at the cell center, you are being carried into the cell in front of you by the velocity at the cell edge. Evaluating the flux at the cell edge as a product of density and velocity, both at the cell edge, is not correct and does not conserve the advected quantity.

Incoming and outgoing fluxes are evaluated as:

$\Phi_{i+1/2} = \dfrac{u_{i+1/2}+|u_{i+1/2}|}{2}\rho_{i} + \dfrac{u_{i+1/2}-|u_{i+1/2}|}{2}\rho_{i+1}$

$\Phi_{i-1/2} = \dfrac{u_{i-1/2}+|u_{i-1/2}|}{2}\rho_{i-1} + \dfrac{u_{i-1/2}-|u_{i-1/2}|}{2}\rho_{i}$

The above treatment of flux differencing ensures upstream-definiteness. In other words, it adjusts the differencing direction according to the sign of velocity.

The Courant-Friedrichs-Lewy (CFL) stability criterion, when doing time differencing with simple first order, forward Euler differencing is given as:

$\mu = \dfrac{u\Delta t}{\Delta x} \leq 1$

Note that in 2 dimensions, the CFL stability criterion is more strict:

$\mu = \dfrac{c\Delta t}{\Delta x} \leq \dfrac{1}{\sqrt{2}}$

where $c$ is velocity magnitude, $\sqrt{u^{2}+v^{2}}$.

Some things to consider. This scheme may or may not be appropriate for your application depending on what kind of process you are simulating. This scheme is highly diffusive, and is appropriate for very smooth flows without sharp gradients. It is also more diffusive for shorter time steps. In the 1-D case, you will obtain an almost exact solution if the gradients are very small, and if $\mu = 1$. In the 2-D case, this is not possible, and diffusion is anisotropic.

If your physical system considers shock waves or high gradients of other sort, you should look into upstream differencing of higher order (e.g. 3rd or 5th order). Also, it may be worthwhile looking into the Flux Corrected Transport family of schemes (Zalesak, 1979, JCP); anti-diffusion correction for the above scheme by Smolarkiewicz (1984, JCP); MPDATA family of schemes by Smolarkiewicz (1998, JCP).

For time differencing, 1st order forward Euler differencing may be satisfactory for your needs. Otherwise, look into higher-order methods such as Runge-Kutta (iterative), or Adams-Bashforth and Adams-Moulton (multi-level).

It would be worthwhile looking into some CFD graduate-level textbook for a summary of above mentioned schemes and many more.

Answer 2

In the 1D case, you don't want to use a forward or central difference scheme for the spatial derivative term $(\frac{d}{dx})$ because they are numerically unstable. Instead, it is better to discretize the equation with an explicit backwards (upwind) finite difference for the spatial derivative:

$\frac{\rho_i^{k+1}-\rho_i^k}{\Delta t} + \frac{\rho_i^k U_i^k-\rho_{i-1}^k U_{i-1}^k}{\Delta x} = 0$.

If the velocities are positive, then this backward scheme is stable. If they are negative, then a forward difference will work. Regardless, there is always a constraint on your choice of $\Delta x$ and $\Delta t$ (courant number) to make the scheme stable.