Issues with 1D Hydrogen Diffusion Simulation Modeling

y5k4d6stjv-crypto · October 31, 2025, 9:55am

Hello everyone, I am trying to reproduce the study from the paper “Sensitivity study of hydrogen Soret transport in yttrium Hydride-Based nuclear fuel.” When calculating the radial hydrogen concentration variation in a 1D fuel rod, I encountered difficulties in setting the boundary conditions.

According to the paper, the boundary condition is not a DirichletBC. Therefore, I attempted to use FluxBC, intending to set the hydrogen flux to zero in this 1D model. However, the final result is a straight line at zero. Does this mean that my boundary condition setup failed? How should I resolve this issue, or what boundary condition should I choose instead?

Best wishes！

Here are the details of my code.

import festim as F
import matplotlib.pyplot as plt
import numpy as np
import os
class ProfileExport(F.VolumeQuantity):
def compute(self):
profile = self.field.solution.x.array[:].copy()
self.data.append(profile)
#model parameters
Rfo=6.35e-3#m
D_0=3.1e-8 #m^2/s
E_D=0.42719#ev
Tfo=684#degrees celsius
Tcl=742#degrees celsius
My=0.088906 #kg/mol
TQ=640 #K
Py=4265 #kg/m^3
Q=0.05514
#set parameter
node_num=1001
time=10
time_step=1
x=1.8
initial_H=5.004e16
cladding_thickness=1e-3

#equilibrium partial pressure of hydrogen gas
#B=225.128np.power(x,5)-1898.266np.power(x,4)+5617.835np.power(x,3)-8264.63np.power(x,2)+6046.559x-1750.929
B=12
d=6.421827e5np.power(x,5)-4.788489e6np.power(x,4)+1.41948e7np.power(x,3)-2.090897e7np.power(x,2)+1.531054e7x-4.435128e6

P_H=np.exp(2.303*B-(d/Tfo))
print(F.version)
my_model = F.Simulation()

my_model.mesh = F.MeshFromVertices(vertices=np.linspace(-Rfo,Rfo, num=node_num))
my_model.materials = F.Material(id=1,D_0=D_0,E_D=E_D,Q=Q)

def temperature_func(x):
return 272+Tfo + (Tcl-Tfo)*(1-np.power((x[0]/Rfo),2))

T_val = lambda x: Tfo + (Tcl-Tfo)*(1-np.power((x/Rfo),2))+272
my_model.T = T_val(F.x)
bc=F.DirichletBC(surfaces=[1, 2], value=1e15, field=0) # H/m3

my_custom_value=5.4e-8*np.exp(-6760/Tfo)*np.sqrt(P_H)/cladding_thickness
bc=F.FluxBC(surfaces=[1,2], value=0,field=0)

#bc=F.DirichletBC(surfaces=[1, 2], value=1e15, field=0) # H/m3
my_model.boundary_conditions = [
bc
]
my_model.initial_conditions=[F.InitialCondition(field=0,value=initial_H)]
#my_model.sources = [F.Source(value=1e16, volume=1, field=0)] # H/m3/s

my_model.settings = F.Settings(
absolute_tolerance=1e10,
relative_tolerance=1e-10,
soret=True
)
derived_quantities = F.DerivedQuantities(
[F.TotalVolume(“solute”, volume=1)],
show_units=True,
)

txt_export = F.TXTExport(field=“solute”, filename=“task12/profile_soret.txt”)

my_model.exports = [derived_quantities, txt_export]

#my_model.dt = F.Stepsize(time_step, milestones=[0.1, 0.2, 0.5, 1]) # s
my_model.dt=None
my_model.settings.transient=False
my_model.initialise()
my_model.run()

my_model.settings.soret = False
txt_export.filename = “task12/profile_no_soret.txt”
my_model.materials[0].Q = None
my_model.initialise()
my_model.run()

import matplotlib.pyplot as plt

data = np.genfromtxt(“task12/profile_soret.txt”, delimiter=“,”, names=True)
data_no_soret = np.genfromtxt(“task12/profile_no_soret.txt”, delimiter=“,”, names=True)

name = data_no_soret.dtype.names[1]

fig, axs = plt.subplots(2, 1, sharex=True)

axs[0].plot(data_no_soret[“x”], data_no_soret[name], label=“No Soret”)

name = data.dtype.names[1]
axs[0].plot(data[“x”], data[name], label=“Soret”)

axs[0].legend()
axs[0].set_ylabel(“Retention [T m$^{-3}$]”)

axs[1].plot(data[“x”], T_val(data[“x”]))
axs[1].set_xlabel(“x [m]”)
axs[1].set_ylabel(“Temperature [K]”)
plt.show()

remidm · October 31, 2025, 3:31pm

Hi @y5k4d6stjv-crypto and welcome to the FESTIM discourse! Do you mind formatting your code with triple quotes?

remidm · October 31, 2025, 3:31pm

What should the boundary condition be in the paper?

y5k4d6stjv-crypto · November 1, 2025, 9:06am

Thank you for your reply!
Here are the details of my code.
'''
import festim as F
import matplotlib.pyplot as plt
import numpy as np
import os
class ProfileExport(F.VolumeQuantity):
    def compute(self):
        profile = self.field.solution.x.array[:].copy()
        self.data.append(profile)
#model parameters
Rfo=6.35e-3#m
D_0=3.1e-8 #m^2/s  
E_D=0.42719#ev
Tfo=684#degrees celsius
Tcl=742#degrees celsius
My=0.088906 #kg/mol
TQ=640 #K
Py=4265 #kg/m^3
Q=0.05514
#set parameter
node_num=1001
time=10
time_step=1
x=1.8
initial_H=5.004e16
cladding_thickness=1e-3


#equilibrium partial pressure of hydrogen gas
#B=225.128*np.power(x,5)-1898.266*np.power(x,4)+5617.835*np.power(x,3)-8264.63*np.power(x,2)+6046.559*x-1750.929
B=12
d=6.421827e5*np.power(x,5)-4.788489e6*np.power(x,4)+1.41948e7*np.power(x,3)-2.090897e7*np.power(x,2)+1.531054e7*x-4.435128e6

P_H=np.exp(2.303*B-(d/Tfo))
print(F.__version__)
my_model = F.Simulation()


my_model.mesh = F.MeshFromVertices(vertices=np.linspace(-Rfo,Rfo, num=node_num))
my_model.materials = F.Material(id=1,D_0=D_0,E_D=E_D,Q=Q)

def temperature_func(x):
    return 272+Tfo + (Tcl-Tfo)*(1-np.power((x[0]/Rfo),2))

T_val = lambda x: Tfo + (Tcl-Tfo)*(1-np.power((x/Rfo),2))+272
my_model.T = T_val(F.x)
bc=F.DirichletBC(surfaces=[1, 2], value=1e15, field=0)  # H/m3

my_custom_value=5.4e-8*np.exp(-6760/Tfo)*np.sqrt(P_H)/cladding_thickness
bc=F.FluxBC(surfaces=[1,2], value=0,field=0)

#bc=F.DirichletBC(surfaces=[1, 2], value=1e15, field=0)  # H/m3
my_model.boundary_conditions = [
    bc
]
my_model.initial_conditions=[F.InitialCondition(field=0,value=initial_H)]
#my_model.sources = [F.Source(value=1e16, volume=1, field=0)]  # H/m3/s



my_model.settings = F.Settings(
    absolute_tolerance=1e10, 
    relative_tolerance=1e-10, 
    soret=True
)
derived_quantities = F.DerivedQuantities(
    [F.TotalVolume("solute", volume=1)],
    show_units=True,
)

txt_export = F.TXTExport(field="solute", filename="task12/profile_soret.txt")

my_model.exports = [derived_quantities, txt_export]

#my_model.dt = F.Stepsize(time_step, milestones=[0.1, 0.2, 0.5, 1])  # s
my_model.dt=None
my_model.settings.transient=False
my_model.initialise()
my_model.run()

my_model.settings.soret = False
txt_export.filename = "task12/profile_no_soret.txt"
my_model.materials[0].Q = None
my_model.initialise()
my_model.run()

import matplotlib.pyplot as plt

data = np.genfromtxt("task12/profile_soret.txt", delimiter=",", names=True)
data_no_soret = np.genfromtxt("task12/profile_no_soret.txt", delimiter=",", names=True)

name = data_no_soret.dtype.names[1]

fig, axs = plt.subplots(2, 1, sharex=True)

axs[0].plot(data_no_soret["x"], data_no_soret[name], label="No Soret")

name = data.dtype.names[1]
axs[0].plot(data["x"], data[name], label="Soret")

axs[0].legend()
axs[0].set_ylabel("Retention [T m$^{-3}$]")

axs[1].plot(data["x"], T_val(data["x"]))
axs[1].set_xlabel("x [m]")
axs[1].set_ylabel("Temperature [K]")
plt.show()
'''

The results from running my code are as follows：

y5k4d6stjv-crypto · November 1, 2025, 9:15am

The boundary conditions are discussed as follows in the paper: Hydrogen atoms dissociate from the metal hydride, combine into gaseous hydrogen molecules at the fuel rod interface, escape into the gas gap, and then permeate from the gas gap through the cladding layer. The condition for the entire model to reach steady state is achieved by setting JH=0 in the stainless steel,

JH=5.4e-8*np.exp(-6760/Tfo)np.sqrt(P_H)/cladding_thickness.Regarding this boundary condition, my code is as follows:

'''
my_custom_value=5.4e-8np.exp(-6760/Tfo)*np.sqrt(P_H)/cladding_thicknessbc=F.FluxBC(surfaces=[1,2], value=0,field=0)
'''

I used FluxBC for the settings, but the final result turned out to be a straight line at zero I don’t know where the problem is.

Thank you again for replying to my message. I would be extremely grateful if I could get your answer.

VVKulagin · November 14, 2025, 5:08am

Hi @y5k4d6stjv-crypto

Have you solved the problem already? If you still have issues with the model, can you provide more details on the problem? A sketch of the geometry might help.

Looking at your code I see some misleading lines (such as initial conditions for the steady-state problem, several overlapping BCs) and, probably, a units issue (e.g. your comments state that Tfo and Tcl are in Celcius, but you use them in expressions for preassure and permiability). Would you mind clearing the code a bit to make it easier to help you?

You might also need to change simulations tolerances, since you work with moles (at least moles appear in eq. (18) of the referenced paper)

y5k4d6stjv-crypto · November 26, 2025, 2:02am

Thanks for your advice,my problem is solved now.

remidm · November 26, 2025, 3:21pm

@y5k4d6stjv-crypto would you mind sharing how you fixed your issues for future users?

Topic		Replies	Views
Understanding depth profiles User Support	6	93	February 20, 2025
Multi-material hydrogen permeation 3D model transient solution not converging User Support	5	89	August 9, 2024
Error with HenrysBC and zero activation energy User Support festim2	8	76	April 8, 2025
Coupled simulation - application of tolerances User Support festim2	4	70	September 19, 2025
FESTIM Capabilities for Modeling Hydrogen Dissolving + Mixing Into Another Material User Support	5	78	February 18, 2025

Issues with 1D Hydrogen Diffusion Simulation Modeling

Related topics