def diffusion_flux_east(D=D_e, c_right=c_E, delta=delta):
""" Compute the diffusion flux to the east cell. """
return docs.diffusion_flux(D, c_P, c_right, delta=delta)
def diffusion_flux_west(D=D_w, c_left=c_W, delta=delta):
""" Compute the diffusion flux to the west cell. """
return docs.diffusion_flux(D, c_left, c_P, delta=delta)2 System discretization
NoteAbout the code
Some code elements are intentionally hidden from the rendered version. For all the details, please check the source code in the GitHub repository and execute the notes by yourself.
To solve the diffusion equation numerically, we employ the finite volume method (FVM). The first step to establish a FVM scheme consists of integrating the governing equation over a control volume \(P\) from \(w\) (west face) to \(e\) (east face) and over a time interval from \(0\) to \(\tau\). The left-hand side is first integrated over time as it represents a storage term, meaning that it measures how much the average concentration changes in the control volume over the time-step. The right-hand side is a flux term, and using similar arguments its natural integration is performed over space. Also, the order of integration can be exchanged due to the smoothness assumptions considered here. For more details, please check Patankar (1980).
\[ \int_{w}^{e}\int_{0}^{\tau} \dfrac{\partial{}c}{\partial{}t}\:dt\:dx= \int_{0}^{\tau}\int_{w}^{e} \dfrac{\partial{}}{\partial{}x} \left(D\dfrac{\partial{}c}{\partial{}x}\right)\:dx\:dt \]
This double integration allows us to convert the partial differential equation into a form suitable for numerical discretization. Evaluating the innermost integrals leads to:
\[ \int_{w}^{e} \left(c_{P}-c_{P}^{0}\right)\:dx= \int_{0}^{\tau} \left[ \left(D\dfrac{\partial{}c}{\partial{}x}\right)_{e}- \left(D\dfrac{\partial{}c}{\partial{}x}\right)_{w} \right] \:dt \]
Here, \(c_P^0\) and \(c_P\) represent the concentration at the center of the control volume \(P\) at times \(0\) and \(\tau\), respectively. A convention of not adding a superscript to time \(\tau\) is adopted here, unless required by the context. The right-hand side represents the net flux into the control volume through its east and west faces. Assuming the concentration is uniform within each control volume and the fluxes are constant over the time interval, we obtain:
\[ \left(c_{P}-c_{P}^{0}\right) \dfrac{\delta}{\tau}= \left(D\dfrac{\partial{}c}{\partial{}x}\right)_{e}- \left(D\dfrac{\partial{}c}{\partial{}x}\right)_{w} \]
where \(\delta\) is the width of the control volume. This equation represents a discrete time-step of the diffusion equation for a single control volume, where one still needs to specify a discretization of the gradients, i.e. how the fluxes are evaluated at a cell face. A simple introductory approach is to use forward differences, where the composition gradients are approximated as first differences between neighbors.
To approximate the aforementioned gradients, the next cell introduces the expressions for the fluxes at east and west faces. These helper functions ensure the order of arguments is consistent with the finite difference definition given above.
For standardizing equations as organized as \(F(x)=0\) we override sympy.Eq with Eq class used in what follows. This allows us to easily extract coefficients and tabulate them in a consistent way across the notes. Below we compose the diffusion equation with help of SymPy and the above extensions. Notice that when compositing the right-hand side rhs we multiply by -1 because we are integrating the negative of the fluxes.
J_w = diffusion_flux_west()
J_e = diffusion_flux_east()
lhs = (c_P - c_0) * delta / tau
rhs = -1 * (J_e - J_w)
eq = docs.Eq(lhs - rhs)
eq\(\displaystyle - \frac{D_{e} c_{E}}{\delta} + \frac{D_{e} c_{P}}{\delta} + \frac{D_{w} c_{P}}{\delta} - \frac{D_{w} c_{W}}{\delta} - \frac{c_{0} \delta}{\tau} + \frac{c_{P} \delta}{\tau} = 0\)
Here, \(c_E\) and \(c_W\) are the concentrations at the eastern and western neighboring control volumes, respectively, and \(D_e\) and \(D_w\) are the effective diffusion coefficients at the east and west faces of the control volume. Patankar (1980) shows that the evaluation of these diffusion coefficients is not arbitrary and for conservation (here the continuity of flux across the interface) to be respected one must take the harmonic mean of neighboring cells:
\[ D_{k}=2\frac{D_{i}D_{j}}{D_{i}+D_{j}} \]
A final step to get a linear form of the problem is multiplying both sides by \(\tau/\delta\) and rearranging terms. This form clearly shows the contribution of each neighboring concentration to the rate of change at point P. The coefficient of \(c_P\) is negative, representing the outflow from the central node, while the coefficients of \(c_E\) and \(c_W\) are positive, representing inflow from the neighboring nodes.
\(\displaystyle - \frac{D_{e} c_{E} \tau}{\delta^{2}} + \frac{D_{e} c_{P} \tau}{\delta^{2}} + \frac{D_{w} c_{P} \tau}{\delta^{2}} - \frac{D_{w} c_{W} \tau}{\delta^{2}} - c_{0} + c_{P} = 0\)
Introducing dimensionless coefficients \(Fo_{k}\), the Fourier number for mass transfer (more generally referred to as \(\mathcal{Fo}_{m}\)), we obtain the final discrete form. It represents the dimensionless diffusion number, which determines the stability and accuracy of the numerical scheme. The subscript \(k\) refers to either the east (\(e\)) or west (\(w\)) face.
\[ Fo_{k}=\dfrac{\tau{}D_{k}}{\delta^2} \]
This definition is provided below for helping with the replacements that follow.
SUBS_FOURIER = [
(docs.fourier_number(D_e, tau, delta), Fo_e),
(docs.fourier_number(D_w, tau, delta), Fo_w),
(docs.fourier_number(D_g, tau, delta), Fo_g),
(docs.fourier_number(D_b, tau, delta), Fo_b),
]Performing the replacements with help of SymPy:
eq = eq.subs(SUBS_FOURIER)
eq\(\displaystyle - Fo_{e} c_{E} + Fo_{e} c_{P} + Fo_{w} c_{P} - Fo_{w} c_{W} - c_{0} + c_{P} = 0\)
Notice that so far nothing has been said with respect to the instant at which the right-hand side is evaluated. The obvious solution would be to take \(c_{0}\equiv{}c_{P}^{0}\) as it is currently known, but for numerical stability reasons we will continue here with an implicit scheme so that all values other than \(c_{0}\) are defined to be evaluated at time \(\tau\), i.e. \(c_{K}\equiv{}c_{K}^{\tau}\), in phase with the notation convention introduced above.
Because it is more convenient for symbolic manipulations with SymPy, the problem is currently written as \(F(c)=0\), i.e. all elements are on the left-hand side. Care must be taken to properly interpret the signs of the coefficients below.
coefs = eq.coefficients(variables=[c_0, c_E, c_P, c_W, c_G, c_B])
coefs\(\displaystyle \left\{ c_{0} : -1, \ c_{E} : - Fo_{e}, \ c_{P} : Fo_{e} + Fo_{w} + 1, \ c_{W} : - Fo_{w}\right\}\)
Because it is more readable, a helper utility is provided to tabulate them:
eq.tabulate(variables=[c_0, c_E, c_P, c_W, c_G, c_B])| Variable | Coefficient | |
|---|---|---|
| 0 | c_0 | -1 |
| 1 | c_E | -Fo_e |
| 2 | c_P | Fo_e + Fo_w + 1 |
| 3 | c_W | -Fo_w |