Advection-Diffusion Problems

• Consider the general PDE:

Or

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

• 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

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

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

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

• 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:

• 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$.

import ipywidgets as wgt
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


def f(Pe):
    ϕw = 2
    ϕe = 1
    ϕp = 1/2*(ϕw*(Pe/2 + 1) - ϕe*(Pe/2 - 1))
    ϕ = np.array([ϕw,ϕp,ϕe])
    x = np.array([0, 1, 2])
    plt.rc("font", size=14)
    plt.plot(x,ϕ,'o-', markersize=10, lw=3)
    plt.ylim([0,3])
    plt.xlabel('x')
    plt.ylabel('ϕ')
    plt.grid()
    plt.show() 

wgt.interact(f, Pe=(-5,5,0.1));     # try adding string and list arguments...

Following Patankar, we can write our above equation as

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.

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

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

Solve for $\phi_P$:

Or,

• 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.