u_g is a profile now

Author: hmonroedd

examples/sparging_standard.py

@@ -39,7 +39,6 @@
     dispersion_on=True,
     constant_profiles=False,
 )
-
 if __name__ == "__main__":
     # my_simulation.exports = ["pressure", "J_T2"]
     output = my_simulation.solve(fast_solve=True)
@@ -49,7 +48,7 @@
     # P = output.exported_fields["pressure"]
     # plt.plot(output.x_ct, P)
     # output.exports_to_csv(FOLDER)
-    output.export = ["pressure, "]
+    # output.export = ["pressure, "]
 
     output.profiles_to_csv(FOLDER)
     output.profiles_to_cdf(FOLDER)

src/sparging/correlations.py

@@ -28,7 +28,7 @@ class CorrelationType(enum.Enum):  # TODO do we really use it ?
     MORTON_NUMBER = "Mo"
     SCHMIDT_NUMBER = "Sc"
     REYNOLDS_NUMBER = "Re"
-    SUPERFICIAL_GAS_VELOCITY = "u_g0"
+    SUPERFICIAL_GAS_VELOCITY = "u_g"
     BUBBLE_VELOCITY = "v_g0"
     GAS_PHASE_DISPERSION = "E_g"
     LIQUID_PHASE_DISPERSION = "E_l"
@@ -308,16 +308,6 @@ def get_list(self, corr_type: CorrelationType) -> list[Correlation]:
 )
 all_correlations.append(v_g0)
 
-# superficial gas velocity at the inlet
-u_g0 = Correlation(
-    identifier="u_g0",
-    function=lambda Vdot_g0, area: (Vdot_g0 / area).to("m/s"),
-    corr_type=CorrelationType.SUPERFICIAL_GAS_VELOCITY,
-    input_units=["m^3/s", "m^2"],
-    output_units="m/s",
-)
-all_correlations.append(u_g0)
-
 
 h_l_higbie = Correlation(
     identifier="h_l_higbie",
@@ -361,27 +351,28 @@ def get_list(self, corr_type: CorrelationType) -> list[Correlation]:
 # liquid phase axial dispersion coefficient
 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"
+    function=lambda tank_diameter, u_g: ureg.Quantity(
+        0.678 * tank_diameter.magnitude**1.4 * u_g(0 * ureg.m).magnitude ** 0.3,
+        "m**2/s",
     ),
     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"],
+    input_units=["m", PROFILE],
     output_units="m**2/s",
 )
 all_correlations.append(E_l)
 
 # gas phase axial dispersion coefficient
 E_g = Correlation(
     identifier="E_g",
-    function=lambda tank_diameter, u_g0: (
-        0.2 * ureg("1/m") * tank_diameter**2 * u_g0
+    function=lambda tank_diameter, u_g: (
+        0.2 * ureg("1/m") * tank_diameter**2 * u_g(0 * ureg.m)
     ),  # gas phase axial dispersion coefficient
     corr_type=CorrelationType.GAS_PHASE_DISPERSION,
     source="Malara 1995",
     description="gas phase axial dispersion coefficient [m2/s], Malara 1995",
-    input_units=["m", "m/s"],
+    input_units=["m", PROFILE],
     output_units="m**2/s",
 )
 all_correlations.append(E_g)
@@ -585,3 +576,15 @@ def get_h_briggs(Re: float, Sc: float, D_l: float, d_b: float) -> float:
     description="specific interfacial area profile",
 )
 all_correlations.append(a)
+
+u_g = Profile(
+    identifier="u_g",
+    function=lambda Vdot_g0, area, P_g: (
+        lambda z: (Vdot_g0 / area).to("m/s") * (P_g(0 * ureg.m) / P_g(z))
+    ),
+    corr_type=CorrelationType.SUPERFICIAL_GAS_VELOCITY,
+    input_units=["m^3/s", "m^2", PROFILE],
+    output_units="m/s",
+    description="superficial gas velocity profile",
+)
+all_correlations.append(u_g)

src/sparging/inputs.py

@@ -65,7 +65,6 @@ def copy(self):
 class SpargingParameters:
     h_l: pint.Quantity | Correlation
     eps_g: pint.Quantity | Correlation | None = None
-    u_g0: pint.Quantity | Correlation | None = None
     d_b: pint.Quantity | Correlation | None = None
     rho_g: pint.Quantity | Correlation | None = None
     E_g: pint.Quantity | Correlation | None = None
@@ -80,7 +79,6 @@ def copy(self):
 class SimulationInput:
     height: pint.Quantity
     area: pint.Quantity
-    u_g0: pint.Quantity
     temperature: pint.Quantity
     h_l: pint.Quantity
     K_s: pint.Quantity
@@ -100,7 +98,6 @@ class SimulationInput:
     required_keys = (
         "height",
         "area",
-        "u_g0",
         "temperature",
         "h_l",
         "K_s",
@@ -114,6 +111,7 @@ class SimulationInput:
         "P_g",
         "eps_g",
         "a",
+        "u_g",
     )  # these parameters will be used to solve the model
     graph: nx.Graph | None = None
     """ Stores the intermediate parameters and their relationships that built the SimulationInput. 
@@ -128,6 +126,7 @@ class SimulationInput:
     P_g: Callable[[pint.Quantity], pint.Quantity] | None = None
     eps_g: Callable[[pint.Quantity], pint.Quantity] | None = None
     a: Callable[[pint.Quantity], pint.Quantity] | None = None
+    u_g: Callable[[pint.Quantity], pint.Quantity] | None = None
 
     @property
     def volume(self):
@@ -138,21 +137,25 @@ def a_0(self):
         return self.a(0 * ureg.m)
 
     @property
-    def eps_g_0(self):
+    def eps_g0(self):
         return self.eps_g(0 * ureg.m)
 
     @property
-    def eps_l_0(self):
-        return 1 - self.eps_g_0
+    def eps_l0(self):
+        return 1 - self.eps_g0
 
     @property
-    def P_l_0(self):
+    def P_l0(self):
         return self.P_l(0 * ureg.m)
 
     @property
-    def P_g_0(self):
+    def P_g0(self):
         return self.P_g(0 * ureg.m)
 
+    @property
+    def u_g0(self):
+        return self.u_g(0 * ureg.m)
+
     def set_S_T(self, val: pint.Quantity):
         self.Q_T = (val.to("molT/m**3/s") * self.volume).to("molT/s")
 
@@ -161,7 +164,7 @@ def get_S_T(self) -> pint.Quantity:
 
     def get_tau(self) -> pint.Quantity:
         """characteristic time of the sparger under the small partial pressure (SPP) approximation"""
-        return (self.eps_l_0 / (self.h_l * self.a_0)).to("seconds")
+        return (self.eps_l0 / (self.h_l * self.a_0)).to("seconds")
 
     def get_c_T2_SS(self) -> pint.Quantity:
         return (self.get_S_T() * 1 / (self.h_l * self.a_0)).to("molT2/m^3")
@@ -177,7 +180,7 @@ def get_Pi_number(self) -> pint.Quantity:
             * self.height
             * self.h_l
             * self.a_0
-            / (self.eps_g_0 * self.graph.nodes["v_g0"]["value"])
+            / (self.eps_g0 * self.graph.nodes["v_g0"]["value"])
         ).to("dimensionless")
 
     def get_dP_dx(self) -> pint.Quantity:
@@ -196,16 +199,13 @@ def get_Bo(self) -> pint.Quantity:
         returns Bodenstein number = ratio of convective to dispersive transport for the gas phase
         corresponds to Peclet number at the scale of the tank
         """
-        # return (self.eps_g * self.u_g0 * self.height / self.E_g).to("dimensionless")
-        return (
-            (self.graph.nodes["Vdot_g0"]["value"] / self.area) * self.height / self.E_g
-        ).to("dimensionless")
+        return (self.u_g0 * self.height / self.E_g).to("dimensionless")
 
     def test_eps_g(
         self,
     ):  # to see if the two definitions of superficial velocity are consistent -> TODO remove
         print(
-            f"{self.eps_g_0 * self.graph.nodes['v_g0']['value']} vs {self.graph.nodes['Vdot_g0']['value'] / self.area} vs {self.graph.nodes['u_g0']['value']}"
+            f"{self.eps_g0 * self.graph.nodes['v_g0']['value']} vs {self.graph.nodes['Vdot_g0']['value'] / self.area} vs {self.u_g0}"
         )
 
     def __post_init__(self):

src/sparging/model.py

@@ -307,7 +307,6 @@ def solve(
         T = self.sim_input.temperature.to("K").magnitude
         E_g = self.sim_input.E_g.to("m**2/s").magnitude
         E_l = self.sim_input.E_l.to("m**2/s").magnitude
-        u_g0 = self.sim_input.u_g0.to("m/s").magnitude
         Q_T2 = self.sim_input.Q_T.to("molT2/s").magnitude
 
         dt = (
@@ -337,40 +336,13 @@ def solve(
         u_n = dolfinx.fem.Function(V)
         v_c, v_y = ufl.TestFunctions(V)
 
-        # set initial concentration
-        c_T2_0_ufl_expr = dolfinx.fem.Constant(
-            mesh, self.sim_input.c_T2_0.to("molT2/m**3").magnitude
-        ) * self.normalize_profile(
-            self.sim_input.profile_c_T2_0, tank_height, mesh, V_profile
-        )
-        u_n.sub(0).interpolate(
-            dolfinx.fem.Expression(
-                c_T2_0_ufl_expr, V.sub(0).element.interpolation_points
-            )
-        )
-
-        c_T2, y_T2 = ufl.split(u)
-        c_T2_n, y_T2_n = ufl.split(u_n)
-
-        vel_x = u_g0  # TODO velocity should vary with hydrostatic pressure
-        vel = dolfinx.fem.Constant(mesh, PETSc.ScalarType([vel_x]))
-
-        h_l_const = dolfinx.fem.Constant(mesh, PETSc.ScalarType(h_l))
-
-        gen_T2_ave = dolfinx.fem.Constant(
-            mesh, Q_T2 / tank_volume * self.sim_input.signal_irr(0 * ureg.s)
-        )  # magnitude of the generation term
-
-        gen_T2 = gen_T2_ave * self.normalize_profile(
-            self.sim_input.profile_source_T, tank_height, mesh, V_profile
-        )
-
         # define spatially varying profiles
         if not self.constant_profiles:
             eps_g = dolfinx.fem.Function(V_profile)
             eps_l = dolfinx.fem.Function(V_profile)
             a = dolfinx.fem.Function(V_profile)
             P_g = dolfinx.fem.Function(V_profile)
+            u_g = dolfinx.fem.Function(V_profile)
 
             eps_g.interpolate(
                 lambda x: (
@@ -387,16 +359,49 @@ def solve(
             P_g.interpolate(
                 lambda x: self.sim_input.P_g(x[0] * ureg.m).to("Pa").magnitude
             )
+            u_g.interpolate(
+                lambda x: self.sim_input.u_g(x[0] * ureg.m).to("m/s").magnitude
+            )
         else:
             # use values at z=0 for constant profiles
-            eps_g_0 = self.sim_input.eps_g_0.to("dimensionless").magnitude
+            eps_g0 = self.sim_input.eps_g0.to("dimensionless").magnitude
             a_0 = self.sim_input.a_0.to("1/m").magnitude
-            P_g_0 = self.sim_input.P_g_0.to("Pa").magnitude
+            P_g0 = self.sim_input.P_g0.to("Pa").magnitude
+            u_g0 = self.sim_input.u_g0.to("m/s").magnitude
 
-            eps_g = dolfinx.fem.Constant(mesh, PETSc.ScalarType(eps_g_0))
-            eps_l = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1 - eps_g_0))
+            eps_g = dolfinx.fem.Constant(mesh, PETSc.ScalarType(eps_g0))
+            eps_l = dolfinx.fem.Constant(mesh, PETSc.ScalarType(1 - eps_g0))
             a = dolfinx.fem.Constant(mesh, PETSc.ScalarType(a_0))
-            P_g = dolfinx.fem.Constant(mesh, PETSc.ScalarType(P_g_0))
+            P_g = dolfinx.fem.Constant(mesh, PETSc.ScalarType(P_g0))
+            u_g = dolfinx.fem.Constant(mesh, PETSc.ScalarType(u_g0))
+
+        # set initial concentration
+        c_T2_0_ufl_expr = dolfinx.fem.Constant(
+            mesh, self.sim_input.c_T2_0.to("molT2/m**3").magnitude
+        ) * self.normalize_profile(
+            self.sim_input.profile_c_T2_0, tank_height, mesh, V_profile
+        )
+        u_n.sub(0).interpolate(
+            dolfinx.fem.Expression(
+                c_T2_0_ufl_expr, V.sub(0).element.interpolation_points
+            )
+        )
+
+        c_T2, y_T2 = ufl.split(u)
+        c_T2_n, y_T2_n = ufl.split(u_n)
+
+        # vel_x = u_g0  # TODO velocity should vary with hydrostatic pressure
+        # vel = dolfinx.fem.Constant(mesh, PETSc.ScalarType([vel_x]))
+
+        h_l_const = dolfinx.fem.Constant(mesh, PETSc.ScalarType(h_l))
+
+        gen_T2_ave = dolfinx.fem.Constant(
+            mesh, Q_T2 / tank_volume * self.sim_input.signal_irr(0 * ureg.s)
+        )  # magnitude of the generation term
+
+        gen_T2 = gen_T2_ave * self.normalize_profile(
+            self.sim_input.profile_source_T, tank_height, mesh, V_profile
+        )
 
         # VARIATIONAL FORMULATION
 
@@ -431,7 +436,9 @@ def solve(
         F += (
             1
             / (const.R * T)
-            * ufl.inner(ufl.dot(ufl.grad(P_g * y_T2), vel), v_y)
+            # * ufl.inner(ufl.dot(ufl.grad(P_g * y_T2), vel), v_y)
+            * ufl.grad(u_g * P_g * y_T2)[0]
+            * v_y
             * ufl.dx
         )
 
@@ -458,13 +465,7 @@ def solve(
 
         # Danckwert BC at gas inlet
         P_T2_inlet = 0
-        F += (
-            1
-            / (const.R * T)
-            * vel_x
-            * ufl.inner((P_g * y_T2 - P_T2_inlet), v_y)
-            * ds(1)
-        )
+        F += 1 / (const.R * T) * u_g * ufl.inner((P_g * y_T2 - P_T2_inlet), v_y) * ds(1)
 
         # n = ufl.FacetNormal(mesh)
 
@@ -527,7 +528,7 @@ def post_process(t):
 
             flux_T2 = dolfinx.fem.assemble_scalar(
                 dolfinx.fem.form(
-                    vel_x * P_g / (const.R * T) * y_T2_post * tank_area * ds(2)
+                    u_g * P_g / (const.R * T) * y_T2_post * tank_area * ds(2)
                 )
             )
             flux_T2_2 = dolfinx.fem.assemble_scalar(