Merge pull request #94 from geo-fluid-dynamics/paper/gallium-melting

sapphire/simulations/enthalpy_porosity.py

@@ -114,6 +114,8 @@ def momentum(self):
 
     def energy(self):
 
+        Re = self.reynolds_number
+
         Pr = self.prandtl_number
 
         Ste = self.stefan_number
@@ -139,7 +141,7 @@ def energy(self):
         dx = fe.dx(degree = self.quadrature_degree)
 
         return (psi_T*(CT_t + 1./Ste*phil_t + dot(u, grad(C*T))) \
-            + 1./Pr*dot(grad(psi_T), k*grad(T)))*dx
+            + 1./(Re*Pr)*dot(grad(psi_T), k*grad(T)))*dx
 
     def extra_time_discrete_terms(self):
 
@@ -303,6 +305,8 @@ def strong_residual(sim, solution):
     r_u += d*u
 
 
+    Re = sim.reynolds_number
+
     Pr = sim.prandtl_number
 
     Ste = sim.stefan_number
@@ -320,6 +324,6 @@ def strong_residual(sim, solution):
     k = phase_dependent_material_property(k_sl)(phil)
 
     r_T = diff(C*T, t) + 1./Ste*diff(phil, t) + dot(u, grad(C*T)) \
-        - 1./Pr*div(k*grad(T))
+        - 1./(Re*Pr)*div(k*grad(T))
 
     return r_p, r_u, r_T

sapphire/simulations/examples/melt_gallium.py

@@ -17,6 +17,17 @@
 and their buoyancy equation for Gallium is $b(T) = Pr Ra T$.
 We chose before $t_r = \nu_l / x^2_r$ which sets $Re = 1$.
 The choice of $t_r = \rho_l c_l x_r^2 / k_l$ in \cite{belhamadia2019adaptive} sets $Re = 1/Pr$.
+
+\cite{belhamadia2019adaptive} uses the temperature scaling $\tilde{T} = (T - T_r) / \Delta T$.
+They chose reference temperature $T_r = 301.3 K$
+and set nondimensional melting temperature $T_f = 0.1525$.
+So their reference temperature was not chosen as the melting temperature.
+They set the initial temperature to $T_c = 0$ and hot wall to $T_h = 1$.
+This means that they chose $T_r = T_c$ and therefore $T_c = 301.3 K$.
+The dimensional melting temperature is $0.1525*(9.7 K) + 301.3 K$ => $T_f = 302.8 K$.
+
+We use the temperature scaling $\tilde{T} = (T - T_f)/ \Delta T$.
+Therefore, $\tilde{T}_f = 0.$ and $\tilde{T}_c = -0.1546$.
 """
 import firedrake as fe
 import sapphire.simulations.navier_stokes_boussinesq
@@ -25,21 +36,23 @@
 
 reference_length = 6.35  # cm
 
-reference_temperature = 301.3  # K
+coldwall_temperature = 301.3  # K
 
 reference_temperature_range = 9.7  # K
 
 reference_time = 292.90  # s
 
+
+
 class Simulation(sapphire.simulations.examples.melt_octadecane.Simulation):
 
     def __init__(self, *args,
             rayleigh_number = 7.e5,
             prandtl_number = 0.0216,
             stefan_number = 0.046,
-            liquidus_temperature = 0.1525,
+            liquidus_temperature = 0.,
             hotwall_temperature = 1.,
-            initial_temperature = 0.,
+            initial_temperature = -0.1546,
             cutoff_length = 0.5,
             element_degrees = (1, 2, 2),
             mesh_dimensions = (20, 40),

sapphire/simulations/navier_stokes_boussinesq.py

@@ -77,6 +77,8 @@ def momentum(self):
 
     def energy(self):
 
+        Re = self.reynolds_number
+
         Pr = self.prandtl_number
 
         _, u, T = fe.split(self.solution)
@@ -85,7 +87,7 @@ def energy(self):
 
         dx = fe.dx(degree = self.quadrature_degree)
 
-        return (psi_T*dot(u, grad(T)) + dot(grad(psi_T), 1./Pr*grad(T)))*dx
+        return (psi_T*dot(u, grad(T)) + dot(grad(psi_T), 1./(Re*Pr)*grad(T)))*dx
 
     def weak_form_residual(self):
 
@@ -152,7 +154,7 @@ def strong_residual(sim, solution):
 
     r_u = grad(u)*u + grad(p) - 2./Re*div(sym(grad(u))) + b
 
-    r_T = dot(u, grad(T)) - 1./Pr*div(grad(T))
+    r_T = dot(u, grad(T)) - 1./(Re*Pr)*div(grad(T))
 
     return r_p, r_u, r_T
 

tests/validation/test__melt_gallium.py

@@ -8,14 +8,14 @@ def test__melt_gallium__regression(tmpdir):
         mesh_dimensions = (20, 40),
         quadrature_degree = 4,
         time_stencil_size = 3,
-        timestep_size = 0.025,
+        timestep_size = 0.0125,
         liquidus_smoothing_factor = 0.05,
         solid_velocity_relaxation_factor = 1.e-10,
         output_directory_path = tmpdir)
 
-    sim.run(endtime = 0.15)
+    sim.run(endtime = 0.3)
 
     print("Liquid area = {}".format(sim.liquid_area))
 
-    assert(round(sim.liquid_area, 2) == 0.43)
+    assert(round(sim.liquid_area, 2) == 0.17)
 

tests/verification/test__unsteady_navier_stokes_boussinesq.py

@@ -91,7 +91,7 @@ def test__verify_second_order_spatial_convergence_via_mms():
         sim_kwargs = sim_kwargs,
         manufactured_solution = space_verification_solution,
         dirichlet_boundary_conditions = dirichlet_boundary_conditions,
-        meshes = [fe.UnitSquareMesh(n, n) for n in (4, 8, 16, 32, 64, 128)],
+        meshes = [fe.UnitSquareMesh(n, n) for n in (8, 16, 32)],
         norms = ("L2", "H1", "H1"),
         expected_orders = (2, 2, 2),
         decimal_places = 1,