Merge pull request #225 from Lee01Atom/ver-1kc

Ver-1kc-1

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

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

doc/content/verification_and_validation/index.md

@@ -33,6 +33,7 @@ TMAP8 also contains [example cases](examples/tmap_index.md), which showcase how
 | ver-1jb | [Radioactive Decay of Mobile Tritium in a Slab with a Distributed Trap Concentration](ver-1jb.md) |
 | ver-1ka | [Simple Volumetric Source](ver-1ka.md)                                                            |
 | ver-1kb | [Henry’s Law Boundaries with No Volumetric Source](ver-1kb.md)                                    |
+| ver-1kc | [Sieverts’ Law Boundaries with No Volumetric Source](ver-1kc.md)                                  |
 
 # List of benchmarking cases
 

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

@@ -0,0 +1,64 @@
+# ver-1kc-1
+
+# Sieverts’ Law Boundaries with No Volumetric Source
+
+## General Case Description
+
+Two enclosures are separated by a membrane that allows diffusion according to Sievert's law, with no volumetric source present. 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$ is modeled across a one-dimensional domain split into two enclosures. 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 pressures are specified as $10^{5}$ Pa for Enclosure 1 and $10^{-10}$ Pa for Enclosure 2.
+
+The diffusion process in each of the two enclosures can be described by
+
+\begin{equation}
+\frac{\partial C_1}{\partial t} = \nabla D \nabla C_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, and $D$ denotes the diffusivity.
+
+This case is similar to the [ver-1kb](ver-1kb.md) case, with the key difference being that sorption here follows Sieverts' law instead of Henry's law.
+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 = \frac{10}{\sqrt{RT}}$, which is expected to lead to $C_1 = 10 \sqrt{C_2}$ at equilibrium.
+As illustrated in [ver-1kc_comparison_time_k10], the pressure jump maintains a ratio of $\frac{C_1}{\sqrt{C_2}} \approx 10$, which is consistent with the relationship $C_1 = K (RT C_2)^n$ for $K = \frac{10}{\sqrt{RT}}$ and $n=0.5$ The concentration ratio between enclosures 1 and 2 in [ver-1kc_concentration_ratio_k10] shows that the results obtained with TMAP8 are consistent with the analytical results derived from the sorption law for $K \sqrt{RT}=10$. As shown in [ver-1kc_mass_conservation_k10], mass is conserved between the two enclosures over time, with a variation in mass of only $0.4$ %. This variation in mass can be further minimized by refining the mesh, i.e., increasing the number of segments in the domain.
+
+!media comparison_ver-1kc.py
+       image_name=ver-1kc_comparison_time_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kc_comparison_time_k10
+       caption=Evolution of species concentration over time governed by Sieverts' law with $K = \frac{10}{\sqrt{RT}}$.
+
+!media comparison_ver-1kc.py
+       image_name=ver-1kc_concentration_ratio_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kc_concentration_ratio_k10
+       caption=Concentrations ratio between enclosures 1 and 2 at the interface for $K = \frac{10}{\sqrt{RT}}$. This verifies TMAP8's ability to apply Sieverts' law across the interface.
+
+!media comparison_ver-1kc.py
+       image_name=ver-1kc_mass_conservation_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kc_mass_conservation_k10
+       caption=Total mass conservation across both enclosures over time for $K = \frac{10}{\sqrt{RT}}$.
+
+## Input files
+
+!style halign=left
+The input file for this case can be found at [/ver-1kc.i]. To limit the computational costs of the test cases, the tests run a version of the file with a coarser mesh and less number of time steps. More information about the changes can be found in the test specification file for this case [/ver-1kc/tests].
+
+!bibtex bibliography

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

@@ -0,0 +1,78 @@
+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-1kc/gold/ver-1kc_out_k10.csv"
+else:                                  # if in test folder
+    csv_folder_k10 = "./gold/ver-1kc_out_k10.csv"
+expt_data_k10 = pd.read_csv(csv_folder_k10)
+TMAP8_time_k10 = expt_data_k10['time']
+TMAP8_pressure_enclosure_1_k10 = expt_data_k10['pressure_enclosure_1']
+TMAP8_pressure_enclosure_2_k10 = expt_data_k10['pressure_enclosure_2']
+concentration_enclosure_1_at_interface_k10 = expt_data_k10['concentration_enclosure_1_at_interface']
+pressure_enclosure_2_at_interface_k10 = expt_data_k10['pressure_enclosure_2_at_interface']
+mass_conservation_sum_encl1_encl2_k10 = expt_data_k10['mass_conservation_sum_encl1_encl2']
+concentration_ratio_k10 = expt_data_k10['concentration_ratio']
+
+# Subplot 1 : Pressure 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, TMAP8_pressure_enclosure_1_k10, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
+ax.plot(TMAP8_time_k10, TMAP8_pressure_enclosure_2_k10, label=r"T$_2$ Enclosure 2", c='tab:blue', linestyle='-')
+ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
+ax.set_xlabel('Time (s)')
+ax.set_ylabel('Pressure (Pa)')
+ax.set_xlim(0, TMAP8_time_k10.max())
+ax.set_ylim(bottom=0)
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+fig.savefig('ver-1kc_comparison_time_k10.png', bbox_inches='tight', dpi=300)
+
+# Subplot 2: Solubility and concentration ratios vs time
+
+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_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_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-1kc_concentration_ratio_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_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.min(mass_conservation_sum_encl1_encl2_k10)*100
+print("Percentage of mass variation: ", mass_variation_percentage)
+fig.savefig('ver-1kc_mass_conservation_k10.png', bbox_inches='tight', dpi=300)
+