Advanced functionality

Advanced functionality#

FESTIM has an internal solver HeatTransferProblem for heat transfer problems. This can be helpful for multiphysics coupling or simulations where temperatures cannot be prescribed, but must be solved for. This tutorial shows how to use FESTIM’s heat-transfer solvers and coupling it to hydrogen transport simulations.

Objectives:

Solve steady-state and transient heat-transfer problems in FESTIM
Couple transient hydrogen transport and heat-transfer simulations

Defining and running a steady-state heat transfer simulation#

Consider the following 2D heat transfer problem, where we want to find the temperature field by solving the steady-state heat transfer equation:

\[ \nabla \cdot (\lambda \ \nabla T) + \dot{Q} = 0 \]

where \(\lambda\) is the thermal conductivity \((\text{W/mK})\) and \(\dot{{Q}}\) is the heat transfer rate \((\text{W})\).

Additionally, we define the following boundary conditions in our problem:

Left boundary (\(x=0\)): Fixed temperature

\[T(x,y) = 350 + 20 \cos(x)\, \sin(y)\]

Top and bottom boundaries (\(y=0\) and \(y=1\)): Convective heat flux

\[q_n = h(x) \left(T_{\text{ext}}(x) - T\right), \quad h(x) = 100 x, \quad T_{\text{ext}}(x) = 300 + 3 y\]

Right boundary (\(x=1\)): Prescribed heat flux

\[q_n = 10 + 3 \cos(x) + \sin(y)\]

Volume source term:

\[Q(x) = 1 + 0.1 x\]

import festim as F

heat_transfer_model = F.HeatTransferProblem()

/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

We first define a thermal conductivity function \( \lambda \) and assign it to our material:

def thermal_cond_function(T):
    return 3 + 0.1 * T

mat = F.Material(D_0=1, E_D=0.1, thermal_conductivity=thermal_cond_function)

We can now add heat sources (HeatSource), fixed temperature (FixedTemperatureBC), and heat flux (HeatFluxBC) boundary conditions. We also prescribe a convective heat flux by defining an external temperature T_ext and heat transfer coefficient h:

import dolfinx
from mpi4py import MPI
import numpy as np

nx = ny = 20
mesh_fenics = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, nx, ny)
heat_transfer_model.mesh = F.Mesh(mesh=mesh_fenics)

volume_subdomain = F.VolumeSubdomain(id=1, material=mat)

top_bot = F.SurfaceSubdomain(id=2, locator=lambda x: np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0)))
left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0.0))
right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1.0))

heat_transfer_model.subdomains = [volume_subdomain, top_bot, left, right]

heat_transfer_model.sources = [
    F.HeatSource(value=lambda x: 1 + 0.1 * x[0], volume=volume_subdomain)
]

import ufl

fixed_temperature_left = F.FixedTemperatureBC(
    subdomain=left, value=lambda x: 350 + 20 * ufl.cos(x[0]) * ufl.sin(x[1])
)

def h_coeff(x):
    return 100 * x[0]

def T_ext(x):
    return 300 + 3 * x[1]

convective_heat_transfer = F.HeatFluxBC(
    subdomain=top_bot, value=lambda x, T: h_coeff(x) * (T_ext(x) - T)
)

heat_flux = F.HeatFluxBC(
    subdomain=right, value=lambda x: 10 + 3 * ufl.cos(x[0]) + ufl.sin(x[1])
)

heat_transfer_model.boundary_conditions = [
    fixed_temperature_left,
    convective_heat_transfer,
    heat_flux,
]

heat_transfer_model.settings = F.Settings(
    transient=False,
    atol=1e-09,
    rtol=1e-09,
)

heat_transfer_model.initialise()
heat_transfer_model.run()

2026-05-22 15:24:00.478 (   0.660s) [    78ECB2B73740]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=

Running a hydrogen transport simulation#

To use this temperature field for a hydrogen transport simulation, we need to define another problem hydrogen_problem using HydrogenTransportProblem. We simply assign the output from the heat transfer simulation to the temperature attribute of our hydrogen simulation:

hydrogen_problem = F.HydrogenTransportProblem()

hydrogen_problem.mesh = heat_transfer_model.mesh
H = F.Species("H")
hydrogen_problem.species = [H]
hydrogen_problem.temperature = heat_transfer_model.u

hydrogen_problem.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=left, value=1, species=H),
    F.FixedConcentrationBC(subdomain=right, value=0, species=H),
]

hydrogen_problem.subdomains = heat_transfer_model.subdomains

hydrogen_problem.settings = F.Settings(
    transient=False,
    atol=1e-09,
    rtol=1e-09,
)

hydrogen_problem.initialise()
hydrogen_problem.run()

The governing equation for transient heat transfer is:

\[\rho \ C_p \frac{\partial T}{\partial t} = \nabla \cdot (\lambda \ \nabla T) + \dot{Q} \]

Therefore, to run a transient heat transfer simulation (with no hydrogen transport coupling), we must add the density and heat_capacity:

mat = F.Material(D_0=1 ,E_D=0.01, thermal_conductivity=thermal_cond_function, density=20, heat_capacity=50)
volume_subdomain = F.VolumeSubdomain(id=1, material=mat)

heat_transfer_model.subdomains = [volume_subdomain, top_bot, left, right]

heat_transfer_model.settings = F.Settings(
    atol=1e-09,
    rtol=1e-09,
    stepsize=0.1,
    final_time=5
)

heat_transfer_model.initialise()
heat_transfer_model.run() 

Temperature from a coupled transient heat transfer simulation#

For a coupled heat-transfer and hydrogen transient simulation, we need to use another problem class CoupledTransientHeatTransferHydrogenTransport to ensure the solution solves each problem at each step.

First, we must define a HydrogenTransportProblem and HeatTransferProblem with each model’s settings. We also define a common mesh:

import festim as F
import dolfinx
from mpi4py import MPI
import numpy as np
import ufl

nx = ny = 20
mesh_fenics = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, nx, ny)
mesh = F.Mesh(mesh=mesh_fenics)

mat = F.Material(D_0=1, E_D=0.01, thermal_conductivity=3, density=2, heat_capacity=5)

volume_subdomain = F.VolumeSubdomain(id=1, material=mat)
top_bot = F.SurfaceSubdomain(id=2, locator=lambda x: np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0)))
left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0.0))
right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1.0))
subdomains = [volume_subdomain, top_bot, left, right]

heat_transfer_model = F.HeatTransferProblem()
hydrogen_problem = F.HydrogenTransportProblem()

heat_transfer_model.mesh = mesh   
hydrogen_problem.mesh = mesh

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

hydrogen_problem.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=top_bot, value=1, species=H),
    F.FixedConcentrationBC(subdomain=left, value=0, species=H),
]

heat_transfer_model.subdomains = subdomains
hydrogen_problem.subdomains = subdomains

hydrogen_problem.settings = F.Settings(
    transient=True,
    atol=1e-09,
    rtol=1e-09,
    stepsize=1,
    final_time=50
)

fixed_temperature_left = F.FixedTemperatureBC(
    subdomain=left, value=lambda x: 350 + 20 * ufl.cos(x[0]) * ufl.sin(x[1])
)

heat_transfer_model.boundary_conditions = [
    fixed_temperature_left,
]

heat_transfer_model.settings = F.Settings(
    transient=True,
    atol=1e-09,
    rtol=1e-09,
    stepsize=1,
    final_time=50
)

Finally, we define and solve a new problem using the CoupledTransientHeatTransferHydrogenTransport class:

problem = F.CoupledTransientHeatTransferHydrogenTransport(
    heat_problem=heat_transfer_model,
    hydrogen_problem=hydrogen_problem
    )
    
problem.initialise()
problem.run()

Advanced functionality

Contents

Advanced functionality#

Defining and running a steady-state heat transfer simulation#

Running a hydrogen transport simulation#

Temperature from a coupled transient heat transfer simulation#