diff --git a/sparging101.py b/sparging101.py new file mode 100644 index 0000000..231559e --- /dev/null +++ b/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) diff --git a/example2.py b/sparging_LIBRA1L.py similarity index 57% rename from example2.py rename to sparging_LIBRA1L.py index 3382f5c..fe29bb7 100644 --- a/example2.py +++ b/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,8 +32,8 @@ 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"), ) @@ -42,35 +41,36 @@ 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) diff --git a/src/sparging/animation.py b/src/sparging/animation.py index c58ef5a..6baacfd 100644 --- a/src/sparging/animation.py +++ b/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) diff --git a/src/sparging/config.py b/src/sparging/config.py index 1ec8f5c..bf9bd10 100644 --- a/src/sparging/config.py +++ b/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 diff --git a/src/sparging/correlations.py b/src/sparging/correlations.py index 9f8f09b..e8f4c0c 100644 --- a/src/sparging/correlations.py +++ b/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 diff --git a/src/sparging/inputs.py b/src/sparging/inputs.py index b61b7e4..28a5597 100644 --- a/src/sparging/inputs.py +++ b/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(): diff --git a/src/sparging/model.py b/src/sparging/model.py index 94fc5bd..49ca119 100644 --- a/src/sparging/model.py +++ b/src/sparging/model.py @@ -1,3 +1,4 @@ +from __future__ import annotations from mpi4py import MPI import dolfinx import basix @@ -16,7 +17,7 @@ import json from pathlib import Path -from sparging.config import ureg +from sparging.config import ureg, const_g from sparging.inputs import SimulationInput @@ -24,6 +25,7 @@ if TYPE_CHECKING: import pint +from collections.abc import Callable hours_to_seconds = 3600 days_to_seconds = 24 * hours_to_seconds @@ -37,20 +39,24 @@ @dataclass -class SimulationResults: +class SimulationResults: # TODO implement pint in this class # TODO change list to np.array times: list c_T2_solutions: list y_T2_solutions: list + J_T2_solutions: list x_ct: np.ndarray x_y: np.ndarray inventories_T2_salt: np.ndarray - source_T2: list + sources_T2: list fluxes_T2: list sim_input: SimulationInput + dt: int | None = None + dx: int | None = None keys_to_ignore_results = [ # TODO do it the other way: keys_to_include_results "c_T2_solutions", "y_T2_solutions", + "J_T2_solutions", "x_ct", "x_y", "inventories_T2_salt", @@ -136,11 +142,25 @@ def profiles_to_csv(self, output_path: Path): class Simulation: sim_input: SimulationInput t_final: pint.Quantity - signal_irr: callable[pint.Quantity] - signal_sparging: callable[pint.Quantity] - profile_source_T: callable[pint.Quantity] | None = None - - def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None): + signal_irr: Callable[[pint.Quantity], float] + signal_sparging: Callable[[pint.Quantity], float] + profile_source_T: Callable[[float], float] | None = None + """callable = f:[0,1] -> R+, it takes a dimensionless coordinate: (z / height)""" + profile_pressure_hydrostatic: bool = False + dispersion_on: bool = True + + def hydrostatic_pressure(self, x: pint.Quantity) -> pint.Quantity: + """returns the hydrostatic pressure at a given height x in the tank given P_bottom""" + rho = self.sim_input.rho_l + g = const_g + return (self.sim_input.P_bottom + rho * g * x).to("Pa") + + def solve( + self, + dt: pint.Quantity | None = None, + dx: pint.Quantity | None = None, + fast_solve: bool = False, + ) -> SimulationResults: # unpack pint.Quantities t_final = self.t_final.to("seconds").magnitude tank_height = self.sim_input.height.to("m").magnitude @@ -153,12 +173,21 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None T = self.sim_input.temperature.to("K").magnitude eps_g = self.sim_input.eps_g.to("dimensionless").magnitude E_g = self.sim_input.E_g.to("m**2/s").magnitude - D_l = self.sim_input.D_l.to("m**2/s").magnitude + E_l = self.sim_input.E_l.to("m**2/s").magnitude + D_l = self.sim_input.D_l.to("m**2/s").magnitude # not needed (included in h_l) u_g0 = self.sim_input.u_g0.to("m/s").magnitude - source_T2 = self.sim_input.source_T.to("molT2/s/m**3").magnitude + Q_T2 = self.sim_input.Q_T.to("molT2/s").magnitude - dt = dt.to("seconds").magnitude if dt is not None else t_final / 1000 - dx = dx.to("m").magnitude if dx is not None else tank_height / 1000 + dt = ( + dt.to("seconds").magnitude + if dt is not None + else (t_final / 1000 if not fast_solve else t_final / 50) + ) + dx = ( + dx.to("m").magnitude + if dx is not None + else (tank_height / 1000 if not fast_solve else tank_height / 50) + ) eps_l = 1 - eps_g # MESH AND FUNCTION SPACES @@ -184,37 +213,62 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None h_l_const = dolfinx.fem.Constant(mesh, PETSc.ScalarType(h_l)) - gen_T2_mag = dolfinx.fem.Constant( - mesh, source_T2 * self.signal_irr(0 * ureg.s) + gen_T2_ave = dolfinx.fem.Constant( + mesh, Q_T2 / tank_volume * self.signal_irr(0 * ureg.s) ) # magnitude of the generation term - gen_T2 = None + if self.profile_source_T is not None: # spatially varying profile is provided - gen_T2_prof = dolfinx.fem.Function(V_profile) - gen_T2_prof.interpolate( - lambda x: x[0] * 0 + self.profile_source_T(x[0] * ureg.m) + arbitrary_profile = dolfinx.fem.Function(V_profile) + arbitrary_profile.interpolate( + lambda x: x[0] * 0 + self.profile_source_T(x[0] / tank_height) ) - gen_T2 = gen_T2_mag * gen_T2_prof + profile_integral = dolfinx.fem.assemble_scalar( + dolfinx.fem.form( + arbitrary_profile * 1 / tank_height * ufl.dx # TODO + ) # dimensionless integral + ) + normalized_profile = ( + arbitrary_profile / profile_integral + ) # normalize profile so that its integral over the dimensionless height is 1 else: # homogeneous generation - gen_T2 = gen_T2_mag + normalized_profile = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1.0)) + + gen_T2 = gen_T2_ave * normalized_profile + + P_prof = dolfinx.fem.Function(V_profile) + if self.profile_pressure_hydrostatic: + P_prof.interpolate( + lambda x: self.hydrostatic_pressure(x[0] * ureg.m).magnitude + ) + else: + P_prof.interpolate(lambda x: x[0] * 0 + P_0) + + P = P_prof # VARIATIONAL FORMULATION # mass transfer rate J_T2 = ( - a * h_l_const * (c_T2 - K_s * (P_0 * y_T2 + EPS)) + a * h_l_const * (c_T2 - K_s * (P * y_T2 + EPS)) ) # TODO pressure shouldn't be a constant (use hydrostatic pressure profile), how to deal with this ? -> use fem.Expression ? F = 0 # variational formulation # transient terms F += eps_l * ((c_T2 - c_T2_n) / dt) * v_c * ufl.dx - F += eps_g * 1 / (const.R * T) * (P_0 * (y_T2 - y_T2_n) / dt) * v_y * ufl.dx - - # diffusion/dispersion terms #TODO shouldn't use D_l, transport of T in liquid is dominated by dispersive effects due to gas sparging, find dispersion coeff for steady liquid in gas bubbles - F += eps_l * D_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx - - # NOTE remove diffusive term for gas for now for mass balance - # F += eps_g * E_g * ufl.dot(ufl.grad(P_0 * y_T2), ufl.grad(v_y)) * ufl.dx + F += eps_g * 1 / (const.R * T) * (P * (y_T2 - y_T2_n) / dt) * v_y * ufl.dx + + # dispersive terms + if self.dispersion_on is True: + F += eps_l * E_l * ufl.dot(ufl.grad(c_T2), ufl.grad(v_c)) * ufl.dx + F += ( + eps_g + * E_g + * 1 + / (const.R * T) + * ufl.dot(ufl.grad(P * y_T2), ufl.grad(v_y)) + * ufl.dx + ) # mass exchange (coupling term) F += J_T2 * v_c * ufl.dx - J_T2 * v_y * ufl.dx @@ -226,7 +280,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None F += ( 1 / (const.R * T) - * ufl.inner(ufl.dot(ufl.grad(P_0 * y_T2), vel), v_y) + * ufl.inner(ufl.dot(ufl.grad(P * y_T2), vel), v_y) * ufl.dx ) @@ -241,12 +295,7 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None dolfinx.fem.Constant(mesh, 0.0), dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, gas_inlet_facets), V.sub(1), - ) - bc2 = dolfinx.fem.dirichletbc( - dolfinx.fem.Constant(mesh, 0.0), - dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, gas_outlet_facets), - V.sub(0), - ) + ) # y_T2 = 0 at gas inlet # Custom measure all_facets = np.concatenate((gas_inlet_facets, gas_outlet_facets)) @@ -260,12 +309,17 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None problem = NonlinearProblem( F, u, - bcs=[bc1, bc2], + bcs=[bc1], # Neumann BCs on c_T2 at inlet and outlet are naturally enforced petsc_options_prefix="librasparge", # petsc_options={"snes_monitor": None}, ) # initialise post processing + + # we define a function and an expression for J_T2 for use in post processing + J_T2_func = dolfinx.fem.Function(V_profile) + J_T2_expr = dolfinx.fem.Expression(J_T2, V_profile.element.interpolation_points) + V0_ct, ct_dofs = u.function_space.sub(0).collapse() coords = V0_ct.tabulate_dof_coordinates()[:, 0] ct_sort_coords = np.argsort(coords) @@ -276,46 +330,48 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None y_sort_coords = np.argsort(coords) x_y = coords[y_sort_coords] - times = [] - c_T2_solutions = [] - y_T2_solutions = [] - sources_T2 = [] - fluxes_T2 = [] - inventories_T2_salt = [] - - # SOLVE - t = 0 - while t < t_final: - gen_T2_mag.value = source_T2 * self.signal_irr(t * ureg.s) - h_l_const.value = h_l * self.signal_sparging(t * ureg.s) - """ utiliser ufl.conditional TODO""" - - problem.solve() - - # update previous solution - u_n.x.array[:] = u.x.array[:] - - # post process + coords_profile = V_profile.tabulate_dof_coordinates()[:, 0] + num_profile_dofs = ( + V_profile.dofmap.index_map.size_local + ) * V_profile.dofmap.index_map_bs + + profile_dofs = np.arange(num_profile_dofs) + profile_sort_coords = np.argsort(coords_profile) + # NOTE currently we don't use x_profile and use another x in the plotting script + x_profile = coords_profile[profile_sort_coords] + + # NOTE maybe we could take this function out and it would take a SimulationResults object as input + u + other things... + def post_process(t): + """ + Post-process the solution at time t. + Extract solution profiles, compute fluxes and inventories, and store them in lists for later analysis. + """ c_T2_post, y_T2_post = u.split() c_T2_vals = u.x.array[ct_dofs][ct_sort_coords] y_T2_vals = u.x.array[y_dofs][y_sort_coords] - - # store time and solution - # TODO give units to the results + J_T2_func.interpolate(J_T2_expr) + J_T2_vals = J_T2_func.x.array[profile_dofs][ct_sort_coords] times.append(t) c_T2_solutions.append(c_T2_vals.copy()) y_T2_solutions.append(y_T2_vals.copy()) + J_T2_solutions.append(J_T2_vals.copy()) sources_T2.append( - source_T2 * self.signal_irr(t * ureg.s) * tank_volume + Q_T2 * self.signal_irr(t * ureg.s) ) # total T generation rate in the tank [mol/s] TODO useless: signal_irr is already given + n = ufl.FacetNormal(mesh) + + # flux_T2 = dolfinx.fem.assemble_scalar( + # dolfinx.fem.form( + # eps_g * vel_x * P / (const.R * T) * y_T2_post * tank_area * ds(2) + # ) + # ) # total T flux at the outlet [mol/s] flux_T2 = dolfinx.fem.assemble_scalar( dolfinx.fem.form( - tank_area * vel_x * P_0 / (const.R * T) * y_T2_post * ds(2) + vel_x * P / (const.R * T) * y_T2_post * tank_area * ds(2) ) - ) # total T flux at the outlet [mol/s] - + ) # TODO replace with integral of J over volume # flux_T_inlet = dolfinx.fem.assemble_scalar( # dolfinx.fem.form( # tank_area @@ -328,15 +384,16 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None # ) # ) # total T dispersive flux at the inlet [mol/s] - n = ufl.FacetNormal(mesh) flux_T2_inlet = dolfinx.fem.assemble_scalar( - dolfinx.fem.form(-E_g * ufl.inner(ufl.grad(P_0 * y_T2_post), n) * ds(1)) + dolfinx.fem.form( + -eps_g * E_g * ufl.inner(ufl.grad(P * y_T2_post), n) * ds(1) + ) ) # total T dispersive flux at the inlet [Pa T2 /s/m2] - flux_T2_inlet *= 1 / (const.R * T) * T2_to_T # convert to mol T/s/m2 - flux_T2_inlet *= tank_area # convert to mol T/s + flux_T2_inlet *= 1 / (const.R * T) # mol T2/s/m2 + flux_T2_inlet *= tank_area # convert to molT2/s - # fluxes_T2.append(flux_T2 + flux_T2_inlet) - fluxes_T2.append(flux_T2) + fluxes_T2.append(flux_T2 + flux_T2_inlet) + # fluxes_T2.append(flux_T2) inventory_T2_salt = dolfinx.fem.assemble_scalar( dolfinx.fem.form(c_T2_post * ufl.dx) @@ -344,6 +401,30 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None inventory_T2_salt *= tank_area # get total amount of T2 in [mol] inventories_T2_salt.append(inventory_T2_salt) + # TODO Initialize results storage with zeros -> there's an inconsistency: times[0] = 0 but inventories[0] != 0 -> screws comparison with analytical solutions + t = 0 + times = [] + c_T2_solutions = [] + y_T2_solutions = [] + J_T2_solutions = [] + sources_T2 = [] + fluxes_T2 = [] + inventories_T2_salt = [] + post_process(t) + + # SOLVE + while t < t_final: + # update time-dependent terms + gen_T2_ave.value = Q_T2 / tank_volume * self.signal_irr(t * ureg.s) + h_l_const.value = h_l * self.signal_sparging(t * ureg.s) + + problem.solve() + + # update previous solution + u_n.x.array[:] = u.x.array[:] + + post_process(t) + # advance time t += dt @@ -355,11 +436,14 @@ def solve(self, dt: pint.Quantity | None = None, dx: pint.Quantity | None = None times=times, c_T2_solutions=c_T2_solutions, y_T2_solutions=y_T2_solutions, + J_T2_solutions=J_T2_solutions, x_ct=x_ct, x_y=x_y, inventories_T2_salt=inventories_T2_salt, - source_T2=sources_T2, + sources_T2=sources_T2, fluxes_T2=fluxes_T2, sim_input=self.sim_input, + dt=dt, + dx=dx, ) return results diff --git a/test/standard_input.json b/test/standard_input.json index 17b74fa..77f4c5f 100644 --- a/test/standard_input.json +++ b/test/standard_input.json @@ -7,8 +7,10 @@ "h_l": "2.750e-05 m / s", "K_s": "6.366e-04 mol / m ** 3 / Pa", "P_bottom": "1.208e+05 Pa", + "rho_l": "1.991e+03 kg / m ** 3", "eps_g": "4.091e-04", "E_g": "1.111e-02 m ** 2 / s", + "E_l": "1.648e-01 m ** 2 / s", "D_l": "2.857e-09 m ** 2 / s", "source_T": "8.303e-16 molT / m ** 3 / s" } \ No newline at end of file diff --git a/test/test_solve.py b/test/test_solve.py index aceb223..963fdd4 100644 --- a/test/test_solve.py +++ b/test/test_solve.py @@ -4,6 +4,7 @@ import pytest import dataclasses from pint import DimensionalityError +import numpy as np standard_input = SimulationInput( @@ -13,12 +14,14 @@ temperature=600 * ureg.celsius, a=0.5 * ureg("1/m"), h_l=3e-5 * ureg("m/s"), + rho_l=2000 * ureg("kg/m^3"), K_s=1e-4 * ureg("mol/m**3/Pa"), P_bottom=1.2 * ureg.bar, eps_g=0.001 * ureg.dimensionless, E_g=1e-2 * ureg("m^2/s"), + E_l=1e-1 * ureg("m^2/s"), D_l=3e-9 * ureg("m^2/s"), - source_T=8e-16 * ureg("molT/m^3/s"), + Q_T=1e8 * ureg("T/s"), ) @@ -83,3 +86,53 @@ def test_model_solve_wrong_argument(standard_simulation): """ with pytest.raises(AttributeError, match="object has no attribute 'to'"): standard_simulation.solve(dt=0.01) + + +@pytest.mark.parametrize("case", ("nominal", "volume_change", "unormalized_profile")) +def test_source_T_normalization(case): + """Tests that the source term is correctly normalized so that the total inventory after a fixed time is always the same + no matter the volume or the given profile function.""" + Q_T = standard_input.Q_T + t_final = 50 * ureg.seconds + # BUILD + match case: + case "nominal": + my_simulation = Simulation( + standard_input, + t_final=t_final, + signal_irr=lambda t: 1, + signal_sparging=lambda t: 0, + ) + case "volume_change": + modified_input = dataclasses.replace( + standard_input, + height=23.2 * standard_input.height, + area=0.123 * standard_input.area, + ) + my_simulation = Simulation( + modified_input, + t_final=t_final, + signal_irr=lambda t: 1, + signal_sparging=lambda t: 0, + ) + case "unormalized_profile": + my_simulation = Simulation( + standard_input, + t_final=t_final, + signal_irr=lambda t: 1, + signal_sparging=lambda t: 0, + profile_source_T=lambda xi: 3 + xi, # not normalized + ) + # RUN + output = my_simulation.solve(fast_solve=True) + + # TEST + n_result = ureg.Quantity(output.inventories_T2_salt[-1], "molT2").magnitude + # n_theory = (Q_T * ureg.Quantity(output.times[-1], "s")).to("molT2").magnitude # TODO change to this + n_theory = (Q_T * t_final).to("molT2").magnitude + # print(len(output.times), len(output.inventories_T2_salt)) + # print(output.times[0], output.inventories_T2_salt[0]) + # print(output.times[-1], output.inventories_T2_salt[-1]) + assert np.isclose(n_result, n_theory, atol=0, rtol=1e-2), print( + f"n_result = {n_result}, n_theory = {n_theory}" + ) # TODO reduce tolerance when vectors initialization with zero is fixed