Add Surface Tension (#387)

Co-authored-by: wilfonba <bwilfong3@gatech.edu> Co-authored-by: Benjamin Wilfong <bwilfong@dt-login01.delta.ncsa.illinois.edu> Co-authored-by: Henry LE BERRE <hberre3@gatech.edu>
MFlowCode · May 30, 2024 · 5b818eb · 5b818eb
1 parent 33984ea
commit 5b818eb
Show file tree

Hide file tree

Showing 35 changed files with 2,280 additions and 334 deletions.
diff --git a/.gitignore b/.gitignore
@@ -52,6 +52,8 @@ examples/*/*.err
 examples/*/viz/
 examples/*.jpg
 examples/*.png
+examples/*/workloads/
+workloads/
 
 benchmarks/*batch/*/
 benchmarks/*/D/*

diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -321,8 +321,9 @@ Additional details on this specification can be found in [The Naca Airfoil Serie
 | `cv`   ** | Real   | Sffened-gas parameter $c_v$ of fluid.          |
 | `qv`   ** | Real   | Stiffened-gas parameter $q$ of fluid.          |
 | `qvp`  ** | Real   | Stiffened-gas parameter $q'$ of fluid.         |
+| `sigma`   | Real   | Surface tension coefficient                    |
 
-Fluid material's parameters. All parameters should be prepended with `fluid_pp(i)` where $i$ is the fluid index.
+Fluid material's parameters. All parameters except for sigma should be prepended with `fluid_pp(i)` where $i$ is the fluid index.
 
 *: Parameters that work only with `model_eqns`=2.
 

diff --git a/examples/2D_laplace_pressure_jump/case.py b/examples/2D_laplace_pressure_jump/case.py
@@ -0,0 +1,144 @@
+#!/usr/bin/env python3
+
+# This case file demonstrates the Laplace pressure jump of a water droplet in air. The laplace pressure jump
+# in 2D is given by delta = sigma / r where delta is the pressure jump, sigma is the surface tension coefficient,
+# and r is the radius of the droplet. The results of this simulation agree with theory to well within 1%
+# relative error.
+
+import math
+import json
+
+l = 0.375
+
+# Numerical setup
+r0 = 0.15
+x0 = 0
+x1 = l
+y0 = 0
+y1 = l
+
+Nx = 49
+Ny = 49
+
+eps = 1e-9
+
+mydt = 5e-6
+
+#Configuration case dictionary
+data = {
+    # Logistics =============================
+        #'case_dir'          : '\'.\'',
+        'run_time_info'     : 'T',
+    # =======================================
+
+    # Computational Domain ==================
+        'x_domain%beg'      : x0,
+        'x_domain%end'      : x1,
+        'y_domain%beg'      : y0,
+        'y_domain%end'      : y1,
+        'm'                 : Nx,
+        'n'                 : Ny,
+        'p'                 : 0,
+        'cyl_coord'        : 'F',
+        'dt'                : mydt,
+        't_step_start'      : 0,
+        't_step_stop'       : 100000,
+        't_step_save'       : 1000,
+        #'t_step_stop'       : 100,
+        #'t_step_save'       : 100,
+    # =======================================
+
+    # Simulation Algorithm ==================
+        'model_eqns'        : 3,
+        'alt_soundspeed'    : 'F',
+        'adv_alphan'        : 'T',
+        'mixture_err'       : 'T',
+        'mpp_lim'           : 'F',
+        'time_stepper'      : 3,
+        #'recon_type'        : 1,
+        #'muscl_order'       : 2,
+        #'muscl_lim'         : 2,
+        'weno_order'        : 5,
+        'avg_state'         : 2,
+        'weno_eps'          : 1e-16,
+        'mapped_weno'       : 'T',
+        'null_weights'      : 'F',
+        'mp_weno'           : 'T',
+        'weno_Re_flux'      : 'F',
+        'riemann_solver'    : 2,
+        'wave_speeds'       : 1,
+        'bc_x%beg'          : -2,
+        'bc_x%end'          : -3,
+        'bc_y%beg'          : -2,
+        'bc_y%end'          : -3,
+        'num_patches'       : 2,
+        'num_fluids'        : 2,
+        'weno_avg'          : 'T',
+    # =======================================
+
+    # Database Structure Parameters =========
+        'format'            : 1,
+        'precision'         : 2,
+        'prim_vars_wrt'     : 'T',
+        'cons_vars_wrt'     : 'T',
+        'cf_wrt'            : 'T',
+        'parallel_io'       : 'T',
+    # =======================================
+
+        'sigma'             : 8,
+        #'flux_lim'          : 2,
+        #'flux_wrt(1)'       : 'T',
+        #'flux_wrt(2)'       : 'T',
+        #'cf_grad_wrt'       : 'T',
+
+    # Fluid Parameters (Water) ==============
+        'fluid_pp(1)%gamma'            : 1.E+00/(2.1E+00-1.E+00),
+        'fluid_pp(1)%pi_inf'           : 2.1E+00*1.E+06/(2.1E+00-1.E+00),
+        #'fluid_pp(1)%Re(1)'             : 1.e3,
+    # =======================================
+
+    # Fluid Parameters (Gas) ================
+        'fluid_pp(2)%gamma'            : 1.E+00/(1.4E+00-1.E+00),
+        'fluid_pp(2)%pi_inf'           : 0.E+00,
+        #'fluid_pp(2)%Re(1)'             : 1.81e5,
+    # =======================================
+
+    # Air Patch ==========================
+        'patch_icpp(1)%geometry'    : 3,
+        'patch_icpp(1)%x_centroid'  : 0,
+        'patch_icpp(1)%y_centroid'  : 0,
+        'patch_icpp(1)%length_x'    : 2,
+        'patch_icpp(1)%length_y'    : 2,
+        'patch_icpp(1)%vel(1)'      : 0.0,
+        'patch_icpp(1)%vel(2)'      : 0.0,
+        'patch_icpp(1)%vel(3)'      : 0.0,
+        'patch_icpp(1)%pres'        : 100000,
+        'patch_icpp(1)%alpha_rho(1)': eps*1000,
+        'patch_icpp(1)%alpha_rho(2)': (1-eps)*1,
+        'patch_icpp(1)%alpha(1)'    : eps,
+        'patch_icpp(1)%alpha(2)'    : 1-eps,
+        'patch_icpp(1)%cf_val'      : 0,
+    # ======================================
+
+    # Water Patch ========================
+        'patch_icpp(2)%alter_patch(1)' : 'T',
+        'patch_icpp(2)%smoothen'    : 'T',
+        'patch_icpp(2)%smooth_patch_id': 1,
+        'patch_icpp(2)%smooth_coeff': 0.95,
+        'patch_icpp(2)%geometry'    : 2,
+        'patch_icpp(2)%x_centroid'  : 0,
+        'patch_icpp(2)%y_centroid'  : 0,
+        'patch_icpp(2)%radius'      : r0,
+        'patch_icpp(2)%vel(1)'      : 0.,
+        'patch_icpp(2)%vel(2)'      : 0.,
+        'patch_icpp(2)%pres'        : 100000,
+        'patch_icpp(2)%alpha_rho(1)': (1-eps)*1000,
+        'patch_icpp(2)%alpha_rho(2)': eps*1,
+        'patch_icpp(2)%alpha(1)'    : (1-eps),
+        'patch_icpp(2)%alpha(2)'    : eps,
+        'patch_icpp(2)%cf_val'      : 1,
+    # ======================================
+
+}
+
+print(json.dumps(data))
diff --git a/examples/3D_recovering_sphere/case.py b/examples/3D_recovering_sphere/case.py
@@ -0,0 +1,146 @@
+#!/usr/bin/env python3
+
+# This simulation shows the early stages of a cubic droplet recovering a spherical shape due to capillary
+# forces. While the relaxation is not complete, it demonstrates the expecteed symmetric behavior.
+
+import math
+import json
+
+l = 0.375
+
+# Numerical setup
+r0 = 0.15
+x0 = 0
+x1 = l
+y0 = 0
+y1 = l
+z0 = 0
+z1 = l
+
+
+Nx = 199
+Ny = 199
+Nz = 199
+
+eps = 1e-9
+
+mydt = 1e-6
+
+#Configuration case dictionary
+data = {
+    # Logistics =============================
+        #'case_dir'          : '\'.\'',
+        'run_time_info'     : 'T',
+    # =======================================
+
+    # Computational Domain ==================
+        'x_domain%beg'      : x0,
+        'x_domain%end'      : x1,
+        'y_domain%beg'      : y0,
+        'y_domain%end'      : y1,
+        'z_domain%beg'      : z0,
+        'z_domain%end'      : z1,
+        'm'                 : Nx,
+        'n'                 : Ny,
+        'p'                 : Nz,
+        'cyl_coord'        : 'F',
+        'dt'                : mydt,
+        't_step_start'      : 400000,
+        't_step_stop'       : 1000000,
+        't_step_save'       : 2000,
+        #'t_step_stop'       : 100,
+        #'t_step_save'       : 100,
+    # =======================================
+
+    # Simulation Algorithm ==================
+        'model_eqns'        : 3,
+        'alt_soundspeed'    : 'F',
+        'adv_alphan'        : 'T',
+        'mixture_err'       : 'T',
+        'mpp_lim'           : 'F',
+        'time_stepper'      : 3,
+        'weno_order'        : 3,
+        'avg_state'         : 2,
+        'weno_eps'          : 1e-16,
+        'mapped_weno'       : 'T',
+        'null_weights'      : 'F',
+        'mp_weno'           : 'F',
+        'weno_Re_flux'      : 'F',
+        'riemann_solver'    : 2,
+        'wave_speeds'       : 1,
+        'bc_x%beg'          : -2,
+        'bc_x%end'          : -3,
+        'bc_y%beg'          : -2,
+        'bc_y%end'          : -3,
+        'bc_z%beg'          : -2,
+        'bc_z%end'          : -3,
+        'num_patches'       : 2,
+        'num_fluids'        : 2,
+        'weno_avg'          : 'T',
+    # =======================================
+
+    # Database Structure Parameters =========
+        'format'            : 1,
+        'precision'         : 2,
+        'alpha_wrt(1)'      : 'T',
+        # 'prim_vars_wrt'     : 'T',
+        'cf_wrt'            : 'T',
+        'parallel_io'       : 'T',
+    # =======================================
+
+        'sigma'             : 8.0,
+
+    # Fluid Parameters (Water) ==============
+        'fluid_pp(1)%gamma'            : 1.E+00/(2.1E+00-1.E+00),
+        'fluid_pp(1)%pi_inf'           : 2.1E+00*1.E+06/(2.1E+00-1.E+00),
+        #'fluid_pp(1)%Re(1)'             : 1.e3,
+    # =======================================
+
+    # Fluid Parameters (Gas) ================
+        'fluid_pp(2)%gamma'            : 1.E+00/(1.4E+00-1.E+00),
+        'fluid_pp(2)%pi_inf'           : 0.E+00,
+        #'fluid_pp(2)%Re(1)'             : 1.81e5,
+    # =======================================
+
+    # Air Patch ==========================
+        'patch_icpp(1)%geometry'    : 9,
+        'patch_icpp(1)%x_centroid'  : 0,
+        'patch_icpp(1)%y_centroid'  : 0,
+        'patch_icpp(1)%z_centroid'  : 0,
+        'patch_icpp(1)%length_x'    : 2,
+        'patch_icpp(1)%length_y'    : 2,
+        'patch_icpp(1)%length_z'    : 2,
+        'patch_icpp(1)%vel(1)'      : 0.0,
+        'patch_icpp(1)%vel(2)'      : 0.0,
+        'patch_icpp(1)%vel(3)'      : 0.0,
+        'patch_icpp(1)%pres'        : 100000,
+        'patch_icpp(1)%alpha_rho(1)': eps*1000,
+        'patch_icpp(1)%alpha_rho(2)': (1-eps)*1,
+        'patch_icpp(1)%alpha(1)'    : eps,
+        'patch_icpp(1)%alpha(2)'    : 1-eps,
+        'patch_icpp(1)%cf_val'      : 0,
+    # ======================================
+
+    # Water Patch ========================
+        'patch_icpp(2)%alter_patch(1)' : 'T',
+        'patch_icpp(2)%geometry'    : 9,
+        'patch_icpp(2)%x_centroid'  : 0,
+        'patch_icpp(2)%y_centroid'  : 0,
+        'patch_icpp(2)%z_centroid'  : 0,
+        'patch_icpp(2)%length_x'    : r0,
+        'patch_icpp(2)%length_y'    : r0,
+        'patch_icpp(2)%length_z'    : r0,
+        'patch_icpp(2)%vel(1)'      : 0.,
+        'patch_icpp(2)%vel(2)'      : 0.,
+        'patch_icpp(2)%vel(3)'      : 0.,
+        'patch_icpp(2)%pres'        : 100000,
+        'patch_icpp(2)%alpha_rho(1)': (1-eps)*1000,
+        'patch_icpp(2)%alpha_rho(2)': eps*1,
+        'patch_icpp(2)%alpha(1)'    : (1-eps),
+        'patch_icpp(2)%alpha(2)'    : eps,
+        'patch_icpp(2)%cf_val'      : 1,
+    # ======================================
+
+}
+
+print(json.dumps(data))
diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp
@@ -22,5 +22,6 @@ module m_constants
     integer, parameter :: num_patches_max = 10
     integer, parameter :: pathlen_max = 400
     integer, parameter :: nnode = 4    !< Number of QBMM nodes
+    real(kind(0d0)), parameter :: capillary_cutoff = 1e-6 !< color function gradient magnitude at which to apply the surface tension fluxes
 
 end module m_constants
diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
@@ -191,6 +191,8 @@ module m_derived_types
         integer :: hcid
         !! id for hard coded initial condition
 
+        real(kind(0d0)) :: cf_val !! color function value
+
     end type ic_patch_parameters
 
     type ib_patch_parameters

diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
@@ -1014,6 +1014,10 @@ contains
                         qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)
                     end do
 
+                    if (sigma /= dflt_real) then
+                        qK_prim_vf(c_idx)%sf(j, k, l) = qK_cons_vf(c_idx)%sf(j, k, l)
+                    end if
+
                 end do
             end do
         end do
@@ -1163,6 +1167,11 @@ contains
                             end if
                         end do
                     end if
+
+                    if (sigma /= dflt_real) then
+                        q_cons_vf(c_idx)%sf(j, k, l) = q_prim_vf(c_idx)%sf(j, k, l)
+                    end if
+
                 end do
             end do
         end do

diff --git a/src/post_process/m_checker.f90 b/src/post_process/m_checker.f90
@@ -473,6 +473,18 @@ subroutine s_check_inputs()
             end if
         end if
 
+        ! Constraints on the surface tension model
+        if (sigma /= dflt_real .and. sigma < 0d0) then
+            call s_mpi_abort('The surface tension coefficient must be'// &
+                             'greater than or equal to zero. Exiting ...')
+        elseif (sigma /= dflt_real .and. model_eqns /= 3) then
+            call s_mpi_abort("The surface tension model requires"// &
+                             'model_eqns=3. Exiting ...')
+        elseif (sigma == dflt_real .and. cf_wrt) then
+            call s_mpi_abort('cf_wrt can only be anabled if the surface'// &
+                             'coefficient is set')
+        end if
+
     end subroutine s_check_inputs
 
 end module m_checker