Skip to content

Commit 681854b

Browse files
lkadzlkadz
committed
1 parent db7a621 commit 681854b

File tree

9 files changed

+1672
-0
lines changed

9 files changed

+1672
-0
lines changed

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
../../../../test/tests/ver-1kc/comparison_ver-1kc.py

‎doc/content/verification_and_validation/index.md

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,8 @@ TMAP8 also contains [example cases](examples/tmap_index.md), which showcase how
3232
| ver-1ja | [Radioactive Decay of Mobile Tritium in a Slab](ver-1ja.md) |
3333
| ver-1jb | [Radioactive Decay of Mobile Tritium in a Slab with a Distributed Trap Concentration](ver-1jb.md) |
3434
| ver-1ka | [Simple Volumetric Source](ver-1ka.md) |
35+
| ver-1kb | [Henry’s Law Boundaries with No Volumetric Source](ver-1kb.md) |
36+
| ver-1kc | [Sieverts’ Law Boundaries with No Volumetric Source](ver-1kb.md) |
3537

3638
# List of benchmarking cases
3739

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

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
# ver-1kc
2+
3+
# Sieverts’ Law Boundaries with No Volumetric Source
4+
5+
## General Case Description
6+
7+
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.
8+
9+
## Case Set Up
10+
11+
This verification problem is taken from [!cite](ambrosek2008verification).
12+
13+
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.
14+
15+
The diffusion process in each of the two enclosures can be described by
16+
17+
\begin{equation}
18+
\frac{\partial C_1}{\partial t} = \nabla D \nabla C_1,
19+
\end{equation}
20+
and
21+
\begin{equation}
22+
\frac{\partial C_2}{\partial t} = \nabla D \nabla C_2,
23+
\end{equation}
24+
25+
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.
26+
27+
The concentration in Enclosure 1 is related to the partial pressure and concentration in Enclosure 2 via the interface sorption law:
28+
29+
\begin{equation}
30+
C_1 = K P_2^n = K \left( C_2 RT \right)^n
31+
\end{equation}
32+
33+
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$.
34+
35+
## Results
36+
37+
Two subcases are considered. In the first subcase, we assume that $K=1/RT$, which is expected to lead to $C_1=C_2^{1/2}$ at equilibrium. In the second, $K=10/RT$, which is expected to lead to $C_1 = 10 C_2^{1/2}$. This second case is added to exercise TMAP8 in a case with a concentration jump.
38+
In the first subcase, the pressure evolution in both enclosures is shown in [ver-1kc_comparison_time] as a function of time. Both pressures find equilibrium, which is consistent with $C_1 = K (RT C_2)^n$ for $K=1/RT$ and $n=0.5$. The concentration ratio between enclosures 1 and 2 in [ver-1kc_concentration_ratio] shows that the results obtained with TMAP8 are consistent with the analytical results derived from the sorption law for $K R T=1$. As shown in [ver-1kc_mass_conservation], mass is conserved between the two enclosures over time, with a variation in mass of only $1.3$ %. This variation in mass can be further minimized by refining the mesh, i.e., increasing the number of segments in the domain.
39+
40+
!media comparison_ver-1kc.py
41+
image_name=ver-1kc_comparison_time.png
42+
style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
43+
id=ver-1kc_comparison_time
44+
caption=Equilibration of species pressures under Sieverts' law for $K = 1/RT$.
45+
46+
!media comparison_ver-1kc.py
47+
image_name=ver-1kc_concentration_ratio.png
48+
style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
49+
id=ver-1kc_concentration_ratio
50+
caption=Concentrations ratio between enclosures 1 and 2 at the interface for $K = 1/RT$. This verifies TMAP8's ability to apply the sorption law across the interface without a concentration jump.
51+
52+
!media comparison_ver-1kc.py
53+
image_name=ver-1kc_mass_conservation.png
54+
style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
55+
id=ver-1kc_mass_conservation
56+
caption=Total mass conservation across both enclosures over time for $K = 1/RT$.
57+
58+
In the second subcase, as illustrated in [ver-1kc_comparison_time_k10], the pressure jump maintains a ratio of $C_1/C_2^{1/2} \approx 10$, which is consistent with the relationship $C_1 = K (RT C_2)^n$ for $K=10/RT$ and $n=1$. 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 RT=10$. Additionally, [ver-1kc_mass_conservation_k10] verifies that mass is conserved between the two enclosures over time, with a variation in mass of only $0.4$ %. As in the previous case, this variation in mass can be further reduced by refining the mesh.
59+
60+
!media comparison_ver-1kc.py
61+
image_name=ver-1kc_comparison_time_k10.png
62+
style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
63+
id=ver-1kc_comparison_time_k10
64+
caption=Pressures jump of species under Sieverts' law for $K = 10/RT$.
65+
66+
!media comparison_ver-1kc.py
67+
image_name=ver-1kc_concentration_ratio_k10.png
68+
style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
69+
id=ver-1kc_concentration_ratio_k10
70+
caption=Concentrations ratio between enclosures 1 and 2 at the interface for $K = 10/RT$. This verifies TMAP8's ability to apply the sorption law across the interface with a concentration jump.
71+
72+
!media comparison_ver-1kc.py
73+
image_name=ver-1kc_mass_conservation_k10.png
74+
style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
75+
id=ver-1kc_mass_conservation_k10
76+
caption=Total mass conservation across both enclosures over time for $K = 10/RT$.
77+
78+
## Input files
79+
80+
!style halign=left
81+
The input file for this case can be found at [/ver-1kc.i], which is also used as tests in TMAP8 at [/ver-1kc/tests].
82+
83+
!bibtex bibliography

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

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
import numpy as np
2+
import pandas as pd
3+
import os
4+
from matplotlib import gridspec
5+
import matplotlib.pyplot as plt
6+
7+
# Changes working directory to script directory (for consistent MooseDocs usage)
8+
script_folder = os.path.dirname(__file__)
9+
os.chdir(script_folder)
10+
11+
# Extract columns for time, pressures, concentration_enclosure_1_at_interface, and concentration_enclosure_2_at_interface
12+
if "/TMAP8/doc/" in script_folder: # if in documentation folder
13+
csv_folder = "../../../../test/tests/ver-1kc/gold/ver-1kc_out_k1.csv"
14+
else: # if in test folder
15+
csv_folder = "./gold/ver-1kc_out_k1.csv"
16+
expt_data = pd.read_csv(csv_folder)
17+
TMAP8_time = expt_data['time']
18+
TMAP8_pressure_enclosure_1 = expt_data['pressure_enclosure_1']
19+
TMAP8_pressure_enclosure_2 = expt_data['pressure_enclosure_2']
20+
concentration_enclosure_1_at_interface = expt_data['concentration_enclosure_1_at_interface']
21+
pressure_enclosure_2_at_interface = expt_data['pressure_enclosure_2_at_interface']
22+
mass_conservation_sum_encl1_encl2 = expt_data['mass_conservation_sum_encl1_encl2'].values
23+
concentration_ratio = expt_data['concentration_ratio']
24+
25+
# Subplot 1: Pressure vs time
26+
fig = plt.figure(figsize=[6.5,5.5])
27+
gs = gridspec.GridSpec(1,1)
28+
ax = fig.add_subplot(gs[0])
29+
30+
ax.plot(TMAP8_time, TMAP8_pressure_enclosure_1, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
31+
ax.plot(TMAP8_time, TMAP8_pressure_enclosure_2, label=r"T$_2$ Enclosure 2", c='tab:blue', linestyle='-')
32+
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
33+
ax.set_xlabel('Time (s)')
34+
ax.set_ylabel('Pressure (Pa)')
35+
ax.set_xlim(0, TMAP8_time.max())
36+
ax.set_ylim(bottom=0)
37+
ax.legend(loc="best")
38+
ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
39+
fig.savefig('ver-1kc_comparison_time.png', bbox_inches='tight', dpi=300)
40+
41+
# Subplot 2: Solubility and concentration ratios vs time
42+
fig = plt.figure(figsize=[6.5,5.5])
43+
gs = gridspec.GridSpec(1,1)
44+
ax = fig.add_subplot(gs[0])
45+
46+
solubility_ratio = [1] * len(TMAP8_time[1:])
47+
ax.plot(TMAP8_time[1:], concentration_ratio[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
48+
ax.plot(TMAP8_time[1:], solubility_ratio, label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
49+
ax.set_yticks(np.arange(0, 3, 1))
50+
ax.set_xlim(0,TMAP8_time.max())
51+
ax.set_xlabel('Time (s)')
52+
ax.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / C_{\text{encl2}}$")
53+
ax.legend(loc="best")
54+
ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
55+
RMSE = np.sqrt(np.mean((concentration_ratio[1:]-solubility_ratio)**2) )
56+
RMSPE = RMSE*100/np.mean(solubility_ratio)
57+
x_pos = TMAP8_time.max() / 7200
58+
y_pos = 0.9 * ax.get_ylim()[1]
59+
ax.text(x_pos, y_pos, 'RMSPE = %.3f ' % RMSPE + '%', fontweight='bold')
60+
fig.savefig('ver-1kc_concentration_ratio.png', bbox_inches='tight', dpi=300)
61+
62+
# Subplot 3: Mass Conservation Sum Encl 1 and 2 vs Time
63+
64+
fig = plt.figure(figsize=[6.5,5.5])
65+
gs = gridspec.GridSpec(1,1)
66+
ax = fig.add_subplot(gs[0])
67+
68+
ax.plot(TMAP8_time, mass_conservation_sum_encl1_encl2, c='tab:blue')
69+
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.3e}'.format(val)))
70+
ax.set_xlabel('Time (s)')
71+
ax.set_ylabel(r"Mass Conservation Sum Encl 1 and 2 (mol/m$^3$)")
72+
ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
73+
mass_variation_percentage = (np.max(mass_conservation_sum_encl1_encl2)-np.min(mass_conservation_sum_encl1_encl2))/np.min(mass_conservation_sum_encl1_encl2)*100
74+
print("Percentage of mass variation: ", mass_variation_percentage)
75+
fig.savefig('ver-1kc_mass_conservation.png', bbox_inches='tight', dpi=300)
76+
77+
# Repeat the same for K=10/RT
78+
79+
if "/TMAP8/doc/" in script_folder: # if in documentation folder
80+
csv_folder_k10 = "../../../../test/tests/ver-1kc/gold/ver-1kc_out_k10.csv"
81+
else: # if in test folder
82+
csv_folder_k10 = "./gold/ver-1kc_out_k10.csv"
83+
expt_data_k10 = pd.read_csv(csv_folder_k10)
84+
TMAP8_time_k10 = expt_data_k10['time']
85+
TMAP8_pressure_enclosure_1_k10 = expt_data_k10['pressure_enclosure_1']
86+
TMAP8_pressure_enclosure_2_k10 = expt_data_k10['pressure_enclosure_2']
87+
concentration_enclosure_1_at_interface_k10 = expt_data_k10['concentration_enclosure_1_at_interface']
88+
pressure_enclosure_2_at_interface_k10 = expt_data_k10['pressure_enclosure_2_at_interface']
89+
mass_conservation_sum_encl1_encl2_k10 = expt_data_k10['mass_conservation_sum_encl1_encl2'].values
90+
concentration_ratio_k10 = expt_data_k10['concentration_ratio']
91+
92+
# Subplot 1 : Pressure vs time
93+
94+
fig = plt.figure(figsize=[6.5,5.5])
95+
gs = gridspec.GridSpec(1,1)
96+
ax = fig.add_subplot(gs[0])
97+
98+
ax.plot(TMAP8_time_k10, TMAP8_pressure_enclosure_1_k10, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
99+
ax.plot(TMAP8_time_k10, TMAP8_pressure_enclosure_2_k10, label=r"T$_2$ Enclosure 2", c='tab:blue', linestyle='-')
100+
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
101+
ax.set_xlabel('Time (s)')
102+
ax.set_ylabel('Pressure (Pa)')
103+
ax.set_xlim(0, TMAP8_time_k10.max())
104+
ax.set_ylim(bottom=0)
105+
ax.legend(loc="best")
106+
ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
107+
fig.savefig('ver-1kc_comparison_time_k10.png', bbox_inches='tight', dpi=300)
108+
109+
# Subplot 2: Solubility and concentration ratios vs time
110+
111+
fig = plt.figure(figsize=[6.5,5.5])
112+
gs = gridspec.GridSpec(1,1)
113+
ax = fig.add_subplot(gs[0])
114+
115+
solubility_ratio = [10] * len(TMAP8_time_k10[1:])
116+
ax.plot(TMAP8_time_k10[1:], concentration_ratio_k10[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
117+
ax.plot(TMAP8_time_k10[1:], solubility_ratio, label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
118+
ax.set_yticks(np.arange(0, 21, 10))
119+
ax.set_xlim(0,TMAP8_time_k10.max())
120+
ax.set_xlabel('Time (s)')
121+
ax.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / C_{\text{encl2}}$")
122+
ax.legend(loc="best")
123+
ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
124+
RMSE = np.sqrt(np.mean((concentration_ratio_k10[1:]-solubility_ratio)**2))
125+
RMSPE = RMSE*100/np.mean(solubility_ratio)
126+
x_pos = TMAP8_time_k10.max() / 7200
127+
y_pos = 0.9 * ax.get_ylim()[1]
128+
ax.text(x_pos, y_pos, 'RMSPE = %.3f ' % RMSPE + '%', fontweight='bold')
129+
fig.savefig('ver-1kc_concentration_ratio_k10.png', bbox_inches='tight', dpi=300)
130+
131+
# Subplot 3 : Mass Conservation Sum Encl 1 and 2 vs Time
132+
133+
fig = plt.figure(figsize=[6.5,5.5])
134+
gs = gridspec.GridSpec(1,1)
135+
ax = fig.add_subplot(gs[0])
136+
137+
ax.plot(TMAP8_time_k10, mass_conservation_sum_encl1_encl2_k10, c='tab:blue')
138+
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.3e}'.format(val)))
139+
ax.set_xlabel('Time (s)')
140+
ax.set_ylabel(r"Mass Conservation Sum Encl 1 and 2 (mol/m$^3$)")
141+
ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
142+
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
143+
print("Percentage of mass variation: ", mass_variation_percentage)
144+
fig.savefig('ver-1kc_mass_conservation_k10.png', bbox_inches='tight', dpi=300)
145+

0 commit comments

Comments
 (0)