class StandardProblem:
""" Standard format of finite volume diffusion problem in 1-D. """
@staticmethod
def get_lhs(**kw):
""" Standard format of inner LHS of problem. """
return (c_P - c_0) * delta / tau
@staticmethod
def get_rhs(**kw):
""" Standard format of inner RHS of problem. """
J_e = kw.get("J_e", docs.diffusion_flux(D_e, c_P, c_E, delta))
J_w = kw.get("J_w", docs.diffusion_flux(D_w, c_W, c_P, delta))
return -1 * (J_e - J_w)
@staticmethod
def full_expression(subs_lhs=[], subs_rhs=[], subs_fourier=False, **kw):
""" Canonical form of problem set as `f(c)=0`. """
lhs = StandardProblem.get_lhs(**kw).subs(subs_lhs)
rhs = StandardProblem.get_rhs(**kw).subs(subs_rhs)
eq = docs.Eq(lhs - rhs) * tau / delta
return eq.subs(SUBS_FOURIER) if subs_fourier else eq3 Boundary conditions
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.
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 in the previous section 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\)
3.1 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(D_b, c_B, 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 Fo_{b} c_{B} + 2 Fo_{b} c_{P} - Fo_{e} c_{E} + Fo_{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 | -Fo_e |
| 2 | c_B | -2*Fo_b |
| 3 | c_P | 2*Fo_b + Fo_e + 1 |
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(D_b, c_B, delta/2)
J_w2 = diffusion_flux_west(D_g, c_G, delta)
docs.Eq(J_w1 - J_w2).tabulate()| Variable | Coefficient | |
|---|---|---|
| 0 | c_G | -D_g/delta |
| 1 | c_P | -2*D_b/delta + 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 = docs.Eq(J_w1 - J_w2).solve_for(c_G)[0]
c_G_sol = c_G_sol.subs(D_b, 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(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_B | -2*Fo_g |
| 2 | c_E | -Fo_e |
| 3 | c_P | Fo_e + 2*Fo_g + 1 |
3.2 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.
3.3 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 = -docs.convection_flux(c_B, h, c_inf)
q_be = +docs.convection_flux(c_B, h, c_inf)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 3.1. 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(D_b, c_B, delta/2)
J_e = diffusion_flux_east(D_b, c_B, delta/2)
c_Bw = docs.Eq(J_w - q_bw).solve_for(c_B)[0]
c_Be = docs.Eq(J_e - q_be).solve_for(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\(\displaystyle - 2 Fo_{b} c_{B} + 2 Fo_{b} c_{P} - Fo_{e} c_{E} + Fo_{e} c_{P} - c_{0} + c_{P} = 0\)
| Variable | Coefficient | |
|---|---|---|
| 0 | c_0 | -1 |
| 1 | c_E | -Fo_e |
| 2 | c_B | -2*Fo_b |
| 3 | c_P | 2*Fo_b + Fo_e + 1 |