Example tritium release models with transient mass transport coefficients

example_freezing.py

from libra_toolbox.tritium.model import Model, ureg

import numpy as np

baby_diameter = 14 * ureg.cm  # TODO confirm with CAD
baby_radius = 0.5 * baby_diameter
baby_volume = 1 * ureg.L
baby_cross_section = np.pi * baby_radius**2
baby_height = baby_volume / baby_cross_section

optimised_ratio = 1.7e-2
k_top = 8.9e-8 * ureg.m * ureg.s**-1
k_wall = optimised_ratio * k_top


def k_wall_fun(t):
    try:
        if 10 * ureg.day < t < 10 * ureg.day + 6 * ureg.h:
            return 0 * k_wall
        else:
            return k_wall
    except ValueError:
        idx = np.where((10 * ureg.day < t) & (t < 10 * ureg.day + 6 * ureg.h))
        k = np.ones_like(t) * k_wall
        k[idx] = 0 * k_wall
        return k


def k_top_fun(t):
    try:
        if 10 * ureg.day < t < 10 * ureg.day + 6 * ureg.h:
            return 0 * k_top
        else:
            return k_top
    except ValueError:
        idx = np.where((10 * ureg.day < t) & (t < 10 * ureg.day + 6 * ureg.h))
        k = np.ones_like(t) * k_top
        k[idx] = 0 * k_top
        return k


# Neutron rate
# calculated from Kevin's activation foil analysis from run 100 mL #7
# TODO replace for values for this run
P383_neutron_rate = 4.95e8 * ureg.neutron * ureg.s**-1
A325_neutron_rate = 2.13e8 * ureg.neutron * ureg.s**-1

neutron_rate_relative_uncertainty = 0.089
neutron_rate = (
    P383_neutron_rate
) / 2  # the neutron rate is divided by two to acount for the double counting (two detectors)

baby_model = Model(
    radius=baby_radius,
    height=baby_height,
    TBR=2e-3,
    neutron_rate=neutron_rate,
    irradiations=[
        (0 * ureg.h, 12 * ureg.h),
        (24 * ureg.h, 36 * ureg.h),
    ],
    k_top=k_top_fun,
    k_wall=k_wall_fun,
)

baby_model.run(40 * ureg.day)

import matplotlib.pyplot as plt
from libra_toolbox.tritium.plotting import (
    plot_irradiation,
    plot_integrated_wall_release,
    plot_integrated_top_release,
    plot_salt_inventory,
    plot_top_release,
    plot_wall_release,
)

fig, axs = plt.subplots(nrows=3, sharex=True, figsize=(10, 10))

plt.sca(axs[0])
plt.title("Salt inventory")
plot_salt_inventory(baby_model)
plot_irradiation(baby_model, color="tab:red", alpha=0.25)

plt.sca(axs[1])
plt.title("Release rates")
plot_wall_release(baby_model, label="wall")
plot_top_release(baby_model, label="top")
plot_irradiation(baby_model, color="tab:red", alpha=0.25)

plt.annotate(
    "Salt is frozen",
    (10 * ureg.day, 0.12),
    (15 * ureg.day, 0.2),
    arrowprops=dict(facecolor="black", width=0.5, headwidth=5),
)

plt.legend()

plt.sca(axs[2])
plt.title("Cumulative release")
plot_integrated_wall_release(baby_model, label="wall")
plot_integrated_top_release(baby_model, label="top")
plot_irradiation(baby_model, color="tab:red", alpha=0.25)
plt.legend()


plt.show()

example_h2.py

from libra_toolbox.tritium.model import Model, ureg

import numpy as np

baby_diameter = 14 * ureg.cm  # TODO confirm with CAD
baby_radius = 0.5 * baby_diameter
baby_volume = 1 * ureg.L
baby_cross_section = np.pi * baby_radius**2
baby_height = baby_volume / baby_cross_section

optimised_ratio = 1.7e-2
k_top = 8.9e-8 * ureg.m * ureg.s**-1
k_wall = optimised_ratio * k_top


time_when_h2_is_added = 10 * ureg.day

increase_factor_with_h2 = 5

def k_wall_fun(t):
    try:
        if t < time_when_h2_is_added:
            return 1 * k_wall
        else:
            return increase_factor_with_h2* k_wall
    except ValueError:
        idx = np.where(t > time_when_h2_is_added)
        k = np.ones_like(t) * k_wall
        k[idx] = increase_factor_with_h2 * k_wall
        return k


def k_top_fun(t):
    try:
        if t < time_when_h2_is_added:
            return 1 * k_top
        else:
            return increase_factor_with_h2 * k_top
    except ValueError:
        idx = np.where(t > time_when_h2_is_added)
        k = np.ones_like(t) * k_top
        k[idx] = increase_factor_with_h2 * k_top
        return k


# Neutron rate
# calculated from Kevin's activation foil analysis from run 100 mL #7
# TODO replace for values for this run
P383_neutron_rate = 4.95e8 * ureg.neutron * ureg.s**-1
A325_neutron_rate = 2.13e8 * ureg.neutron * ureg.s**-1

neutron_rate_relative_uncertainty = 0.089
neutron_rate = (
    P383_neutron_rate
) / 2  # the neutron rate is divided by two to acount for the double counting (two detectors)

baby_model = Model(
    radius=baby_radius,
    height=baby_height,
    TBR=2e-3,
    neutron_rate=neutron_rate,
    irradiations=[
        (0 * ureg.h, 12 * ureg.h),
        (24 * ureg.h, 36 * ureg.h),
    ],
    k_top=k_top_fun,
    k_wall=k_wall_fun,
)

baby_model.run(40 * ureg.day)

import matplotlib.pyplot as plt
from libra_toolbox.tritium.plotting import (
    plot_irradiation,
    plot_integrated_wall_release,
    plot_integrated_top_release,
    plot_salt_inventory,
    plot_top_release,
    plot_wall_release,
)

fig, axs = plt.subplots(nrows=3, sharex=True, figsize=(10, 10))

plt.sca(axs[0])
plt.title("Salt inventory")
plot_salt_inventory(baby_model)
plot_irradiation(baby_model, color="tab:red", alpha=0.25)

plt.sca(axs[1])
plt.title("Release rates")
plot_wall_release(baby_model, label="wall")
plot_top_release(baby_model, label="top")
plot_irradiation(baby_model, color="tab:red", alpha=0.25)

plt.annotate(
    f"Sweep gas turned to H2 \n mass transports x{increase_factor_with_h2}",
    (time_when_h2_is_added, 0.12),
    (15 * ureg.day, 0.2),
    arrowprops=dict(facecolor="black", width=0.5, headwidth=5),
)

plt.legend()

plt.sca(axs[2])
plt.title("Cumulative release")
plot_integrated_wall_release(baby_model, label="wall")
plot_integrated_top_release(baby_model, label="top")
plot_irradiation(baby_model, color="tab:red", alpha=0.25)
plt.legend()


plt.show()