Basic functionality

Basic functionality#

Materials are key components of hydrogen transport simulations. They hold the properties like diffusivity, solubility and even thermal properties like thermal conductivity or heat capacity. Read more about the Materials class and syntax at Materials.

Objectives:

Learn how to define material properties (thermal, solubility, diffusivity)
Learn how to define materials on different subdomains
Solve a multi-material diffusion problem

Defining material properties#

We can define a material using the diffusion exponential pre-factor \(D_0\) and activation energy \(E_D\). By default, FESTIM assumes they follow an Arrhenius law, which is of the following form:

\[ D = D_0 \exp{(-E_D/k_B T)} \]

where \(k_B\) is the Boltzmann constant in eV/K and \(T\) is the temperature in K. To define a material using these two properties:

import festim as F

mat = F.Material(D_0=1.11e-6, E_D=0.4)  # m2/s, eV

/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

When considering chemical potential conservation at material interfaces, hydrogen solubility can be defined using the solubility coeffeicient prefactor K_S_0, solubility activation energy E_K_S, and solubility law solubility_law (either "henry" or "sievert"):

mat.K_S_0 = 1.0
mat.E_K_S = 3.0
mat.solubility_law = "sievert"

Defining thermal properties#

Users can define thermal properties, such as thermal conductivity, heat capacity, and density for their materials. The simplest way to define them:

mat.thermal_conductivity = 10.0  # W/m/K
mat.heat_capacity = 500.0  # J/kg/K
mat.density = 8000.0  # kg/m3

Users can also define these as a function of temperature:

import ufl

mat.thermal_conductivity = lambda T: 3*T + 2*ufl.exp(-20*T)
mat.heat_capacity = lambda T: 4*T + 8
mat.density = lambda T: 7*T + 5

Defining materials on different subdomains#

Volume subdomains are used to assign different materials or define regions with specific physical properties. Each volume subdomain must be associated with a festim.Material object. Read more about subdomains here.

Consider the following volume with two subdomains separated halfway through the mesh:

To define one material on each subdomain:

mat = F.Material(D_0=1.11e-6, E_D=0.4)  # m2/s, eV

top = F.VolumeSubdomain(id=1, material=mat, locator=lambda x: x[0] >= 0.5)
bottom = F.VolumeSubdomain(id=2, material=mat, locator=lambda x: x[0] < 0.5)

Similarly, for two materials:

mat1 = F.Material(D_0=1.11e-6, E_D=0.4)  # m2/s, eV
mat2 = F.Material(D_0=2e-6, E_D=0.3)  # m2/s, eV

top = F.VolumeSubdomain(id=1, material=mat1, locator=lambda x: x[0] >= 0.5)
bottom = F.VolumeSubdomain(id=2, material=mat2, locator=lambda x: x[0] < 0.5)

Multi-material example#

Considering the following 2D example, where hydrogen diffuses through a 2D domain composed of two materials with different diffusion and solubility properties. The top half (Material A) has a lower diffusion coefficient and solubility than the bottom half (Material B). The interface at \(𝑦=0.5\) clearly separates the two materials, and the steady-state hydrogen distribution illustrates how material properties impact transport.

First, we create the mesh for our discontinuous (materials have different solubility properties) problem. Note that we use HydrogenTransportProblemDiscontinuous to define our multi-material problem:

import festim as F
import numpy as np
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import pyvista
from dolfinx import plot

fenics_mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
festim_mesh = F.Mesh(fenics_mesh)

my_model = F.HydrogenTransportProblemDiscontinuous()

Note

Be sure to use an even number of cells in each direction when creating the mesh in a discontinuous problem.

Next, we define solubility properties. If the solubilties were the same, we’d expect to see a smooth “continous” concentration profile. For this problem, we have different solubilities:

# Top material (Material A)
material_top = F.Material(
    name="Material_A",
    D_0=1,    # m²/s
    E_D=0,     # eV
    K_S_0=2,    # mol/m³/Pa  (solubility prefactor)
    E_K_S=0,     # eV (activation energy for solubility)
)

# Bottom material (Material B)
material_bottom = F.Material(
    name="Material_B",
    D_0=2,    # m²/s
    E_D=0,     # eV
    K_S_0=3,    # mol/m³/Pa
    E_K_S=0,     # eV
)

Now we must assemble the subdomains for our volumes, surfaces, interfaces, and species. A penalty term is also used when defining the interface to ensure a well-behaved solution at the interface between both materials:

top_volume = F.VolumeSubdomain(id=3, material=material_top, locator=lambda x: x[1] >= 0.5)
bottom_volume = F.VolumeSubdomain(id=4, material=material_bottom, locator=lambda x: x[1] <= 0.5)

top_surface = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1.0))
bottom_surface = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[1], 0.0))

my_model.mesh = festim_mesh
my_model.subdomains = [top_surface, bottom_surface, top_volume, bottom_volume]

my_model.interfaces = [F.Interface(5, (bottom_volume, top_volume), penalty_term=1000)]
my_model.surface_to_volume = {
    top_surface: top_volume,
    bottom_surface: bottom_volume,
}

H = F.Species("H")
my_model.species = [H]

for species in my_model.species:
    species.subdomains = [bottom_volume, top_volume]

Finally, we define the temperature, boundary conditions, and then solve:

my_model.temperature = 400

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=top_surface, value=1.0, species=H),
    F.FixedConcentrationBC(subdomain=bottom_surface, value=0.0, species=H),
]

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

my_model.initialise()
my_model.run()

To visualize the results of the multi-material problem, we need to look at each subdomain separately. We can use H.subdomain_to_post_processing_solution to define a plotting function make_ugrid that visualizes each domain:

def make_ugrid(solution):
    topology, cell_types, geometry = plot.vtk_mesh(solution.function_space)
    u_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    u_grid.point_data["c"] = solution.x.array.real
    u_grid.set_active_scalars("c")
    return u_grid

pyvista.set_jupyter_backend("html")

u_plotter = pyvista.Plotter()
u_grid_top = make_ugrid(H.subdomain_to_post_processing_solution[top_volume])
u_grid_bottom = make_ugrid(H.subdomain_to_post_processing_solution[bottom_volume])
u_plotter.add_mesh(u_grid_top, cmap="magma", show_edges=False)
u_plotter.add_mesh(u_grid_bottom, cmap="magma", show_edges=False)
u_plotter.view_xy()
u_plotter.add_text("Hydrogen concentration in multi-material problem", font_size=12)

if not pyvista.OFF_SCREEN:
    u_plotter.show()
else:
    figure = u_plotter.screenshot("concentration.png")

2026-05-22 15:21:35.501 (   2.118s) [    7F1B4205E740]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=

Success! We see diffusion from the top surface downwards, with a discontinuity at the interface!