Diffusion from a Depleting Source

🏷 Tags: 1D MES transient

This verification case [3] consists of an enclosure containing an initial quantity of hydrogen gas (at an initial pressure \(P_0\)). The hydrogen can dissociate on the inner surface of the enclosure and then permeate through the walls. As hydrogen particles escape the enclosure, the pressure decreases and as a result, so does the surface concentration on the inner walls. This 1D problem is, therefore, coupled. At each time step, the flux of particles escaping the enclosure is computed, and the internal pressure is updated.

An analytical solution to this problem can be obtained. The analytical expression of the concentration in the wall is given by Ambrosek et al [3]:

(24)\[\begin{align} c(x=0, t) &= 0\\ c(x=l, t) &= K_H \ P(t) \\ c(x,t) &= 2 L P_{0} K_H \sum_{n=1}^{\infty} \frac{e^{- D t \alpha_n^2} \sin{\left(x \alpha_n \right)}}{\left(L + l \left(L^{2} + \alpha_n^2\right)\right) \sin{\left(l \alpha_n \right)}} \end{align}\]

where \(L = \dfrac{K_H \ T \ A \ k_B}{V}\), \(K_H\) is the solubility of the material, \(T\) is the temperature, \(A\) is the enclosure surface area, \(V\) is the volume of the enclosure, \(l\) is the thickness of the wall, and \(\alpha_n\) are the roots of:

(25)\[\begin{equation} \alpha_n = \frac{L}{\tan(\alpha_n \ l)} \end{equation}\]

The flux at the non-exposed surface can be expressed as:

(26)\[\begin{align} J(x=0, t) &= D \nabla c |_{x=0} \nonumber \\ &= 2 L P_0 D K_H \sum_{n=1}^{\infty} \frac{e^{- D t \alpha_n^2} \alpha_n}{\left(L + l \left(L^{2} + \alpha_n^2\right)\right) \sin{\left(l \alpha_n \right)}} \end{align}\]

The cumulative release \(R\) is then:

(27)\[\begin{align} \mathrm{R} &= \int_0^t J(x=0, t) dt\nonumber \\ &= 2 L P_0 S \sum_{n=1}^{\infty} \frac{1 - e^{- D t \alpha_n^{2}}}{\left(L + l \left(L^{2} + \alpha_n^2\right)\right) \alpha_n\sin{\left(l \alpha_n \right)}} \end{align}\]

Multiplying by the area and normalising by the initial quantity in the enclosure, one can obtain the fractional release:

(28)\[\begin{align} \mathrm{FR} &= \frac{\mathrm{R} \ A}{P_0 \ V / (k_B T)} \nonumber\\ &= \frac{2 L S A}{V / (k_B T)} \sum_{n=1}^{\infty} \frac{1 - e^{- D t \alpha_n^{2}}}{\left(L + l \left(L^{2} + \alpha_n^2\right)\right) \alpha_n\sin{\left(l \alpha_n \right)}} \end{align}\]

FESTIM Code

import festim as F
import numpy as np

encl_vol = 5.20e-11  # m^3
encl_surf = 2.16e-6  # m^2
l = 3.3e-5  # m
R = 8.314  #
avogadro = 6.022e23  # mol-1
temperature = 2373  # K
initial_pressure = 1e6  # Pa
solubility = 7.244e22 / temperature  # H/m3/Pa
diffusivity = 2.6237e-11  # m2/s


class PressureExport(F.SurfaceQuantity):
    def __init__(self):
        super().__init__(field=F.Species(), surface=left_surface)

    def compute(self):
        self.value = left_bc.value / solubility
        self.data.append(self.value)


class CustomHydrogenTransportProblem(F.HydrogenTransportProblem):
    def iterate(self):
        super().iterate()

        # Update pressure based on flux
        left_flux_val = -left_flux.value
        old_pressure = left_bc.value / solubility
        new_pressure = (
            old_pressure
            - (left_flux_val * encl_surf / encl_vol * R * self.temperature / avogadro)
            * self.dt.value
        )

        left_bc.value = new_pressure * solubility

        self.bc_forms[0] = self.create_dirichletbc_form(left_bc)


model = CustomHydrogenTransportProblem()

vertices = np.linspace(0, l, 150)

model.mesh = F.Mesh1D(vertices)

material = F.Material(D_0=diffusivity, E_D=0.0)

left_surface = F.SurfaceSubdomain1D(id=1, x=0)
right_surface = F.SurfaceSubdomain1D(id=2, x=l)
volume = F.VolumeSubdomain1D(id=3, borders=[0, l], material=material)
model.subdomains = [left_surface, volume, right_surface]

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

left_bc = F.DirichletBC(
    subdomain=left_surface, value=initial_pressure * solubility, species=H
)
model.boundary_conditions = [
    left_bc,
    F.DirichletBC(subdomain=right_surface, value=0, species=H),
]

model.temperature = temperature

left_flux = F.SurfaceFlux(field=H, surface=left_surface)
right_flux = F.SurfaceFlux(field=H, surface=right_surface)
pressure_export = PressureExport()

tot_vol = F.TotalVolume(field=H, volume=volume)
model.exports = [left_flux, right_flux, pressure_export]

model.settings = F.Settings(atol=1e8, rtol=1e-10, final_time=140)

model.settings.stepsize = F.Stepsize(0.05)

model.initialise()
model.run()

/home/docs/checkouts/readthedocs.org/user_builds/festim-vv-report/conda/festim-2/lib/python3.11/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

Comparison with exact solution

Show code cell source

Hide code cell source

import matplotlib.pyplot as plt

# https://stackoverflow.com/questions/28766692/how-to-find-the-intersection-of-two-graphs/28766902#28766902

k = 1.38065e-23  # J/mol Boltzmann constant


def get_roots(L, l, alpha_max, step=0.0001):
    """Gets the roots of alpha = L / tan(alpha * l)

    Args:
        L (float): parameter L
        l (float): parameter l
        alpha_max (float): the maximum alpha to consider
        step (float, optional): the step discretising alphas.
            The smaller the step, the more accurate the roots.
            Defaults to 0.0001.

    Returns:
        np.array: array of roots
    """
    alphas = np.arange(0, alpha_max, step=step)[1:]

    f = alphas

    g = L / np.tan(alphas * l)

    idx = np.argwhere(np.diff(np.sign(f - g))).flatten()

    # remove one every other idx
    idx = idx[::2]

    roots = alphas[idx]
    return roots


def analytical_expression_fractional_release_TMAP7(t, P_0, D, S, V, T, A, l):
    """
    FR = 1 - P(t) / P_0
    where P(t) is the pressure at time t and P_0 is the initial pressure

    Reference: https://doi.org/10.13182/FST05-A967 (Equation 4)
    Note: in the report, the expression of FR is given as P(T)/P_0, but it shown as 1 - P(t)/P_0 in the graph (Figure 1)

    Args:
        t (float, ndarray): time (s)
        P_0 (float): initial presure (Pa)
        D (float): diffusivity (m2/s)
        S (float): solubility (H/m3/Pa)
        V (float): enclosure volume (m3)
        T (float): temperature (K)
        A (float): enclosure surface area (m2)
        l (float): slab length (m)
    """
    L = S * T * A * k / V

    roots = get_roots(L=L, l=l, alpha_max=1e6, step=1)
    roots = roots[:, np.newaxis]
    summation = np.exp(-(roots**2) * D * t) / (l * (roots**2 + L**2) + L)
    summation = np.sum(summation, axis=0)

    pressure = 2 * P_0 * L * summation
    fractional_release = 1 - pressure / P_0
    return fractional_release


def RMSPE(x_FESTIM, x_analytical):
    error = np.sqrt(np.mean((x_FESTIM - x_analytical) ** 2)) / x_analytical.mean()
    return error


# ------------ post processing ----------------

t = np.array(right_flux.t)
right_flux_data = np.array(right_flux.data)
pressures = np.array(pressure_export.data)
fractional_release = 1 - pressures / initial_pressure

plt.figure()
plt.plot(t, fractional_release, label="FESTIM", lw=2, color="tab:blue")

analytical = analytical_expression_fractional_release_TMAP7(
    t=t,
    P_0=initial_pressure,
    D=diffusivity,
    S=solubility,
    V=encl_vol,
    T=temperature,
    A=encl_surf,
    l=l,
)
print(f"RMSPE = {RMSPE(fractional_release, analytical)*100:.2f}%")

plt.plot(
    t,
    analytical,
    label="analytical",
    color="tab:orange",
    linestyle="--",
    dashes=(3, 3),
    lw=2,
)
plt.legend(frameon=False)
plt.xlabel("Time, s")
plt.ylabel("Fractional release")
plt.gca().spines[["right", "top"]].set_visible(False)

plt.show()

RMSPE = 0.22%

Show code cell source

Hide code cell source

def analytical_expression_flux(t, P_0, D, S, V, T, A, l):
    """
    value of the flux at the external surface (not in contact with pressure)
    J = -D * dc/dx

    Args:
        t (float, ndarray): time (s)
        P_0 (float): initial presure (Pa)
        D (float): diffusivity (m2/s)
        S (float): solubility (H/m3/Pa)
        V (float): enclosure volume (m3)
        T (float): temperature (K)
        A (float): enclosure surface area (m2)
        l (float): slab length (m)
    """
    L = S * T * A * k / V

    roots = get_roots(L=L, l=l, alpha_max=1e7, step=1)
    roots = roots[:, np.newaxis]

    summation = (np.exp(-(roots**2) * D * t) * roots) / (
        (l * (roots**2 + L**2) + L) * np.sin(roots * l)
    )
    summation = np.sum(summation, axis=0)

    flux = 2 * S * P_0 * L * D * summation
    flux[0] = 0
    return flux


analytical_flux = analytical_expression_flux(
    t=t,
    P_0=initial_pressure,
    D=diffusivity,
    S=solubility,
    V=encl_vol,
    T=temperature,
    A=encl_surf,
    l=l,
)

print(f"RMSPE = {RMSPE(right_flux_data, analytical_flux)*100:.2f}%")

plt.plot(t, right_flux_data, label="FESTIM", lw=2, color="tab:blue")
plt.plot(
    t,
    analytical_flux,
    color="tab:orange",
    linestyle="--",
    dashes=(3, 3),
    label="analytical",
    lw=2,
)

plt.legend(frameon=False)
plt.ylim(bottom=0)
plt.ylabel(r"Flux at outer surface, H m$^{-2}$s$^{-1}$")
plt.xlabel("Time, s")
plt.gca().spines[["right", "top"]].set_visible(False)
plt.ylim(bottom=0)

plt.show()

RMSPE = 0.29%