Temperature dependence

Temperature dependence#

A lot of hydrogen transport processes are thermally activated, meaning they are goverened by a property following an Arrhenius law. For instance, the diffusivity \(D\) can be written as

\[ D = D_0 \exp{\left(\frac{-E_D}{k_B ~T}\right)} \]

where \(D_0\) is called the pre-exponential factor, \(E_D\) the activation energy, \(T\) the temperature, and \(k_B\) the Boltzmann constant.

Non-constant temperature#

We will take the same example as above but setting a temperature dependence for \(D\), \(k\), and \(p\) and a non-constant temperature.

(24)#\[\begin{equation} T(t) = \begin{cases} 400 ~\mathrm{K} \quad \text{for } t < 15\\ 800 ~\mathrm{K} \quad \text{otherwise} \end{cases} \end{equation}\]

import dolfinx
import numpy as np
import ufl
from dolfinx.fem.petsc import NonlinearProblem
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
import basix
from scifem import assemble_scalar

nx = ny = 20

domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.quadrilateral)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)

cg_element = basix.ufl.element("Lagrange", domain.basix_cell(), degree=1)

mixed_element = basix.ufl.mixed_element([cg_element, cg_element])

V = fem.functionspace(domain, mixed_element)

u = fem.Function(V)
u_n = fem.Function(V)

cm, ct = ufl.split(u)
cm_n, ct_n = ufl.split(u_n)

v_cm, v_ct = ufl.TestFunctions(V)

def inlet(x):
    return np.logical_and(np.isclose(x[0], 0), x[1] <= 0.5)

def outlet(x):
    return np.logical_and(np.isclose(x[0], 1), x[1] >= 0.5)

V0, submap = V.sub(0).collapse()

dofs_outlet = fem.locate_dofs_geometrical((V.sub(0), V0), outlet)
dofs_inlet = fem.locate_dofs_geometrical((V.sub(0), V0), inlet)

c_inlet = fem.Constant(domain, 1.0)
c_outlet = fem.Constant(domain, 0.0)

bc_outlet = fem.dirichletbc(c_outlet, dofs_outlet[0], V.sub(0))
bc_inlet = fem.dirichletbc(c_inlet, dofs_inlet[0], V.sub(0))

T = dolfinx.fem.Constant(domain, 400.0)

def arrhenius_rate(k0, Ek, T, mod=ufl):
    boltzmann_constant = 8.617333262145e-5  # eV/K
    return k0 * mod.exp(-Ek/(T*boltzmann_constant))

k = arrhenius_rate(1e2, 0.2, T) # trapping rate
p = arrhenius_rate(5e9, 1.5, T)  # detrapping rate
n = 0.5  # total trapping sites
D = arrhenius_rate(1e3, 0.2, T) # diffusion coefficient

dt = dolfinx.fem.Constant(domain, 0.3)

F_mobile_transient = (cm - cm_n)/dt* v_cm * ufl.dx
F_trapped_transient = (ct - ct_n)/dt * v_ct * ufl.dx

trapping = k * cm * (n - ct)
detrapping = p * ct

F_mobile = (
    D*ufl.dot(ufl.grad(cm), ufl.grad(v_cm)) * ufl.dx
    + trapping * v_cm * ufl.dx
    - detrapping * v_cm * ufl.dx
)
F_trapped = -trapping * v_ct * ufl.dx + detrapping * v_ct * ufl.dx

F = F_mobile_transient + F_trapped_transient + F_mobile + F_trapped

# taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28
use_superlu = PETSc.IntType == np.int64  # or PETSc.ScalarType == np.complex64
sys = PETSc.Sys()  # type: ignore
if sys.hasExternalPackage("mumps") and not use_superlu:
    linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
    linear_solver = "superlu_dist"
else:
    linear_solver = "petsc"

petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2,
    "snes_atol": 1e-10,
    "snes_rtol": 1e-10,
    "snes_max_it": 100,
    "snes_divergence_tolerance": "PETSC_UNLIMITED",
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": linear_solver,
    # "snes_monitor": None,
}

problem = NonlinearProblem(
    F,
    u,
    bcs=[bc_outlet, bc_inlet],
    petsc_options=petsc_options,
    petsc_options_prefix="poisson_transient_temp",
)

import matplotlib as mpl
import pyvista
from dolfinx import plot

c_m_post = u.split()[0].collapse()
c_t_post = u.split()[1].collapse()

grid_c_m = pyvista.UnstructuredGrid(*plot.vtk_mesh(c_m_post.function_space))
grid_c_t = pyvista.UnstructuredGrid(*plot.vtk_mesh(c_t_post.function_space))

grid_c_m.point_data["c_m"] = c_m_post.x.array
grid_c_t.point_data["c_t"] = c_t_post.x.array

viridis = mpl.colormaps.get_cmap("viridis").resampled(50)
sargs = dict(
    title_font_size=25,
    label_font_size=20,
    fmt="%.2e",
    color="black",
    position_x=0.1,
    position_y=0.8,
    width=0.8,
    height=0.1,
)

plotter = pyvista.Plotter(shape=(1, 2))
plotter.open_gif("transient_temperature.gif", fps=7)

plotter.subplot(0, 0)
plotter.view_xy(bounds=[0, 1, 0, 1, 0, 0])
_ = plotter.add_mesh(
    grid_c_m,
    show_edges=False,
    lighting=False,
    cmap=viridis,
    scalar_bar_args=sargs,
    clim=[0, 1],
)

plotter.subplot(0, 1)
plotter.view_xy(bounds=[0, 1, 0, 1, 0, 0])

_ = plotter.add_mesh(
    grid_c_t,
    show_edges=False,
    lighting=False,
    cmap=viridis,
    scalar_bar_args=sargs,
    clim=[0, 0.2],
)

2026-05-22 15:22:38.465 (   0.479s) [    7311AB727740]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=

inventories_cm = []
inventories_ct = []
times = []

t = 0.0
t_final = 20
n_it = 0

while t < t_final:
    t += dt.value
    n_it += 1
    times.append(t)

    # solve the problem with the current u_n as previous solution
    problem.solve()
    converged = problem.solver.getConvergedReason()
    num_iter = problem.solver.getIterationNumber()
    assert converged > 0, f"Solver did not converge, got {converged}."
    print(
        f"Time: {t:.2f} ({n_it=}). \n Solver converged after {num_iter} iterations with converged reason {converged}."
    )

    # update u_n with the current solution u
    u_n.x.array[:] = u.x.array[:]

    # update inlet value to show transient response
    c_inlet.value = 1.0 if t < 5 else 0.0
    T.value = 300.0 if t < 15 else 800.0

    # post processing
    c_m_post = u.split()[0].collapse()
    c_t_post = u.split()[1].collapse()

    # Update plot
    grid_c_m.point_data["c_m"][:] = c_m_post.x.array
    grid_c_t.point_data["c_t"][:] = c_t_post.x.array
    plotter.write_frame()

    # compute inventory
    inventories_cm.append(assemble_scalar(c_m_post * ufl.dx))
    inventories_ct.append(assemble_scalar(c_t_post * ufl.dx))

plotter.close()

Time: 0.30 (n_it=1). 
 Solver converged after 3 iterations with converged reason 2.

Time: 0.60 (n_it=2). 
 Solver converged after 2 iterations with converged reason 2.
Time: 0.90 (n_it=3). 
 Solver converged after 2 iterations with converged reason 2.
Time: 1.20 (n_it=4). 
 Solver converged after 2 iterations with converged reason 2.
Time: 1.50 (n_it=5). 
 Solver converged after 2 iterations with converged reason 2.
Time: 1.80 (n_it=6). 
 Solver converged after 2 iterations with converged reason 2.

Time: 2.10 (n_it=7). 
 Solver converged after 2 iterations with converged reason 2.
Time: 2.40 (n_it=8). 
 Solver converged after 2 iterations with converged reason 2.
Time: 2.70 (n_it=9). 
 Solver converged after 2 iterations with converged reason 2.
Time: 3.00 (n_it=10). 
 Solver converged after 2 iterations with converged reason 2.
Time: 3.30 (n_it=11). 
 Solver converged after 2 iterations with converged reason 2.

Time: 3.60 (n_it=12). 
 Solver converged after 2 iterations with converged reason 2.
Time: 3.90 (n_it=13). 
 Solver converged after 2 iterations with converged reason 2.
Time: 4.20 (n_it=14). 
 Solver converged after 2 iterations with converged reason 2.
Time: 4.50 (n_it=15). 
 Solver converged after 2 iterations with converged reason 2.
Time: 4.80 (n_it=16). 
 Solver converged after 2 iterations with converged reason 2.

Time: 5.10 (n_it=17). 
 Solver converged after 2 iterations with converged reason 2.
Time: 5.40 (n_it=18). 
 Solver converged after 2 iterations with converged reason 2.
Time: 5.70 (n_it=19). 
 Solver converged after 2 iterations with converged reason 2.
Time: 6.00 (n_it=20). 
 Solver converged after 2 iterations with converged reason 2.
Time: 6.30 (n_it=21). 
 Solver converged after 2 iterations with converged reason 2.

Time: 6.60 (n_it=22). 
 Solver converged after 2 iterations with converged reason 2.
Time: 6.90 (n_it=23). 
 Solver converged after 2 iterations with converged reason 2.
Time: 7.20 (n_it=24). 
 Solver converged after 2 iterations with converged reason 2.
Time: 7.50 (n_it=25). 
 Solver converged after 2 iterations with converged reason 2.
Time: 7.80 (n_it=26). 
 Solver converged after 2 iterations with converged reason 2.

Time: 8.10 (n_it=27). 
 Solver converged after 2 iterations with converged reason 2.
Time: 8.40 (n_it=28). 
 Solver converged after 1 iterations with converged reason 2.
Time: 8.70 (n_it=29). 
 Solver converged after 1 iterations with converged reason 2.
Time: 9.00 (n_it=30). 
 Solver converged after 1 iterations with converged reason 2.
Time: 9.30 (n_it=31). 
 Solver converged after 1 iterations with converged reason 2.

Time: 9.60 (n_it=32). 
 Solver converged after 1 iterations with converged reason 2.
Time: 9.90 (n_it=33). 
 Solver converged after 1 iterations with converged reason 2.
Time: 10.20 (n_it=34). 
 Solver converged after 1 iterations with converged reason 2.
Time: 10.50 (n_it=35). 
 Solver converged after 1 iterations with converged reason 2.

Time: 10.80 (n_it=36). 
 Solver converged after 1 iterations with converged reason 2.
Time: 11.10 (n_it=37). 
 Solver converged after 1 iterations with converged reason 2.
Time: 11.40 (n_it=38). 
 Solver converged after 1 iterations with converged reason 2.
Time: 11.70 (n_it=39). 
 Solver converged after 1 iterations with converged reason 2.

Time: 12.00 (n_it=40). 
 Solver converged after 1 iterations with converged reason 2.
Time: 12.30 (n_it=41). 
 Solver converged after 1 iterations with converged reason 2.
Time: 12.60 (n_it=42). 
 Solver converged after 1 iterations with converged reason 2.
Time: 12.90 (n_it=43). 
 Solver converged after 1 iterations with converged reason 2.

Time: 13.20 (n_it=44). 
 Solver converged after 1 iterations with converged reason 2.
Time: 13.50 (n_it=45). 
 Solver converged after 1 iterations with converged reason 2.
Time: 13.80 (n_it=46). 
 Solver converged after 1 iterations with converged reason 2.
Time: 14.10 (n_it=47). 
 Solver converged after 1 iterations with converged reason 2.

Time: 14.40 (n_it=48). 
 Solver converged after 1 iterations with converged reason 2.
Time: 14.70 (n_it=49). 
 Solver converged after 1 iterations with converged reason 2.
Time: 15.00 (n_it=50). 
 Solver converged after 1 iterations with converged reason 2.
Time: 15.30 (n_it=51). 
 Solver converged after 2 iterations with converged reason 2.

Time: 15.60 (n_it=52). 
 Solver converged after 2 iterations with converged reason 2.
Time: 15.90 (n_it=53). 
 Solver converged after 2 iterations with converged reason 2.
Time: 16.20 (n_it=54). 
 Solver converged after 2 iterations with converged reason 2.

Time: 16.50 (n_it=55). 
 Solver converged after 2 iterations with converged reason 2.
Time: 16.80 (n_it=56). 
 Solver converged after 2 iterations with converged reason 2.
Time: 17.10 (n_it=57). 
 Solver converged after 2 iterations with converged reason 2.
Time: 17.40 (n_it=58). 
 Solver converged after 2 iterations with converged reason 2.

Time: 17.70 (n_it=59). 
 Solver converged after 2 iterations with converged reason 2.
Time: 18.00 (n_it=60). 
 Solver converged after 2 iterations with converged reason 2.
Time: 18.30 (n_it=61). 
 Solver converged after 2 iterations with converged reason 2.
Time: 18.60 (n_it=62). 
 Solver converged after 2 iterations with converged reason 2.

Time: 18.90 (n_it=63). 
 Solver converged after 1 iterations with converged reason 2.
Time: 19.20 (n_it=64). 
 Solver converged after 1 iterations with converged reason 2.
Time: 19.50 (n_it=65). 
 Solver converged after 1 iterations with converged reason 2.
Time: 19.80 (n_it=66). 
 Solver converged after 1 iterations with converged reason 2.

Time: 20.10 (n_it=67). 
 Solver converged after 1 iterations with converged reason 2.

transient_temp

When the inlet concentration drops, the mobile inventory quickly drops to zero. But the trapped inventory doesn’t. That’s because the detrapping rate \(p\) is too low at this temperature.

It’s only when the temperature is increased at \(t=15\) that the trapped hydrogen starts detrapping and leave.

import matplotlib.pyplot as plt
plt.stackplot(times, inventories_cm, inventories_ct, labels=["mobile", "trapped"])
plt.ylabel("Inventory")
plt.xlabel("Time")
plt.legend(reverse=True)
plt.show()

../../_images/9573e588465b9e2f581eca7bd1ca5b7f0ac7495e62faf32b577d952d1012cfc5.png

Non-homogeneous temperature#

Now we will do the same thing but with a non-homogeneous temperature (ie. varying in space)

(25)#\[\begin{equation} T(x,y) = 300 + 300~x + 200~y \end{equation}\]

nx = ny = 30

domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.quadrilateral)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)

cg_element = basix.ufl.element("Lagrange", domain.basix_cell(), degree=1)

mixed_element = basix.ufl.mixed_element([cg_element, cg_element])

V = fem.functionspace(domain, mixed_element)

u = fem.Function(V)
u_n = fem.Function(V)

cm, ct = ufl.split(u)
cm_n, ct_n = ufl.split(u_n)

v_cm, v_ct = ufl.TestFunctions(V)

def inlet(x):
    return np.logical_and(np.isclose(x[0], 0), x[1] <= 0.5)

def outlet(x):
    return np.logical_and(np.isclose(x[0], 1), x[1] >= 0.5)

V0, submap = V.sub(0).collapse()

dofs_outlet = fem.locate_dofs_geometrical((V.sub(0), V0), outlet)
dofs_inlet = fem.locate_dofs_geometrical((V.sub(0), V0), inlet)

c_inlet = fem.Constant(domain, 1.0)
c_outlet = fem.Constant(domain, 0.0)

bc_outlet = fem.dirichletbc(c_outlet, dofs_outlet[0], V.sub(0))
bc_inlet = fem.dirichletbc(c_inlet, dofs_inlet[0], V.sub(0))

Because \(T\) is now a function of space, it needs to become a dolfinx.fem.Function. We create a new functionspace for T, and then create a Function from it.

Then, we call the .interpolate() method with an appropriate lambda function of x.

Note

In this context, x[0] is the first coordinate (\(x\)) and x[1] is the second one (\(y\))

V_cg = dolfinx.fem.functionspace(domain, ("CG", 1))
T = dolfinx.fem.Function(V_cg)

T.interpolate(lambda x: 300.0 + 300.0*x[0] + 200.0*x[1])

def arrhenius_rate(k0, Ek, T, mod=ufl):
    boltzmann_constant = 8.617333262145e-5  # eV/K
    return k0 * mod.exp(-Ek/(T*boltzmann_constant))

k = arrhenius_rate(1e2, 0.2, T) # trapping rate
p = arrhenius_rate(5e9, 1.5, T)  # detrapping rate
n = 0.5  # total trapping sites
D = arrhenius_rate(1e3, 0.2, T) # diffusion coefficient

dt = dolfinx.fem.Constant(domain, 0.6)

F_mobile_transient = (cm - cm_n)/dt* v_cm * ufl.dx
F_trapped_transient = (ct - ct_n)/dt * v_ct * ufl.dx

trapping = k * cm * (n - ct)
detrapping = p * ct

F_mobile = (
    D*ufl.dot(ufl.grad(cm), ufl.grad(v_cm)) * ufl.dx
    + trapping * v_cm * ufl.dx
    - detrapping * v_cm * ufl.dx
)
F_trapped = -trapping * v_ct * ufl.dx + detrapping * v_ct * ufl.dx

F = F_mobile_transient + F_trapped_transient + F_mobile + F_trapped

# taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28
use_superlu = PETSc.IntType == np.int64  # or PETSc.ScalarType == np.complex64
sys = PETSc.Sys()  # type: ignore
if sys.hasExternalPackage("mumps") and not use_superlu:
    linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
    linear_solver = "superlu_dist"
else:
    linear_solver = "petsc"

petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2,
    "snes_atol": 1e-10,
    "snes_rtol": 1e-10,
    "snes_max_it": 100,
    "snes_divergence_tolerance": "PETSC_UNLIMITED",
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": linear_solver,
    # "snes_monitor": None,
}

problem = NonlinearProblem(
    F,
    u,
    bcs=[bc_outlet, bc_inlet],
    petsc_options=petsc_options,
    petsc_options_prefix="non_homogeneous_temperature",
)

import matplotlib as mpl
c_m_post = u.split()[0].collapse()
c_t_post = u.split()[1].collapse()

grid_c_m = pyvista.UnstructuredGrid(*plot.vtk_mesh(c_m_post.function_space))
grid_c_t = pyvista.UnstructuredGrid(*plot.vtk_mesh(c_t_post.function_space))
grid_T = pyvista.UnstructuredGrid(*plot.vtk_mesh(T.function_space))

grid_c_m.point_data["c_m"] = c_m_post.x.array
grid_c_t.point_data["c_t"] = c_t_post.x.array
grid_T.point_data["T"] = T.x.array

viridis = mpl.colormaps.get_cmap("viridis").resampled(50)
sargs = dict(
    title_font_size=25,
    label_font_size=15,
    fmt="%.2e",
    color="black",
    position_x=0.1,
    position_y=0.8,
    width=0.8,
    height=0.1,
)

plotter = pyvista.Plotter(shape=(1, 3))
plotter.open_gif("non_homogeneous_temperature.gif", fps=7)

plotter.subplot(0, 0)
plotter.view_xy(bounds=[0, 1, 0, 1, 0, 0])
_ = plotter.add_mesh(
    grid_c_m,
    show_edges=False,
    lighting=False,
    cmap=viridis,
    scalar_bar_args=sargs,
    clim=[0, 1],
)

plotter.subplot(0, 1)
plotter.view_xy(bounds=[0, 1, 0, 1, 0, 0])

_ = plotter.add_mesh(
    grid_c_t,
    show_edges=False,
    lighting=False,
    cmap=viridis,
    scalar_bar_args=sargs,
    clim=[0, n],
)

plotter.subplot(0, 2)
plotter.view_xy(bounds=[0, 1, 0, 1, 0, 0])

_ = plotter.add_mesh(
    grid_T,
    show_edges=False,
    lighting=False,
    cmap="inferno",
    scalar_bar_args=sargs,
)

t = 0.0
t_final = 10
n_it = 0

while t < t_final:
    t += dt.value
    n_it += 1

    # solve the problem with the current u_n as previous solution
    problem.solve()
    converged = problem.solver.getConvergedReason()
    num_iter = problem.solver.getIterationNumber()
    assert converged > 0, f"Solver did not converge, got {converged}."
    print(
        f"Time: {t:.2f} ({n_it=}). \n Solver converged after {num_iter} iterations with converged reason {converged}."
    )

    # update u_n with the current solution u
    u_n.x.array[:] = u.x.array[:]

    # update inlet value to show transient response
    c_inlet.value = 1.0 if t < 5 else 0.0

    # post processing
    c_m_post = u.split()[0].collapse()
    c_t_post = u.split()[1].collapse()

    # Update plot
    grid_c_m.point_data["c_m"][:] = c_m_post.x.array
    grid_c_t.point_data["c_t"][:] = c_t_post.x.array
    plotter.write_frame()

plotter.close()

Time: 0.60 (n_it=1). 
 Solver converged after 3 iterations with converged reason 2.
Time: 1.20 (n_it=2). 
 Solver converged after 2 iterations with converged reason 2.
Time: 1.80 (n_it=3). 
 Solver converged after 2 iterations with converged reason 2.

Time: 2.40 (n_it=4). 
 Solver converged after 2 iterations with converged reason 2.
Time: 3.00 (n_it=5). 
 Solver converged after 2 iterations with converged reason 2.
Time: 3.60 (n_it=6). 
 Solver converged after 2 iterations with converged reason 2.
Time: 4.20 (n_it=7). 
 Solver converged after 2 iterations with converged reason 2.

Time: 4.80 (n_it=8). 
 Solver converged after 2 iterations with converged reason 2.
Time: 5.40 (n_it=9). 
 Solver converged after 2 iterations with converged reason 2.
Time: 6.00 (n_it=10). 
 Solver converged after 2 iterations with converged reason 2.
Time: 6.60 (n_it=11). 
 Solver converged after 2 iterations with converged reason 2.

Time: 7.20 (n_it=12). 
 Solver converged after 2 iterations with converged reason 2.
Time: 7.80 (n_it=13). 
 Solver converged after 2 iterations with converged reason 2.
Time: 8.40 (n_it=14). 
 Solver converged after 2 iterations with converged reason 2.
Time: 9.00 (n_it=15). 
 Solver converged after 2 iterations with converged reason 2.

Time: 9.60 (n_it=16). 
 Solver converged after 2 iterations with converged reason 2.
Time: 10.20 (n_it=17). 
 Solver converged after 2 iterations with converged reason 2.

conc_non_hom

Temperature dependence

Contents

Temperature dependence#

Non-constant temperature#

Non-homogeneous temperature#