15  Diffusion in 1-D

15.1 Import toolbox

from abc import ABC
from sympy import symbols, init_printing
from diffusion import (
    FvmVars,
    BcVars,
    CellVars,
    FaceVars,
    StandardEquation as Eq,
    StandardProblem,
    harmonic_mean,
    diffusion_flux,
    convection_flux,
    diffusion_flux_east,
    diffusion_flux_west,
    SUBS_FOURIER,
)
import casadi as cs
import numpy as np

init_printing(use_latex="mathjax")

15.2 Introductory Fickean diffusion

Assume a simple concentration-dependent diffusion in the absence of external driving forces; in this introduction we chose to model the diffusive flux through Fick’s constitutive (first) law stated below, where \(D\) is the diffusion coefficient that may depend on position \(x\) and/or concentration \(c\), i.e. \(D\equiv{}D(x, c)\). It states that diffusive flux of an species has the composition gradient as a driver force; the negative sign indicates that diffusion occurs from high to low concentration regions.

\[ J=-D\dfrac{\partial{}c}{\partial{}x} \]

We can anticipate here that this can be the case for simple diluted systems, but interactions between species must be considered in concentrated solutions, what is generalized by the Irreversible Thermodynamics by Onsager (1931a) and Onsager (1931b). Here we emphasize again the word model: empirical observation has led scientists to represent reality through such an approximation; different models can be used for the same physics with their applicability limited to certain scenarios. Properly selecting an applicable model is the role of the researcher/engineer, what is outside our scope here; in what follows we will focus on how to solve balance equations for a given model.

Simply put, diffusion equation solves for a local mass balance; and mathematically, a local balance leads to a divergence operation. Following this idea, it can be shown that the time derivative of concentration \(c\) at one location is given by the divergence - which collapses to the spatial derivative - of the negative of mass flux given by Fick’s (first) law (or in general any other constitutive law/model describing the mass flux); sometimes this form is called Fick’s second law in the literature. Given the possible non-linearity introduced by diffusion coefficient \(D\), the derivative in the right-hand side is not expanded. In Cartesian coordinate system, the governing diffusion equation in one dimension writes:

\[ \dfrac{\partial{}c}{\partial{}t}= \dfrac{\partial{}}{\partial{}x} \left(D\dfrac{\partial{}c}{\partial{}x}\right) \]

This is the form of the diffusion equation that is discussed in what follows.

15.3 Internal discretization

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.

Below we compose this expression with help of SymPy. 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 = (CellVars.c_P - CellVars.c_0) * FvmVars.delta / FvmVars.tau
rhs = -1 * (J_e - J_w)

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

eq = eq * FvmVars.tau / FvmVars.delta
eq

\(\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 \(\alpha_{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.

\[ \alpha_{k}=\dfrac{\tau{}D_{k}}{\delta^2} \]

Performing the replacements with help of SymPy:

eq = eq.subs(SUBS_FOURIER)
eq

\(\displaystyle - \alpha_{e} c_{E} + \alpha_{e} c_{P} + \alpha_{w} c_{P} - \alpha_{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()
coefs

\(\displaystyle \left\{ c_{0} : -1, \ c_{E} : - \alpha_{e}, \ c_{P} : \alpha_{e} + \alpha_{w} + 1, \ c_{W} : - \alpha_{w}\right\}\)

Because it is more readable, a helper utility is provided to tabulate them:

eq.tabulate()

	Variable	Coefficient
0	c_0	-1
1	c_E	-alpha_e
2	c_P	alpha_e + alpha_w + 1
3	c_W	-alpha_w

15.4 Boundary conditions

So far all we have derived are the internal balances of the system, we still lack proper boundary conditions. It should be evident by now that we have formulated a tridiagonal matrix for representing the 1-D system: there are 3 unknowns at step \(\tau\) (\(c_W\), \(c_P\), and \(c_E\)) and a currently known field \(c_0\).

In what follows we will workout the three families of boundary conditions: Dirichlet (type 1), Neumann (type 2), and the more general Robin (type 3). Following this sequence will create the required intuition for deriving the progressively more complex conditions. In practice, for diffusion in solids often only the implementation of Robin conditions is required, as it will be discussed later. Our goal at presenting types 1 and 2 is simply to introduce the tooling required for working out type 3 boundary condition; they will not be treated in full detail, leaving the full discussion for later.

The implementation of boundary conditions requires special treatment of the control volumes at the domain boundaries. The general form for a boundary control volume is:

\[ \left(c_{P}-c_{P}^{0}\right) \dfrac{\delta}{\tau}= \text{flux}_{\text{interior}} - \text{flux}_{\text{boundary}} \]

It is worth here repeating the partial development derived above for an analogy:

\[ \left(c_{P}^{\tau}-c_{P}^{0}\right) \dfrac{\delta}{\tau}= D_{e}\dfrac{c_{E}-c_{P}}{\delta}- D_{w}\dfrac{c_{P}-c_{W}}{\delta} \]

Please, notice that from now on we make use of class StandardProblem for the end of representing a canonical discrete equation. It encodes the results obtained in the previous section and is also repeated below:

StandardProblem.full_expression()

\(\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\)

15.5 Type 1 BC - Dirichlet

Let’s start by the Dirichlet, or fixed concentration, boundary condition (type 1), for which the concentration at the boundary face is prescribed, i.e. \(c_{\text{B}} = c_{\text{prescribed}}\). Notice that before using this we need to make a choice regarding the computational grid (mesh); you could have half-volumes on boundaries, fictitious volumes (sometimes called ghost cells), or immersed nodes over the boundary.

Partial control volumes on boundaries are not desirable as they introduce implementation complexities in code and lack conservation (especially in higher dimensions), see a discussion by Maliska (2004) (section 3.7) for details. In what follows we discuss immersed nodes and extended domains with ghost cells. This last solution implies considerably more unknowns in 3-D problems, for which different variants of FVM will be discussed at a more advanced point.

Immersed boundary: Assume the last cell on west boundary is indexed as \(P\); over its west boundary we find a node where prescribed composition \(c_{B}\) is enforced - and the corresponding diffusion coefficient \(D_{b}\). Notice that this node does not belong to the cell, nor represents a half-cell, but lies over its boundary. To keep it simple, let’s assume \(c_{B}\) is constant (time-dependency would slightly complicate the development and it is worth postponing it to the algorithm development discussion), so keeping our fully implicit scheme leads to We start by computing the mass flux from the boundary node to the first cell center; this flux is computed using half-cell length \(\delta/2\) for the derivative approximation.

J_w = diffusion_flux_west(FaceVars.D_b, CellVars.c_B, FvmVars.delta/2)
J_w

\(\displaystyle - \frac{2 D_{b} \left(- c_{B} + c_{P}\right)}{\delta}\)

This expression of flux is used as replacement for the west boundary in the full equation:

eq = StandardProblem.full_expression(J_w=J_w, subs_fourier=True)
eq

\(\displaystyle - 2 \alpha_{b} c_{B} + 2 \alpha_{b} c_{P} - \alpha_{e} c_{E} + \alpha_{e} c_{P} - c_{0} + c_{P} = 0\)

The coefficients for the first row of the matrix problem are summarized below:

eq.tabulate()

	Variable	Coefficient
0	c_0	-1
1	c_E	-alpha_e
2	c_P	2*alpha_b + alpha_e + 1
3	c_B	-2*alpha_b

There are two main features in this equation: the coefficient of cell \(P\) now depends on the prescribed composition through \(\alpha_{b}\) and a modification is applied to the right-hand side at boundaries where solution is prescribed. That’s a good reminder that the array of coefficients must correspond to the number of interfaces, not cells - more on this later.

The immersed boundary approach presented above is simple. Nonetheless, it requires a different calculation of boundary nodes with respect to the volume, but that comes with the advantage that no additional equation has been added to the system. Let’s now explore the extended domain alternative.

Ghost cell: As above, assume the last cell on west boundary is indexed as \(P\). Consider another ghost cell to its left being denoted \(G\) (fictitious cell in place of \(W\)). The node \(c_{B}\) where composition is enforced and the same details presented above still apply here. Here we repeat the immersed node flux computation and do the same with respect to ghost cell \(G\).

Our goal is to identify an equation for \(c_{G}\) so that it is solved along with the interior points and ensures the same flux as the immersed node. To identify the compatible ghost composition we equate the fluxes by both definitions and solve for \(c_{G}\). The following steps repeat what we have seem so far and should be self-explanatory.

J_w1 = diffusion_flux_west(FaceVars.D_b, CellVars.c_B, FvmVars.delta/2)
J_w2 = diffusion_flux_west(FaceVars.D_g, CellVars.c_G, FvmVars.delta)

Eq(J_w1 - J_w2).tabulate()

	Variable	Coefficient
0	c_P	-2*D_b/delta + D_g/delta
1	c_G	-D_g/delta
2	c_B	2*D_b/delta

Care must be taken for interpreting the above result; here \(c_{G}\) takes the diagonal position and \(c_{P}\) is the first east cell in the first row of the problem’s matrix. There is no previous state (no coefficient on \(c_{0}\)), but a source term based on \(c_{B}\) that must be placed on the right-hand side for solution. This is the equation we are looking for to extend the system.

Since this equation has \(D_{g}\) as a coefficient, an initial guess for \(c_{G}\) is required. The boundary \(B\) being the same interface for evaluation of \(D_{g}\), for both fluxes to equate one needs to make sure the same diffusivity is evaluated there. Thus we make \(D_{b}\equiv{}D_{g}\), what seems physically reasonable but so far not supported by any mathematical argument:

c_G_sol = Eq(J_w1 - J_w2).solve_for(CellVars.c_G)[0]
c_G_sol = c_G_sol.subs(FaceVars.D_b, FaceVars.D_g)
c_G_sol

\(\displaystyle \frac{2 D_{g} c_{B} - D_{g} c_{P}}{D_{g}}\)

This simply says that the composition is linearly extrapolated, with the distances from cell centers to the boundary being used as weighting factors. It ensures the problem is initialized the same as the immersed boundary as we made \(D_g=D_b\).

CautionCommon mistake

One could think about solving for \(c_G\) and using the result in the standard balance equation, but doing so falls back to the immersed boundary approach. Replacing the solution in J_w2 simply gets back to the initial point and does not produce one extra equation!

J_w = J_w2.subs(CellVars.c_G, c_G_sol)

eq = StandardProblem.full_expression(J_w=J_w, subs_fourier=True)
eq.tabulate()

	Variable	Coefficient
0	c_0	-1
1	c_E	-alpha_e
2	c_P	alpha_e + 2*alpha_g + 1
3	c_B	-2*alpha_g

15.6 Type 2 BC - Neumann

Next, let’s handle Neumann boundary condition (type 2). For a fixed flux boundary condition, the flux at the boundary is prescribed at \(q_{B}\). A special case arises when one imposes zero flux (impermeable boundary for mass transfer, adiabatic boundary for heat transfer) where \(q_{B}=0\). It is worthless formulating the problem for this case as the equation is trivial, so we will assume a general \(q_{B}\) in what follows.

\[ -D\dfrac{\partial{}c}{\partial{}x}\bigg|_{B} = q_{B} \]

When dealing with Neumann boundary conditions one must take care of the sign of fluxes, especially in 1-D. In higher dimensions, when applying Gauss theorem to convert the volume integral into a surface integral one normally works with vector calculus formulations; here we have dropped the vector notation when integrating over \(dx\), but the concept of surface orientation must be recalled at this point. At the west boundary the flux into domain is positive and conversely, at the east boundary the flux out of domain is positive.

\[ \begin{aligned} \text{(west)}\qquad & \left(c_{P}^{\tau}-c_{P}^{0}\right)\dfrac{\delta}{\tau}= D_{e}\dfrac{c_{E}-c_{P}}{\delta} + q_{B} \\[12pt] \text{(east)}\qquad & \left(c_{P}^{\tau}-c_{P}^{0}\right)\dfrac{\delta}{\tau}= -q_{B} - D_{w}\dfrac{c_{P}-c_{W}}{\delta} \end{aligned} \]

We will not elaborate more as the most relevant condition represented by this class, the impermeable condition, can be represented with a type 3 boundary condition by setting the mass transfer coefficient \(h\) to zero.

15.7 Type 3 BC - Robin

Finally, we generalize the above discussion for a Robin boundary condition (type 3). Also known as a convective or mixed boundary condition. In this case the flux is proportional to the concentration difference between the boundary and a reference external composition \(c_{\infty}\). This often models mass transfer between a reacting fluid and a solid, but the discussion of boundary layer formation is out of the present scope and we summarize it simply by the mass transfer coefficient \(h\).

\[ -D\dfrac{\partial c}{\partial x}\bigg|_{B} = h(c_{B} - c_{\infty}) \]

NoteFallback to other kinds

In fact, this is not the canonical form of such a boundary condition, but the reasonable engineering statement. Below we transform the equation to a format more commonly found the in mathematical literature, where all elements should be interpreted as being evaluated at boundary \(B\).

\[ \alpha{}c + \beta\dfrac{\partial{}c}{\partial{}x} = g \]

where

\[ \alpha = 1 \quad\text{and}\quad \beta = \dfrac{D}{h} \quad\text{and}\quad g = c_{\infty} \]

Under this form we can easily jump to conclusions with respect to its relationship to the other types of boundary conditions. First, simply by setting \(\beta=0\) we find ourselves with a type 1 condition where \(c=\alpha^{-1}g\) - or in engineering notation \(c_{B}=c_{\infty}\); this happens as mass transfer efficiency is ideal, as \(\beta\to{}0\) when \(h\to{}\infty\). In computational practice it is enough to enforce that \(h\gg{}D\) to approach a Dirichlet boundary condition.

A type 2 condition would be established with \(\alpha=0\), but that is not possible straight from the way we have posed the equation. By multiplying everything by \(h\) and setting its value to zero we get to an impermeable Neumann condition, the only one we consider here for the mass transfer context (because in general it is not simple to enforce a mass uptake by a material).

That said, type 3 condition is all we need for a first solid state diffusion solver. As for the type 2 condition, we need to take care of signs here.

The (scalar) convention we used so far implies a convective flux that is positive outwards the domain, \(J=h(c_{b}-c_{\infty})>0\) if \(c_{b}>c_{\infty}\), i.e., the chemical potential in the solid is higher than the environment. The west boundary being oriented towards the \(-y\) axis, we need to account for its sign; in later chapters we will see that the boundary condition reasoning is actually simpler to follow in higher dimensions as we implicitly treat the face normal vectors through the inner product.

q_bw = -convection_flux(CellVars.c_B)
q_be = +convection_flux(CellVars.c_B)

In what follows we derive the immersed boundary equation for type 3 condition. The extension to ghost cells is trivial and can be worked out with the ideas presented in section Section 15.5. We solve for \(c_{B}\) at both sides of the domain by equating the convective flux to the one established internally with the immersed boundary node:

J_w = diffusion_flux_west(FaceVars.D_b, CellVars.c_B, FvmVars.delta/2)
J_e = diffusion_flux_east(FaceVars.D_b, CellVars.c_B, FvmVars.delta/2)

c_Bw = Eq(J_w - q_bw).solve_for(CellVars.c_B)[0]
c_Be = Eq(J_e - q_be).solve_for(CellVars.c_B)[0]

This leads to

c_Bw

\(\displaystyle \frac{2 D_{b} c_{P} - c_{\infty} \delta h}{2 D_{b} - \delta h}\)

and

c_Be

\(\displaystyle \frac{2 D_{b} c_{P} - c_{\infty} \delta h}{2 D_{b} - \delta h}\)

If we take, e.g., the west boundary equation with the previous definition of flux evaluated from the immersed boundary node, we retrieve the corresponding equation for the system. Notice that \(c_{B}\) is a known value at the current iteration, so in solution it can be moved to the right-hand side and act as a flux term. This gives a grasp of the solution algorithm that makes the subject for the next section.

eq = StandardProblem.full_expression(J_w=J_w, subs_fourier=True)
eq.tabulate()

	Variable	Coefficient
0	c_0	-1
1	c_E	-alpha_e
2	c_P	2*alpha_b + alpha_e + 1
3	c_B	-2*alpha_b

15.8 Arbitrary grid spacing

For establishing the equations for an arbitrary grid spacing in 1-D, we move forward to a more general derivation. Instead of expressing the problem over a single coordinate, we introduce the control volume \(V\) that will be revisited again when solving for 3-D geometries. For convenience the minus sign is moved outside the integral in what follows.

\[ \int_{V}\int_{0}^{\tau} \dfrac{\partial{}c}{\partial{}t}\:dt\:dV= -\int_{0}^{\tau}\int_{V} \nabla\cdot\vec{J}\:dV\:dt \]

The inner integration on the right-hand side is treated first. Using the divergence theorem, we convert the volume integral of the divergence of flux into a surface integral over the control volume’s faces. Since in FVM we are dealing with a finite number of faces, we can express the surface integral as a sum over all faces \(f\) of the control volume, where \(J_{f}\) is the flux through face \(f\) and \(A_{f}\) is the area of that face.

\[ \int_{V}\nabla\cdot\vec{J}\:dV= \int_{S}\vec{J}\cdot\vec{n}\:dS= \sum_{f}J_{f}A_{f} \]

Substituting this result back into the original equation and integrating the left-hand side over time gives:

\[ \int_{V}\left(c_{P}-c_{P}^{0}\right)\:dV= -\int_{0}^{\tau}\sum_{f}J_{f}A_{f}\:dt \]

And finally integrating the last level:

\[ \left(c_{P}-c_{P}^{0}\right)\dfrac{V}{\tau}= -\vec{J}_{f}\cdotp\vec{A}_{f} \]

In 1-D, there are only two faces: the west face (\(w\)) and the east face (\(e\)). Thus, the summation over faces simplifies to just these two contributions:

\[ \left(c_{P}-c_{P}^{0}\right)\dfrac{V}{\tau}= -J_{w}(-\hat{i})A_{w}(-\hat{i}) -J_{e}(+\hat{i})A_{e}(+\hat{i}) \]

Notice that as we are dealing with the actual surfaces, the signs of the terms will be naturally handled by the flux definitions. The final step consists of substituting the expressions for the diffusive fluxes at each face, which can be approximated using finite differences. Below we account already for the face orientation (sign). For the west face:

\[ -J_{w}(-\hat{i})A_{w}(-\hat{i})=D_{w}\dfrac{c_{W}-c_{P}}{\delta_{w}}A_{w} \]

And for the east face:

\[ -J_{e}(+\hat{i})A_{e}(+\hat{i})=D_{e}\dfrac{c_{E}-c_{P}}{\delta_{e}}A_{e} \]

Substituting these flux expressions back into the integrated equation gives:

\[ \left(c_{P}-c_{P}^{0}\right)\dfrac{V}{\tau}= D_{w}\dfrac{c_{W}-c_{P}}{\delta_{w}}A_{w} D_{e}\dfrac{c_{E}-c_{P}}{\delta_{e}}A_{e} \]

Rearranging terms leads to the final discrete equation for an arbitrary grid spacing in 1-D implicit formulation:

\[ -\dfrac{D_{w}A_{w}}{\delta_{w}}c_{W} +\left( \dfrac{V}{\tau}+ \dfrac{D_{w}A_{w}}{\delta_{w}}+ \dfrac{D_{e}A_{e}}{\delta_{e}} \right)c_{P} -\dfrac{D_{e}A_{e}}{\delta_{e}}c_{E}= \dfrac{V}{\tau}c_{P}^{0} \]

or assuming a unit cross-sectional area \(A\) so that \(V=A\delta_{P}\), we have:

\[ -\dfrac{D_{w}\tau}{\delta_{w}\delta_{P}}c_{W} +\left( 1 + \dfrac{D_{w}\tau}{\delta_{w}\delta_{P}}+ \dfrac{D_{e}\tau}{\delta_{e}\delta_{P}} \right)c_{P} -\dfrac{D_{e}\tau}{\delta_{e}\delta_{P}}c_{E}= c_{P}^{0} \]

For the boundary conditions we apply the same approach as before:

\[ -J_{b}(-\hat{i})A_{b}(-\hat{i})= D_{b}\dfrac{c_{B}-c_{P}}{\frac{1}{2}\delta_{P}}A_{b}= hA_{b}(c_{\infty}-c_{B}) \]

Before simplifying this expression to find out \(c_{B}\) it is worth introducing the Sherwood number \(\mathcal{Sh}\) defined as follows, where \(L\) represents a characteristic length of the problem and \(D\) the diffusion coefficient; here \(L=\delta_{P}/2\) as half-cell length was used in flux computations and

\[ \mathcal{Sh}=\dfrac{hL}{D} \quad\text{and}\quad \mathcal{Sh}_{b}=\dfrac{h\delta_{P}}{2D_{b}} \]

leading to the final boundary node concentration

\[ c_{B}= \dfrac{c_{P}+\mathcal{Sh}_{b}\:c_{\infty}}{1 + \mathcal{Sh}_{b}} \]

15.9 Solution algorithm

Using a segregated solution of non-linear boundaries or immersed nodes:

\[ \begin{bmatrix} A_{(0)} & A_{E} & 0 & & & \dots & 0 & 0 \vphantom{\dfrac{1}{1}}\\ A_{W} & A_{P} & A_{E} & & & \dots & 0 & 0 \vphantom{\dfrac{1}{1}}\\ 0 & A_{W} & A_{P} & A_{E} & 0 & \dots & 0 & 0 \vphantom{\dfrac{1}{1}}\\ 0 & 0 & \ddots & \ddots & \ddots & \ddots & \vdots & \vdots \vphantom{\dfrac{1}{1}}\\ \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & 0 & 0 \vphantom{\dfrac{1}{1}}\\ 0 & 0 & \dots & 0 & A_{W} & A_{P} & A_{E} & 0 \vphantom{\dfrac{1}{1}}\\ 0 & 0 & \dots & 0 & 0 & A_{W} & A_{P} & A_{E} \vphantom{\dfrac{1}{1}}\\ 0 & 0 & \dots & 0 & 0 & 0 & A_{W} & A_{(n)} \vphantom{\dfrac{1}{1}}\\ \end{bmatrix}\cdotp %-------------------------------------- \begin{bmatrix} c_{(0)}^{\tau} \vphantom{\dfrac{1}{1}}\\ c_{(1)}^{\tau} \vphantom{\dfrac{1}{1}}\\ c_{(2)}^{\tau} \vphantom{\dfrac{1}{1}}\\ c_{(3)}^{\tau} \vphantom{\dfrac{1}{1}}\\ \vdots \vphantom{\dfrac{1}{1}}\\ c_{(N-3)}^{\tau} \vphantom{\dfrac{1}{1}}\\ c_{(N-2)}^{\tau} \vphantom{\dfrac{1}{1}}\\ c_{(N-1)}^{\tau} \vphantom{\dfrac{1}{1}}\\ \end{bmatrix}= %-------------------------------------- \begin{bmatrix} c_{(0)}^{0} \vphantom{\dfrac{1}{1}}\\ c_{(1)}^{0} \vphantom{\dfrac{1}{1}}\\ c_{(2)}^{0} \vphantom{\dfrac{1}{1}}\\ c_{(3)}^{0} \vphantom{\dfrac{1}{1}}\\ \vdots \vphantom{\dfrac{1}{1}}\\ c_{(N-3)}^{0} \vphantom{\dfrac{1}{1}}\\ c_{(N-2)}^{0} \vphantom{\dfrac{1}{1}}\\ c_{(N-1)}^{0} \vphantom{\dfrac{1}{1}}\\ \end{bmatrix}+ %-------------------------------------- \begin{bmatrix} B_{(0)} \vphantom{\dfrac{1}{1}}\\ 0 \vphantom{\dfrac{1}{1}}\\ 0 \vphantom{\dfrac{1}{1}}\\ 0 \vphantom{\dfrac{1}{1}}\\ \vdots \vphantom{\dfrac{1}{1}}\\ 0 \vphantom{\dfrac{1}{1}}\\ 0 \vphantom{\dfrac{1}{1}}\\ B_{(N-1)} \vphantom{\dfrac{1}{1}}\\ \end{bmatrix} \]

15.10 Python prototype

NoteDraft 1 - cell-node approach

class Cell:
    __slots__ = ("_x",)

    def __init__(self, center) -> None:
        self._x = center

class NodeCell(Cell):
    def __repr__(self) -> str:
        return f"<NodeCell x={self._x} />"

class LineCell(Cell):
    __slots__ = ("_d",)

    def __init__(self, center, length) -> None:
        super().__init__(center)
        self._d = length

    def __repr__(self) -> str:
        return f"<LineCell x={self._x} d={self._d} />"

class CellFace:
    __slots__ = ("_center", "_normal",)

    def __init__(self, center, normal) -> None:
        self._center = center
        self._normal = normal

class CellFace1D(CellFace):
    def __init__(self, center, normal) -> None:
        super().__init__(center, normal)

    def __repr__(self) -> str:
        return f"<CellFace1D center={self._center} normal={self._normal} />"

class Grid1D:
    __slots__ = ["_cells", "_faces"]

    def __init__(self, x):
        if len(x) < 3:
            raise ValueError("At least 3 cells are required!")

        # TODO create an networkx unstructured graph.

        # Coordinates of cell faces (first and last coordinate are
        # implicitly the boundaries of the system).
        f = np.hstack((x[0], (x[2:-1] + x[1:-2]) / 2, x[-1]))

        # Cell lengths are simply the distance between faces:
        L = f[1:] - f[:-1]

        # Create inner cells (LineCell objects):
        inner = map(lambda a: LineCell(*a), zip(x[1:-1], L[:]))

        # Stack boundaries and inner cells:
        self._cells = [NodeCell(x[0]), *inner, NodeCell(x[-1])]

        # Assembly faces:
        self._faces = [CellFace1D(g, +1) for g in f[1:]]
        self._faces.insert(0, CellFace1D(f[0], -1))

Lx = 1.00
dx = 0.1

xc1 = dx/2
xcn = Lx - dx/2 + 0.01*dx

x = np.hstack((0, np.arange(xc1, xcn, dx), Lx))

grid = Grid1D(x)
for c in grid._faces:
    print(c)

<CellFace1D center=0.0 normal=-1 />
<CellFace1D center=0.1 normal=1 />
<CellFace1D center=0.20000000000000004 normal=1 />
<CellFace1D center=0.30000000000000004 normal=1 />
<CellFace1D center=0.4 normal=1 />
<CellFace1D center=0.5000000000000001 normal=1 />
<CellFace1D center=0.6000000000000001 normal=1 />
<CellFace1D center=0.7000000000000002 normal=1 />
<CellFace1D center=0.8000000000000003 normal=1 />
<CellFace1D center=0.9000000000000001 normal=1 />
<CellFace1D center=1.0 normal=1 />

x = np.array(symbols("x:10"))

grid = Grid1D(x)
for c in grid._cells:
    print(c)

for c in grid._faces:
    print(c)

<NodeCell x=x0 />
<LineCell x=x1 d=-x0 + x1/2 + x2/2 />
<LineCell x=x2 d=-x1/2 + x3/2 />
<LineCell x=x3 d=-x2/2 + x4/2 />
<LineCell x=x4 d=-x3/2 + x5/2 />
<LineCell x=x5 d=-x4/2 + x6/2 />
<LineCell x=x6 d=-x5/2 + x7/2 />
<LineCell x=x7 d=-x6/2 + x8/2 />
<LineCell x=x8 d=-x7/2 - x8/2 + x9 />
<NodeCell x=x9 />
<CellFace1D center=x0 normal=-1 />
<CellFace1D center=x1/2 + x2/2 normal=1 />
<CellFace1D center=x2/2 + x3/2 normal=1 />
<CellFace1D center=x3/2 + x4/2 normal=1 />
<CellFace1D center=x4/2 + x5/2 normal=1 />
<CellFace1D center=x5/2 + x6/2 normal=1 />
<CellFace1D center=x6/2 + x7/2 normal=1 />
<CellFace1D center=x7/2 + x8/2 normal=1 />
<CellFace1D center=x9 normal=1 />

NoteDraft 2 - OpenFOAM approach

Coordinate = tuple[float, ...]


class GridGraph(ABC):
    """ A computational mesh similar to the one employed by OpenFOAM.

    References:
    - https://doc.cfd.direct/notes/cfd-general-principles/computational-mesh
    - https://doc.cfd.direct/openfoam/user-guide-v13/mesh-files
    """
    __slots__ = (
        "vertex", # a list of every vertex used to define all faces;
        "faces",  # a list of faces, each defined by a sequence of vertex indices;
        "owner",  # a list of owner cell indices associated with the faces;
        "neigh",  # a list of neighbour cell indices associated with faces;
        "patches" # grouping boundary faces into a set of patches.
    )
    def __init__(self, *args, **kwargs) -> None:
        self.vertex: list[Coordinate] = []
        self.faces: list[int] = []
        self.owner: list[int] = []
        self.neigh: list[int] = []
        self.patches: dict[int, list[int]] = []


class GridGraph1D(GridGraph):
    """ One dimensional implementation of GridGraph. """
    def __init__(self, *args, **kwargs) -> None:
        super().__init__(*args, **kwargs)

        (xf,) = args

        self.vertex = xf

First we define the \(x\) coordinates of faces in grid:

Lx = 1.00
dx = 0.1

xf = np.arange(0, Lx+dx/2, dx)

For a pseudo-one-dimensional mesh we need each point over \(x\) to actually be a face of unit area; this can be accomplished, for instance, by surrounding the axis with a unit square such as provided in the pseudo-code below:

0: (x, -1/2, -1/2)
1: (x, +1/2, -1/2)
2: (x, -1/2, +1/2)
3: (x, +1/2, +1/2)

square = [
    (-1/2, -1/2),
    (+1/2, -1/2),
    (-1/2, +1/2),
    (+1/2, +1/2),
]

x =  [(xf_i, y, z) for xf_i in xf for (y, z) in square]
x

\(\displaystyle \left[ \left( 0.0, \ -0.5, \ -0.5\right), \ \left( 0.0, \ 0.5, \ -0.5\right), \ \left( 0.0, \ -0.5, \ 0.5\right), \ \left( 0.0, \ 0.5, \ 0.5\right), \ \left( 0.1, \ -0.5, \ -0.5\right), \ \left( 0.1, \ 0.5, \ -0.5\right), \ \left( 0.1, \ -0.5, \ 0.5\right), \ \left( 0.1, \ 0.5, \ 0.5\right), \ \left( 0.2, \ -0.5, \ -0.5\right), \ \left( 0.2, \ 0.5, \ -0.5\right), \ \left( 0.2, \ -0.5, \ 0.5\right), \ \left( 0.2, \ 0.5, \ 0.5\right), \ \left( 0.3, \ -0.5, \ -0.5\right), \ \left( 0.3, \ 0.5, \ -0.5\right), \ \left( 0.3, \ -0.5, \ 0.5\right), \ \left( 0.3, \ 0.5, \ 0.5\right), \ \left( 0.4, \ -0.5, \ -0.5\right), \ \left( 0.4, \ 0.5, \ -0.5\right), \ \left( 0.4, \ -0.5, \ 0.5\right), \ \left( 0.4, \ 0.5, \ 0.5\right), \ \left( 0.5, \ -0.5, \ -0.5\right), \ \left( 0.5, \ 0.5, \ -0.5\right), \ \left( 0.5, \ -0.5, \ 0.5\right), \ \left( 0.5, \ 0.5, \ 0.5\right), \ \left( 0.6, \ -0.5, \ -0.5\right), \ \left( 0.6, \ 0.5, \ -0.5\right), \ \left( 0.6, \ -0.5, \ 0.5\right), \ \left( 0.6, \ 0.5, \ 0.5\right), \ \left( 0.7, \ -0.5, \ -0.5\right), \ \left( 0.7, \ 0.5, \ -0.5\right), \ \left( 0.7, \ -0.5, \ 0.5\right), \ \left( 0.7, \ 0.5, \ 0.5\right), \ \left( 0.8, \ -0.5, \ -0.5\right), \ \left( 0.8, \ 0.5, \ -0.5\right), \ \left( 0.8, \ -0.5, \ 0.5\right), \ \left( 0.8, \ 0.5, \ 0.5\right), \ \left( 0.9, \ -0.5, \ -0.5\right), \ \left( 0.9, \ 0.5, \ -0.5\right), \ \left( 0.9, \ -0.5, \ 0.5\right), \ \left( 0.9, \ 0.5, \ 0.5\right), \ \left( 1.0, \ -0.5, \ -0.5\right), \ \left( 1.0, \ 0.5, \ -0.5\right), \ \left( 1.0, \ -0.5, \ 0.5\right), \ \left( 1.0, \ 0.5, \ 0.5\right)\right]\)

# grid = GridGraph1D(xf)
# grid.vertex

Maliska, Clovis R. 2004. Transferência de Calor e Mecânica Dos Fluidos Computacional. 2nd ed. Rio de Janeiro: LTC.

Onsager, Lars. 1931a. “Reciprocal Relations in Irreversible Processes. i.” Physical Review 37: 405–26.

———. 1931b. “Reciprocal relations in irreversible processes. II.” Physical Review 38: 2265–79.

Patankar, Suhas V. 1980. Numerical Heat Transfer and Fluid Flow. Taylor & Francis. https://doi.org/10.1201/9781482234213.