Initial condition, dummy export and new test case

Validation/DWO/exp.csv

@@ -0,0 +1,66 @@
+x, y
+304.6498408177643, 17622949651356056
+315.59076058557696, 10689316241137326
+315.59076058557696, 10689316241137326
+322.42883544045986, 24556583061574784
+334.0535626937608, 33801427608533092
+337.47260012120216, 49979905565710130
+337.47260012120216, 49979905565710130
+345.67828994706167, 100826550573980820
+355.2515947438977, 130872295351595310
+360.03824714231575, 144739562172032770
+364.8248995407338, 163229251265949380
+364.8248995407338, 163229251265949380
+375.7658193085464, 190963784906824300
+375.7658193085464, 190963784906824300
+380.55247170696447, 237188007641615800
+384.65531661989417, 260300119009011600
+390.80958398928874, 285723441513146940
+395.5962363877068, 331947664247938400
+400.38288878612485, 371238253572511200
+405.1695411845428, 401283998350125760
+411.3238085539375, 429018531991000700
+422.2647283217501, 442885798811438140
+422.2647283217501, 442885798811438140
+429.7866106621213, 431329743127740200
+435.9408780315159, 405906420623604900
+440.04372294444556, 366615831299032100
+448.24941277030507, 352748564478594600
+457.13891008165285, 329636453111198900
+468.7636373349537, 292657074923365600
+468.7636373349537, 292657074923365600
+475.6017121898367, 276478596966188640
+487.2264394431377, 251055274462053280
+498.16735921095017, 232565585368136670
+505.6892415513215, 218698318547699200
+521.416813717552, 190963784906824300
+529.6225035434117, 167851673539428540
+543.2986532531773, 133183506488334880
+554.23957302099, 107760183984199540
+562.4452628468493, 84648072616803780
+575.4376050711271, 82336861480064200
+588.4299472954045, 47668694428970550
+588.4299472954045, 47668694428970550
+599.3708670632172, 33801427608533092
+609.6279793455415, 24556583061574784
+617.833669171401, 19934160788095630
+630.8260113956785, 13000527377876904
+630.8260113956785, 13000527377876904
+647.2373910473974, 1444471694179021
+653.3916584167919, 24556583061574784
+660.9135407571631, 13000527377876904
+669.8030380685109, 3755682830918597.5
+678.6925353798588, -10111583989518862
+686.8982252057183, -3177950579300132
+697.8391449735307, -3177950579300132
+706.0448347993904, -3177950579300132
+713.5667171397615, -3177950579300132
+729.2942893059922, -5489161716039709
+729.2942893059922, -5489161716039709
+738.1837866173398, -7800372852779285
+746.3894764431994, -7800372852779285
+755.2789737545471, -7800372852779285
+763.4846635804066, -7800372852779285
+772.3741608917544, -7800372852779285
+792.8883854564031, -866739442560555.6
+792.8883854564031, -866739442560555.6

Validation/DWO/test.py

@@ -0,0 +1,187 @@
+import festim as F
+import fenics as f
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+import sympy as sp
+import warnings
+
+warnings.filterwarnings("ignore")
+
+################### PARAMETERS ###################
+N_A_const = 6.022e23  # Avogadro, mol^-1
+e = 1.602e-19
+M_D2 = 4.028e-3 / N_A_const  # the H2 mass, kg mol^-1
+
+# Exposure conditions
+P_D2 = 2e-5  # Pa
+T0 = 300  # K
+ramp = 5  # K/s
+t_imp = 3000  # exposure duration, s
+t_storage = 3600  # storage time, s
+t_TDS = 100  # s
+L = 2e-3 / 2  # half thickness, m
+
+# W properties
+rho_W = 6.3382e28  # W atomic concentration, m^-3
+n_IS = 6 * rho_W  # concentration of interstitial sites, m^-3
+n_surf = 1.416e19  # concentration of adsorption sites, m^-2
+lamda_des = 1 / np.sqrt(n_surf)  # average distance between adsorption sites, m
+lambda_abs = n_surf / n_IS  # distance between subsurface IS and adsorption sites, m
+nu0 = 1e13  # attempt frequency, s^-1
+
+D0 = 1.9e-7 / np.sqrt(2)  # diffusivity pre-factor, m^2 s^-1
+E_diff = 0.2  # diffusion activation energy, eV
+
+# Transitions
+E_bs = E_diff  # energy barrier from bulk to surface, eV
+E_diss = 0  # energy barrier for D2 dissociation, eV
+Q_sol = 1  # heat of solution, eV
+
+
+################### FUNCTIONS ###################
+def S0(T):
+    # the capturing coefficient
+    return f.exp(-E_diss / F.k_B / T)
+
+
+def E_des(surf_conc):
+    theta_D = surf_conc / n_surf
+    theta_D0 = 0.253
+    del_theta_D = 0.180
+    E0 = 1.142
+    dE = 0.346
+    alpha = 0.303
+    beta = 8.902
+
+    E_FD = E0 + dE / (1 + f.exp((theta_D - theta_D0) / del_theta_D))
+    E_des = E_FD * (1 - alpha * f.exp(-beta * (1 - theta_D)))
+
+    return E_des
+
+
+def E_sb(surf_conc):
+    # energy barrier from surface to bulk, eV
+    return (E_des(surf_conc) - E_diss) / 2 + E_bs + Q_sol
+
+
+def J_vs(surf_conc, T, t):
+    J_ads = (
+        2
+        * S0(T)
+        * (1 - surf_conc / n_surf) ** 2
+        * f.conditional(t <= t_imp, 1.52e18, 0)
+    )
+
+    J_des = (
+        2 * nu0 * (lamda_des * surf_conc) ** 2 * f.exp(-E_des(surf_conc) / F.k_B / T)
+    )
+    return J_ads - J_des
+
+
+################### MODEL ###################
+W_model = F.Simulation(log_level=20)
+
+# Mesh
+vertices = np.linspace(0, L, num=500)
+W_model.mesh = F.MeshFromVertices(np.sort(vertices))
+
+# Materials
+tungsten = F.Material(id=1, D_0=D0, E_D=E_diff)
+W_model.materials = tungsten
+
+surf_conc = F.SurfaceConcentration(
+    k_sb=nu0,
+    E_sb=E_sb,
+    k_bs=nu0,
+    E_bs=E_bs,
+    l_abs=lambda_abs,
+    N_s=n_surf,
+    N_b=n_IS,
+    J_vs=J_vs,
+    surfaces=1,
+    initial_condition=F.InitialCondition(field="adsorbed", value=0),
+    t=F.t,
+)
+
+W_model.surface_concentrations = [surf_conc]
+
+traps = F.Traps(
+    [
+        F.Trap(
+            k_0=nu0 / n_IS,
+            E_k=E_diff,
+            p_0=nu0,
+            E_p=0.912,
+            density=7.7e-6,
+            materials=tungsten,
+        ),
+        F.Trap(
+            k_0=nu0 / n_IS,
+            E_k=E_diff,
+            p_0=nu0,
+            E_p=1.176,
+            density=7.9e-7,
+            materials=tungsten,
+        ),
+    ]
+)
+# W_model.traps = traps
+
+W_model.T = F.Temperature(
+    value=sp.Piecewise(
+        (T0, F.t <= t_imp + t_storage), (T0 + ramp * (F.t - t_imp - t_storage), True)
+    )
+)
+
+W_model.dt = F.Stepsize(
+    initial_value=1e-7,
+    stepsize_change_ratio=1.25,
+    max_stepsize=lambda t: 10 if t < t_imp + t_storage else 0.1,
+    dt_min=1e-9,
+)
+
+W_model.settings = F.Settings(
+    absolute_tolerance=1e8,
+    relative_tolerance=1e-7,
+    maximum_iterations=50,
+    final_time=t_imp + t_storage + t_TDS,
+)
+
+# Exports
+results_folder = "test/"
+
+derived_quantities = F.DerivedQuantities(
+    [F.AdsorbedHydrogen(surface=1), F.TotalSurface(field="T", surface=1)],
+    nb_iterations_between_compute=1,
+    filename=results_folder + "Clean_W.csv",
+)
+
+W_model.exports = [derived_quantities]
+
+W_model.initialise()
+W_model.run()
+
+results = pd.read_csv(results_folder + "Clean_W.csv", header=0)
+results["Flux"] = results["Concentration of adsorbed H on surface 1"]
+for i in range(len(results)):
+    surf_conc = results["Concentration of adsorbed H on surface 1"][i]
+    T = results["Total T surface 1"][i]
+    results["Flux"][i] = (
+        2 * nu0 * (lamda_des * surf_conc) ** 2 * np.exp(-E_des(surf_conc) / F.k_B / T)
+    )
+
+
+plt.plot(
+    results["Total T surface 1"][results["t(s)"] >= t_imp + t_storage],
+    results["Flux"][results["t(s)"] >= t_imp + t_storage],
+    linewidth=3,
+)
+
+exp = pd.read_csv(f"./exp.csv", header=None, skiprows=1, sep=",")
+plt.scatter(exp[0], exp[1], marker="x", s=75, linewidths=1.2)
+
+plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
+plt.xlabel(r"Time (s)")
+plt.xlim(300, 800)
+plt.show()