Add ver-1kd

doc/content/verification_and_validation/figures/comparison_ver-1kd.py

@@ -0,0 +1 @@
+../../../../test/tests/ver-1kd/comparison_ver-1kd.py

doc/content/verification_and_validation/index.md

@@ -35,6 +35,7 @@ TMAP8 also contains [example cases](examples/tmap_index.md), which showcase how
 | ver-1kb | [Henry’s Law Boundaries with No Volumetric Source](ver-1kb.md)                                    |
 | ver-1kc-1 | [Sieverts’ Law Boundaries with No Volumetric Source](ver-1kc-1.md)                              |
 | ver-1kc-2 | [Sieverts’ Law Boundaries with Chemical Reaction and No Volumetric Source](ver-1kc-2.md)        |
+| ver-1kd | [Sieverts’ Law Boundaries with Chemical Reaction and Volumetric Source](ver-1kd.md)               |
 
 # List of benchmarking cases
 

doc/content/verification_and_validation/ver-1kd.md

@@ -0,0 +1,129 @@
+# ver-1kd
+
+# Sieverts’ Law Boundaries with Chemical Reactions and Volumetric Source
+
+## General Case Description
+
+Two enclosures are separated by a membrane that allows diffusion according to Sieverts' law and chemical reactions, with a volumetric source in Enclosure 1. Enclosure 2 has twice the volume of Enclosure 1.
+
+## Case Set Up
+
+This verification problem is taken from [!cite](ambrosek2008verification).
+
+This setup describes a diffusion system in which tritium T$_2$, dihydrogen H$_2$ and HT are modeled across a one-dimensional domain split into two enclosures. Compared to the [ver-1kc-2](ver-1kc-2.md) case, we now incorporate a T$_2$ tritium volumetric source in Enclosure 1. The volumetric source rate is set to $S_{\text{T}_2} = 10^{23} mol/m^3/s$.
+The total system length is $2.5 \times 10^{-4}$ m, divided into 100 segments. The system operates at a constant temperature of 500 Kelvin. Initial tritium T$_2$ and dihydrogen H$_2$ pressures are specified as $10^{5}$ Pa for Enclosure 1 and $10^{-10}$ Pa for Enclosure 2. Initially, there is no HT in either enclosure.
+
+The reaction between the species is described as follows
+
+\begin{equation}
+\text{H}_2 + \text{T}_2 \leftrightarrow 2\text{HT}
+\end{equation}
+
+The kinematic evolutions of the species are given by the following equations
+
+\begin{equation}
+\frac{d C_{\text{HT}}}{dt} = 2K_1 C_{\text{H}_2} C_{\text{T}_2} - K_2 C_{\text{HT}}^2
+\end{equation}
+
+\begin{equation}
+\frac{d C_{\text{H}_2}}{dt} = -K_1 C_{\text{H}_2} C_{\text{T}_2} + \frac{1}{2} K_2 C_{\text{HT}}^2
+\end{equation}
+
+\begin{equation}
+\frac{d C_{\text{T}_2}}{dt} = -K_1 c_{\text{H}_2} C_{\text{T}_2} + \frac{1}{2} K_2 C_{\text{HT}}^2
+\end{equation}
+
+where $K_1$ and $K_2$ represent the reaction rates for the forward and reverse reactions, respectively.
+
+At equilibrium, the time derivatives are zero
+
+\begin{equation}
+2K_1 C_{\text{H}_2} C_{\text{T}_2} - K_2 C_{\text{HT}}^2 = 0
+\end{equation}
+
+From this, we can derive the same equilibrium condition as used in TMAP7:
+
+\begin{equation}
+P_{\text{HT}} = \eta \sqrt{P_{\text{H}_2} P_{\text{T}_2}}
+\end{equation}
+
+where the equilibrium constant $\eta$ is defined as
+
+\begin{equation} \label{eq:eta}
+\eta = \sqrt{\frac{2K_1}{K_2}}
+\end{equation}
+
+Similarly to TMAP7, the equilibrium constant $\eta$ has been set to a fixed value of $\eta = 2$.
+
+The diffusion and generation processes for each species in the two enclosures can be expressed by
+
+\begin{equation}
+\frac{\partial C_1}{\partial t} = \nabla D \nabla C_1 + S_1,
+\end{equation}
+and
+\begin{equation}
+\frac{\partial C_2}{\partial t} = \nabla D \nabla C_2,
+\end{equation}
+
+where $C_1$ and $C_2$ represent the concentration fields in enclosures 1 and 2 respectively, $t$ is the time, $D$ denotes the diffusivity, $S$ the volumetric source rate, $V_1$ the volume of enclosure 1, $k$ the Boltzmann constant and $T$ the temperature.
+Note that the diffusivity may vary across different species and enclosures. However, in this case, it is assumed to be identical for all.
+
+The concentration in Enclosure 1 is related to the partial pressure and concentration in Enclosure 2 via the interface sorption law:
+
+\begin{equation}
+C_1 = K P_2^n = K \left( C_2 RT \right)^n
+\end{equation}
+
+where $R$ is the ideal gas constant in J/mol/K, $T$ is the temperature in K, $K$ is the solubility, and $n$ is the exponent of the sorption law. For Sieverts' law, $n=0.5$.
+
+## Results
+
+We assume that $K = 10/\sqrt{RT}$, which is expected to result in $C_1 = 10 \sqrt{C_2}$ at equilibrium.
+
+As illustrated in [ver-1kd_comparison_time_k10], since the chemical reactions occur immediately, an initial quantity of HT is present in Enclosure 1, while T$_2$ and H$_2$ initially drop to half their original amounts.
+
+The pressures of T$_2$ and H$_2$ in Enclosure 1 decrease due to diffusion into Enclosure 2. The constant flow rate of the T$_2$ source slows down the pressure drop of T$_2$ in Enclosure 2. Moreover, the volumetric source increases the amount of T$_2$ in Enclosure 1, thereby enhancing its chemical reactions with H$_2$. Consequently, in Enclosure 1, H$_2$ pressure gradually decreases over time, while HT pressure rises.
+
+Similarly, in Enclosure 2, the pressures of T$_2$ and H$_2$ increase, with the rise being more pronounced for T$_2$ due to the continuous supply from the source. As a result, more H$_2$ reacts with T$_2$, further contributing to the increase in HT pressure.
+
+For verification purposes, it is crucial to ensure that the chemical equilibrium between HT, T$_2$ and H$_2$ is achieved. This can be verified in both enclosures by examining the ratio between $P_{\text{HT}}$ and $\sqrt{P_{\text{H}_2} P_{\text{T}_2}}$, which must equal $\eta=2$.
+As shown in [ver-1kd_equilibrium_constant_k10], this ratio approaches $\eta=2$ for both enclosures at equilibrium. To reach this equilibrium, the ratio of $K_1$ and $K_2$ must respect [eq:eta]. The values of $K_1$ and $K_2$ must also be large enough to ensure that the kinetics of chemical reactions are faster than diffusion or surface permeation to be closer to the equilibrium assumption imposed in TMAP7. Here, the equilibrium in enclosure 1 is achieved rapidly. Increasing $K_1$ and $K_2$ would also enable a quicker attainment of equilibrium in enclosure 2. However, using very high values for $K_1$ and $K_2$ would lead to an unnecessary increase in computational costs.
+
+The concentration ratios for T$_2$, H$_2$, and HT between enclosures 1 and 2, shown in [ver-1kd_concentration_ratio_T2_k10], [ver-1kd_concentration_ratio_H2_k10], and [ver-1kd_concentration_ratio_HT_k10], demonstrate that the results obtained with TMAP8 are consistent with the analytical results derived from the sorption law for $K \sqrt{RT} = 10$.
+
+!media comparison_ver-1kd.py
+       image_name=ver-1kd_comparison_time_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kd_comparison_time_k10
+       caption=Evolution of species concentration over time governed by Sieverts' law with $K = 10/\sqrt{RT}$ and $\eta = \sqrt{2K_1/K_2}$.
+
+!media comparison_ver-1kd.py
+       image_name=ver-1kd_equilibrium_constant_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kd_equilibrium_constant_k10
+       caption=Equilibrium constant as a function of time for $\eta = \sqrt{2K_1/K_2}=2$.
+
+!media comparison_ver-1kd.py
+       image_name=ver-1kd_concentration_ratio_T2_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kd_concentration_ratio_T2_k10
+       caption=T$_2$ concentration ratio between enclosures 1 and 2 at the interface for $K = 10/\sqrt{RT}$ and $\eta = \sqrt{2K_1/K_2}$. This verifies TMAP8's ability to apply Sieverts' law across the interface.
+
+!media comparison_ver-1kd.py
+       image_name=ver-1kd_concentration_ratio_H2_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kd_concentration_ratio_H2_k10
+       caption=H$_2$ concentration ratio between enclosures 1 and 2 at the interface for $K = 10/\sqrt{RT}$ and $\eta = \sqrt{2K_1/K_2}$. This verifies TMAP8's ability to apply Sieverts' law across the interface.
+
+!media comparison_ver-1kd.py
+       image_name=ver-1kd_concentration_ratio_HT_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kd_concentration_ratio_HT_k10
+       caption=HT concentration ratio between enclosures 1 and 2 at the interface for $K = 10/\sqrt{RT}$ and $\eta = \sqrt{2K_1/K_2}$. This verifies TMAP8's ability to apply Sieverts' law across the interface.
+
+## Input files
+
+!style halign=left
+The input file for this case can be found at [/ver-1kd.i].
+
+!bibtex bibliography

test/tests/ver-1kd/comparison_ver-1kd.py

@@ -0,0 +1,171 @@
+import numpy as np
+import pandas as pd
+import os
+from matplotlib import gridspec
+import matplotlib.pyplot as plt
+
+# Changes working directory to script directory (for consistent MooseDocs usage)
+script_folder = os.path.dirname(__file__)
+os.chdir(script_folder)
+
+# Load experimental data
+if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+    csv_folder_k10 = "../../../../test/tests/ver-1kd/gold/ver-1kd_out_k10.csv"
+else:                                  # if in test folder
+    csv_folder_k10 = "./gold/ver-1kd_out_k10.csv"
+expt_data_k10 = pd.read_csv(csv_folder_k10)
+TMAP8_time_k10 = expt_data_k10['time']
+TMAP8_pressure_H2_enclosure_1_k10 = expt_data_k10['pressure_H2_enclosure_1']
+TMAP8_pressure_H2_enclosure_2_k10 = expt_data_k10['pressure_H2_enclosure_2']
+TMAP8_pressure_T2_enclosure_1_k10 = expt_data_k10['pressure_T2_enclosure_1']
+TMAP8_pressure_T2_enclosure_2_k10 = expt_data_k10['pressure_T2_enclosure_2']
+TMAP8_pressure_HT_enclosure_1_k10 = expt_data_k10['pressure_HT_enclosure_1']
+TMAP8_pressure_HT_enclosure_2_k10 = expt_data_k10['pressure_HT_enclosure_2']
+mass_conservation_sum_encl1_encl2_k10 = expt_data_k10['mass_conservation_sum_encl1_encl2']
+concentration_ratio_H2_k10 = expt_data_k10['concentration_ratio_H2']
+concentration_ratio_T2_k10 = expt_data_k10['concentration_ratio_T2']
+concentration_ratio_HT_k10 = expt_data_k10['concentration_ratio_HT']
+equilibrium_constant_encl_1 = expt_data_k10['equilibrium_constant_encl_1']
+equilibrium_constant_encl_2 = expt_data_k10['equilibrium_constant_encl_2']
+
+
+# Subplot 1 : Pressure vs time
+
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+# Plot for Enclosure 1
+line1 = ax.plot(TMAP8_time_k10, TMAP8_pressure_H2_enclosure_1_k10, label=r"H$_2$ Enclosure 1", c='tab:red', linestyle='-')
+line2 = ax.plot(TMAP8_time_k10, TMAP8_pressure_T2_enclosure_1_k10, label=r"T$_2$ Enclosure 1", c='tab:orange', linestyle='--')
+line3 = ax.plot(TMAP8_time_k10, TMAP8_pressure_HT_enclosure_1_k10, label=r"HT Enclosure 1", c='tab:purple', linestyle='-')
+
+ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
+ax.set_xlabel('Time (s)')
+ax.set_ylabel('Pressure Enclosure 1 (Pa)')
+ax.set_ylim(bottom=0)
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+
+# Plot for Enclosure 2
+ax2 = ax.twinx()
+line4 = ax2.plot(TMAP8_time_k10, TMAP8_pressure_H2_enclosure_2_k10, label=r"H$_2$ Enclosure 2", c='tab:blue', linestyle='-')
+line5 = ax2.plot(TMAP8_time_k10, TMAP8_pressure_T2_enclosure_2_k10, label=r"T$_2$ Enclosure 2", c='tab:green', linestyle='--')
+line6 = ax2.plot(TMAP8_time_k10, TMAP8_pressure_HT_enclosure_2_k10, label=r"HT Enclosure 2", c='tab:cyan', linestyle='-')
+
+ax2.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
+ax2.set_ylabel('Pressure Enclosure 2 (Pa)')
+ax2.set_ylim(bottom=0)
+ax.set_xlim(0,TMAP8_time_k10.max())
+
+# Combine legends
+lines_left  = line1 + line2 + line3
+lines_right = line4 + line5 + line6
+all_lines = lines_left + lines_right
+all_labels = [l.get_label() for l in all_lines]
+
+ax.legend(all_lines, all_labels, loc='upper left')
+fig.savefig('ver-1kd_comparison_time_k10.png', bbox_inches='tight', dpi=300)
+
+# Subplot 2: Solubility and concentration ratios vs time
+
+## Subplot 2.1: for H2
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+solubility_ratio = [10] * len(TMAP8_time_k10[1:])
+ax.plot(TMAP8_time_k10[1:], concentration_ratio_H2_k10[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
+ax.plot(TMAP8_time_k10[1:], solubility_ratio, label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
+ax.set_yticks(np.arange(0, 21, 10))
+ax.set_xlim(0,TMAP8_time_k10.max())
+ax.set_xlabel('Time (s)')
+ax.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / \sqrt{C_{\text{encl2}}}$")
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+RMSE = np.sqrt(np.mean((concentration_ratio_H2_k10[1:]-solubility_ratio)**2))
+RMSPE = RMSE*100/np.mean(solubility_ratio)
+x_pos = TMAP8_time_k10.max() / 7200
+y_pos = 0.9 * ax.get_ylim()[1]
+ax.text(x_pos, y_pos, 'RMSPE = %.3f ' % RMSPE + '%', fontweight='bold')
+fig.savefig('ver-1kd_concentration_ratio_H2_k10.png', bbox_inches='tight', dpi=300)
+
+## Subplot 2.2: for T2
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+solubility_ratio = [10] * len(TMAP8_time_k10[1:])
+ax.plot(TMAP8_time_k10[1:], concentration_ratio_T2_k10[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
+ax.plot(TMAP8_time_k10[1:], solubility_ratio, label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
+ax.set_yticks(np.arange(0, 21, 10))
+ax.set_xlim(0,TMAP8_time_k10.max())
+ax.set_xlabel('Time (s)')
+ax.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / \sqrt{C_{\text{encl2}}}$")
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+RMSE = np.sqrt(np.mean((concentration_ratio_T2_k10[1:]-solubility_ratio)**2))
+RMSPE = RMSE*100/np.mean(solubility_ratio)
+x_pos = TMAP8_time_k10.max() / 7200
+y_pos = 0.9 * ax.get_ylim()[1]
+ax.text(x_pos, y_pos, 'RMSPE = %.3f ' % RMSPE + '%', fontweight='bold')
+fig.savefig('ver-1kd_concentration_ratio_T2_k10.png', bbox_inches='tight', dpi=300)
+
+## Subplot 2.3: for HT
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+solubility_ratio = [10] * len(TMAP8_time_k10[1:])
+ax.plot(TMAP8_time_k10[1:], concentration_ratio_HT_k10[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
+ax.plot(TMAP8_time_k10[1:], solubility_ratio, label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
+ax.set_yticks(np.arange(0, 21, 10))
+ax.set_xlim(0,TMAP8_time_k10.max())
+ax.set_xlabel('Time (s)')
+ax.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / \sqrt{C_{\text{encl2}}}$")
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+RMSE = np.sqrt(np.mean((concentration_ratio_HT_k10[1:]-solubility_ratio)**2))
+RMSPE = RMSE*100/np.mean(solubility_ratio)
+x_pos = TMAP8_time_k10.max() / 7200
+y_pos = 0.9 * ax.get_ylim()[1]
+ax.text(x_pos, y_pos, 'RMSPE = %.3f ' % RMSPE + '%', fontweight='bold')
+fig.savefig('ver-1kd_concentration_ratio_HT_k10.png', bbox_inches='tight', dpi=300)
+
+
+# Subplot 3 : Mass Conservation Sum Encl 1 and 2 vs Time
+
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+ax.plot(TMAP8_time_k10, mass_conservation_sum_encl1_encl2_k10, c='tab:blue')
+ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.3e}'.format(val)))
+ax.set_xlabel('Time (s)')
+ax.set_xlim(-TMAP8_time_k10.max()/100,TMAP8_time_k10.max())
+ax.set_ylabel(r"Mass Conservation Sum Encl 1 and 2 (mol/m$^3$)")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+mass_variation_percentage = (np.max(mass_conservation_sum_encl1_encl2_k10)-np.min(mass_conservation_sum_encl1_encl2_k10))/np.max(mass_conservation_sum_encl1_encl2_k10)*100
+print("Percentage of mass variation: ", mass_variation_percentage)
+fig.savefig('ver-1kd_mass_conservation_k10.png', bbox_inches='tight', dpi=300)
+
+# Subplot 4 : Equilibrium constant in enclosures 1 and 2 vs Time
+
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+ax.plot(TMAP8_time_k10, equilibrium_constant_encl_2, label = r"Enclosure 2", c='tab:red')
+ax.plot(TMAP8_time_k10, equilibrium_constant_encl_1, label = r"Enclosure 1", c='tab:blue')
+ax.axhline(y=2, color='tab:green', linestyle='--', label='TMAP7 Equilibrium Constant')
+ax.set_xlabel('Time (s)')
+ax.set_xlim(0,TMAP8_time_k10.max())
+ax.set_ylabel(r"Equilibrium constant $P_{\text{HT}} / \sqrt{P_{\text{H}_2} P_{\text{T}_2}}$")
+ax.set_ylim(bottom=0)
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+print("Relative variation to equilibrium constant in enclosure 1", abs(equilibrium_constant_encl_1[len(equilibrium_constant_encl_1)-1]-2)/2 * 100)
+print("Relative variation to equilibrium constant in enclosure 2", abs(equilibrium_constant_encl_2[len(equilibrium_constant_encl_2)-1]-2)/2 * 100)
+fig.savefig('ver-1kd_equilibrium_constant_k10.png', bbox_inches='tight', dpi=300)
+
+
+