Implantation flux as a function of desorption flux

[KulaginVladimir]

Hi, All! I wonder, is it possible to implement a time-dependent implanatation flux, the value of which will depend on the desorption flux or on the diffusion flux at x = 0? In the case of zero H concentration at the surface it would look like: $S(x) = \alpha\cdot\left(\left.D(T)\frac{\partial c}{\partial x}\right)\right|_{x=0}\cdot f(x)$, where $\alpha < 1$, $f(x)~-$ implantation profile.

[RemDelaporteMathurin]

Hi @KulaginVladimir and welcome 🎉

There is no native support for adding a surface flux in the governing equations directly.
@jhdark do you know of a way to do this in fenics?

However, it is possible to have a time (and space) dependent implantation flux. Consider:

import festim as F

my_source = F.ImplantationFlux(
    flux=2*festim.t,
    imp_depth=5e-9,
    width=5e-9,
    volume=1,
)

Depending on your application case, it might be possible to evaluate the surface flux at each time step and correct the value of the implantation flux.

What is your application case?

[KulaginVladimir]

@RemDelaporteMathurin, thank you for the fast response. I intend to study the H retention under transient heat/particle fluxes. As I've already tested, the problem can be sucessfully solved with FESTIM and the results agree well with TMAP7.

As a next step, I'd like to analyse whether an additional flux from desorbed particles would affect the retention (assume that desorbed particles somehow dissociate/become ionised/thermalise in the near-surface plasma, e.g. divertor, and come back to the target). In the test-code snippet below (I think I cannot post the full code now), I apply the Dirichlet BCs for the mobile H concentration, therfore the desorption flux should be ($D\dfrac{\partial c}{\partial x}$).

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

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+1.e-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):
    tau = 250e-6
    q_max = 10e6
    return q_max*sp.exp(-(tau/t)**2)
def cooling(T, mobile):
    return -1.98e7 * (3.35e-3 * (T-373) + 1.2e-1 * (1 - f.exp(-6.9e-2*(T-373))))
def flux(t):
    q = 1.6e-19 #C
    q_max = 10e6
    E = 115
    return q_max/(E+13.6)/q*sp.exp(-(tau/t)**2)

w_atom_density = 6.31e28  # atom/m3

my_model = F.Simulation()

vertices = np.concatenate([
    np.linspace(0, 5e-8, num=500),
    np.linspace(5e-8, 1e-5, num=200),
    np.linspace(1e-5, 1e-4, num=50),
    np.linspace(1e-4, 6e-3, num=20),])

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

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

my_model.materials = tungsten

trap = F.Trap(
        k_0=1.97e-7/(1.1e-10**2*w_atom_density),
        E_k=0.2,
        p_0=1e13,
        E_p=1.5,
        density=1e-5*w_atom_density,
        materials=tungsten
    )

my_model.traps = [trap]

impl_source = F.ImplantationFlux(
    flux=flux(F.t),  # H/m2/s
    imp_depth=25e-10,  # m
    width=15e-10,  # m
    volume=1
)

my_model.sources = [impl_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="solute")
]

absolute_tolerance = 1e-4
relative_tolerance = 1e-6

my_model.T = F.HeatTransferProblem(
    absolute_tolerance=absolute_tolerance,
    relative_tolerance=relative_tolerance,
    transient=True,
    initial_value=373,
    maximum_iterations=200
)

my_model.dt = F.Stepsize(
    initial_value=1e-5,
    stepsize_change_ratio=1.01,
    t_stop = 0.,
    stepsize_stop_max = 1e-4,
    dt_min=1e-7
)

my_model.settings = F.Settings(
    absolute_tolerance=absolute_tolerance,
    relative_tolerance=relative_tolerance,
    transient=True,
    final_time=100,
)

results_folder = "./"
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 = 10,
                                         filename=results_folder + "derived_quantities_"+str(sys.argv[1])+".csv")
my_model.exports = [derived_quantities]
my_model.initialise()
my_model.run()

[SamueleMeschini]

In Fenics you might try to add the term related to your boundary condition in the variational formulation, i.e. your bilinear form would look like (with the BC applied only to the left boundary):

$a(u,v) = \int_\Omega \nabla(u) \cdot \nabla(v) dx - \int_{\partial \Gamma_{left}} \alpha D \nabla (u)v f(0)\ ds$

So the code for the bilinear form would be:

a = inner(grad(u), grad(v))*dx - \grad(u)*v*f*ds(1)

where you want to compute the second term on one of the boundary only. It should work but I am not 100% sure. The linear form should remain the same with the source term defined as $f(x)$.

Alternatively, it may be possible to do something similar to the convective or mass transfer flux BC if one can compute the gradient of $c$ within the Flux BC class - in this case there would be no need to modify the variational formulation.

[RemDelaporteMathurin]

As a next step, I'd like to analyse whether an additional flux from desorbed particles would affect the retention

Oh I think I understand. So:

1. particles are implanted
2. the retro-desorb at the surface, form a molecule
3. they are re-ionised and reimplanted (what energy? i suppose that would be fairly low compared to the initial incident energy?)

Is that correct?

Your implantation depth is very small (less than a nanometer). Under this regime, everything equilibrates very quickly and the retro-desorbed flux is almost equal to the implantation flux (see Eq 2.9 to 2.11 of this thesis).
So maybe you can use this to your advantage and do something like:

impl_source = F.ImplantationFlux(
    flux=flux(F.t),  # H/m2/s
    imp_depth=25e-10,  # m
    width=15e-10,  # m
    volume=1
)

alpha = 0.05  # 5% of the desorbed particles are reimplanted
re_implanted_source = F.ImplantationFlux(
    flux=alpha * flux(F.t),  # H/m2/s
    imp_depth=...,  # depends on the energy
    width=...,  # depends on the energy
    volume=1
)

my_model.sources = [impl_source, re_implanted_source]

Now this assumes that the stage 3 of this recirculation process is negligible (ie the particles that are reimplanted can also desorb but the contribution of these is negligible)

Hope this helps

[RemDelaporteMathurin]

So the code for the bilinear form would be:

a = inner(grad(u), grad(v))*dx - \grad(u)*v*f*ds(1)

@SamueleMeschini this adds a flux on surface 1 based on the source term f, correct?

I think we are looking for the opposite: adding a volumetric source, which value is based on D grad(u)

Something like:
F = inner(grad(u), grad(v))*dx - assemble(-D * \grad(u) * n * ds(1) )*f*v*dx
with n the facet normal?

But I don't know if this is possible. Maybe a question for the FEniCS discourse!

[SamueleMeschini]

Correct, I misunderstood the question.

Something like: F = inner(grad(u), grad(v))*dx - assemble(-D * \grad(u) * n * ds(1) )*f*v*dx with n the facet normal?

This should definitely work!

[KulaginVladimir]

@RemDelaporteMathurin

Oh I think I understand. So:

1. particles are implanted
2. the retro-desorb at the surface, form a molecule
3. they are re-ionised and reimplanted (what energy? i suppose that would be fairly low compared to the initial incident energy?)

Is that correct?

Yeah, you are right. Though I don't know exactly what energy such particles could have, but definitely low.

Your implantation depth is very small (less than a nanometer). Under this regime, everything equilibrates very quickly and the retro-desorbed flux is almost equal to the implantation flux (see Eq 2.9 to 2.11 of this thesis).

Right, this is a kind of solution. I assume it can be generalised for two implantation fluxes (low energy: ~100 eV, high energy: ~1keV). Thus, I will try the suggestion. However, it is not the case for small irradiation times and it does not allow to precisely estimate time-dependece of the recycling coefficient ($\Gamma_{des}/\Gamma_{impl}$).

[RemDelaporteMathurin]

@KulaginVladimir You are correct this approximation is valid once the fluxes have reached equilibrium.

In that case, maybe something along these lines would work:

import festim as F
import fenics as f
import sympy as sp
import numpy as np


class CustomImplantationFlux(F.Source):
    def __init__(self, imp_depth, width, volume):
        self.imp_depth = imp_depth
        self.width = width

        # gaussian distribution
        distribution = (
            1
            / (self.width * (2 * np.pi) ** 0.5)
            * sp.exp(-0.5 * ((F.x - self.imp_depth) / self.width) ** 2)
        )

        # flux is a new symbol
        flux = sp.Symbol("flux")

        # create full sympy expression for value
        value = flux * distribution

        # call parent class constructor
        super().__init__(0, volume, field="0")

        # override value attribute
        self.value = f.Expression(sp.printing.ccode(value), t=0, flux=0, degree=2)


class RecyclingSimulation(F.Simulation):
    def iterate(self):
        # call parent class iterate method
        super().iterate()

        # compute outgassing flux
        outgassing_flux_value = outgassing_flux_export.compute()

        # new flux value is 5% of outgassing flux
        new_flux_value = 0.05 * np.abs(outgassing_flux_value)

        # update implanatation flux value
        source_from_recycling.value.flux = new_flux_value


my_model = RecyclingSimulation()

implantation_flux = F.ImplantationFlux(
    flux=flux(F.t), imp_depth=5e-9, width=2e-9, volume=1
)
source_from_recycling = F.CustomImplantationFlux(imp_depth=5e-9, width=2e-9, volume=1)

my_model.sources = [implantation_flux, source_from_recycling]

outgassing_flux_export = F.SurfaceFlux("solute", 1)
my_model.exports = [F.DerivedQuantities(outgassing_flux_export)]

...

my_model.initialise()
my_model.run()

What we do here is:

• Create a custom source class inheriting from ImplantationFlux which builds the C++ code with a flux variable.
• Create a custom simulation class inheriting from Simulation which behaves exactly the same way except that after each iteration, it evaluates the outgassing flux and then corrects the flux value in the source_from_recycling

I haven't tested it but I don't see why this wouldn't work.

However, this explicit way of updating the source requires that the timestep remains small enough.

Let us know!

EDIT: made a few fixes to make sure it works

[KulaginVladimir]

@RemDelaporteMathurin, it definetely works! Thanks for the operative feedback.

Some comparison for the case of cyclic loads: