Damaged induced traps#
Objectives
Learn how to model dynamic extrinsic traps in a FESTIM model
Research indicates that neutron-like damage in materials creates defects that trap hydrogen, potentially increasing hydrogen retention, which eventually saturates after a certain damage level. A proposed model aims to predict the temporal evolution of these damage-induced trap densities, considering the competition between hydrogen migration and trapping properties.
This equation describes the temporal evolution of trap density \(n_{i}\). The creation term, \(\Phi \cdot K \left[1-\frac{n_{i}}{n_{\mathrm{max},\Phi}}\right]\), describes the number of traps formed by damage to a given saturation point \(n_{\mathrm{max},\Phi}\). \(\Phi\) is the damage rate in dpa \(\text{s}^{-1}\), and can vary spatially. \(K\) is the trap creation factor in traps \(\text{m}^{-3} \text{dpa}^{-1}\). \(A\) is the trap annealing factor in \(\text{s}^{-1}\) and is described as an Arrhenius law where \(A = A_{0} \cdot \exp(-E_{A}/(k_{B} \cdot T))\).
We can use festim to show how such a trap may influence hydrogen retention in tungsten.
We will consider 2 models of 2mm of tungsten, one with a standard intrinsic trap, and another with a neutron-damage induced trap.
FESTIM models#
Let’s start with a mesh and material properties. We will consider the diffusion transport properties from Frauenfelder.
import numpy as np
import festim as F
import h_transport_materials as htm
# define mesh
my_mesh = F.MeshFromVertices(np.linspace(0, 2e-03, num=1000))
# define materials
D = htm.diffusivities.filter(material="tungsten").filter(author="frauenfelder")[0]
tungsten = F.Material(D_0=D.pre_exp.magnitude, E_D=D.act_energy.magnitude, id=1)
my_materials = F.Materials(
[F.Material(D_0=D.pre_exp.magnitude, E_D=D.act_energy.magnitude, id=1)]
)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 3
1 import numpy as np
2 import festim as F
----> 3 import h_transport_materials as htm
4
5 # define mesh
6 my_mesh = F.MeshFromVertices(np.linspace(0, 2e-03, num=1000))
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/h_transport_materials/__init__.py:25
22 Rg = 8.314 * ureg.Pa * ureg.m**3 * ureg.mol**-1 * ureg.K**-1
23 avogadro_nb = 6.022e23 * ureg.particle * ureg.mol**-1
---> 25 from pybtex.database import parse_file
26 from pathlib import Path
28 bib_database = parse_file(str(Path(__file__).parent) + "/references.bib")
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/pybtex/database/__init__.py:44
42 from pybtex.errors import report_error
43 from pybtex.py3compat import fix_unicode_literals_in_doctest, python_2_unicode_compatible
---> 44 from pybtex.plugin import find_plugin
47 # for python2 compatibility
48 def indent(text, prefix):
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/pybtex/plugin/__init__.py:26
2 # Copyright (c) 2006-2021 Andrey Golovizin
3 # Copyright (c) 2014 Matthias C. M. Troffaes
4 #
(...) 21 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
22 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
25 import os.path # splitext
---> 26 import pkg_resources
28 from pybtex.exceptions import PybtexError
31 class Plugin(object):
ModuleNotFoundError: No module named 'pkg_resources'
Next we can define the intrinsic and damage-induced extrinsic trap. We will consider the tungsten is being damaged at a rate of 5 dpa/fpy.
# define traps
lambda_val = 1.1e-10 # m
atom_density_W = 6.3222e28
trap_W = F.Trap(
k_0=tungsten.D_0 / (lambda_val**2 * 6 * atom_density_W),
E_k=tungsten.E_D,
p_0=1e13,
E_p=1.04,
density=2.4e22,
materials=tungsten,
)
fpy = 3600 * 24 * 365
trap_W_damage = F.NeutronInducedTrap(
k_0=tungsten.D_0 / (lambda_val**2 * 6 * atom_density_W),
E_k=tungsten.E_D,
p_0=1e13,
E_p=1.85,
A_0=6.1838e-03,
E_A=0.30,
phi=5 / fpy,
K=5.0e26,
n_max=4.7e25,
materials=tungsten,
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 6
2 lambda_val = 1.1e-10 # m
3 atom_density_W = 6.3222e28
4
5 trap_W = F.Trap(
----> 6 k_0=tungsten.D_0 / (lambda_val**2 * 6 * atom_density_W),
7 E_k=tungsten.E_D,
8 p_0=1e13,
9 E_p=1.04,
NameError: name 'tungsten' is not defined
Next we can define the other necessary parameters. The tungsten will be modelled at constant temperature of 800K and subjected to a source of hydrogen implanting on one surface at a flux of \(10^{20}\) with an associated implantation depth of 3nm.
The simulation will be a transient simulating exposure in a full power year (fpy).
# define temperature
temperature = 800
# define boundary conditions
my_boundary_conditions = [
F.ImplantationDirichlet(
surfaces=1,
phi=1e20,
R_p=3e-09,
D_0=tungsten.D_0,
E_D=tungsten.E_D,
),
]
# define exports
my_exports = [
F.DerivedQuantities(
[F.TotalVolume("solute", volume=1)],
show_units=True,
),
]
# define temporal settings
my_dt = F.Stepsize(
initial_value=0.1,
stepsize_change_ratio=1.1,
dt_min=1e-1,
)
# define settings
my_settings = F.Settings(
transient=True,
final_time=fpy,
absolute_tolerance=1e10,
relative_tolerance=1e-10,
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 10
6 F.ImplantationDirichlet(
7 surfaces=1,
8 phi=1e20,
9 R_p=3e-09,
---> 10 D_0=tungsten.D_0,
11 E_D=tungsten.E_D,
12 ),
13 ]
NameError: name 'tungsten' is not defined
Let’s now create a F.Simulation() object and fill it with parameters for a standard case and run the model:
results_folder = "task17/"
my_exports_standard = my_exports + [
F.TXTExport(
field="retention",
filename=results_folder + "retention_profile_standard.txt",
)
]
my_standard_model = F.Simulation(
log_level=40,
mesh=my_mesh,
materials=my_materials,
boundary_conditions=my_boundary_conditions,
traps=F.Traps([trap_W]),
temperature=temperature,
exports=my_exports_standard,
dt=my_dt,
settings=my_settings
)
my_standard_model.initialise()
my_standard_model.run()
Let’s run another simulation but accounting for the damage-induced trap:
my_exports_damaged = my_exports + [
F.TXTExport(
field="retention",
filename=results_folder + "retention_profile_damaged.txt",
)
]
my_damaged_model = F.Simulation(
log_level=40,
mesh=my_mesh,
materials=my_materials,
boundary_conditions=my_boundary_conditions,
traps=F.Traps([trap_W, trap_W_damage]),
temperature=temperature,
exports=my_exports_damaged,
dt=my_dt,
settings=my_settings
)
my_damaged_model.initialise()
my_damaged_model.run()
Results#
We can plot the concentration profiles with or without the damage included over time in a GIF. With very high levels of damage, it is clear that damage induced traps could have a huge impact on the results.
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
# read results data
results_folder = "task17/"
data_standard = np.genfromtxt(
results_folder + "retention_profile_standard.txt", delimiter=",", skip_header=1
)
data_damaged = np.genfromtxt(
results_folder + "retention_profile_damaged.txt", delimiter=",", skip_header=1
)
times = my_damaged_model.exports[0].t
# Separate the x values and y values for each time step
x_standard, y_standard = data_standard[:, 0] * 1e3, data_standard[:, 1:]
x_damaged, y_damaged = data_damaged[:, 0] * 1e3, data_damaged[:, 1:]
# Set up the figure and axis
fig, ax = plt.subplots()
(line1,) = ax.plot([], [], lw=2, color="black", label="standard")
(line2,) = ax.plot([], [], lw=2, color="red", label="damaged")
# Initialization function to set up the background of each frame
def init():
line1.set_data([], [])
line2.set_data([], [])
ax.set_xlim(0, 2)
ax.set_ylim(1e14, 1e26)
ax.set_yscale("log")
ax.legend()
ax.set_xlabel(r"x (mm)")
ax.set_ylabel(r"T retention (m$^{-3}$)")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
return line1, line2
def animate(i):
y1 = y_standard[:, i]
y2 = y_damaged[:, i]
line1.set_data(x_standard, y1)
line2.set_data(x_damaged, y2)
one_min, one_hour, one_day = 60, 3600, 3600 * 24
t = times[i] # NOTE this only works because both models have the same timesteps
if t < one_min:
ax.set_title(f"Time: {t:.1f} seconds")
elif t < one_hour:
ax.set_title(f"Time: {t/one_min:.1f} minutes")
elif t < one_day:
ax.set_title(f"Time: {t/one_hour:.1f} hours")
else:
ax.set_title(f"Time: {t/one_day:.1f} days")
return line1, line2
# Call the animator
ani = animation.FuncAnimation(
fig, animate, init_func=init, frames=y_damaged.shape[1], interval=50, blit=True
)
plt.close(fig)
# Display the animation in the Jupyter Notebook
HTML(ani.to_jshtml())
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[6], line 7
3 from IPython.display import HTML
4
5 # read results data
6 results_folder = "task17/"
----> 7 data_standard = np.genfromtxt(
8 results_folder + "retention_profile_standard.txt", delimiter=",", skip_header=1
9 )
10 data_damaged = np.genfromtxt(
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_npyio_impl.py:1978, in genfromtxt(fname, dtype, comments, delimiter, skip_header, skip_footer, converters, missing_values, filling_values, usecols, names, excludelist, deletechars, replace_space, autostrip, case_sensitive, defaultfmt, unpack, usemask, loose, invalid_raise, max_rows, encoding, ndmin, like)
1976 fname = os.fspath(fname)
1977 if isinstance(fname, str):
-> 1978 fid = np.lib._datasource.open(fname, 'rt', encoding=encoding)
1979 fid_ctx = contextlib.closing(fid)
1980 else:
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_datasource.py:192, in open(path, mode, destpath, encoding, newline)
155 """
156 Open `path` with `mode` and return the file object.
157
(...) 188
189 """
191 ds = DataSource(destpath)
--> 192 return ds.open(path, mode, encoding=encoding, newline=newline)
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/numpy/lib/_datasource.py:529, in DataSource.open(self, path, mode, encoding, newline)
526 return _file_openers[ext](found, mode=mode,
527 encoding=encoding, newline=newline)
528 else:
--> 529 raise FileNotFoundError(f"{path} not found.")
FileNotFoundError: task17/retention_profile_standard.txt not found.