4  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} \]

The interpolation of the diffusivities is no longer a simply harmonic mean as in the case of uniform grid spacing, but it can be computed using the diffusivity values at the neighboring nodes and the distances to the faces. For interface \(ij\) we have:

\[ D_{ij}= \dfrac{\Delta_{ij}}{D_{i}\delta_{i}^{-1}/2 + D_{j}\delta_{j}^{-1}/2} = \dfrac{(\delta_{i}+\delta_{j})/2}{D_{i}\delta_{i}^{-1}/2 + D_{j}\delta_{j}^{-1}/2} = \dfrac{(\delta_{i}+\delta_{j})}{D_{i}/\delta_{i} + D_{j}/\delta_{j}} \]

The above expression can be easily derived by considering the flux through the face as a series of two resistances (one for each neighboring node) and applying the concept of equivalent resistance in series. It falls back to the harmonic mean in the case of uniform grid spacing.

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}) \]

Applying this expression to the full problem, we have that the length squared to be used in the Fourier number for the boundary is \(\delta_{P}^{2}/2\). 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}} \]

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} \]