two-bulb-imp.tex

\documentclass[11]{Report}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{float}
\makeatletter
\newcommand*\bigcdot{\mathpalette\bigcdot@{.4}}
\newcommand*\bigcdot@[2]{\mathbin{\vcenter{\hbox{\scalebox{#2}{$\m@th#1\bullet$}}}}}
\makeatother
\usepackage{siunitx} % Required for alignment
\usepackage{titlesec}

\sisetup{
  round-mode          = places, % Rounds numbers
  round-precision     = 2, % to 2 places
}

\begin{document}
\author{Derek W. Harrison}
\title{Application of the Finite Volume Method to the Maxwell-Stefan model for diffusion}
\maketitle


\section*{Introduction}
The Maxwell-Stefan equations for multi-component diffusion are solved using the finite volume method. Time discretization is fully implicit. To validate the method the Maxwell-Stefan equations are used to model the three-component twin-bulb experiment and the results are compared with results obtained from an explicit numerical scheme as well as experimental observations.

\section*{Twin-bulb model}
The twin-bulb experiment [\ref{eqn:my_ref_2}] consists of two small compartments (bulbs) connected by a tube through which the components can diffuse. The compartments are filled with $H_2$, $N_2$ and $CO_2$. Diffusion through the tube can be modeled using the Maxwell-Stefan equations [\ref{eqn:my_ref_1}]:
\begin{equation}
\label{eqn:eqn_2}
-\left( \frac{\partial \ln{\gamma_i}}{\partial \ln{x_i}} + 1 \right) \nabla x_i = \sum_{j \neq i} \frac{x_j \textbf{J}_i - x_i \textbf{J}_j}{c_t D_{ij}}
\end{equation}
For ideal systems the activity coefficient $\gamma_i$ of component $i$ is equal to unity. The left side of (\ref{eqn:eqn_2}) then simplifies, resulting in:
\begin{equation}
\label{eqn:eqn_3}
- \nabla x_i = \sum_{j \neq i} \frac{x_j \textbf{J}_i - x_i \textbf{J}_j}{c_t D_{ij}}
\end{equation}
From a mass balance follows that the change in local composition at any given time is:
\begin{equation}
\label{eqn:eqn_4}
c_t \frac{\partial x_i}{\partial t} = - \nabla \bigcdot \textbf{J}_i
\end{equation}
Diffusion occurs at constant pressure. To preserve the total concentration the fluxes of the different components sum up to zero:
\begin{equation}
\label{eqn:eqn_5}
\sum_{i} \textbf{J}_i = 0
\end{equation}
\section*{Method}
To compute the composition in the bulbs the model equations (\ref{eqn:eqn_3}) - (\ref{eqn:eqn_5}) are solved using the finite volume method. Time discretization is fully implicit, achieved by eliminating the flux-components from the model equations. Central differencing is used for the diffusion terms.

An explicit time discretization scheme is used as a reference to validate and obtain a measure for the stability of the implicit approach.

The mole fractions of $H_2$, $N_2$ and $CO_2$ in the first compartment (bulb $1$) are initially 0.501, 0.499 and 0.0, respectively. In the second compartment (bulb $2$) the mole fractions of $H_2$, $N_2$ and $CO_2$ are initially 0.0, 0.501 and 0.499, respectively. The diffusivities are $D_{12} = 8.33e-5$  $(m^2/s)$, $D_{13} = 6.8e-5$ $(m^2/s)$ and $D_{23} = 1.68e-5$ $(m^2/s)$. The volumes of the compartments are $5e-4$ $(m^3)$ and the tube connecting the compartments has a length of $1e-2$ $(m)$ and a diameter of $2e-3$ $(m)$. 

\section*{Results} 
The results obtained from solving the model equations are shown in figure \ref{fig:fig_1}. The results are in excellent agreement with the results obtained from the explicit method (results are not shown) and agree well with experimental observations [\ref{eqn:my_ref_3}].

\begin{figure}
\includegraphics[width=\linewidth]{graph_of_results.png}
\caption{The mole fraction as a function of time (h).}
\label{fig:fig_1}
\end{figure}

\section*{Discussion}
Explicit time discretization was only stable when small timesteps were used, whereas implicit discretization was stable even when large timesteps were used. Analysis of the coefficients of the linear systems obtained from implicit discretization of the model equations showed the coefficients satisfy the Scarborough criterion and have the same sign. Therefore, the coefficients meet the requirements for convergence and boundedness. The implicit scheme also had the lowest running time, suggesting that the efficiency gained by using the implicit scheme weighs more heavily than the cost of solving a linear system at each timestep.

\section*{Conclusion}
The Maxwell-Stefan equations were solved using the finite volume method and the results were found to coincide with those of the explicit method and agree well with experimental observations, thereby validating the method. Regarding the discretization strategy, it was found that implicit time discretization is significantly more stable than explicit discretization and led to lower running times. 

\section*{Appendix}
Here the one-dimensional case of the Maxwell-Stefan equations modeling the three-component twin-bulb experiment is elaborated on. The one-dimensional case of (\ref{eqn:eqn_3}) for component $1$ is:
\begin{equation}
\label{eqn:eqn_6}
-c_t \frac{\partial x_1}{\partial z} = \frac{x_2 J_1 - x_1 J_2}{D_{12}} + \frac{x_3 J_1 - x_1 J_3}{D_{13}}
\end{equation}
The equation for component $2$ is:
\begin{equation}
\label{eqn:eqn_7}
-c_t \frac{\partial x_2}{\partial z} = \frac{x_1 J_2 - x_2 J_1}{D_{12}} + \frac{x_3 J_2 - x_2 J_3}{D_{23}}
\end{equation}
The change in local mole fraction of component $1$ is:
\begin{equation}
\label{eqn:eqn_8}
c_t \frac{\partial x_1}{\partial t} = - \frac{\partial J_1}{\partial z}
\end{equation}
The change in local mole fraction of component $2$ is:
\begin{equation}
\label{eqn:eqn_9}
c_t \frac{\partial x_2}{\partial t} = - \frac{\partial J_2}{\partial z}
\end{equation}
To facilitate the elimination of the fluxes from the equations above equation (\ref{eqn:eqn_6}) and (\ref{eqn:eqn_7}) are rewritten:
\begin{equation}
\label{eqn:eqn_10}
c_t \frac{\partial x_1}{\partial z} = a_1 J_1 + a_2 J_2
\end{equation}
\begin{equation}
\label{eqn:eqn_11}
c_t \frac{\partial x_2}{\partial z} = b_1 J_1 + b_2 J_2
\end{equation}
With $a_1$, $a_2$, $b_1$ and $b_2$ given by:
\begin{equation}
\label{eqn:eqn_12}
a_1 = \left( \frac{1}{D_{12}} - \frac{1}{D_{13}} \right) x_2 + \frac{1}{D_{13}}
\end{equation}
\begin{equation}
\label{eqn:eqn_13}
a_2 = \left( \frac{1}{D_{13}} - \frac{1}{D_{12}} \right) x_1
\end{equation}
\begin{equation}
\label{eqn:eqn_14}
b_1 = \left( \frac{1}{D_{23}} - \frac{1}{D_{12}} \right) x_2
\end{equation}
\begin{equation}
\label{eqn:eqn_15}
b_2 = \left( \frac{1}{D_{12}} - \frac{1}{D_{23}} \right) x_1 + \frac{1}{D_{23}}
\end{equation}
The fluxes can now be written in terms of the composition gradients:
\begin{equation}
\label{eqn:eqn_16}
J_1 = \beta_1 \frac{\partial x_1}{\partial z} + \beta_2 \frac{\partial x_2}{\partial z}
\end{equation}
\begin{equation}
\label{eqn:eqn_17}
J_2 = \alpha_1 \frac{\partial x_1}{\partial z} + \alpha_2 \frac{\partial x_2}{\partial z}
\end{equation}
With $\beta_1$, $\beta_2$, $\alpha_1$ and $\alpha_2$ given by:
\begin{equation}
\label{eqn:eqn_18}
\beta_1 = -\frac{c_t}{a_1} - \frac{a_2 \alpha_1}{a_1}
\end{equation}
\begin{equation}
\label{eqn:eqn_19}
\beta_2 = - \frac{a_2 \alpha_2}{a_1}
\end{equation}
\begin{equation}
\label{eqn:eqn_20}
\alpha_1 = - \frac{c_t}{a_2 - \frac{a_1 b_2}{b_1}}
\end{equation}
\begin{equation}
\label{eqn:eqn_21}
\alpha_2 = \frac{a_1 c_t}{a_2 b_1 - a_1 b_2}
\end{equation}
Equations (\ref{eqn:eqn_8}), (\ref{eqn:eqn_9}), (\ref{eqn:eqn_16}) and (\ref{eqn:eqn_17}) are the set of equations modeling the twin-bulb experiment.
\subsubsection*{Nomenclature}
\begin{table}[H]
    \begin{tabular}{c l}  
      $a_i$ & Proportionality factor ($ m^{-2} \cdot s $) \\     
      $b_i$ & Proportionality factor ($ m^{-2} \cdot s $) \\     
      $c_t$ & Concentration ($mol \cdot m^{-3}$) \\
      $D_{ij}$ & Diffusivity ($m^2 \cdot s^{-1}$) \\
      $\textbf{J}_i$ & Flux vector ($mol \cdot m^{-2} \cdot s^{-1}$) \\
      $J_i$ & Flux component ($mol \cdot m^{-2} \cdot s^{-1}$) \\         
      $x_i$ & Mole fraction (-)\\
      $z$ & Axial coordinate ($m$) \\
    \end{tabular}
\end{table}
\subsubsection*{Greek}
\begin{table}[H]
    \begin{tabular}{c l} 
      $\alpha_i$ & Proportionality factor ($mol \cdot m^{-1} \cdot s^{-1}$) \\
      $\beta_i$ & Proportionality factor ($mol \cdot m^{-1} \cdot s^{-1}$) \\
      $\gamma_i$ & Activity coefficient (-) \\
    \end{tabular}
\end{table}
\subsubsection*{Subscripts}
\begin{table}[H]
    \begin{tabular}{c l} 
      $1$ & Index for $H_2$ (-)\\
      $2$ & Index for $N_2$ (-)\\
      $3$ & Index for $CO_2$ (-)\\
      $i$ & Component index (-)\\
      $j$ & Component index (-) \\
      $t$ & Total (-) \\
    \end{tabular}
\end{table}

%\titleformat*{\section}{\normalfont}
\renewcommand{\bibname}{References}
\begin{thebibliography}{5}
\bibitem{Duncan}
\label{eqn:my_ref_2}
Duncan, J.B., Toor, H.L. \textit{AIChE J.}, 1962, \textbf{8}, 38–41.
\bibitem{Taylor} 
\label{eqn:my_ref_1}
Taylor, R., Krishna, R. \textit{Multicomponent Mass Transfer}. New York: Wiley, 1993.

\bibitem{Taylor} 
\label{eqn:my_ref_3}
Krishna, R. Uphill diffusion in multicomponent mixtures. \textit{Chem. Soc. Rev.}, 2015, \textbf{44}, 2812.

\end{thebibliography}

\end{document}