Merge pull request #37 from thermotools/soret_error

Fix stupid Soret coefficient bug

cpp/KineticGas_properties.cpp

@@ -175,13 +175,14 @@ Eigen::MatrixXd KineticGas::interdiffusion(double T, double Vm, const vector1d&
     if (dependent_idx < 0 && frame_of_reference == FrameOfReference::solvent) dependent_idx = solvent_idx;
     while (dependent_idx < 0) dependent_idx += Ncomps;
     if (frame_of_reference == FrameOfReference::zarate_x){
+        if (!is_idealgas) throw std::runtime_error("Zarate relation for diffusion only implemented for ideal gas!");
         Eigen::MatrixXd D = interdiffusion(T, Vm, x, N, FrameOfReference::CoN, dependent_idx, dependent_idx, true);
         return D;
     }
     else if (frame_of_reference == FrameOfReference::zarate){
         Eigen::MatrixXd X = get_zarate_X_matr(x, dependent_idx);
         Eigen::MatrixXd Dx = interdiffusion(T, Vm, x, N, FrameOfReference::zarate_x, dependent_idx);
-        return X * Dx * X.inverse();
+        return X.inverse() * Dx * X;
     }
     else if (frame_of_reference == FrameOfReference::zarate_w){
         Eigen::MatrixXd W = get_zarate_W_matr(x, dependent_idx);
@@ -257,6 +258,7 @@ Eigen::VectorXd KineticGas::thermal_diffusion_coeff(double T, double Vm, const v
     DT /= AVOGADRO;
 
     if (use_zarate){
+        if (!is_idealgas) throw std::runtime_error("Zarate relation for thermal diffusion only implemented for ideal gas!");
         Eigen::VectorXd DT_indep(Ncomps - 1);
         Eigen::MatrixXd D_indep = interdiffusion(T, Vm, x, N, frame_of_reference, dependent_idx, solvent_idx);
 
@@ -319,6 +321,7 @@ Eigen::MatrixXd KineticGas::thermal_diffusion_factor(double T, double Vm, const
 }
 
 Eigen::VectorXd KineticGas::soret_coefficient(double T, double Vm, const std::vector<double>& x, int N, int dependent_idx){
+    if (!is_idealgas) throw std::runtime_error("Zarate relation for Soret coefficient only implemented for ideal gas!");
     Eigen::MatrixXd D = interdiffusion(T, Vm, x, N, FrameOfReference::zarate, dependent_idx, dependent_idx, false);
     Eigen::VectorXd DT = thermal_diffusion_coeff(T, Vm, x, N, FrameOfReference::zarate, dependent_idx);
     return D.partialPivLu().solve(DT);

docs/memo/binary_limit/binary_limit/diffusion.tex

@@ -1,7 +1,7 @@
 \section{Diffusion}
 This section discusses the binary limit of the diffusion coefficients in a ternary system. The purpose of the section is to give insight into how the force-flux in a binary system relate to those in a ternary, how to interpret the diffusion coefficients in a ternary system, and raise awareness regarding potential pitfalls when attempting to model a ternary system as a pseudo-binary system.
 
-\subsection{Notation}
+\subsection{Notation}\label{sec:diffusion_notation}
 
 In the following text, $D_{ii}^{(b)}$ is used to denote the diffusion coefficient of a binary system, where component $i$ is the independent component. $D_{ii}^{(k)}$ denotes the diagonal elements of the Fick diffusion matrix in a ternary system, where component $k$ is the dependent component. $D_{ij}$ denotes the non-diagonal elements of the Fick diffusion matrix in a ternary system, where components $i$ and $j$ are the independent components.
 

docs/memo/binary_limit/binary_limit/main.tex

@@ -20,6 +20,7 @@ \section{Introduction}
 \input{diffusion}
 \clearpage
 \input{thermal_diffusion}
+\input{soret_coefficient}
 
 \clearpage
 

docs/memo/binary_limit/binary_limit/soret_coefficient.tex

@@ -0,0 +1,66 @@
+\section{Soret coefficients}
+
+The Soret coefficient is a measure of the steady state separation in a mixture induced by a temperature gradient. Ortiz de Zárate \cite{ortiz2019definition} defines the Soret coefficients of a multicomponent mixture ($s$ components) through
+\begin{equation}
+    \begin{bmatrix}
+        x_1 (1 - x1) & x_1 x_2 & \hdots & x_1 x_{s - 1} \\
+        x_2 x_1 & x_2 (1 - x_2) & \hdots & x_2 x_{s - 1} \\
+        \vdots & & \ddots & \vdots \\
+        x_{s - 1} x_1 & & & x_{s - 1}(1 - x_{s - 1}) 
+    \end{bmatrix}
+    \begin{pmatrix}
+        S_{T, 1} \\
+        S_{T, 2} \\
+        \vdots \\
+        S_{T, s - 1}
+    \end{pmatrix}
+    = 
+    - \frac{1}{\nabla T}
+    \begin{pmatrix}
+        \nabla x_1 \\
+        \nabla x_2 \\
+        \vdots \\
+        \nabla x_{s - 1}
+    \end{pmatrix},
+\end{equation}
+or more compactly,
+\begin{equation}
+    \Mat{X} \Vec{S}_T = - \frac{\nabla \Vec{x}}{\nabla T}.
+\end{equation}
+Similarly to the thermal diffusion coefficients, this definition carries the advantage that 
+\begin{equation}
+    \Mat{X} \Vec{S}_T = - \frac{\nabla \Vec{x}}{\nabla T} \quad \iff \quad \Mat{W} \Vec{S}_T = - \frac{\nabla \Vec{w}}{\nabla T},
+\end{equation}
+such that mole- and mass fractions can be used interchangably with the same Soret coefficients.
+
+In the state with vanishing mass fluxes ($\Vec{J} = \Vec{0}$). From the condition of vanishing mass fluxes, we find that we can compute the Soret coefficients as\footnote{See the memo on definitions of the diffusion coefficient for notes on $\Mat{D}^{(z)}$.}
+\begin{equation}
+    \begin{split}
+        \Vec{J} = - c \left( \Mat{X} \Vec{D}_{T}^{(z)} \nabla T + \Mat{X} \Mat{D}^{(z)} \Mat{X}^{-1} \nabla \Vec{x} \right) &= \Vec{0} \\
+        - \frac{\nabla \Vec{x}}{\nabla T} &= \Mat{X} \left(\Mat{D}^{(z)}\right)^{-1}  \Vec{D}_{T}^{(z)} \\
+        \Vec{S}_T &= \left(\Mat{D}^{(z)}\right)^{-1}  \Vec{D}_{T}^{(z)}.
+    \end{split}
+\end{equation}
+
+For a binary system (1, 2), this definition thus yields
+\begin{equation}
+    S_{T, 1}^{(b)} = \frac{D_{T, 1}^{(z, b)}}{D_{11}^{(z, b)}}.
+\end{equation}
+
+From the preceding relations it is possible to show that when using this definition of the Soret coefficient \cite{ortiz2019definition},
+\begin{equation}
+    \begin{split}
+        \lim_{x_1 \to 0} S_{T, 2} = S_{T, 2}^{(b, 3)} \\
+        \lim_{x_2 \to 0} S_{T, 1} = S_{T, 1}^{(b, 3)} \\
+        \lim_{x_3 \to 0} S_{T, 1} - S_{T, 2} = S_{T, 1}^{(b, 2)},
+    \end{split}
+    \label{eq:soret_binary_limit}
+\end{equation}
+where $S_{T, i}^{b, j}$ denotes the Soret coefficient of component $i$ in a binary mixture with $j$, with species $j$ being the dependent species. Figure \ref{fig:soret_binary_limit} shows how this convergence behaviour is obeyed.
+
+\begin{figure}
+    \centering
+    \includegraphics[width=.85\textwidth]{soret_limit.pdf}
+    \caption{The convergence of the ternary Soret coefficient to the corresponding binaries as indicated in Eq. \eqref{eq:soret_binary_limit}.}
+    \label{fig:soret_binary_limit}
+\end{figure}

docs/memo/binary_limit/binary_limit/thermal_diffusion.tex

@@ -10,6 +10,8 @@ \subsection{Notation}
 
 The notation $D_{T,i}^{(z,tj)}$ denotes the thermal diffusion coefficient of species $i$ in a ternary mixture with species $j$ taken to be the dependent species, as defined in ref. \cite{ortiz2019definition} while $D_{T,i}^{(z,bj)}$ denotes the thermal diffusion coefficient of species $i$ in a binary mixture, with species $j$ taken to be the dependent species.
 
+Boldface roman font ($\Vec{v}$) is used to indicate vectors, and slanted underlined boldface figures ($\Mat{D}$) are used to indicate matrices. The notation $\nabla \Vec{v}$, should be understood not as a dot product, but as a vector of gradients (i.e. $\nabla \Vec{v} \equiv (\nabla v_1, \nabla v_2, ...)^\top$, in which case the equations should be understood to hold componentwise for these gradients.
+
 \subsection{Definitions of the thermal diffusion coefficient}
 
 The default definition of the thermal diffusion coefficient in a multicomponent mixture in the KineticGas package is
@@ -62,7 +64,7 @@ \subsection{Definitions of the thermal diffusion coefficient}
 where $D_{T,1}^{(z,t)}$ is the thermal diffusion coefficient of species 1 in a ternary mixture, and $D_{T,1}^{(z,b3)}$ is the thermal diffusion coefficient of species 1 in a binary mixture with species 3, both as defined by Eq. \eqref{eq:zarate_def}.
 that is: For a ternary system (1, 2, 3), when the mole fraction of species 2 tends to zero, the thermal diffusion coefficient of species 1 approaches the thermal diffusion coefficient of species 1 in the binary mixture (1, 3). This relation follows from an argument analogous to that in section \ref{sec:diff_indep}.
 
-We can relate the thermal diffusion coefficients $\Vec{D}_T^{(n)}$ to the $\Vec{D}_T^{(z)}$, by rewriting Eq. \eqref{eq:DTn_def} as 
+For an ideal gas can relate the thermal diffusion coefficients $\Vec{D}_T^{(n)}$ to the $\Vec{D}_T^{(z)}$, by rewriting Eq. \eqref{eq:DTn_def} as 
 \begin{equation}
     \begin{split}
         \Vec{J}^{(n, n)} &= - \Mat{D}^{(n)} ( c \nabla \Vec{x} + \Vec{x} \nabla c) + \Vec{D}_T^{(n)} \nabla \ln T \\
@@ -72,11 +74,31 @@ \subsection{Definitions of the thermal diffusion coefficient}
 \end{equation}
 such that $\Vec{D}_T^{(z)}$ is given by the solution to
 \begin{equation}
-    - c \Mat{X} \Vec{D}_T^{(z)} = \frac{1}{T} \left( \Vec{D}_T^{(n)} + c \Mat{D}^{(n)}\Vec{x}\right).
+    - c \Mat{X} \Vec{D}_T^{(z)} = \frac{1}{T} \left( \Vec{D}_T^{(n)} + c \Mat{D}^{(n)}\Vec{x}\right), \quad \text{Ideal gas}
     \label{eq:DTz_DTn_relate}
 \end{equation}
+where we have made use of $\nabla c = c \nabla \ln T$ for an ideal gas. For a non-ideal gas the corresponding expression is
+\begin{equation}
+    \begin{split}
+        \Vec{J}^{(n, n)} &= - \Mat{D}^{(n)} \left[\ppder{\vec{c}}{T}_{\vec{x}, p} \nabla T + \Mat{\Gamma}_c \nabla \vec{x}\right] + \frac{1}{T}\Vec{D}_T^{(n)} \nabla T \\
+        &= - \Mat{D}^{(n)} \Mat{\Gamma}_c \nabla \vec{x} + \left[\frac{1}{T}\Vec{D}_T^{(n)} - \Mat{D}^{(n)} \ppder{\vec{c}}{T}_{\vec{x}, p}\right] \nabla T,
+    \end{split}
+\end{equation}
+where 
+\begin{equation}
+    \left[\Mat{\Gamma}_c\right]_{ij} = \ppder{c_i}{x_j}_{T, p, x_{k \neq j}},
+\end{equation}
+and it can be useful to note that 
+\begin{equation}
+    \begin{split}
+        c_i &= \frac{x_i}{v} \implies \\
+        \ppder{c_i}{T}_{p, \vec{x}} &= - \frac{x_i}{v^2} \ppder{v}{T}_{p, \vec{x}}, \\
+        \ppder{c_i}{x_j}_{T, p, x_{k \neq j}} &= \frac{\delta_{ij}}{v} - \frac{c_i v_j}{v},
+    \end{split}
+\end{equation}
+and to keep in mind that these matrices and vectors are in $\mathbb{R}^{s - 1}$, because the dependent component is excluded. \textbf{Note:} At the time of writing, the expressions for non-ideal gas have not been implemented.
 
-By the same manipulation, and noting that we can write $\nabla \vec{x} = \Mat{T}^{x \leftmapsto w} \nabla \Vec{w}$, we can rewrite equation $\eqref{eq:Dtm_def}$ as 
+Returning to an ideal gas, using $\nabla c = c \nabla \ln T$, and noting that we can write $\nabla \vec{x} = \Mat{T}^{x \leftmapsto w} \nabla \Vec{w}$, we can rewrite equation $\eqref{eq:Dtm_def}$ as 
 \begin{equation}
     \begin{split}
         \Vec{J}^{(n, m)} &= - \Mat{D}^{(m)} \nabla \Vec{c} + \Vec{D}_T^{(m)} \nabla \ln T \\
@@ -96,24 +118,30 @@ \subsection{Definitions of the thermal diffusion coefficient}
     \begin{split}
         \frac{1}{c} \Mat{X}^{-1}\left( \Vec{D}_T^{(n)} + c \Mat{D}^{(n)}\Vec{x}\right) &= \frac{1}{\rho} \Mat{W}^{-1} \diag(\Vec{M}) \left( \Vec{D}_T^{(m)} + c \Mat{D}^{(m)}\Vec{x}\right) \\
         &= \frac{1}{\rho} \Mat{W}^{-1} \diag(\Vec{M}) \Mat{\Psi}^{m \leftmapsto n} \left( \Vec{D}_T^{(n)} + c \Mat{D}^{(n)}\Vec{x}\right) \\
-        \frac{1}{c} \Mat{X}^{-1} \Mat{\Psi}^{n \leftmapsto m} \left( \Vec{D}_T^{(m)} + c \Mat{D}^{(m)}\Vec{x}\right) &= \frac{1}{\rho} \Mat{W}^{-1} \diag(\Vec{M}) \left( \Vec{D}_T^{(m)} + c \Mat{D}^{(m)}\Vec{x}\right)
+        \Mat{I} &= \frac{c}{\rho} \Mat{X} \Mat{W}^{-1} \diag(\Vec{M}) \Mat{\Psi}^{m \leftmapsto n}, \quad \text{and} \\
+        \frac{1}{c} \Mat{X}^{-1} \Mat{\Psi}^{n \leftmapsto m} \left( \Vec{D}_T^{(m)} + c \Mat{D}^{(m)}\Vec{x}\right) &= \frac{1}{\rho} \Mat{W}^{-1} \diag(\Vec{M}) \left( \Vec{D}_T^{(m)} + c \Mat{D}^{(m)}\Vec{x}\right) \\
+        \frac{\rho}{c} \diag(\Vec{M}^{-1}) \Mat{W} \Mat{X}^{-1} \Mat{\Psi}^{n \leftmapsto m} &= \Mat{I}
     \end{split}
 \end{equation}
 
-As a second side-effect of this procedure, we can relate the diffusion matrices, as defined by Eqs. \eqref{eq:zarate_def} to $\Mat{D}^{(m)}$ and $\Mat{D}^{(n)}$ as
+As a second side-effect of this procedure, we can relate the diffusion matrices (for an ideal gas), as defined by Eqs. \eqref{eq:zarate_def} to $\Mat{D}^{(m)}$ and $\Mat{D}^{(n)}$ as
 \begin{equation}
     \begin{split}
         \Mat{D}^{(x)} &= \Mat{D}^{(n)} \\
         \Mat{D}^{(w)} &= \frac{c}{\rho} \diag(\Vec{M}) \Mat{D}^{(m)} \Mat{T}^{x \leftmapsto w}
     \end{split}
+    \label{eq:zarate_D_impl}
 \end{equation}
 and remark that Ortiz de Zárate gives the relation
 \begin{equation}
-    \Mat{W}^{-1} \Mat{D}^{(w)} \Mat{W} = \Mat{X}^{-1} \Mat{D}^{(x)} \Mat{X},
+    \Mat{W}^{-1} \Mat{D}^{(w)} \Mat{W} = \Mat{X}^{-1} \Mat{D}^{(x)} \Mat{X} \equiv \Mat{D}^{(z)},
+    \label{eq:zarate_D_equiv}
 \end{equation}
 which can serve as a second consistency check.
 
-All the consistency checks mentioned above have been carried out numerically and been found to pass.
+All the consistency checks mentioned above have been carried out numerically and been found to pass. 
+
+In light of Eq. \eqref{eq:zarate_D_impl}, the simplest way to implement these diffusion matrices is likely to implement $\Mat{D}^{(n)}$, and to compute $\Mat{D}^{(w)}$ and $\Mat{D}^{(z)}$ using the relations in Eq. \eqref{eq:zarate_D_equiv}. This is the approach taken in the KineticGas package.
 
 \subsubsection{The binary case}
 
@@ -147,4 +175,4 @@ \subsection{The binary limit}
     \includegraphics[width=.85\textwidth]{ternary_DT_large.pdf}
     \caption{The same thermal diffusion coefficients as in Fig. \ref{fig:ternary_DT}, with deviations and for a larger composition span.}
     \label{fig:ternary_DT_large}
-\end{figure}
+\end{figure}

tests/test_binary_limit.py

@@ -117,5 +117,63 @@ def test_th_diffusion(model, N, silent=True):
 
     assert all((diff / DT_b) < 1e-3)
 
+@pytest.mark.parametrize('model', models)
+@pytest.mark.parametrize('N', [2, 3])
+def test_soret(model, N, silent=True):
+    T = 300
+    Vm = 1 / 10
+    x_bin = [0.3, 0.7]
+    bin_12 = model('H2O,O2', is_idealgas=True)
+    ternary = model('H2O,O2,N2', is_idealgas=True)
+
+    x_bin = [0.6, 0.4]
+    ST_1b2 = bin_12.soret_coefficient(T, Vm, x_bin, N=N, use_zarate=True)[0]
+    ST_1t, ST_2t = ternary.soret_coefficient(T, Vm, get_ternary_x(x_bin, 0.5), N=N, use_zarate=True)
+    diff = abs((ST_1t - ST_2t) - ST_1b2)
+    x3_lst = np.logspace(-1, -5, 10)
+
+    print(ST_1t - ST_2t, np.diff([ST_1t, ST_2t]), ST_1b2)
+
+    if silent is False:
+        plt.plot(x3_lst, [- np.diff(ternary.soret_coefficient(T, Vm, get_ternary_x(x_bin, x3), N=N, use_zarate=True)) for x3 in x3_lst])
+        plt.plot(x3_lst, ST_1b2 * np.ones_like(x3_lst))
+        plt.xscale('log')
+        plt.show()
+
+    for x3 in x3_lst:
+        ST_1t, ST_2t = ternary.soret_coefficient(T, Vm, get_ternary_x(x_bin, x3), N=N, use_zarate=True)
+        new_diff = abs((ST_1t - ST_2t) - ST_1b2)
+        assert new_diff < diff
+        diff = new_diff
+
+    bin_23 = model('O2,N2', is_idealgas=True)
+    def get_ternary_x1(x1):
+        return [x1, x_bin[0] * (1 - x1), x_bin[1] * (1 - x1)]
+
+    ST_2b3 = bin_23.soret_coefficient(T, Vm, x_bin, N=N, use_zarate=True)[0]
+    _, ST_2t = ternary.soret_coefficient(T, Vm, get_ternary_x1(0.11), N=N, use_zarate=True)
+    diff = abs(ST_2t - ST_2b3)
+    x1_lst = np.logspace(-1, -5, 10)
+    for x1 in x1_lst:
+        _, ST_2t = ternary.soret_coefficient(T, Vm, get_ternary_x1(x1), N=N, use_zarate=True)
+        new_diff = abs(ST_2t - ST_2b3)
+        assert new_diff < diff
+        diff = new_diff
+
+    bin_13 = model('H2O,N2', is_idealgas=True)
+    def get_ternary_x2(x2):
+        return [x_bin[0] * (1 - x2), x2, x_bin[1] * (1 - x2)]
+
+    ST_1b3 = bin_13.soret_coefficient(T, Vm, x_bin, N=N, use_zarate=True)[0]
+    ST_1t, _ = ternary.soret_coefficient(T, Vm, get_ternary_x2(0.5), N=N, use_zarate=True)
+    diff = abs(ST_1t - ST_1b3)
+    x2_lst = np.logspace(-1, -5, 10)
+    for x2 in x2_lst:
+        ST_1t, _ = ternary.soret_coefficient(T, Vm, get_ternary_x2(x2), N=N, use_zarate=True)
+        new_diff = abs(ST_1t - ST_1b3)
+        assert new_diff < diff
+        diff = new_diff
+
+
 if __name__ == "__main__":
-    test_diffusion(models[0], 3)
+    test_soret(models[0], 2, silent=True)