Conservative finite-difference expression for the advection equation

Question

Following on from the earlier question I am trying to derive a finite-difference scheme for the advection equation which is conservative. It was suggested that for advection equation with variable velocity,

$$ \frac{\partial u}{\partial x} + \frac{\partial}{\partial x}\left( \boldsymbol{v}(x)u\right) = 0 $$ that the chain rule should not be applied if conservation is required because the flux of $u$ is the fundamental quantity.

The product $\boldsymbol{v}(x)⋅∇u$ is not conservative transport if $\boldsymbol{v}(x)$ is not divergence-free (i.e., constant in 1D). You should stick with the conservative form and resist the urge to apply the chain rule when evaluating conservation properties.

Question

Is the following a reasonable approach for discretization?

• I have kept the flux together (i.e the product $\boldsymbol{v}_ju_j$) while applying the discretization.
• However to write as a matrix equation where $u_j$ is the unknown they eventually must be separated.
• I have a slight concern that this is incorrect because the only difference between the constant velocity case and the varying velocity case is that that $\boldsymbol{v}$ becomes a vector.
• If we applied the chain rule we could get exactly the same equation but with the additional $u \frac{\partial\boldsymbol{v(x)}}{\partial x}$ term which is treated as a source term.

For example, if we apply Crank-Nicolson scheme to the advection equation,

$$ \frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} + \left[ \frac{1-\beta}{2\Delta x} \left( \boldsymbol{v}_{j+1}^{n} u_{j+1}^{n} - \boldsymbol{v}_{j-1}^{n} u_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \boldsymbol{v}_{j+1}^{n+1} u_{j+1}^{n+1} - \boldsymbol{v}_{j-1}^{n+1} u_{j-1}^{n+1} \right) \right] = 0 $$

Using the following substitution, $$ \boldsymbol{r}_j^n = \frac{\boldsymbol{v}_j^n}{2}\frac{\Delta t}{\Delta x} $$ and move the $n+1$ terms to the l.h.s gives,

$$ u_{j}^{n+1} + \beta\boldsymbol{r}_{j+1}^{n+1} u_{j+1}^{n+1} - \beta\boldsymbol{r}_{j-1}^{n+1} u_{j-1}^{n+1} = u_{j}^{n} - (1- \beta)\boldsymbol{r}_{j+1}^{n} u_{j+1}^{n} + (1- \beta)\boldsymbol{r}_{j-1}^{n} u_{j-1}^{n} $$

We can write this as the matrix equation, $$ \begin{pmatrix} 1 & \beta r_2^{n+1} & & 0 \\ -\beta r_1^{n+1} & 1 & \beta r_3^{n+1} & \\ & \ddots & \ddots & \ddots \\ & -\beta r_{J-2}^{n+1} & 1 & \beta r_J^{n+1} \\ 0 & & -\beta r_{J-1}^{n+1} & 1 \\ \end{pmatrix} \begin{pmatrix} u_{1}^{n+1} \\ u_{2}^{n+1} \\ \vdots \\ u_{J-1}^{n+1} \\ u_{J}^{n+1} \end{pmatrix} = \begin{pmatrix} 1 & -(1 - \beta)r_2^n & & 0 \\ (1 - \beta) r_1^n & 1 & -(1 - \beta)r_3^n & \\ & \ddots & \ddots & \ddots \\ & (1 - \beta) r_{J-2}^n & 1 & -(1 - \beta)r_J^n \\ 0 & & (1 - \beta) r_{J-1}^n & 1 \\ \end{pmatrix} \begin{pmatrix} u_{1}^{n} \\ u_{2}^{n} \\ \vdots \\ u_{J-1}^{n} \\ u_{J}^{n} \end{pmatrix} $$

---

Update I Corrected a sign error suggested by @Geoff-Oxberry.

Answer 1

In 1-D, the first discretization you've presented is correct. The matrix equation does not look right, though. For starters, it's not clear what your boundary conditions would be or how you would incorporate them.

In $N$ dimensions, the advection equation looks like

\begin{align} \frac{\partial{u}}{\partial{t}} + \sum_{i=1}^{N}\frac{\partial}{\partial{x^{i}}}(v^{i} u) &= 0, \end{align}

adopting the convention that superscripts refer to components of vectors; you would discretize this equation in similar fashion to the 1-D case.

If $v$ (or in $N$-D, $\mathbf{v}$) is a known function of $x$ (respectively, $\mathbf{x}$), then the values $v_{j}$ (resp., $\mathbf{v}_{j}$) become coefficients for the various $u_{j}$ terms in the discrete equation.