or,
$$\frac{d\phi}{dx} = \frac{\Gamma}{u}\frac{d^2\phi}{dx^2}.$$import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
xL = np.linspace(0,1,1000)
def phihat(xL, Pe):
return (np.exp(Pe*xL)-1)/(np.exp(Pe)-1)
plt.rc('font', size=14)
plt.plot(xL, phihat(xL,-20), label=r'$P_e=-20$')
plt.plot(xL, phihat(xL,-5), label=r'$P_e=-5$')
plt.plot(xL, phihat(xL,1E-12), label=r'$P_e=0$')
plt.plot(xL, phihat(xL,5), label=r'$P_e=5$')
plt.plot(xL, phihat(xL,20), label=r'$P_e=20$')
plt.legend(frameon=False, bbox_to_anchor=(1.0, 0.8))
plt.xlabel('x/L')
plt.ylabel(r'$\hat{\phi}$');
For mixed problems, we can do better by using the exact solution to construct an exact finite difference method.
$$\frac{d\phi}{dx} = \frac{\Gamma}{u}\frac{d^2\phi}{dx^2} \rightarrow \frac{d}{dx}\underbrace{\left(u\phi - \Gamma\frac{d\phi}{dx}\right)}_{\mbox{Flux: F}} = 0.$$Hence, our equation is
$$\frac{dF}{dx} = 0.$$For finite volume we have
$$F_e - F_w = 0.$$$$F_e = u_e\phi_e - \Gamma_e\left(\frac{d\phi}{dx}\right)_e,$$$$F_w = u_w\phi_w - \Gamma_w\left(\frac{d\phi}{dx}\right)_w.$$Here, we are using $(x/\Delta x)_e=1/2$, that is Practice A, where the face is midway between grid points. The above expression simplifies nicely (by cancellation of the exponential terms) to
$$F_e = u_e(\phi_P - \gamma),$$Now, we insert $\gamma$, and repeat all of this for $F_w$ to obtain
$$F_e = u_e\left(\phi_P + \frac{\phi_P-\phi_E}{e^{P_e}-1}\right).$$ $$F_w = u_w\left(\phi_W + \frac{\phi_W-\phi_P}{e^{P_w}-1}\right).$$Finally, for $F_e-F_w=0$ we can write our scheme in terms of coefficients of $\phi_W$, $\phi_P$ and $\phi_E$:
$$ \phi_W\left[-u_w-\frac{u_w}{e^{P_w}-1}\right] + \phi_P\left[u_e+\frac{u_e}{e^{P_e}-1} + \frac{u_w}{e^{P_w}-1}\right] + \phi_E\left[\frac{-u_e}{e^{P_e}-1}\right] = 0. $$If we know the exact solution, why bother going through the effort to use it in a finite volume method using a grid? Why not just use the solution?
We can use this method to treat the combined advection and diffusion fluxes $F$ for more complicated situations that we don't know (or don't want to find) the exact solution for. For example:
Exponentials are computationally expensive to compute.
Replace the exponential terms in the above coefficients with polynomial approximations.
$$\frac{u}{e^P-1} = \frac{\Gamma}{\Delta x}\frac{P}{e^P-1}.$$
Consider the last exponential term, call it $T$:
P = np.linspace(-5,8,1001) # 1001 to avoid 0 in the list
T = P/(np.exp(P)-1)
Tfit = np.zeros(len(P))
Tfit[P<=2] = 1-0.5*P[P<=2]
Tfit[P<-2] = -P[P<-2]
fig, ax = plt.subplots()
ax.plot(P,T, 'k-', label='exp')
ax.plot(P,Tfit, 'r-', label='fit')
plt.ylim([-2,6])
plt.xlim([-6,9])
ax.axvline(x=0, color='grey', ls=':')
ax.axhline(y=0, color='grey', ls=':')
plt.ylabel(r'$\frac{P}{e^P-1}$', fontsize='18')
plt.xlabel('P')
plt.legend(frameon=False);