Hydrogen transport: advanced

This section discusses how to implement advanced boundary conditions (BCs) for hydrogen transport problems.

Objectives:

• Choose solubility laws (Henry’s or Siervert’s)

• Implement surface reaction boundary conditions

• Model isotopic exchange as a recombination flux

Understanding background for solubility laws and surface reactions

Solubility laws

A solubility law prescribes the equilibrium relationship between the hydrogen concentration in a solid and the surrounding environment. Two common solubility laws for hydrogen are Henry’s law and Sieverts’ law:

Henry’s law

Henry’s law assumes a linear relationship between the hydrogen concentration in the solid and the partial pressure of hydrogen in the gas:

\[ c_s = k_H \, P_H \]

where:

• \(c_s\) is the surface concentration,

• \(P_H\) is the hydrogen partial pressure in the surrounding environment,

• \(k_H\) is the Henry’s law constant (material and temperature-dependent).

Tip

Users interested in molten salt interactions with hydrogen will usually need to use Henry’s Law.

Sieverts’ law

Sieverts’ law applies when hydrogen dissolves in metals as atomic hydrogen. The solubility is proportional to the square root of the hydrogen partial pressure (due to diatomic hydrogen dissociating into atoms in the solid):

\[ c_s = k_S \, \sqrt{P_H} \]

where:

• \(k_S\) is the Sieverts’ constant (which depends on temperature and the metal)

Tip

Users interested in liquid metal interactions with hydrogen will usually need to use Sievert’s Law.

Imposing a solubility law

Users can define the surface concentration using either Sieverts’ law. For Sieverts’ law of solubility, we can use SievertsBC:

import festim as F

boundary = F.SurfaceSubdomain(id=1)
H = F.Species(name="Hydrogen")

custom_pressure_value = lambda t: 2 + t
my_sieverts_bc = F.SievertsBC(subdomain=3, S_0=2, E_S=0.1, species=H, pressure=custom_pressure_value)

/home/docs/checkouts/readthedocs.org/user_builds/festim-workshop/conda/latest/lib/python3.12/site-packages/festim/coupled_heat_hydrogen_problem.py:1: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  import tqdm.autonotebook

Similarly, for Henry’s law of solubility, we can use HenrysBC:

pressure_value = lambda t: 5 * t
my_henrys_bc = F.HenrysBC(subdomain=3, H_0=1.5, E_H=0.2, species=H, pressure=pressure_value)

To use these BCs in your simulation, add them to boundary_conditions:

my_model = F.HydrogenTransportProblem()
my_model.boundary_conditions = [my_sieverts_bc, my_henrys_bc]

Note

Using HenrysBC or SievertsBC is just a convenient way for users to define a FixedConcentrationBC, using the solubility law to calculate the concentration.

---

Simple surface reaction

Surface reactions between adsorbed species and gas-phase products can be represented using the SurfaceReactionBC class. This boundary condition defines a reversible reaction of the form:

\[ A + B \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad C \]

where:

• \(\mathrm{A}, \mathrm{B}\) are surface reactants

• \(\mathrm{C}\) is the product species in the gas phase

• \(\mathrm{K}_r\) and \(\mathrm{K}_d\) are forward and backward reaction rates, respectively

Mathematical formulation

Let:

• \(c_A, c_B\) be the surface concentrations of \(\mathrm{A}\) and \(\mathrm{B}\)

• \(P_C\) be the partial pressure of \(\mathrm{C}\)

• \(T\) the surface temperature

• \(k_B\) the Boltzmann constant

The forward and backward reaction rates follow Arrhenius laws:

\[ K_r = k_{r0} \exp\!\left(-\frac{E_{kr}}{k_B T}\right), \qquad K_d = k_{d0} \exp\!\left(-\frac{E_{kd}}{k_B T}\right) \]

The net surface reaction rate is:

\[ K = K_r \, c_A c_B - K_d \, P_C \]

The resulting flux of species \(\mathrm{A}\) into the surface is:

\[ \mathbf{J}_A \cdot \mathbf{n} = K \]

If \(\mathrm{A=B}\), the total flux becomes:

\[ \mathbf{J}_A \cdot \mathbf{n} = 2 K \]

where \(\mathbf{n}\) is the outward normal vector at the surface.

Recombination and dissociation

If \(\text{Species A = Species B}\), recombination or dissociation can also be modeled with SurfaceReactionBC, e.g.:

\[ \mathrm{H + H} \quad \overset{K_r}{\underset{K_d}{\rightleftharpoons}} \quad \mathrm{H_2} \]

where

\[ K = K_r \, c_H^2 - K_d \, P_{H_2} \]

and corresponding flux of atomic hydrogen is:

\[ \mathbf{J}_H \cdot \mathbf{n} = - 2 K \]

Implementing surface reaction boundary conditions

Users can impose a surface reaction on a boundary using SurfaceReactionBC.

First, create a boundary to impose the BC and the reactant species:

import festim as F

boundary = F.SurfaceSubdomain(id=1)
A = F.Species("A")
B = F.Species("B")

Now, we add these as arguments to the SurfaceReactionBC class. We’ll also need to assign a gas_pressure (corresponding to the partial pressure of the product species), and the forward/backward rate (k_r0, k_d0) and energy (E_kr, E_kd) coefficients (see above to learn more about these parameters):

my_bc = F.SurfaceReactionBC(
    reactant=[A, B],
    gas_pressure=1e5,
    k_r0=1,
    E_kr=0.1,
    k_d0=0,
    E_kd=0,
    subdomain=boundary,
)

Finally, add the BC to your problem’s boundary_conditions attribute using a list:

my_model = F.HydrogenTransportProblem()
my_model.boundary_conditions = [my_bc]

Note

Using SurfaceReactionBC is just a convenient way for users to define surface fluxes, which can be done manually (as shown below in the isotopic exchange example).

---

Isotopic exchange

Isotopic exchange occurs when hydrogenic isotopes swap positions, for example:

(1)\[\mathrm{T + T} \rightleftharpoons \mathrm{T_2}\]

(2)\[\mathrm{T + H_2} \rightleftharpoons \mathrm{H} + \mathrm{HT}\]

Mathematical formulation

Let:

• \(c_T\) be the surface concentrations of tritium

• \(P_\mathrm{H_2}, P_\mathrm{HT}\) the partial pressures of \(\mathrm{H_2}\) and \(\mathrm{HT}\)

• \(K_{r}\) the rate coefficient of reaction (1)

• \(K_{r}^\ast\) the rate coefficient of reaction (2)

\[ K_{r} = K_{r0} \exp\!\left(-\frac{E_{Kr}}{k_B T}\right) \]

\[ K_{r}^\ast = K_{r0}^\ast \exp\!\left(-\frac{E_{Kr}^\ast}{k_B T}\right) \]

If \(c_{H_2} \gg c_{HT}\), the flux of tritium is approximated as:

\[ \phi_T = \phi_\mathrm{recombination} + \phi_\mathrm{exchange} \]

with $\( \phi_\mathrm{recombination} = - K_{r} c_T^2 \\ \phi_\mathrm{exchange} = - K_{r}^\ast P_{H_2} c_T \)$

These fluxes can be implemented in FESTIM using ParticleFluxBC with user-defined expressions for each reaction term, as shown below.

Note

\(K_{r0}^\ast\) and \(K_{r0}\) have different dimensions since \(P_\mathrm{H_2}\) is a pressure.

Note

This example is not tracking the dissolution and transport of protium in the metal. And assumes that the partial presure of \(\mathrm{H_2}\) is constant.

Note

For more complex isotopic exchanges, we can also use SurfaceReactionBC add other reactions. See modeling isotopic exchange for multiple species to learn more.

Implementation

Let’s work through a 1D example to illustrate the effect of isotopic exchange. We’ll consider tritium diffusion from left to right, with a large partial pressure of \(H_2\) at the right boundary (where isotopic exchange can occur).

First, we’ll set up our mesh and materials:

import numpy as np

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

right_subdomain = F.SurfaceSubdomain1D(id=2, x=1)
left_subdomain = F.SurfaceSubdomain1D(id=3, x=0)

material = F.Material(D_0=1e-2, E_D=0)

vol = F.VolumeSubdomain1D(id=5, material=material, borders=[0, 1])
my_model.subdomains = [vol, left_subdomain, right_subdomain]

Isotopic exchange reactions can be modeled as user-defined expressions using ufl.

import ufl
import festim as F

Kr_0_custom = 1
E_Kr_custom = 0.0  # eV
P_H2 = 10  # assumed constant H2 concentration in

def isotopic_exchange_flux_func(c, T):
    return - (Kr_0_custom * ufl.exp(-E_Kr_custom / (F.k_B * T))) * P_H2 * c

Now we’ll add our species tritium to our model and define a custom flux to represent isotopic exchange. We use ParticleFluxBC and the custom flux function we defined above to define this BC on the right boundary:

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

isotopic_exchange_flux = F.ParticleFluxBC(
    value=isotopic_exchange_flux_func,
    subdomain=right_subdomain,
    species_dependent_value={"c": tritium},
    species=tritium,
)

t2_recomb = F.SurfaceReactionBC(
    reactant=[tritium, tritium],
    gas_pressure=0,
    k_r0=1e-22,
    E_kr=0,
    k_d0=0,
    E_kd=0,
    subdomain=right_subdomain,
)

For diffusion across the mesh, we also define a fixed concentration on the left surface:

my_model.boundary_conditions = [
    isotopic_exchange_flux,
    t2_recomb,
    F.FixedConcentrationBC(subdomain=left_subdomain, value=1e20, species=tritium),
]

We’ll export the flux using SurfaceFlux (see exporting derived quantities to learn more) and save this data to a variable data_with_H2. We also create a Profile1DExport to be able to plot the concentration profile later:

surface_flux = F.SurfaceFlux(field=tritium, surface=right_subdomain)

profile = F.Profile1DExport(field=tritium, subdomain=vol)

my_model.temperature = 300
my_model.settings = F.Settings(atol=1e-10, rtol=1e-8, stepsize=1, final_time=120)

my_model.exports = [surface_flux, profile]
my_model.initialise()
my_model.run()

data_with_H2 = surface_flux.data.copy()

Let’s plot the evolution of the tritium concentration profile:

Show code cell source

Hide code cell source

import matplotlib.pyplot as plt

import matplotlib.animation as animation
from IPython.display import HTML

# Turn off interactive plotting to prevent static display
plt.ioff()

# Create the figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# Set up the plot limits and labels
ax.set_ylabel("Tritium concentration")
ax.set_xlabel("Position")

# Initialize empty line objects for the animation
line_current, = ax.plot([], [], color="black", linewidth=2, zorder=1000)
lines_previous = []

def animate(frame):
    # Clear previous iteration lines
    for line in lines_previous:
        line.remove()
    lines_previous.clear()
    
    # Plot all previous iterations in light grey
    for i in range(frame):
        line_prev, = ax.plot(profile.x, profile.data[i], color="lightgrey", linewidth=0.8, alpha=0.6)
        lines_previous.append(line_prev)
    
    # Plot current iteration in black
    line_current.set_data(profile.x, profile.data[frame])
    
    ax.set_title(f"Time: {profile.t[frame]:.2f}")

    return [line_current] + lines_previous

# Create animation
anim = animation.FuncAnimation(
    fig, animate, frames=len(profile.data), interval=100, blit=False, repeat=True
)

# Close the figure to prevent static display
plt.close(fig)

# Display only the animation
HTML(anim.to_jshtml())

Once Loop Reflect

We can compare the surface flux to the case where we have no isotopic exchange by removing the isotopic exchange boundary condition and saving the results to data_no_H2:

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=left_subdomain, value=1e20, species=tritium),
    t2_recomb,
]

profile = F.Profile1DExport(field=tritium, subdomain=vol)
my_model.exports = [surface_flux, profile]
my_model.initialise()
my_model.run()
data_no_H2 = surface_flux.data.copy()

Without isotopic exchange the gradient is much lower, indicating that isotopic exchange helps enhance surface outgassing of tritium!

Show code cell source

Hide code cell source

# Create animation
anim = animation.FuncAnimation(
    fig, animate, frames=len(profile.data), interval=100, blit=False, repeat=True
)

# Close the figure to prevent static display
plt.close(fig)

# Display only the animation
HTML(anim.to_jshtml())

Once Loop Reflect

Let’s plot both data_with_H2 and data_no_H2 to see the temporal evolution of the surface flux for each case:

Show code cell source

Hide code cell source

import matplotlib.pyplot as plt

x = my_model.mesh.mesh.geometry.x[:, 0]
t = surface_flux.t
plt.plot(t, data_no_H2, label="No H_2")
plt.plot(t, data_with_H2,label="With H_2")

plt.xlabel("Time")
plt.ylabel("Surface Flux (right)")
plt.legend(reverse=True)
plt.show()

We see that the flux is higher with isotopic exchange, as we’d expect.