Merge pull request #75 from rekomodo/tmap_1db

initial implementation of tmap_1db

report/_toc.yml

@@ -28,6 +28,7 @@ parts:
         - file: verification/mes/tmap_1b
         - file: verification/mes/tmap_1c
         - file: verification/mes/tmap_1da
+        - file: verification/mes/tmap_1db
   - caption: Validation
     chapters:
       - file: validation/plasma-driven-permeation/plasma-driven-permeation

report/verification/mes/tmap_1da.md

@@ -28,7 +28,7 @@ $$
 where
 
 $\lambda \ \mathrm{(m)}$ is the lattice parameter, \
-$\nu \ (\mathrm{s}^{-1})$, \
+$\nu \ (\mathrm{s}^{-1})$, the Debye frequency \
 $\rho$ is the trapping site fraction, \
 $c_m (\text{atom} \ \mathrm{m}^{-3})$ is the mobile atom concentration,
 

report/verification/mes/tmap_1db.md

@@ -0,0 +1,138 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.4
+kernelspec:
+  display_name: festim-report
+  language: python
+  name: python3
+---
+
+# Strong Trapping Regime
+
+```{tags} 1D, MES, transient, trapping
+```
+
+This verification case consists of a slab of depth $l = 1 \times 10^{-3} \ \mathrm{m}$ with one trap under the strong trapping regime.
+
+A trapping parameter $\zeta$ is defined by
+
+$$
+    \zeta = \frac{\lambda ^ 2 \nu}{D_0 \rho} \exp \left(\frac{E_k - E_p}{k_B T}\right) + \frac{c_m}{n}
+$$
+
+where
+
+$\lambda \ \mathrm{(m)}$ is the lattice parameter, \
+$\nu \ (\mathrm{s}^{-1})$, the Debye frequency \
+$n$ is the trapping site fraction, \
+$c_m (\text{atom} \ \mathrm{m}^{-3})$ is the mobile atom concentration,
+
+and the effective diffusivity $\mathrm{D_\text{eff}}$ is defined by
+
+$$
+    D_\text{eff} = \frac{D}{1 + \frac{1}{\zeta}}
+$$
+
+Then the breakthrough time is given by
+
+$$
+    \tau = \frac{l^2}{2\pi^2 D_\text{eff}}
+$$
+
+This analytical solution was obtained from TMAP7's V&V report {cite}`ambrosek_verification_2008`, case Val-1db.
+
++++
+
+## FESTIM Code
+
+```{code-cell} ipython3
+:tags: [hide-cell]
+
+import festim as F
+import numpy as np
+import scipy.constants as const
+import matplotlib.pyplot as plt
+
+l = 0.1e-3
+D_0 = 1.9e-7
+T = 1000
+n_t = 1e-3
+w_atom_density = 6.3382e28
+N_A = const.N_A
+
+my_model = F.Simulation()
+vertices = np.concatenate(
+    [np.linspace(0, 0.9 * l, num=200), np.linspace(0.9 * l, l, num=100)]
+)
+my_model.mesh = F.MeshFromVertices(vertices)
+my_model.materials = F.Material(id=1, D_0=D_0, E_D=0.2)
+my_model.T = T
+
+my_model.boundary_conditions = [
+    F.DirichletBC(surfaces=1, value=0.0088 * N_A, field=0),
+    F.DirichletBC(surfaces=2, value=0, field=0),
+]
+
+trap = F.Trap(
+    k_0=1.58e7 / N_A,
+    E_k=0.2,
+    p_0=1e13,
+    E_p=2.5,
+    density=n_t * w_atom_density,
+    materials=my_model.materials[0],
+)
+
+my_model.traps = [trap]
+my_model.settings = F.Settings(
+    absolute_tolerance=1e10,
+    relative_tolerance=1e-10,
+    final_time=4000,  # s
+)
+
+my_model.dt = F.Stepsize(
+    initial_value=1e-3,
+    dt_min=1e-3,
+    stepsize_change_ratio=1.1,
+    max_stepsize=lambda t: 1 if 3250 >= t >= 3100 else None,
+)
+
+derived_quantities = F.DerivedQuantities([F.HydrogenFlux(surface=2)], show_units=True)
+
+my_model.exports = [derived_quantities, F.XDMFExport("1", checkpoint=False)]
+
+my_model.initialise()
+my_model.run()
+```
+
+## Comparison with exact solution
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+times = derived_quantities.t
+computed_flux = derived_quantities.filter(surfaces=2).data
+
+plt.plot(times, np.abs(computed_flux) / 2, label="FESTIM")
+
+D_eff = D_0 * np.exp(-0.2 / F.k_B / T)
+S = 2.9e-5 * np.exp(-1 / F.k_B / T)
+c0m = np.sqrt(1e5) * S * 1.0525e5
+
+time_exact = l**2 * n_t / (2 * c0m * D_eff) * w_atom_density / N_A
+
+plt.axvline(x=time_exact, color="r", label="analytical")
+
+plt.ylim(bottom=0)
+plt.xlabel("Time (s)")
+plt.ylabel("Downstream flux (H2/m2/s)")
+# plt.ylim([1e-6, 0.5e17])
+plt.xlim([3000, 3500])
+
+plt.legend()
+plt.show()
+```