Change VOF density flux form in mdot to cancel out buoyancy (#1192)

* Change density form in mdot to cancel out buoyancy * Fix LHS and re-bless * Format text --------- Co-authored-by: whorne <[email protected]> Co-authored-by: psakievich <[email protected]>
Exawind · Jun 21, 2023 · 923bc26 · 923bc26
1 parent cbcb4d7
commit 923bc26
Show file tree

Hide file tree

Showing 4 changed files with 16 additions and 20 deletions.
diff --git a/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold b/reg_tests/test_files/VOFDroplet/VOFDroplet.norm.gold
@@ -1,3 +1,3 @@
-0.0001417205047386008 11 0.0001
-4.686504088256447e-05 12 0.0002
-1.667681063887448e-05 13 0.0003
+2.67343413675555e-05 11 0.0001
+6.859619665424383e-06 12 0.0002
+3.313464463487743e-06 13 0.0003
diff --git a/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml b/reg_tests/test_files/VOFDroplet/VOFDroplet.yaml
@@ -29,7 +29,7 @@ linear_solvers:
 realms:
 
   - name: realm_1
-    mesh: "generated:150x150x4|bbox:-0.5,-0.5,-0.5,0.5,0.5,0.5|sideset:xXyYzZ|show" 
+    mesh: "generated:75x75x75|bbox:-0.5,-0.5,-0.5,0.5,0.5,0.5|sideset:xXyYzZ|show" 
     use_edges: yes
     automatic_decomposition_type: rcb
 

diff --git a/src/edge_kernels/ContinuityEdgeSolverAlg.C b/src/edge_kernels/ContinuityEdgeSolverAlg.C
@@ -88,22 +88,15 @@ ContinuityEdgeSolverAlg::execute()
       const DblType pressureL = pressure.get(nodeL, 0);
       const DblType pressureR = pressure.get(nodeR, 0);
 
-      const DblType densityL =
-        om_solveIncompressibleEqn * density.get(nodeL, 0) +
-        solveIncompressibleEqn;
-      const DblType densityR =
-        om_solveIncompressibleEqn * density.get(nodeR, 0) +
-        solveIncompressibleEqn;
-
-      const DblType udiagL =
-        udiag.get(nodeL, 0) * (om_solveIncompressibleEqn +
-                               solveIncompressibleEqn * density.get(nodeL, 0));
-      const DblType udiagR =
-        udiag.get(nodeR, 0) * (om_solveIncompressibleEqn +
-                               solveIncompressibleEqn * density.get(nodeR, 0));
+      const DblType densityL = density.get(nodeL, 0);
+      const DblType densityR = density.get(nodeR, 0);
 
+      const DblType udiagL = udiag.get(nodeL, 0);
+      const DblType udiagR = udiag.get(nodeR, 0);
       const DblType projTimeScale = 0.5 * (1.0 / udiagL + 1.0 / udiagR);
       const DblType rhoIp = 0.5 * (densityL + densityR);
+      const DblType denScale =
+        (1.0 / rhoIp) * solveIncompressibleEqn + om_solveIncompressibleEqn;
 
       DblType axdx = 0.0;
       DblType asq = 0.0;
@@ -137,7 +130,9 @@ ContinuityEdgeSolverAlg::execute()
           kxj * GjIp * nocFac;
       }
       tmdot /= tauScale;
-      const DblType lhsfac = -asq * inv_axdx * projTimeScale / tauScale;
+      tmdot *= denScale;
+      const DblType lhsfac =
+        -asq * inv_axdx * projTimeScale * denScale / tauScale;
 
       // Left node entries
       smdata.lhs(0, 0) = -lhsfac;

diff --git a/src/user_functions/DropletVOFAuxFunction.C b/src/user_functions/DropletVOFAuxFunction.C
@@ -38,12 +38,13 @@ DropletVOFAuxFunction::do_evaluate(
 
     const double x = coords[0];
     const double y = coords[1];
-    const double interface_thickness = 0.005;
+    const double z = coords[2];
+    const double interface_thickness = 0.025;
 
     fieldPtr[0] = 0.0;
     fieldPtr[0] += -0.5 * (std::erf(y / interface_thickness) + 1.0) + 1.0;
 
-    auto radius = std::sqrt(x * x + (y - 0.15) * (y - 0.15)) - 0.075;
+    auto radius = std::sqrt(x * x + (y - 0.25) * (y - 0.25) + z * z) - 0.075;
     fieldPtr[0] += -0.5 * (std::erf(radius / interface_thickness) + 1.0) + 1.0;
 
     fieldPtr += fieldSize;