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

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

interactive(children=(FloatSlider(value=0.0, description='Pe', max=5.0, min=-5.0), Output()), _dom_classes=('w…

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.