Advection-Diffusion Problems¶

• Consider the general PDE:

$$\phi_t + u\phi_x = \Gamma\phi_{xx},$$

Or

$$\frac{\partial\phi}{\partial t} + u\frac{\partial\phi}{\partial x} = \Gamma \frac{\partial^2\phi}{\partial x^2},$$

where $u$ is a velocity, $\phi$ is our scalar variable, and $\Gamma$ is a diffusion coefficient (like $\alpha$).

• The second term in the equation is a advection term. This is a hyperbolic-like term, a wave-like term.
• The equation as a whole is parabolic.

• Example. Let $\phi=\rho Y_i$, $\Gamma=D_i$ with constant $\rho$, $D_i$: $$\frac{\partial\rho Y_i}{\partial t} + \vec{v}\cdot\nabla(\rho Y_i) = \rho D_i\nabla^2Y_i.$$
  • Unsteady species transport equation.

FTCS method¶

$$\frac{\partial\phi}{\partial t} + u\frac{\partial\phi}{\partial x} = \Gamma \frac{\partial^2\phi}{\partial x^2},$$

$$\phi_i^{n+1} = \phi_i^n - \frac{1}{2}\underbrace{\frac{u\Delta t}{\Delta x}}_{c}(\phi_{i+1}^n - \phi_{i-1}^n) + \underbrace{\frac{\Gamma\Delta t}{\Delta x^2}}_{d}(\phi_{i-1}^n - 2\phi_i^n + \phi_{i+1}^n).$$

• FTCS properties:
  • consistent,
  • conditionally stable,
  • $\mathcal{O}(\Delta t)$,
  • $\mathcal{O}(\Delta x^2)$.

• A stability analysis yeilds the following constraints: $$d\le \frac{1}{2},$$ $$c^2 \le 2d,$$ Or, $$c^2\le2d\le 1.$$

• For $d=1/2$, $c<1$.

• Both $c$ and $d$ are time ratios.

  • $d$ is the ratio of our timestep size to the grid-scale diffusion time.
    • $\tau_D = \Delta x^2/\Gamma.$
    • This is the nominal time it takes to diffuse across a cell.
  • $c$ is the ratio of our timestep size to the grid-scale advection time.
    • $\tau_c = \Delta x/u.$
    • This is the nominal time it takes to convect across a cell.
  • Don't take a longer timestep than the physical time of advection or diffusion.

• Note that $c$ and $d$ are related.

BTCS method¶

• BTCS properties:
  • consistent,
  • unconditionally stable,
  • $\mathcal{O}(\Delta t)$,
  • $\mathcal{O}(\Delta x^2)$.

Trouble with central differences¶

• Consider the steady state equation

$$u\frac{\partial\phi}{\partial x} = \Gamma \frac{\partial^2\phi}{\partial x^2}.$$

• Apply a finite volume solution by integrating over a control volume:

$$(u\phi)_e - (u\phi)_w = \Gamma\left(\frac{d\phi}{dx}\right)_e - \Gamma\left(\frac{d\phi}{dx}\right)_w.$$

$$\frac{1}{2}u_e(\phi_E + \phi_P) - \frac{1}{2}u_w(\phi_P + \phi_W) = \frac{\Gamma}{\Delta x}(\phi_E - \phi_P) - \frac{\Gamma}{\Delta x}(\phi_P - \phi_W).$$

• For constant $\rho$, $u$ is constant. Also assume $\Gamma$ is constant. Rearrange the expression:

$$\phi_W\left(-\frac{1}{2}u - \frac{\Gamma}{\Delta x}\right) + \phi_P\left(\frac{\Gamma}{\Delta x}\right) + \phi_E\left(\frac{1}{2}u - \frac{\Gamma}{\Delta x}\right) = 0.$$

Divide through by $\frac{\Gamma}{\Delta x}$:

$$\phi_W\left(-\frac{Pe}{2} -1 \right) + \phi_P(2) + \phi_E\left(\frac{Pe}{2} - 1\right) = 0.$$

• Here, $Pe=u\Delta x/\Gamma$ is a grid Peclet number. If $\Gamma=\nu$, the kinematic viscosity, this would be a grid Reynolds number.
  • $Pe = c/d = \tau_D/\tau_c$.
  • $Pe$ is a advection rate over a diffusion rate.

Solve the above equation for $\phi_P$ in terms of its neighbors:

$$\phi_P = \frac{1}{2}\left[\phi_W\left(\frac{Pe}{2} + 1\right) - \phi_E\left(\frac{Pe}{2}-1\right)\right].$$

• Suppose we have a decreasing profile where $\phi_W>\phi_E$.
  • For a steady state problem without a source term, we expect $\phi_P$ to be between its neighbors $\phi_E$ and $\phi_W$.
  • But if our velocity and grid spacing are set so that $Pe > 2$, then the sign of the $\phi_E$ term in the above equations changes and $\phi_P > \phi_W$, contrary to our expectation.
  • For $Pe>2$, we get unphysical behavior.
  • For given $u$ and $\Gamma$, $Pe$ can be reduced by decreasing $\Delta x$.

Following Patankar, we can write our above equation as

$$C_p \phi_P = C_E\phi_E + C_W\phi_W,$$

where $C$'s are coefficients of the $\phi$'s.

• Here, $C_P=1$, $C_E = -Pe/4 + 1/2$, and $C_W = Pe/4 + 1/2.$
• All coefficients in this form should have the same sign so that increases in the neighbors result in increases in $\phi_P$, and vice-versa.
• If $Pe>2$, then $C_E$ is negitive while the other $C$'s are positive.

Remedy: upwinding¶

The problem above is that we applied a central difference approximation when evaluating the convective term.

• (In the FV approach, we applied a linear interpolation when evaluating the face properties, which is analogous to using a central difference.)
• These approaches effectively give equal weighting to the neighbors when computing a property.
• But advection is directional.
  • You can smell much better when something is upwind of you then when it is downwind of you.
• In evaluating face properties, don't weight each neighbor equally.
• Insteady, advection blows the upstream value to the face.

$$\phi_e = \phi_P,\phantom{xx}\mbox{if}\phantom{xx} u_e>0,$$ $$\phi_e = \phi_E,\phantom{xx}\mbox{if}\phantom{xx} u_e<0,$$

$$\phi_w = \phi_W,\phantom{xx}\mbox{if}\phantom{xx} u_w>0,$$ $$\phi_w = \phi_P,\phantom{xx}\mbox{if}\phantom{xx} u_w<0.$$

There is no change to the diffusion term. Diffusion has no directionality. Use central differences.

• Now, take $u>0$ and solve our PDE:

$$u\phi_e - u\phi_w = \Gamma\left(\frac{d\phi}{dx}\right)_e - \Gamma\left(\frac{d\phi}{dx}\right)_w.$$

$$u\phi_P - u\phi_W = \frac{\Gamma}{\Delta x}(\phi_E - \phi_P) - \frac{\Gamma}{\Delta x}(\phi_P - \phi_W).$$

Solve for $\phi_P$:

$$\phi_P = \phi_E\left[\frac{\frac{\Gamma}{\Delta x}}{u + 2\frac{\Gamma}{\Delta x}}\right] + \phi_W\left[\frac{u+\frac{\Gamma}{\Delta x}}{u + 2\frac{\Gamma}{\Delta x}}\right], $$

Or,

$$\phi_P = \phi_E\left[\frac{1}{Pe+2}\right] + \phi_W\left[\frac{Pe+1}{Pe+2}\right]. $$

• Here, the coefficients are always positive for any value of $Pe$.

• Note, this method is $\mathcal{O}(\Delta t)$, but only $\mathcal{O}(\Delta x)$. The upwinding is first order spatially.