diff --git a/libra_toolbox/tritium/model.py b/libra_toolbox/tritium/model.py index e5f6bc5..ff9b882 100644 --- a/libra_toolbox/tritium/model.py +++ b/libra_toolbox/tritium/model.py @@ -5,6 +5,8 @@ from libra_toolbox.tritium import ureg +from typing import Callable + SPECIFIC_ACT = 3.57e14 * ureg.Bq * ureg.g**-1 MOLAR_MASS = 6.032 / 2 * ureg.g * ureg.mol**-1 @@ -70,8 +72,8 @@ def __init__( height: pint.Quantity, TBR: pint.Quantity, neutron_rate: pint.Quantity, - k_wall: pint.Quantity, - k_top: pint.Quantity, + k_wall: pint.Quantity | Callable[[pint.Quantity], pint.Quantity], + k_top: pint.Quantity | Callable[[pint.Quantity], pint.Quantity], irradiations: list, ) -> None: @@ -154,7 +156,7 @@ def source(self, t): return self.TBR * self.neutron_rate return 0 * self.TBR * self.neutron_rate - def Q_wall(self, c_salt): + def Q_wall(self, t, c_salt): """ Calculate the release rate of tritium through the wall. @@ -167,10 +169,13 @@ def Q_wall(self, c_salt): Returns: pint.Quantity: The release rate of tritium through the wall. """ + if callable(self.k_wall): + k_wall = self.k_wall(t) + else: + k_wall = self.k_wall + return self.A_wall * k_wall * c_salt - return self.A_wall * self.k_wall * c_salt - - def Q_top(self, c_salt): + def Q_top(self, t, c_salt): """ Calculate the release rate of tritium through the top surface of the salt. @@ -183,7 +188,11 @@ def Q_top(self, c_salt): Returns: pint.Quantity: The release rate of tritium through the top. """ - return self.A_top * self.k_top * c_salt + if callable(self.k_top): + k_top = self.k_top(t) + else: + k_top = self.k_top + return self.A_top * k_top * c_salt def rhs(self, t, c): """ @@ -201,8 +210,8 @@ def rhs(self, t, c): return self.volume.to(ureg.m**3) ** -1 * ( self.source(t).to(ureg.particle * ureg.s**-1) - - self.Q_wall(c).to(ureg.particle * ureg.s**-1) - - self.Q_top(c).to(ureg.particle * ureg.s**-1) + - self.Q_wall(t, c).to(ureg.particle * ureg.s**-1) + - self.Q_top(t, c).to(ureg.particle * ureg.s**-1) ) def _generate_time_intervals(self, t_final): @@ -222,7 +231,7 @@ def _generate_time_intervals(self, t_final): for irr in self.irradiations: t0 = irr[0] - tf = irr[1] + tf = min(irr[1], t_final) if previous_tf is not None: time_intervals.append((previous_tf, t0)) @@ -287,7 +296,7 @@ def integrated_release_top(self): ndarray: array with units same size as the number of time steps, the integrated release of tritium through the top surface. """ - top_release = self.Q_top(self.concentrations) + top_release = self.Q_top(t=self.times, c_salt=self.concentrations) integrated_top = cumulative_trapezoid( top_release.to(ureg.particle * ureg.h**-1).magnitude, self.times.to(ureg.h).magnitude, @@ -304,7 +313,7 @@ def integrated_release_wall(self): ndarray: array with units same size as the number of time steps, the integrated release of tritium through the walls. """ - wall_release = self.Q_wall(self.concentrations) + wall_release = self.Q_wall(t=self.times, c_salt=self.concentrations) integrated_wall = cumulative_trapezoid( wall_release.to(ureg.particle * ureg.h**-1).magnitude, self.times.to(ureg.h).magnitude, diff --git a/libra_toolbox/tritium/plotting.py b/libra_toolbox/tritium/plotting.py index b958444..f5933d0 100644 --- a/libra_toolbox/tritium/plotting.py +++ b/libra_toolbox/tritium/plotting.py @@ -285,14 +285,14 @@ def plot_salt_inventory(model: Model, **kwargs): def plot_top_release(model: Model, **kwargs): - top_release = model.Q_top(model.concentrations) + top_release = model.Q_top(t=model.times, c_salt=model.concentrations) top_release = quantity_to_activity(top_release).to(ureg.Bq * ureg.h**-1) l = plt.plot(model.times.to(ureg.day), top_release, **kwargs) return l def plot_wall_release(model: Model, **kwargs): - wall_release = model.Q_wall(model.concentrations) + wall_release = model.Q_wall(t=model.times, c_salt=model.concentrations) wall_release = quantity_to_activity(wall_release).to(ureg.Bq * ureg.h**-1) l = plt.plot(model.times.to(ureg.day), wall_release, **kwargs) return l diff --git a/test/tritium/test_mass_conservation.py b/test/tritium/test_model.py similarity index 67% rename from test/tritium/test_mass_conservation.py rename to test/tritium/test_model.py index d4e62d4..6a2f58e 100644 --- a/test/tritium/test_mass_conservation.py +++ b/test/tritium/test_model.py @@ -32,3 +32,25 @@ def test_simple_case(TBR): assert np.isclose(model.concentrations[-1], 0) assert np.isclose(expected_total_production, computed_total_production, rtol=1e-2) + + +def test_mass_transport_coeffs_time_dependent(): + + def k_top(t): + return 2 / (1 + t.magnitude) * ureg.m * ureg.s**-1 + + def k_wall(t): + return 3 / (1 + t.magnitude) * ureg.m * ureg.s**-1 + + model = Model( + radius=2 * ureg.m, + height=4 * ureg.m, + TBR=1 * ureg.particle * ureg.neutron**-1, + k_top=k_top, + k_wall=k_wall, + irradiations=[(0 * ureg.s, 10 * ureg.s), (60 * ureg.s, 70 * ureg.s)], + neutron_rate=30 * ureg.neutron * ureg.s**-1, + ) + + # run + model.run(100 * ureg.s)