Boundary conditions for HT in FESTIM2

User Support
1

Hi,

I am trying to figure out how to apply boundary conditions in case of mixed isotopologues. If we have H2 and HT in gas, each with own partial pressure - is there a way to apply boundary conditions to find H and T concentrations under the surface of metal? As far as I remember, concentration of T should be proportional to P(HT), not square root. I could not find this in readthedocs, but this might be the usual case of me missing something. :slight_smile:

Best wishes,

Mikhail

2

Yes this is definitely possible.

You can use the SurfaceReactionBC to model this, where you can set the gas partial pressure for H2 and HT separately.

Here is an example in action:

import festim as F
import numpy as np

my_model = F.HydrogenTransportProblem()

protium = F.Species("H")
tritium = F.Species("T")
my_model.species = [protium, tritium]

my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100))

left_surf = F.SurfaceSubdomain1D(id=1, x=0)
right_surf = F.SurfaceSubdomain1D(id=2, x=1)

# assumes the same diffusivity for all species
material = F.Material(D_0={protium: 1e-02, tritium: 8e-03}, E_D={protium: 0.0, tritium: 0.0})

vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material)

my_model.subdomains = [vol, left_surf, right_surf]

my_model.boundary_conditions = [
    # Protium BCs
    F.SurfaceReactionBC(
        reactant=[protium, protium],
        gas_pressure=100,
        k_r0=0,
        E_kr=0,
        k_d0=1,
        E_kd=0,
        subdomain=left_surf,
    ),
    F.SurfaceReactionBC(
        reactant=[protium, tritium],
        gas_pressure=100,
        k_r0=0,
        E_kr=0,
        k_d0=0.8,
        E_kd=0,
        subdomain=left_surf,
    ),
]

my_model.temperature = 300

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=100)

my_model.settings.stepsize = F.Stepsize(1)

my_model.exports = [
    F.VTXSpeciesExport(filename="protium_profile.bp", field=protium),
    F.VTXSpeciesExport(filename="tritium_profile.bp", field=tritium),
]

my_model.initialise()
my_model.run()

Obviously, in this example, H comes from 2 sources: H2 and HT, and T comes only from HT. We also have different diffusion coefficients. However, a limitation is that it assumes the gas pressures remain constant, even if recombination values are included.

Hope this helps!

3

I believe the gas pressure can be a function of space and time by giving a lambda function no?

4

Yes, you’re right, sorry, just meant to say that, any species leaving the system due to recombination, for instance, will not affect the pressure value

5

Hi James,

I modified your example a bit to output steady state total gas content:

import festim as F
import numpy as np

my_model = F.HydrogenTransportProblem()

protium = F.Species(“H”)
tritium = F.Species(“T”)
my_model.species = [protium, tritium]

my_model.mesh = F.Mesh1D(np.linspace(0, 1, 100))

left_surf = F.SurfaceSubdomain1D(id=1, x=0)
right_surf = F.SurfaceSubdomain1D(id=2, x=1)

assumes the same diffusivity for all species

material = F.Material(D_0={protium: 1e-02, tritium: 8e-03}, E_D={protium: 0.0, tritium: 0.0})
my_model.materials = [material]
vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=material)

my_model.subdomains = [vol, left_surf, right_surf]

my_model.boundary_conditions = [

Protium BCs

F.SurfaceReactionBC(
reactant=[protium, protium],
gas_pressure=1,
k_r0=0,
E_kr=0,
k_d0=1,
E_kd=0,
subdomain=left_surf,
),
F.SurfaceReactionBC(
reactant=[protium, tritium],
gas_pressure=0,
k_r0=0,
E_kr=0,
k_d0=0.8,
E_kd=0,
subdomain=left_surf,
),
F.FixedConcentrationBC(subdomain = right_surf, value = 0, species = protium),
F.FixedConcentrationBC(subdomain = right_surf, value = 0, species = tritium),
]

my_model.temperature = 300

my_model.settings = F.Settings(transient=False, atol=1e-10, rtol=1e-10)

my_model.settings.stepsize = F.Stepsize(1)

total_H_conc = F.TotalVolume(
field = protium,
volume = vol,
filename = “./total_H_conc.csv”
)
total_T_conc = F.TotalVolume(
field = tritium,
volume = vol,
filename = “./total_T_conc.csv”
)

my_model.exports = [total_H_conc, total_T_conc]

my_model.initialise()
my_model.run()

Now with H2 pressure=1 and HT pressure=0, total H content is 100; if H2 pressure is increased to 100, total H content also increases 100 times to 10000 - is this correct? :crying_cat:

6

A very gentle reminder. :slight_smile:

7

hi @miklavrent the code block is weirdly formatted
The flux on the surface with SurfaceReactionBC for H (with no HT) is \varphi = -D \nabla c \cdot n= 2 K_d P_\mathrm{H_2} and the concentration is set to zero on the other side.

In steady state, the concentration profile should be linear in the volume c(x) = ax+b.

You can show analytically that c = a (x-1) in your case with a=\frac{-2K_d P_\mathrm{H_2}}{D}
So c \propto -a \propto P_\mathrm{H_2}.

The quantity you are computing is nothing but the space integral \int_\Omega c(x) dx so it makes sense that multiplying the pressure by 100 also mulitplies the inventory by 100.

8

Hi Remi, sorry for the formatting!

Indeed, you are absolutely correct: in the limiting case K_r = 0, Sieverts’ law does not hold. Not sure what this means for reaction with two isotopes…

9

Sievert’s law doesn’t hold if K_r = 0 no. The definition of Sievert’s constant is the square root of the ratio K_d/K_r. It makes even less sense for multiple isotopes and more isotopologues. For this sort of cases, one should only rely on fluxes and surface reactions