tug/doc/ADI_scheme.org

Finite Difference Schemes for the numerical solution of heterogeneous diffusion equation in 2D

Homogeneous diffusion in 1D

Finite differences with nodes as cells' centres

The 1D diffusion equation for spatially constant diffusion coefficients $\alpha$ is:

\begin{align} \frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha \frac{\partial C }{\partial x} \right) \nonumber \\ & = \alpha \frac{\partial^2 C}{\partial x^2} \end{align}

We aim at numerically solving eqn:1 on a spatial grid such as:

/max/tug/src/commit/ee77b5f7f3cbf99c4bbe4ef779dd532b72a3f8f3/doc/images/grid_pqc.pdf

The left boundary is defined on $x=0$ while the center of the first cell - which are the points constituting the finite difference nodes - is in $x=dx/2$, with $dx=L/n$.

The explicit FTCS scheme (as in PHREEQC)

We start by discretizing eqn:1 following an explicit Euler scheme and specifically a Forward Time, Centered Space finite difference.

For each cell index $i \in 1, \dots, n-1$ and assuming constant $\alpha$, we can write:

\begin{equation}\displaystyle \frac{C_i^t+1 -C_i^t}{Δ t} = α\frac{\frac{C^t_i+1-C^t_i}{Δ x}-\frac{C^t_i-C^t_i-1}{Δ x}}{Δ x}

\end{equation}

In practice, we evaluate the first derivatives of $C$ w.r.t. $x$ on the boundaries of each cell (i.e., $(C_{i+1}-C_i)/\Delta x$ on the right boundary of the i-th cell and $(C_{i}-C_{i-1})/\Delta x$ on its left cell boundary) and then repeat the differentiation to get the second derivative of $C$ on the the cell centre $i$.

This discretization works for all internal cells, but not for the domain boundaries ($i=0$ and $i=n$). To properly treat them, we need to account for the discrepancy in the discretization.

For the first (left) cell, whose center is at $x=dx/2$, we can evaluate the left gradient with the left boundary using such distance, calling $l$ the numerical value of a constant boundary condition:

\begin{equation}\displaystyle \frac{C_0^t+1 -C_0^t}{Δ t} = α\frac{\frac{C^t₁-C^t₀}{Δ x}- \frac{C^t₀-l}{\frac{Δ x}{2}}}{Δ x}

\end{equation}

This expression, once developed, yields:

\begin{align}\displaystyle C_0^t+1 & = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( C^t₁-C^t₀- 2 C^t₀+2l \right) \nonumber
& = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( C^t₁- 3 C^t₀ +2l \right)

\end{align}

In case of constant right boundary, the finite difference of point $C_n$ - calling $r$ the right boundary value - is:

\begin{equation}\displaystyle \frac{C_n^t+1 -C_n^t}{Δ t} = α\frac{\frac{r - C^t_n}{\frac{Δ x}{2}}- \frac{C^t_n-C^t_n-1}{Δ x}}{Δ x}

\end{equation}

Which, developed, gives

\begin{align}\displaystyle C_n^t+1 & = C_n^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( 2 r - 2 C^t_n -C^t_n + C^t_n-1 \right) \nonumber
& = C_n^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( 2 r - 3 C^t_n + C^t_n-1 \right)

\end{align}

If on the right boundary we have closed or Neumann condition, the left derivative in eq. eqn:5 becomes zero and we are left with:

\begin{equation}\displaystyle C_n^t+1 = C_n^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ (C^t_n-1 - C^t_n)

\end{equation}

A similar treatment can be applied to the BTCS implicit scheme.

Implicit BTCS scheme

First, we define the Backward Time difference:

\begin{equation} \frac{\partial C^{t+1} }{\partial t} = \frac{C^{t+1}_i - C^{t}_i}{\Delta t} \end{equation}

Second the spatial derivative approximation, evaluated at time level $t+1$:

\begin{equation} \frac{\partial^2 C^{t+1} }{\partial x^2} = \frac{\frac{C^{t+1}_{i+1}-C^{t+1}_{i}}{\Delta x}-\frac{C^{t+1}_{i}-C^{t+1}_{i-1}}{\Delta x}}{\Delta x} \end{equation}

Taking the 1D diffusion equation from eqn:1 and substituting each term by the equations given above leads to the following equation:

\begin{align}\displaystyle \frac{C_i^t+1 - C_i^t}{Δ t} & = α\frac{\frac{C^t+1_i+1-C^t+1_i}{Δ x}-\frac{C^t+1_i-C^t+1_i-1}{Δ x}}{Δ x} \nonumber
& = α\frac{C^t+1_i-1 - 2C^t+1_i + C^t+1_i+1}{Δ x^2}

\end{align}

This only applies to inlet cells with no ghost node as neighbor. For the left cell with its center at $\frac{dx}{2}$ and the constant concentration on the left ghost node called $l$ the equation goes as followed:

\begin{equation}\displaystyle \frac{C_0^t+1 -C_0^t}{Δ t} = α\frac{\frac{C^t+1₁-C^t+1₀}{Δ x}- \frac{C^t+1₀-l}{\frac{Δ x}{2}}}{Δ x}

\end{equation}

This expression, once developed, yields:

\begin{align}\displaystyle C_0^t+1 & = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( C^t+1₁-C^t+1₀- 2 C^t+1₀+2l \right) \nonumber
& = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( C^t+1₁- 3 C^t+1₀ +2l \right)

\end{align}

Now we define variable $s_x$ as followed:

\begin{equation} s_x = \frac{\alpha \cdot \Delta t}{\Delta x^2} \end{equation}

Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model:

\begin{equation}\displaystyle -C^t_0 = (2s_x) ⋅ l + (-1 - 3s_x) ⋅ C^t+1_0 + s_x ⋅ C^t+1_1

\end{equation}

The right boundary follows the same scheme. We now want to show the equation for the rightmost inlet cell $C_n$ with right boundary value $r$:

\begin{equation}\displaystyle \frac{C_n^t+1 -C_n^t}{Δ t} = α\frac{\frac{r-C^t+1_n}{\frac{Δ x}{2}}- \frac{C^t+1_n-C^t+1_n-1}{Δ x}}{Δ x}

\end{equation}

This expression, once developed, yields:

\begin{align}\displaystyle C_n^t+1 & = C_n^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( 2r - 2C^t+1_n - C^t+1_n + C^t+1_n-1 \right) \nonumber
& = C_0^t + \frac{α ⋅ Δ t}{Δ x^2} ⋅ ≤ft( 2r - 3C^t+1_n + C^t+1_n-1 \right)

\end{align}

Now rearrange terms and substituting with $s_x$ leads to:

\begin{equation}\displaystyle -C^t_n = s_x ⋅ C^t+1_n-1 + (-1 - 3s_x) ⋅ C^t+1_n + (2s_x) ⋅ r

\end{equation}

TODO

• Tridiagonal matrix filling

Diffusion in 2D: the Alternating Direction Implicit scheme

In 2D, the diffusion equation (in absence of source terms) and assuming homogeneous but anisotropic diffusion coefficient $\alpha=(\alpha_x,\alpha_y)$ becomes:

\begin{equation} \displaystyle \frac{\partial C}{\partial t} = \alpha_x \frac{\partial^2 C}{\partial x^2} + \alpha_y\frac{\partial^2 C}{\partial y^2} \end{equation}

2D ADI using BTCS scheme

The Alternating Direction Implicit method consists in splitting the integration of eq. eqn:2d in two half-steps, each of which represents implicitly the derivatives in one direction, and explicitly in the other. Therefore we make use of second derivative operator defined in equation eqn:secondderivative in both $x$ and $y$ direction for half the time step $\Delta t$.

Denoting $i$ the grid cell index along $x$ direction, $j$ the index in $y$ direction and $t$ as the time level, the spatially centered second derivatives can be written as:

\begin{align}\displaystyle \frac{∂^2 C^t_i,j}{∂ x^2} &= \frac{C^t_i-1,j - 2C^t_i,j + C^t_i+1,j}{Δ x^2}
\frac{∂^2 C^t_i,j}{∂ y^2} &= \frac{C^t_i,j-1 - 2C^t_i,j + C^t_i,j+1}{Δ y^2}

\end{align}

The ADI scheme is formally defined by the equations:

\begin{equation} \systeme{ \displaystyle \frac{C^{t+1/2}_{i,j}-C^t_{i,j}}{\Delta t/2} = \displaystyle \alpha_x \frac{\partial^2 C^{t+1/2}_{i,j}}{\partial x^2} + \alpha_y \frac{\partial^2 C^{t}_{i,j}}{\partial y^2}, \displaystyle \frac{C^{t+1}_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = \displaystyle \alpha_x \frac{\partial^2 C^{t+1/2}_{i,j}}{\partial x^2} + \alpha_y \frac{\partial^2 C^{t+1}_{i,j}}{\partial y^2} } \end{equation}

\noindent The first of equations eqn:genADI, which writes implicitly the spatial derivatives in $x$ direction, after bringing the $\Delta t / 2$ terms on the right hand side and substituting $s_x=\frac{\alpha_x \cdot \Delta t}{2 \Delta x^2}$ and $s_y=\frac{\alpha_y\cdot\Delta t}{2\Delta y^2}$ reads:

\begin{equation}\displaystyle C^t+1/2_i,j-C^t_i,j = s_x (C^t+1/2_i-1,j - 2C^t+1/2_i,j + C^t+1/2_i+1,j) + s_y (C^t_i,j-1 - 2C^t_i,j + C^t_i,j+1)

\end{equation}

\noindent Separating the known terms (at time level $t$) on the left hand side and the implicit terms (at time level $t+1/2$) on the right hand side, we get:

\begin{equation}\displaystyle -C^t_i,j - s_y (C^t_i,j-1 - 2C^t_i,j + C^t_i,j+1) = - C^t+1/2_i,j + s_x (C^t+1/2_i-1,j - 2C^t+1/2_i,j + C^t+1/2_i+1,j)

\end{equation}

\noindent Equation eqn:sweepX can be solved with a BTCS scheme since it corresponds to a tridiagonal system of equations, and resolves the $C^{t+1/2}$ at each inner grid cell.

The second of equations eqn:genADI can be treated the same way and yields:

\begin{equation}\displaystyle -C^{t + 1/2}_i,j - s_x (C^{t + 1/2}_i-1,j - 2C^{t + 1/2}_i,j + C^{t + 1/2}_i+1,j) = - C^t+1_i,j + s_y (C^t+1_i,j-1 - 2C^t+1_i,j + C^t+1_i,j+1)

\end{equation}

This scheme only applies to inner cells, or else $\forall i,j \in [1, n-1] \times [1, n-1]$. Following an analogous treatment as for the 1D case, and noting $l_x$ and $l_y$ the constant left boundary values and $r_x$ and $r_y$ the right ones for each direction $x$ and $y$, we can modify equations eqn:sweepX for $i=0, j \in [1, n-1]$

\begin{equation}\displaystyle -C^t_0,j - s_y (C^t_0,j-1 - 2C^t_0,j + C^t_0,j+1) = - C^t+1/2_0,j + s_x (C^t+1/2_1,j - 3C^t+1/2_0,j + 2 l_x)

\end{equation}

\noindent Similarly for $i=n, j \in [1, n-1]$:

\begin{equation}\displaystyle -C^t_n,j - s_y (C^t_n,j-1 - 2C^t_n,j + C^t_n,j+1) = - C^t+1/2_n,j + s_x (C^t+1/2_n-1,j - 3C^t+1/2_n,j + 2 r_x)

\end{equation}

\noindent For $i=j=0$:

\begin{equation}\displaystyle -C^t_0,0 - s_y (C^t_0,1 - 3C^t_0,0 + 2l_y) = - C^t+1/2_0,0 + s_x (C^t+1/2_1,0 - 3C^t+1/2_0,0 + 2 l_x)

\end{equation}

Analogous expressions are readily derived for all possible combinations of $i,j \in 0\times n$. In practice, wherever an index $i$ or $j$ is $0$ or $n$, the centered spatial derivatives in $x$ or $y$ directions must be substituted in relevant parts of the sweeping equations \textbf{in both the implicit or the explicit sides} of equations eqn:sweepX and eqn:sweepY by a term

\begin{equation}\displaystyle s(C_forw - 3C + 2 bc)

\end{equation} \noindent where $bc$ is the boundary condition in the given direction, $s$ is either $s_x$ or $s_y$, and $C_{forw}$ indicates the contiguous cell opposite to the boundary. Alternatively, noting the second derivative operator as $\partial_{dir}^2$, we can write in compact form:

\begin{equation} \systeme{ \displaystyle \partial_x^2 C_{0,j} = 2l_x - 3C_{0,j} + C_{1,j} , \displaystyle \partial_x^2 C_{n,j} = 2r_x - 3C_{n,j} + C_{n-1,j} , \displaystyle \partial_y^2 C_{i,0} = 2l_y - 3C_{i,0} + C_{i,1} , \displaystyle \partial_y^2 C_{i,n} = 2r_y - 3C_{i,n} + C_{i,n-1} } \end{equation}

Heterogeneous diffusion

If the diffusion coefficient $\alpha$ is spatially variable, equation eqn:1 can be rewritten as:

\begin{align} \frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x} \right) \end{align}

Discretization of the equation using chain rule

\noindent From the product rule for derivatives we obtain:

\begin{equation} \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right) = \frac{\partial \alpha}{\partial x} \cdot \frac{\partial C}{\partial x} + \alpha \frac{\partial^2 C }{\partial x^2} \end{equation}

\noindent Using a spatially centred second order finite difference approximation at $x=x_i$ for both $\alpha$ and $C$, we have

\begin{align} \frac{\partial \alpha}{\partial x} \cdot \frac{\partial C}{\partial x} + \alpha \frac{\partial^2 C }{\partial x^2} & \simeq \frac{\alpha_{i+1} - \alpha_{i-1}}{2\Delta x}\cdot\frac{C_{i+1} - C_{i-1}}{2\Delta x} + \alpha_i \frac{C_{i+1} - 2 C_{i} + C_{i-1}}{\Delta x^2} \\ \nonumber & = \frac{1}{\Delta x^2} \frac{\alpha_{i+1} - \alpha_{i-1}}{4}(C_{i+1} - C_{i-1}) + \frac{\alpha_{i}}{\Delta x^2}(C_{i+1}-2C_i+C_{i-1})\\ \nonumber & = \frac{1}{\Delta x^2} \left\{A C_{i+1} -2\alpha_i C_i + AC_{i-1})\right\} \end{align}

\noindent having set \[ A = \frac{\alpha_{i+1}-\alpha_{i-1}}{4} + \alpha_i \]

\noindent In 2D the ADI scheme eqn:genADI with heterogeneous diffusion coefficients can thus be written:

\begin{equation} \systeme{ \displaystyle \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = \displaystyle \frac{\partial}{\partial x} \left( \alpha^x_{i,j} \frac{\partial C^{t+1/2}_{i,j}}{\partial x}\right) + \frac{\partial}{\partial y} \left( \alpha^y_{i,j} \frac{\partial C^{t }_{i,j}}{\partial y}\right), \displaystyle \frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = \displaystyle \frac{\partial}{\partial x} \left( \alpha^x_{i,j} \frac{\partial C^{t+1/2}_{i,j}}{\partial x}\right) + \frac{\partial}{\partial y} \left( \alpha^y_{i,j} \frac{\partial C^{t+1}_{i,j}}{\partial y}\right) } \end{equation}

\noindent We define for compactness $S_x=\frac{\Delta t}{2\Delta x^2}$ and $S_y=\frac{\Delta t}{2\Delta y^2}$ and

\begin{equation} \systeme{ \displaystyle A_{i,j} = \displaystyle \frac{\alpha^x_{i+1,j} -\alpha^x_{i-1,j }}{4} + \alpha^x_{i,j}, \displaystyle B_{i,j} = \displaystyle \frac{\alpha^y_{i, j+1}-\alpha^y_{i ,j-1}}{4} + \alpha^y_{i,j} } \end{equation}

\noindent Plugging eq. (eqn:hetdiff_fd) into the first of equations (eqn:genADI_het) - so called "sweep by x" - and putting all implicit terms (at time level $t+1/2$) on the left hand side we obtain:

\begin{equation}\displaystyle

\begin{split} -S_x A_{i,j} C^{t+1/2}_{i+1,j} + (1 + 2S_x\alpha^x_{i,j})C^{t+1/2}_{i,j} - S_x A_{i,j}C^{t+1/2}_{i-1,j} = \\ S_y B_{i,j} C^{t }_{i,j+1} + (1 - 2S_y\alpha^y_{i,j})C^{t }_{i,j} + S_y B_{i,j}C^{t }_{i,j-1} \end{split}

\end{equation}

\noindent In the same way for the second of eq. eqn:genADI_het we have:

\begin{equation}\displaystyle

\begin{split} -S_y B_{i,j} C^{t+1 }_{i,j+1} + (1 + 2S_y\alpha^y_{i,j})C^{t+1 }_{i,j} - S_y B_{i,j}C^{t+1 }_{i,j-1} = \\ S_x A_{i,j} C^{t+1/2}_{i+1,j} + (1 - 2S_x\alpha^x_{i,j})C^{t+1/2}_{i,j} + S_x A_{i,j}C^{t+1/2}_{i-1,j} \end{split}

\end{equation}

\noindent If the diffusion coefficients are constant, $A_{i,j}=B_{i,j}=\alpha$ and the scheme reverts to the homogeneous case. Problem with this discretization is that the terms in $A_{ij}$ and $B_{ij}$ can be negative depending on the derivative of the diffusion coefficient, resulting in unphysical values for the concentrations.

Direct discretization

As noted in literature (LeVeque and Numerical Recipes) a better way is to discretize directly the physical problem (eq. eqn:hetdiff) at points halfway between grid points:

\begin{align*} \begin{cases} \displaystyle \alpha(x_{i+1/2}) \frac{\partial C }{\partial x}(x_{i+1/2}) & \displaystyle = \alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) \\ \displaystyle \alpha(x_{i-1/2}) \frac{\partial C }{\partial x}(x_{i-1/2}) & \displaystyle = \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \end{cases} \end{align*}

\noindent A further differentiation gives us the spatially centered approximation of $\frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)$:

\begin{equation} \begin{aligned} \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)(x_i) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\ &\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+\alpha_{i-1/2}) C_{i} + \alpha_{i-1/2}C_{i-1}\right] \end{aligned} \end{equation}

\noindent The ADI scheme with this approach becomes:

\begin{equation} \left\{ \begin{aligned} \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = & \frac{1}{\Delta x^2} \left[ \alpha_{i+1/2,j} C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^{t+1/2}_{i,j} + \alpha_{i-1/2,j} C^{t+1/2}_{i-1,j}\right] + \\ & \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^t_{i,j} + \alpha_{i,j-1/2}C^{t}_{i,j-1}\right]\\ \frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = & \frac{1}{\Delta y^2} \left[ \alpha_{i+1/2,j}C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^t_{i,j} + \alpha_{i-1/2,j}C^{t+1/2}_{i-1,j}\right] + \\ & \frac{1}{\Delta x^2} \left[ \alpha_{i,j+1/2}C^{t+1}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^{t+1}_{i,j} + \alpha_{i,j-1/2}C^{t+1}_{i,j-1}\right] \end{aligned} \right. \end{equation}

\noindent Doing the usual algebra and separating implicit from explicit terms, the two sweeps become:

\begin{equation} \left\{ \begin{aligned} -S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 + S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} - S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} = \quad & \\ S_y \alpha^y_{i,j+1/2} C^{t }_{i,j+1} + (1 - S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t }_{i,j} + S_y \alpha^y_{i,j-1/2} C^{t }_{i,j-1} & \\[1em] -S_y \alpha^y_{i,j+1/2} C^{t+1 }_{i,j+1} + (1 + S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t+1 }_{i,j} - S_y \alpha^y_{i,j-1/2} C^{t+1 }_{i,j-1} = \qquad & \\ S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 - S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} + S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} \end{aligned} \right. \end{equation}

\bigskip

\noindent The "interblock" diffusion coefficients $\alpha_{i+1/2,j}$ can be arithmetic mean:

\[ \displaystyle \alpha_{i+1/2, j} = \displaystyle \frac{\alpha_{i+1, j} + \alpha_{i, j}}{2} \]

\noindent or the harmonic mean:

\[ \displaystyle \alpha_{i+1/2, j} = \displaystyle \frac{2}{\frac{1}{\alpha_{i+1, j}} + \frac{1}{\alpha_{i, j}}} \]

\pagebreak

Explicit scheme for 2D heterogeneous diffusion

A classical explicit FTCS scheme (forward in time, central in space) for 2D heterogeneous diffusion can be expressed simply leveraging the discretization of equation eqn:CS_het:

\begin{equation} \begin{aligned} \frac{C_{i,j}^{t+1} - C_{i,j}^{t}}{\Delta t} = & \frac{1}{\Delta x^2} \left[ \alpha^x_{i+1/2, j}C^t_{i+1, j} - (\alpha^x_{i+1/2, j} + \alpha^x_{i-1/2, j}) C^t_{i,j} + \alpha^x_{i-1/2,j}C^t_{i-1,j}\right] + \\ & \frac{1}{\Delta y^2} \left[ \alpha^y_{i, j+1/2}C^t_{i, j+1} - (\alpha^y_{i, j+1/2} + \alpha^y_{i, j-1/2}) C^t_{i,j} + \alpha^y_{i,j-1/2}C^t_{i,j-1}\right] \end{aligned} \end{equation}

\noindent where in the RHS only the known concentrations at time $t$ appear. Rearranging the terms, we get:

\begin{equation} \begin{aligned} C_{i,j}^{t+1} = & C^t_{i,j} +\\ & \frac{\Delta t}{\Delta x^2} \left[ \alpha^x_{i+1/2, j}C^t_{i+1, j} - (\alpha^x_{i+1/2, j} + \alpha^x_{i-1/2, j}) C^t_{i,j} + \alpha^x_{i-1/2,j}C^t_{i-1,j}\right] + \\ & \frac{\Delta t}{\Delta y^2} \left[ \alpha^y_{i, j+1/2}C^t_{i, j+1} - (\alpha^y_{i, j+1/2} + \alpha^y_{i, j-1/2}) C^t_{i,j} + \alpha^y_{i,j-1/2}C^t_{i,j-1}\right] \end{aligned} \end{equation}

The Courant-Friedrichs-Lewy stability criterion (cfr Lee, 2017) for this scheme reads:

\begin{equation} \Delta t \leq \frac{1}{2 \max(\alpha_{i,j})} \cdot \frac{1}{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}} \end{equation}

Note that other derivations for the CFL condition are found in literature. For example, the sources cited by Wikipedia solution give:

\begin{equation} \displaystyle \Delta t\leq {\frac {1}{4 \max(\alpha) \left({\frac {1}{\Delta x^{2}}}+{\frac {1}{\Delta y^{2}}}\right)}} \end{equation}

We can produce a more restrictive condition than equation eqn:CFL2DFTCS_Lee by considering the min of the $\Delta x$ and $\Delta y$:

\begin{equation} \Delta t \leq \frac{\min(\Delta x, \Delta y)^2}{4 \max(\alpha_{i,j})} \end{equation}

In practice for the implementation it is advantageous to specify an optional parameter $C$, $C \in [0, 1]$ so that the user can restrict the "inner time stepping":

\begin{equation} \Delta t \leq C \cdot \frac{\min(\Delta x, \Delta y)^2}{4 \max(\alpha_{i,j})} \end{equation}

Boundary conditions

In analogy to the treatment of the 1D homogeneous FTCS scheme (cfr section 1), we need to differentiate the domain boundaries ($i=0$ and $i=n_x$; the same applies to $j$ of course) accounting for the discrepancy in the discretization.

For the zero-th (left) cell, whose center is at $x=dx/2$, we can evaluate the left gradient with the left boundary using such distance, calling $l$ the numerical value of a constant boundary condition, equation eqn:CS_het becomes:

\begin{equation} \begin{aligned} \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)(x_0) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i} \left( \frac{C_{i} - l }{\frac{\Delta x}{2}} \right) \right]\\ &\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+ 2\alpha_i) C_{i} + 2 \alpha_{i}\cdot l\right] \end{aligned} \end{equation}

\noindent Similarly, for $i=n_x$,

\begin{equation} \begin{aligned} \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)(x_n) & \simeq \frac{1}{\Delta x}\left[\alpha_{i} \left( \frac{r -C_{i}}{\frac{\Delta x}{2}} \right) - \alpha_{i-1/2} \left( \frac{C_{i} - C_{i-1}} {\Delta x} \right) \right]\\ &\displaystyle =\frac{1}{\Delta x^2} \left[ 2 \alpha_{i} r - (\alpha_{i+1/2}+ 2\alpha_i) C_{i} + \alpha_{i-1/2}\cdot C_{i-1}\right] \end{aligned} \end{equation}

If on the right boundary we have closed or Neumann condition, the left derivative becomes zero and we are left with:

\begin{equation} \begin{aligned} \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right)(x_n) & \simeq \frac{1}{\Delta x}\left[\cancel{\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right)} - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\ &\displaystyle =\frac{\alpha_{i-1/2}}{\Delta x^2} (C_{i-1} - C_i) \end{aligned} \end{equation}