How to calculate and export recycling coefficient during simulations?

User Support
1

Hello, everyone. I wonder, if it is possible to calculate and export the hydrogen recycling coefficient (R = Desorbed flux / Implanted flux) during simulations? Here is my code snippet:

import festim as F
import fenics as f
import numpy as np
import sympy as sp
from scipy import special

f_transient = 100
r_transient = 2.98e-1
r_en_transient = 1.40e-1
X_transient= 7.34e-8
sigma_transient= 3.82e-8
r_stat = 6.54e-1
r_en_stat = 4.47e-1
X_stat = 3.76e-9
sigma_stat = 2.03e-9
q_stat = 10e6
Gamma_stat = 8e22
q_transient_max = 4e8
T0 = 373
w_atom_density = 6.31e28
E = 9400
q = 1.6e-19
tau = 250e-6

def thermal_cond_function(T):
    return 149.441-45.466e-3*T+13.193e-6*T**2-1.484e-9*T**3+3.866e6/(T+1e-4)**2
def heat_capacity_function(T):
    return (21.868372+8.068661e-3*T-3.756196e-6*T**2+1.075862e-9*T**3+1.406637e4/(T+1.)**2) / 183.84e-3   
def q_tot(t):
    t = t % (1 / f_transient)
    return q_stat + q_transient_max*(1+(tau/t)**2)*(tau/t)**2*sp.exp(-(tau/t)**2)
def cooling(T, mobile):
    return -(q_stat + q_transient_max) * (1.179e-1 + 2.78e-3 * (T-T0) - 1.136e-1 * f.exp(np.log(0.94)*(T-T0)))
def flux_transient(t):
    t = t % (1 / f_transient)
    return q_transient_max*(tau/t)**2*sp.exp(-(tau/t)**2) / E / (1-r_en_transient) / q

my_model = F.Simulation(log_level = 40)

def create_x(list_x, n, dx):
    for i in range(n):
        list_x.append(x[-1]+dx)
    return list_x
x = [0]  
mesh = [
    [50, 1e-10],
    [95, 1e-9],
    [30, 1e-8],
    [36, 1e-7],
    [96, 1e-6],
    [20, 5e-6],
    [20, 1e-5],
    [32, 5e-5],
    [40, 1e-4]
] 

for i in mesh:
    x = create_x(x, i[0], i[1])
vertices = np.array(x)

my_model.mesh = F.MeshFromVertices(vertices=vertices)

tungsten = F.Material(id=1, 
                    D_0=1.97e-7/np.sqrt(2), 
                    E_D=0.2,
                    thermal_cond=thermal_cond_function, 
                    heat_capacity=heat_capacity_function,
                    rho=19250)

my_model.materials = tungsten

my_model.traps = [F.Trap(
        k_0=1.97e-7/(1.1e-10**2*6*w_atom_density)/np.sqrt(2),
        E_k=0.2,
        p_0=1e13,
        E_p=1.0,
        density=1e-4*w_atom_density,
        materials=tungsten
)]

stat_source = F.ImplantationFlux(
    flux=Gamma_stat, 
    imp_depth=X_stat, 
    width=sigma_stat, 
    volume=1
)

transient_source = F.ImplantationFlux(
    flux=flux_transient(F.t),
    imp_depth=X_transient,  # m
    width=sigma_transient,  # m
    volume=1
)

my_model.sources = [stat_source, transient_source]

my_model.boundary_conditions = [
    F.FluxBC(value=q_tot(F.t), field="T", surfaces=1),
    F.CustomFlux(function=cooling, field="T", surfaces=2),
    F.DirichletBC(surfaces=[1,2], value=0, field=0)
]

absolute_tolerance_T = 1.0
relative_tolerance_T = 1e-5
absolute_tolerance_c = 1e10
relative_tolerance_c = 1e-10
max_iter = 500

my_model.T = F.HeatTransferProblem(
    absolute_tolerance=absolute_tolerance_T,
    relative_tolerance=relative_tolerance_T,
    initial_value=T0,
    maximum_iterations=max_iter,
)

in_ts = 1e-5

my_model.dt = F.Stepsize(
    initial_value=in_ts,
    stepsize_change_ratio=1.05,
    t_stop = 1e-5,
    stepsize_stop_max = 5e-5,
    dt_min=1e-8
)

my_model.settings = F.Settings(
    absolute_tolerance=absolute_tolerance_c,
    relative_tolerance=relative_tolerance_c,
    final_time=100,
    soret=True,
    maximum_iterations=max_iter,
    
)

results_folder = "./results"
derived_quantities = [F.DerivedQuantities([F.HydrogenFlux(surface=1), 
                                          F.TotalSurface(field='T', surface=1), 
                                          F.TotalVolume(field='retention', volume=1),
                                          F.TotalVolume(field='solute', volume=1),
                                          F.TotalVolume(field="1", volume=1)],
                                          nb_iterations_between_compute = 1,
                                          filename=results_folder + "derived_quantities.csv")]

my_model.exports = derived_quantities
my_model.initialise()
my_model.run()

Thus, there are two implantation fluxes (stationaty and transient) in the 1D model with Dirichlet boundary conditions on two surfaces.

I thought of creating a custom F.DerivedQuantity to export, based on the F.HydrogenFlux class in the following way:

class RecyclingCoefficient(F.HydrogenFlux):
    def __init__(self, surface) -> None:
        super().__init__(surface)
        self.title = "Recycling coefficient surface {}".format(surface)

    def compute(self):
        des_flux = np.abs(super().compute())
        impl_flux = 0
        for source in my_model.sources:
            impl_flux += source.flux
        return des_flux / impl_flux

but I am not sure. Maybe, there is a better way to do it?

2

Hi @VVKulagin I think this can very well work!

It might need adapting though since source.flux may not be a float. In your case, I think it is a sympy.Expression. You would need to evaluate this expression by replacing the symbol t by the value of the current time. Consider this.

To do so, in addition to this custom RecyclingCoefficient class, you could rewrite the DerivedQuantities.compute method:

class CustomDerivedQuantities(F.DerivedQuantities):
    def compute(self, t):
        super().compute(t)
        for quantity in self.exports:
            if isinstance(quantity, RecyclingCoefficient):
                value = quantity.compute(t)
        self.data[-1].append(value)  # append to the global data list
1 Like
3

Much obliged. It works :slightly_smiling_face:

1 Like