7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You can easily check whether a scheme is conservative by summing over the grid to see if all the fluxes cancel. $\endgroup$ – David Ketcheson May 21 '13 at 14:54
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you for the comments. I think I have fixed the matrix equation, all terms should have the correct sign now. For boundary conditions I was going to specify zero-flux type (Neumann conditions) i.e. $\frac{\partial u}{\partial x} = 0$ at the left and right boundary of the system. I would implement this in the usual way; by introducing ghost points etc. Does this all seem reasonable? $\endgroup$ – boyfarrell May 20 '13 at 4:46
  • $\begingroup$ @boyfarrell No, that is not a valid boundary condition for advection. I suggest checking out a book like LeVeque. $\endgroup$ – Jed Brown May 20 '13 at 12:35
  • $\begingroup$ @JedBrown I see. I think I need to specify an inflow condition on one end (according the LeVeque). I think this is probably a question in it's own right. $\endgroup$ – boyfarrell May 20 '13 at 14:17
  • $\begingroup$ @boyfarrell Yes, you need a valid inflow condition and you don't get to impose an outflow boundary condition (unless you add diffusion). $\endgroup$ – Jed Brown May 20 '13 at 14:23
  • $\begingroup$ Do you mean if the flux has a diffusion and an advection component then, then it's proper to do the following? scicomp.stackexchange.com/questions/5434/… $\endgroup$ – boyfarrell May 20 '13 at 14:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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