Hello! I’m not very familiar with python. Using sp.Piecewise, how can I set the trap distribution with these conditions: layer of oxide (10 nm), layer of steel (2 mm), layer of oxide (10 nm), and the trap is located on both sides of the steel with 1 μm trap layer thickness? I think trap distribution should be in the coordinates described by these inequalities:
(10 nm < x < 10 nm + 1 μm) & (10 nm + 1 μm + 2 mm < x < 10 nm + 1 μm + 2 mm + 1 μm),
if in FESTIM the coordinate x zero point is considered to be on the left side of the first (oxide) layer. Unfortunately, I don’t know how set this inequalities in sp.Piecewise
I know that I can compute only half of the problem due to symmetry but I’m interested in the general principle of setting such complicated trap distribution
If I understand your problem you have 3 materials (2 oxide layers and a steel substrate) and you have a trap that is located only in the steel layer but in the first micron (on each side).
Consider the following steady state example:
import festim as F
import numpy as np
import sympy as sp
my_model = F.Simulation()
L_oxide_1 = 10e-9
L_steel = 2e-3
L_oxide_2 = 10e-9
# Define the mesh with local refinement around the trap layers
my_model.mesh = F.MeshFromVertices(
np.concatenate(
[
np.linspace(0, L_oxide_1, 20),
np.linspace(L_oxide_1, L_oxide_1 + 1.5e-6, 300),
np.linspace(L_oxide_1 + 1.5e-6, L_oxide_1 + L_steel - 1.5e-6, 100),
np.linspace(L_oxide_1 + L_steel - 1.5e-6, L_oxide_1 + L_steel, 300),
np.linspace(L_oxide_1 + L_steel, L_oxide_1 + L_steel + L_oxide_2, 20),
]
)
)
oxide_1 = F.Material(id=1, D_0=1, E_D=0, borders=[0, L_oxide_1])
steel = F.Material(id=2, D_0=1, E_D=0, borders=[L_oxide_1, L_oxide_1 + L_steel])
oxide_2 = F.Material(
id=3, D_0=1, E_D=0, borders=[L_oxide_1 + L_steel, L_oxide_1 + L_steel + L_oxide_2]
)
my_model.materials = [oxide_1, steel, oxide_2]
trap_density = 1 # trap/m3
distribution = sp.Piecewise(
(1, sp.And(L_oxide_1 < F.x, F.x < L_oxide_1 + 1e-6)), # 1 if L_oxide_1 < x < L_oxide_1 + 1e-6
(1, sp.And(L_oxide_1 + L_steel - 1e-6 < F.x, F.x < L_oxide_1 + L_steel)), # 1 if L_oxide_1 + L_steel - 1e-6 < x < L_oxide_1 + L_steel
(0, True), # 0 otherwise
)
my_model.traps = F.Trap(
k_0=1,
E_k=0,
p_0=1e-10,
E_p=0,
materials=[steel], # trap in steel only
density=trap_density * distribution,
)
my_model.boundary_conditions = [F.DirichletBC(value=1, surfaces=[1, 2], field="solute")]
my_model.T = 500 # ignored here since activation energies are zero
my_model.settings = F.Settings(
absolute_tolerance=1e-10,
relative_tolerance=1e-10,
traps_element_type="DG",
transient=False,
)
my_model.exports = [F.XDMFExport(field="1", checkpoint=False)]
my_model.initialise()
my_model.run()
Reading the xdmf file in paraview you obtain this on the left side, you have a trapped concentration of 1 in the first micron of the steel:
Now I understand how it works! Thank you!
Hello again! I tried to run my code with complicated trap distribunion and specific parameters (such as oxide diffusicon energy, steel duffusion energy, etc.) but it seems like the solution doesn’t converge (the concentration oscillations is near to the right side of trap layer)
If you try to run code without traps you will get the same TDS-spectrum that I calculate in TMAP7.
With near-surface traps it looks like all the deuterium go out immmediately.
This is what a TDS spectrum with near-surface traps shoud look like (calculated in TMAP, results in TMAP7 and TMAP8 are the same):
Yes, I set traps_element_type="DG".
I would be grateful if you could point out the mistakes I made in FESTIM.
My FESTIM code
import festim as F
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
my_model = F.Simulation()
oxide_thick = 10e-9
steel_thick = 2e-3
oxide_left = F.Material(
id=1,
D_0=2e-8,
E_D=0.63,
borders=[0, oxide_thick]
)
rusfer = F.Material(
id=2,
D_0=0.85e-7,
E_D=0.165,
borders=[oxide_thick, steel_thick + oxide_thick]
)
oxide_right = F.Material(
id=3,
D_0=2e-8,
E_D=0.63,
borders=[steel_thick + oxide_thick, steel_thick + 2 * oxide_thick]
)
my_model.materials = [oxide_left, rusfer, oxide_right]
my_model.mesh = F.MeshFromVertices(
np.concatenate(
[
np.linspace(0, oxide_thick, 500),
np.linspace(oxide_thick, oxide_thick + 5e-6, 500),
np.linspace(oxide_thick + 5e-6, oxide_thick + steel_thick - 5e-6, 500),
np.linspace(oxide_thick + steel_thick - 5e-6, oxide_thick + steel_thick, 500),
np.linspace(oxide_thick + steel_thick, oxide_thick + steel_thick + oxide_thick, 500),
]
)
)
rusfer_atom_density = 8.21e28 # atom/m3
trap_density = 1.75e-2
distribution = sp.Piecewise(
(1, sp.And(oxide_thick < F.x, F.x < oxide_thick + 1e-6)), # 1 if L_oxide_1 < x < L_oxide_1 + 1e-6
(1, sp.And(oxide_thick + steel_thick - 1e-6 < F.x, F.x < oxide_thick + steel_thick)), # 1 if L_oxide_1 + L_steel - 1e-6 < x < L_oxide_1 + L_steel
(0, True), # 0 otherwise
)
trap_1 = F.Trap(
k_0=0.85e-7/(2.87e-10**2),
E_k=0.165,
p_0=1e13,
E_p=0.75,
# (10 nm < x < 10 nm + 1 mkm) и (10 nm + 1 mkm + 2 mm < x < 10 nm + 1 mkm + 2 mm + 1 mkm)
density= rusfer_atom_density * trap_density * distribution,
materials=rusfer
)
my_model.traps = trap_1
fraction_filled = 0.01
my_model.initial_conditions = [
F.InitialCondition(field="solute", value=0.85e24),
F.InitialCondition(field = "1", value = rusfer_atom_density* trap_density * distribution * fraction_filled),
]
my_model.boundary_conditions = [
F.RecombinationFlux(1e-29, -0.35, 2, surfaces = [1,2])
]
temperature_ramp = 2 # K/s
my_model.T = 293 + F.t*temperature_ramp
my_model.dt = F.Stepsize(
initial_value=1e-3,
stepsize_change_ratio=1.1,
max_stepsize = 2,
dt_min=1e-2,
milestones = [0, 753]
)
my_model.settings = F.Settings(
absolute_tolerance=1e10,
relative_tolerance=1e-09,
final_time=753,
traps_element_type="DG",
)
list_of_derived_quantities = [
F.TotalVolume("solute", volume=1),
F.TotalVolume("retention", volume=1),
F.TotalVolume("1", volume=1),
F.HydrogenFlux(surface=1),
F.HydrogenFlux(surface=2),
F.AverageVolume(field="T", volume=1)
]
derived_quantities = F.DerivedQuantities(
list_of_derived_quantities,
filename="TDS/EK-181.csv", # optional set a filename to export the data to csv
show_units=True,
)
my_model.exports = [derived_quantities]
#my_model.exports = [F.XDMFExport(field="1", checkpoint=False)]
my_model.initialise()
my_model.run()
t = derived_quantities.t
flux_left = derived_quantities.filter(fields="solute", surfaces=1).data
flux_right = derived_quantities.filter(fields="solute", surfaces=2).data
Temp = derived_quantities.filter(fields="T").data
flux_total = -np.array(flux_left) - np.array(flux_right)
plt.plot(Temp, flux_total, '.', linewidth=3)
plt.xlim(293, 1800)
# plt.ylim(bottom=-1.25e18, top=0.6e19)
plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Temperature, K")
plt.grid()
plt.show()
Hi @cruiseraurora so if I understand well you have a good agreement between TMAP7, TMAP8 and FESTIM but you want to add additional traps, and that’s not working (ie. oscillations in the concentration profile). Is this correct?
I have reproduced the issue locally. I believe this is a parametrisation issue, the trapping rate seems very high. Could you provide the parametrisation for TMAP7 with this near-surface trap please?
You can consider replacing the recombination BC by am homogeneous dirichlet BC (ie. c_mobile = 0) for simplicity (while debugging).
Hi @cruiseraurora!
Please, consider this code for a simplified problem (a half of the sample, a simpler mesh, and the updated value of the trapping rate).
Code
import festim as F
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
my_model = F.Simulation()
oxide_thick = 10e-9
steel_thick = 2e-3
oxide_left = F.Material(id=1, D_0=2e-8, E_D=0.63, borders=[0, oxide_thick])
rusfer = F.Material(
id=2, D_0=0.85e-7, E_D=0.165, borders=[oxide_thick, steel_thick / 2 + oxide_thick]
)
my_model.materials = [oxide_left, rusfer]
my_model.mesh = F.MeshFromVertices(
np.concatenate(
[
np.linspace(0, 5 * oxide_thick, 2000),
np.linspace(5 * oxide_thick, steel_thick / 2 + oxide_thick, num=2000),
]
)
)
rusfer_atom_density = 8.21e28 # atom/m3
trap_density = 1.75e-2
distribution = sp.Piecewise(
(
1,
sp.And(oxide_thick < F.x, F.x < oxide_thick + 1e-6),
), # 1 if L_oxide_1 < x < L_oxide_1 + 1e-6
(0, True), # 0 otherwise
)
trap_1 = F.Trap(
k_0=0.85e-7 / (2.87e-10**2 * rusfer_atom_density),
E_k=0.165,
p_0=1e13,
E_p=0.75,
# (10 nm < x < 10 nm + 1 mkm) и (10 nm + 1 mkm + 2 mm < x < 10 nm + 1 mkm + 2 mm + 1 mkm)
density=rusfer_atom_density * trap_density * distribution,
materials=rusfer,
)
my_model.traps = trap_1
fraction_filled = 0.01
my_model.initial_conditions = [
F.InitialCondition(field="solute", value=0.85e24),
F.InitialCondition(
field="1",
value=rusfer_atom_density * trap_density * distribution * fraction_filled,
),
]
my_model.boundary_conditions = [F.RecombinationFlux(1e-29, -0.35, 2, surfaces=1)]
temperature_ramp = 2 # K/s
my_model.T = 293 + F.t * temperature_ramp
my_model.dt = F.Stepsize(
initial_value=1e-3,
stepsize_change_ratio=1.1,
max_stepsize=5,
dt_min=1e-2,
milestones=[0, 753],
)
my_model.settings = F.Settings(
absolute_tolerance=1e12,
relative_tolerance=1e-8,
final_time=753,
traps_element_type="DG",
)
list_of_derived_quantities = [
F.TotalVolume("solute", volume=1),
F.TotalVolume("retention", volume=1),
F.TotalVolume("1", volume=1),
F.HydrogenFlux(surface=1),
F.HydrogenFlux(surface=2),
F.AverageVolume(field="T", volume=1),
]
derived_quantities = F.DerivedQuantities(
list_of_derived_quantities,
filename="TDS/EK-181.csv", # optional set a filename to export the data to csv
show_units=True,
)
TXT = [
F.TXTExport(
field="1",
filename="./trap1.txt",
write_at_last=True,
times=[753],
)
]
my_model.exports = [derived_quantities]
# my_model.exports = [F.XDMFExport(field="1", checkpoint=False)]
my_model.initialise()
my_model.run()
t = derived_quantities.t
flux_left = derived_quantities.filter(fields="solute", surfaces=1).data
# flux_right = derived_quantities.filter(fields="solute", surfaces=2).data
Temp = derived_quantities.filter(fields="T").data
flux_total = -2 * np.array(flux_left) # - np.array(flux_right)
plt.plot(Temp, flux_total, linewidth=3)
plt.xlim(293, 1800)
# plt.ylim(bottom=-1.25e18, top=0.6e19)
plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Temperature, K")
plt.grid()
plt.show()
It produces the following TDS spectrum:
1 Like
Thank you very much! I apologize for not noticing my mistake in the expression for trapping rate
1 Like