Stabilizing Upwinding of Transport Equation with Varying Velocity and Large Gradients

Question

I am attempting to solve the following transport equation using a 2D finite-difference scheme

$$c_t+\mathbf{v}(c)\cdot\nabla c=\frac{1}{Pe}\nabla^2 c,$$

where $c$ is modeling a local concentration, $t$ is time, and $\mathbf{v}$ is the local velocity which depends rather strongly on $c$. For the convection term, $\mathbf{v}\cdot\nabla c$, I write it in the expanded form

$$\mathbf{v}\cdot\nabla c=v_xc_x+v_zc_z$$

and upwind each term separately. By upwinding, I mean that if $v_x>0$, then I use a backwards finite difference to approximate $c_x$, and similarly for $c_z$. The issue that arises is there are sometimes very sharp gradients in $c$ in the $z$ direction that are difficult to properly resolve with a finite number of nodes (the number of nodes required is infeasible). Because of this, the forward and backward finite difference calculations may have large discrepancies. So, when $v_z$ switches sign and the differencing scheme for $c_z$ changes between cells, there is a large change in the value of $c_z$ across the cell border. This seems to be causing unstable oscillations to creep into the solution centered at the point that $v_z$ changes sign.

I have had a few ideas on how to address this issue, but I cannot find anything rigorous. One idea is to somehow interpolate across the cells where $v_z$ changes sign so the change in $c_z$ isn't so abrupt, but I'm not sure what interpolation scheme would be appropriate. I could also use a different discretization method in the $z$-direction (one that is not so sensitive to the large gradients), but again, I'm not sure what would be appropriate. I appreciate any suggestions you may have.

Answer 1

Under the assumption that your code is bug-free I suspect that you violate the CFL condition / stability requirement of the upwind scheme. In this lecture notes (chapter 2.4) it is shown that you need that $$ \max \{ v_x, v_z \} \frac{\Delta t}{\min \{\Delta x, \Delta z\}} \leq 1.$$ I suspect that this requirement is violated, thus you might need to decrease your timestep $\Delta x$ or decrease the spatial resolution $\Delta x, \Delta z.$

Apart from this, your code should damp oscillations because of

1. The diffusive term on the right smears out differences (For $Pe$ not too large).
2. The upwind scheme can also be reformulated (1D case for instruction purposes) as (see again the notes) $$ c_t + v \frac{c_{j+1} - c_j }{2 \Delta x} = \frac{ v}{2 \Delta x} \underbrace{\Big(c_{j+1} - 2 c_j + c_{j-1}\Big)}_{\approx \partial_{xx} c} $$ Thus, it also carries a diffusive term.