Enriched model

sparging101.py

@@ -0,0 +1,108 @@
+from __future__ import annotations
+from sparging.config import ureg
+from sparging import all_correlations
+from sparging import animation
+from sparging.model import Simulation
+from sparging.inputs import (
+    ColumnGeometry,
+    BreederMaterial,
+    OperatingParameters,
+    SpargingParameters,
+    SimulationInput,
+)
+import logging
+from typing import TYPE_CHECKING
+import numpy as np
+
+if TYPE_CHECKING:
+    import pint
+
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.WARNING)
+
+
+geom = ColumnGeometry(
+    area=0.2 * ureg.m**2,
+    height=1 * ureg.m,
+    nozzle_diameter=0.001 * ureg.m,
+    nb_nozzle=10 * ureg.dimensionless,
+)
+
+flibe = BreederMaterial(
+    name="FLiBe",
+)
+
+operating_params = OperatingParameters(
+    temperature=600 * ureg.celsius,
+    P_top=1 * ureg.atm,
+    flow_g_mol=400 * ureg.sccm,
+    tbr=0.1 * ureg("triton / neutron"),
+    n_gen_rate=1e9 * ureg("neutron / s"),
+)
+
+sparging_params = SpargingParameters(
+    h_l=all_correlations("h_l_briggs"),
+)
+
+# class method from_parameters that takes in objects like ColumnGeometry, BreederMaterial, OperatingParameters and returns a SimulationInput object with the appropriate correlations for the given parameters. This method should be able to handle cases where some of the parameters are already provided as correlations and should not overwrite them.
+my_input = SimulationInput.from_parameters(
+    geom, flibe, operating_params, sparging_params
+)
+logger.info(my_input)
+
+print(my_input.get_S_T())
+print(f"{my_input.Q_T.to('molT/s')} = {my_input.Q_T.to('molT2/hour')}")
+print(my_input.volume.to("m^3"))
+print(
+    f"Concentration at steady state (no dispersion, no PP limited): {
+        my_input.get_c_T2_SS().to('molT2/m^3')
+    }"
+)
+T_99 = (np.log(100) * my_input.get_tau()).to("seconds")
+print(f"T_99% = {T_99.to('hours')}")
+print(f"Partial pressure number PP = {my_input.get_PP_number()}")
+
+
+def profile_source_T(z: pint.Quantity | list[float], height: pint.Quantity = None):
+    import numpy as np
+
+    if isinstance(z, (float, np.ndarray, list)):  # non-dimensional height (0 to 1)
+        # return np.pi / 2 * np.sin(np.pi * z)  # normalized
+        return 1 + 1 * np.sin(np.pi * z)  # not normalized
+    if isinstance(z, ureg.Quantity):
+        assert False
+        if height is None:
+            raise ValueError("Must provide height if z is a dimensional quantity")
+        return np.pi / 2 * np.sin(np.pi / height * z)  # normalized
+    else:
+        raise NotImplementedError("z must be either a float or a pint.Quantity")
+    # return 0.5 * (1 + np.cos(0.5 * np.pi / (1 * ureg.m) * z))
+
+
+T_99 = my_input.get_tau()
+my_simulation = Simulation(
+    my_input,
+    t_final=2 * T_99,
+    signal_irr=lambda t: 1 if t < T_99 else 0,
+    signal_sparging=lambda t: 1,
+    profile_pressure_hydrostatic=False,
+    profile_source_T=profile_source_T,
+)
+
+if __name__ == "__main__":
+    # my_simulation.sim_input.E_g *= 1e5
+    # my_simulation.sim_input.E_l *= 1e-5
+    output = my_simulation.solve()
+
+    # # save output to file
+    # output.profiles_to_csv(f"output_{tank_height}m.csv")
+
+    # # plot results
+    # from sparging import plotting
+    # plotting.plot_animation(output)
+
+    # import matplotlib.pyplot as plt
+
+    # plt.plot(output.times, output.inventories_T2_salt)
+    # plt.show()
+    animation.create_animation(output, show_activity=False)

example2.py → sparging_LIBRA1L.py

@@ -18,12 +18,11 @@
 logger = logging.getLogger(__name__)
 logging.basicConfig(level=logging.WARNING)
 
-
 geom = ColumnGeometry(
-    area=0.2 * ureg.m**2,
-    height=1.0 * ureg.m,
-    nozzle_diameter=0.001 * ureg.m,
-    nb_nozzle=10 * ureg.dimensionless,
+    area=170 * ureg.cm**2,
+    height=7 * ureg.cm,
+    nozzle_diameter=1.4 * ureg.mm,
+    nb_nozzle=1 * ureg.dimensionless,
 )
 
 flibe = BreederMaterial(
@@ -33,44 +32,45 @@
 operating_params = OperatingParameters(
     temperature=600 * ureg.celsius,
     P_top=1 * ureg.atm,
-    flow_g_mol=400 * ureg.sccm,
-    tbr=0.1 * ureg("triton / neutron"),
+    flow_g_mol=40 * ureg.sccm,
+    tbr=2e-3 * ureg("triton / neutron"),
     n_gen_rate=1e9 * ureg("neutron / s"),
 )
 
 sparging_params = SpargingParameters(
     h_l=all_correlations("h_l_briggs"),
 )
 
-
 # class method from_parameters that takes in objects like ColumnGeometry, BreederMaterial, OperatingParameters and returns a SimulationInput object with the appropriate correlations for the given parameters. This method should be able to handle cases where some of the parameters are already provided as correlations and should not overwrite them.
 my_input = SimulationInput.from_parameters(
     geom, flibe, operating_params, sparging_params
 )
 logger.info(my_input)
+print(my_input.get_S_T())
+print(my_input.Q_T.to("molT/s"))
 
-
-def profile_source_T(z: pint.Quantity):
-    import numpy as np
-
-    # return np.sin(np.pi / (1 * ureg.m) * z)
-    return 0.5 * (1 + np.cos(0.5 * np.pi / (1 * ureg.m) * z))
-
-
+n_fluence = 2.5e13 * ureg("neutron")
+n_gen_rate = operating_params.n_gen_rate
+t_irr = n_fluence / n_gen_rate
+print(f"t_irr = {t_irr.to('seconds')}")
 my_simulation = Simulation(
     my_input,
-    t_final=3 * ureg.days,
-    signal_irr=lambda t: 1 if t < 12 * ureg.hour else 0,
-    signal_sparging=lambda t: 1,
+    t_final=50 * ureg.days,
+    signal_irr=lambda t: 1 if t < t_irr else 0,
+    signal_sparging=lambda t: 0 if t < t_irr else 1,
+    # signal_sparging=lambda t: 0,
+    profile_pressure_hydrostatic=False,
+    profile_source_T=lambda z: 1,
 )
-output = my_simulation.solve()
 
-# # save output to file
-# output.profiles_to_csv(f"output_{tank_height}m.csv")
+if __name__ == "__main__":
+    output = my_simulation.solve()
 
-# # plot results
-# from sparging import plotting
-# plotting.plot_animation(output)
+    # # save output to file
+    # output.profiles_to_csv(f"output_{tank_height}m.csv")
 
+    # # plot results
+    # from sparging import plotting
+    # plotting.plot_animation(output)
 
-animation.create_animation(output, show_activity=True)
+    animation.create_animation(output, show_activity=False)

src/sparging/animation.py

@@ -43,7 +43,7 @@ def __init__(
         self.x_y = results.x_y
         self.inventories_T2_salt = results.inventories_T2_salt
         self.source_T2 = (
-            None if results.source_T2 is None else np.array(results.source_T2)
+            None if results.sources_T2 is None else np.array(results.sources_T2)
         )
         self.fluxes_T2 = (
             None if results.fluxes_T2 is None else np.array(results.fluxes_T2)

src/sparging/config.py

@@ -9,7 +9,7 @@
 ureg.define(f"molT = {const.N_A} * triton")
 ureg.define(f"molT2 = 2 * {const.N_A} * triton")
 ureg.define("neutron = [neutron] = n")
-ureg.define("sccm = 7.44e-7 mol/s")
+ureg.define("sccm = 7.44e-7 mol/s")  # holds for an ideal gas
 
 
 const_R = const.R * ureg("J/K/mol")  # ideal gas constant

src/sparging/correlations.py

@@ -30,6 +30,7 @@ class CorrelationType(enum.Enum):  # TODO do we really use it ?
     REYNOLDS_NUMBER = "Re"
     BUBBLE_VELOCITY = "u_g0"
     GAS_PHASE_DISPERSION = "E_g"
+    LIQUID_PHASE_DISPERSION = "E_l"
     PRESSURE = "P"
     FLOW_RATE = "flow_g_mol"
     INTERFACIAL_AREA = "a"
@@ -306,10 +307,23 @@ def __contains__(self, key: str | Correlation):
 )
 all_correlations.append(h_l_briggs)
 
+E_l = Correlation(
+    identifier="E_l",
+    function=lambda tank_diameter, u_g0: ureg.Quantity(
+        0.678 * tank_diameter.magnitude**1.4 * u_g0.magnitude**0.3, "m**2/s"
+    ),  # liquid phase axial dispersion coefficient
+    corr_type=CorrelationType.LIQUID_PHASE_DISPERSION,
+    source="Deckwer 1974",
+    description="liquid phase axial dispersion coefficient, assumed equal to diffusivity of tritium in liquid FLiBe",
+    input_units=["m", "m/s"],
+    output_units="m**2/s",
+)
+all_correlations.append(E_l)
+
 E_g = Correlation(
     identifier="E_g",
-    function=lambda tank_diameter, u_g0: get_E_g(
-        diameter=tank_diameter, u_g=u_g0
+    function=lambda tank_diameter, u_g0: (
+        0.2 * ureg("1/m") * tank_diameter**2 * u_g0
     ),  # gas phase axial dispersion coefficient
     corr_type=CorrelationType.GAS_PHASE_DISPERSION,
     source="Malara 1995",
@@ -379,16 +393,16 @@ def __contains__(self, key: str | Correlation):
 )
 all_correlations.append(specific_interfacial_area)
 
-source_T_from_tbr = Correlation(
-    identifier="source_T",
-    function=lambda tbr, n_gen_rate, tank_volume: (
-        tbr * n_gen_rate / tank_volume
+source_T_integral = Correlation(
+    identifier="Q_T",
+    function=lambda tbr, n_gen_rate: (
+        tbr * n_gen_rate
     ),  # source term for tritium generation calculated from TBR and neutron generation rate
     corr_type=CorrelationType.TRITIUM_SOURCE,
-    input_units=["triton/neutron", "neutron/s", "m**3"],
-    output_units="molT/m**3/s",
+    input_units=["triton/neutron", "neutron/s"],
+    output_units="molT/s",
 )
-all_correlations.append(source_T_from_tbr)
+all_correlations.append(source_T_integral)
 
 
 def get_d_b(
@@ -474,10 +488,3 @@ def get_h_briggs(Re: float, Sc: float, D_l: float, d_b: float) -> float:
     Sh = 0.089 * Re**0.69 * Sc**0.33  # Sherwood number
     h_l = Sh * D_l / d_b
     return h_l
-
-
-def get_E_g(diameter: float, u_g: float) -> float:
-    """gas phase axial dispersion coefficient [m2/s], Malara 1995 correlation
-    models dispersion of the gas velocity distribution around the mean bubble velocity"""
-    E_g = 0.2 * ureg("1/m") * diameter**2 * u_g
-    return E_g

src/sparging/inputs.py

@@ -5,6 +5,7 @@
 import inspect
 import numpy as np
 import logging
+from sparging.config import const_R
 
 logger = logging.getLogger(__name__)
 
@@ -46,7 +47,7 @@ class OperatingParameters:
     P_bottom: pint.Quantity | Correlation | None = None
     tbr: pint.Quantity | None = None
     n_gen_rate: pint.Quantity | None = None
-    source_T: pint.Quantity | Correlation | None = None
+    Q_T: pint.Quantity | None = None
 
 
 @dataclass
@@ -57,6 +58,7 @@ class SpargingParameters:
     d_b: pint.Quantity | Correlation | None = None
     rho_g: pint.Quantity | Correlation | None = None
     E_g: pint.Quantity | Correlation | None = None
+    E_l: pint.Quantity | Correlation | None = None
     a: pint.Quantity | Correlation | None = None
 
 
@@ -70,15 +72,47 @@ class SimulationInput:
     h_l: pint.Quantity
     K_s: pint.Quantity
     P_bottom: pint.Quantity
+    rho_l: pint.Quantity
     eps_g: pint.Quantity
     E_g: pint.Quantity
+    E_l: pint.Quantity
     D_l: pint.Quantity
-    source_T: pint.Quantity
+    Q_T: pint.Quantity
+    # normalize_source_T: bool = True
+    # TODO refactoring, put signals in this class
 
     @property
     def volume(self):
         return self.area * self.height
 
+    @property
+    def eps_l(self):
+        return 1 - self.eps_g
+
+    def set_S_T(self, val: pint.Quantity):
+        self.Q_T = (val.to("molT/m**3/s") * self.volume).to("molT/s")
+
+    def get_S_T(self) -> pint.Quantity:
+        return (self.Q_T / self.volume).to("molT/m**3/s")
+
+    def get_tau(self) -> pint.Quantity:
+        """characteristic time of the sparger under the small partial pressure (SPP) approximation"""
+        return (self.eps_l / (self.h_l * self.a)).to("seconds")
+
+    def get_c_T2_SS(self) -> pint.Quantity:
+        return (self.get_S_T() * 1 / (self.h_l * self.a)).to("molT2/m^3")
+
+    def get_PP_number(self) -> pint.Quantity:
+        """Partial pressure number, ratio of the equivalent T concentration at liquid boundary to the bulk liquid concentration.
+        If PP << 1, then we are in the small partial pressure (SPP) regime, and if PP ~= 1, then we are in the partial pressure limited (PPL) regime.
+        """
+        return (
+            self.K_s
+            * (const_R * self.temperature)
+            * self.height
+            / (self.get_tau() * self.u_g0)
+        ).to("dimensionless")
+
     def __post_init__(self):
         # make sure there are only pint.Quantity or callables in the input, otherwise raise an error
         for key in self.__dataclass_fields__.keys():