added options and drag models

wilfonba · wilfonba · commit 8ffc20cf0e28 · 2025-10-07T11:32:44.000-04:00
diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp
@@ -22,7 +22,7 @@ module m_constants
     integer, parameter :: fourier_rings = 5                       !< Fourier filter ring limit
     integer, parameter :: num_fluids_max = 10                     !< Maximum number of fluids in the simulation
     integer, parameter :: num_probes_max = 10                     !< Maximum number of flow probes in the simulation
-    integer, parameter :: num_patches_max = 10
+    integer, parameter :: num_patches_max = 20
     integer, parameter :: num_bc_patches_max = 10
     integer, parameter :: pathlen_max = 400
     integer, parameter :: nnode = 4    !< Number of QBMM nodes
diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
@@ -434,6 +434,9 @@ module m_derived_types
         integer :: nBubs_glb                !< Global number of bubbles
         integer :: vel_model                !< Particle velocity model
         integer :: drag_model               !< Particle drag model
+        logical :: pressure_force           !< Include pressure force translational motion
+        logical :: gravity_force            !< Include gravity force in translational motion
+        logical :: momentum_transfer_force  !< Include momentum transfer from radial dynamics in translational motion
         real(wp) :: c_d                     !< Drag coefficient
         real(wp) :: epsilonb         !< Standard deviation scaling for the gaussian function
         real(wp) :: charwidth        !< Domain virtual depth (z direction, for 2D simulations)
diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp
@@ -479,22 +479,19 @@ contains
         !!  @param var_loc holds the local value to be reduced among
         !!      all the processors in communicator. On output, the variable holds
         !!      the sum, reduced amongst all of the local values.
-    subroutine s_mpi_reduce_int_sum(var_loc)
+    subroutine s_mpi_reduce_int_sum(var_loc, sum)
 
         integer, intent(inout) :: var_loc
+        integer, intent(inout) :: sum
 
 #ifdef MFC_MPI
         integer :: ierr !< Generic flag used to identify and report MPI errors
 
-        ! Temporary storage variable that holds the reduced sum value
-        integer :: var_glb
-
         ! Performing reduction procedure and eventually storing its result
         ! into the variable that was initially inputted into the subroutine
-        call MPI_REDUCE(var_loc, var_glb, 1, MPI_INTEGER, &
+        call MPI_REDUCE(var_loc, sum, 1, MPI_INTEGER, &
                         MPI_SUM, 0, MPI_COMM_WORLD, ierr)
 
-        var_loc = var_glb
 #endif
 
     end subroutine s_mpi_reduce_int_sum
diff --git a/src/simulation/m_body_forces.fpp b/src/simulation/m_body_forces.fpp
@@ -56,16 +56,13 @@ contains
     subroutine s_compute_acceleration(t)
 
         real(wp), intent(in) :: t
+        accel_bf(:) = 0._wp
 
-        if (m > 0) then
-            accel_bf(1) = g_x + k_x*sin(w_x*t - p_x)
-            if (n > 0) then
-                accel_bf(2) = g_y + k_y*sin(w_y*t - p_y)
-                if (p > 0) then
-                    accel_bf(3) = g_z + k_z*sin(w_z*t - p_z)
-                end if
+        #:for DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')]
+            if (bf_${XYZ}$) then
+                accel_bf(${DIR}$) = g_${XYZ}$ + k_${XYZ}$*sin(w_${XYZ}$*t - p_${XYZ}$)
             end if
-        end if
+        #:endfor
 
         $:GPU_UPDATE(device='[accel_bf]')
 
diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp
@@ -450,7 +450,7 @@ contains
         !!  @param fRho Current density
         !!  @param fP Current driving pressure
         !!  @param fR Current bubble radius
-        !!  @param fV Current bubble velocity
+        !!  @param fV Current bubble radial velocity
         !!  @param fR0 Equilibrium bubble radius
         !!  @param fpb Internal bubble pressure
         !!  @param fpbdot Time-derivative of internal bubble pressure
@@ -490,7 +490,7 @@ contains
         real(wp) :: h !< Time step size
         real(wp), dimension(4) :: myR_tmp1, myV_tmp1, myR_tmp2, myV_tmp2 !< Bubble radius, radial velocity, and radial acceleration for the inner loop
         real(wp), dimension(4) :: myPb_tmp1, myMv_tmp1, myPb_tmp2, myMv_tmp2 !< Gas pressure and vapor mass for the inner loop (EL)
-        real(wp) :: fR2, fV2, fpb2, fmass_v2, vTemp, aTemp
+        real(wp) :: fR2, fV2, fpb2, fmass_v2, vTemp, aTemp, f_bTemp
         integer :: l, iter_count
 
         call s_initial_substep_h(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
@@ -571,10 +571,19 @@ contains
                             end do
                         elseif (lag_params%vel_model == 2) then
                             do l = 1, num_dims
-                                aTemp = f_get_acceleration(fPos(l), fR, fVel(l), fmass_n, fmass_v, &
+                                f_bTemp = f_get_bubble_force(fPos(l), fR, fV, fVel, fmass_n, fmass_v, &
+                                                           fRe, fRho, cell, l, q_prim_vf)
+                                aTemp = f_bTemp/(fmass_n + fmass_v)
+                                fPos(l) = fPos(l) + h * fVel(l)
+                                fVel(l) = fVel(l) + h * aTemp
+                            end do
+                        elseif (lag_params%vel_model == 3) then
+                            do l = 1, num_dims
+                                f_bTemp = f_get_bubble_force(fPos(l), fR, fV, fVel, fmass_n, fmass_v, &
                                                            fRe, fRho, cell, l, q_prim_vf)
-                                fPos(l) = fPos(l) + h*fVel(l)
-                                fVel(l) = fVel(l) + h*aTemp
+                                aTemp = 2._wp * f_bTemp / (fmass_n + fmass_v) - 3 * fV * fVel(l) / fR
+                                fPos(l) = fPos(l) + h * fVel(l)
+                                fVel(l) = fVel(l) + h * aTemp
                             end do
                         end if
                     end if
diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
@@ -279,7 +279,7 @@ contains
 
         print *, " Lagrange bubbles running, in proc", proc_rank, "number:", bub_id, "/", id
 
-        call s_mpi_reduce_int_sum(bub_id)
+        call s_mpi_reduce_int_sum(bub_id, bub_id)
 
         if (proc_rank == 0) then
             if (bub_id == 0) call s_mpi_abort('No bubbles in the domain. Check input/lag_bubbles.dat')
@@ -621,7 +621,7 @@ contains
         real(wp) :: myR, myV, myBeta_c, myBeta_t, myR0, myPbdot, myMvdot
         real(wp) :: myPinf, aux1, aux2, myCson, myRho
         real(wp), dimension(3) :: myPos, myVel
-        real(wp) :: gamma, pi_inf, qv
+        real(wp) :: gamma, pi_inf, qv, f_b
         real(wp), dimension(contxe) :: myalpha_rho, myalpha
         real(wp), dimension(2) :: Re
         integer, dimension(3) :: cell
@@ -734,10 +734,18 @@ contains
                         mtn_dveldt(k, l, stage) = 0._wp
                     elseif (lag_params%vel_model == 2) then
                         mtn_dposdt(k, l, stage) = myVel(l)
-                        mtn_dveldt(k, l, stage) = f_get_acceleration(myPos(l), &
-                                                                     myR, myVel(l), &
-                                                                     myMass_n, myMass_v, &
-                                                                     Re(1), myRho, cell, l, q_prim_vf)
+                        f_b = f_get_bubble_force(myPos(l), &
+                                                 myR, myV, myVel, &
+                                                 myMass_n, myMass_v, &
+                                                 Re(1), myRho, cell, l, q_prim_vf)
+                        mtn_dveldt(k, l, stage) = f_b / (myMass_n + myMass_v)
+                    elseif (lag_params%vel_model == 3) then
+                        mtn_dposdt(k, l, stage) = myVel(l)
+                        f_b = f_get_bubble_force(myPos(l), &
+                                                 myR, myV, myVel, &
+                                                 myMass_n, myMass_v, &
+                                                 Re(1), myRho, cell, l, q_prim_vf)
+                        mtn_dveldt(k, l, stage) = 2._wp * f_b / (myMass_n + myMass_v) - 3._wp * myV * myVel(l) / myR
                     else
                         mtn_dposdt(k, l, stage) = 0._wp
                         mtn_dveldt(k, l, stage) = 0._wp
@@ -776,6 +784,7 @@ contains
 
         if (lag_params%solver_approach == 2) then
 
+            ! (q / (1 - beta)) * d(beta)/dt source
             if (p == 0) then
                 $:GPU_PARALLEL_LOOP(collapse=4)
                 do k = 0, p
@@ -813,6 +822,7 @@ contains
 
                 call s_gradient_dir(q_prim_vf(E_idx), q_beta%vf(3), l)
 
+                ! (beta / (1 - beta)) * dP/dl source
                 $:GPU_PARALLEL_LOOP(collapse=3)
                 do k = 0, p
                     do j = 0, n
@@ -839,6 +849,7 @@ contains
 
                 call s_gradient_dir(q_beta%vf(3), q_beta%vf(4), l)
 
+                ! (beta / (1 - beta)) * d(Pu)/dl source
                 $:GPU_PARALLEL_LOOP(collapse=3)
                 do k = 0, p
                     do j = 0, n
@@ -852,7 +863,6 @@ contains
                     end do
                 end do
             end do
-
         end if
 
     end subroutine s_compute_bubbles_EL_source
@@ -1541,6 +1551,13 @@ contains
             nBubs = nBubs + newBubs
         end if
 
+        if (run_time_info) then
+            call s_mpi_reduce_int_sum(nBubs, active_bubs)
+            if (proc_rank == 0 .and. active_bubs == 0) then
+                call s_mpi_abort('No bubbles remain in the domain. Simulation ending.')
+            end if
+        end if
+
         ! Handle MPI transfer of bubbles going to another processor's local domain
         if (num_procs > 1) then
             call nvtxStartRange("LAG-BC-TRANSFER-LIST")
@@ -1745,36 +1762,34 @@ contains
                     end do
                 end do
             end do
-        else
-            if (dir == 2) then
-                ! Gradient in y dir.
-                $:GPU_PARALLEL_LOOP(collapse=3)
-                do k = 0, p
-                    do j = 0, n
-                        do i = 0, m
-                            dq%sf(i, j, k) = q%sf(i, j, k)*(dy(j + 1) - dy(j - 1)) &
-                                             + q%sf(i, j + 1, k)*(dy(j) + dy(j - 1)) &
-                                             - q%sf(i, j - 1, k)*(dy(j) + dy(j + 1))
-                            dq%sf(i, j, k) = dq%sf(i, j, k)/ &
-                                             ((dy(j) + dy(j - 1))*(dy(j) + dy(j + 1)))
-                        end do
+        elseif (dir == 2) then
+            ! Gradient in y dir.
+            $:GPU_PARALLEL_LOOP(collapse=3)
+            do k = 0, p
+                do j = 0, n
+                    do i = 0, m
+                        dq%sf(i, j, k) = q%sf(i, j, k)*(dy(j + 1) - dy(j - 1)) &
+                                         + q%sf(i, j + 1, k)*(dy(j) + dy(j - 1)) &
+                                         - q%sf(i, j - 1, k)*(dy(j) + dy(j + 1))
+                        dq%sf(i, j, k) = dq%sf(i, j, k)/ &
+                                         ((dy(j) + dy(j - 1))*(dy(j) + dy(j + 1)))
                     end do
                 end do
-            else
-                ! Gradient in z dir.
-                $:GPU_PARALLEL_LOOP(collapse=3)
-                do k = 0, p
-                    do j = 0, n
-                        do i = 0, m
-                            dq%sf(i, j, k) = q%sf(i, j, k)*(dz(k + 1) - dz(k - 1)) &
-                                             + q%sf(i, j, k + 1)*(dz(k) + dz(k - 1)) &
-                                             - q%sf(i, j, k - 1)*(dz(k) + dz(k + 1))
-                            dq%sf(i, j, k) = dq%sf(i, j, k)/ &
-                                             ((dz(k) + dz(k - 1))*(dz(k) + dz(k + 1)))
-                        end do
+            end do
+        elseif (dir == 3) then
+            ! Gradient in z dir.
+            $:GPU_PARALLEL_LOOP(collapse=3)
+            do k = 0, p
+                do j = 0, n
+                    do i = 0, m
+                        dq%sf(i, j, k) = q%sf(i, j, k)*(dz(k + 1) - dz(k - 1)) &
+                                         + q%sf(i, j, k + 1)*(dz(k) + dz(k - 1)) &
+                                         - q%sf(i, j, k - 1)*(dz(k) + dz(k + 1))
+                        dq%sf(i, j, k) = dq%sf(i, j, k)/ &
+                                         ((dz(k) + dz(k - 1))*(dz(k) + dz(k + 1)))
                     end do
                 end do
-            end if
+            end do
         end if
 
     end subroutine s_gradient_dir
diff --git a/src/simulation/m_bubbles_EL_kernels.fpp b/src/simulation/m_bubbles_EL_kernels.fpp
@@ -523,10 +523,11 @@ contains
 
     end function f_interpolate_velocity
 
-    !! This function calculates the acceleration of the bubble
+    !! This function calculates the force on a bubble
             !!      based on the pressure gradient, velocity, and drag model.
             !! @param pos Position of the bubble in direction i
             !! @param rad Radius of the bubble
+            !! @param rdot Radial velocity of the bubble
             !! @param vel Velocity of the bubble
             !! @param mg Mass of the gas in the bubble
             !! @param mv Mass of the liquid in the bubble
@@ -536,16 +537,18 @@ contains
             !! @param i Direction of the velocity (1: x, 2: y, 3: z)
             !! @param q_prim_vf Eulerian field with primitive variables
             !! @return a Acceleration of the bubble in direction i
-    pure function f_get_acceleration(pos, rad, vel, mg, mv, Re, rho, cell, i, q_prim_vf) result(a)
+    pure function f_get_bubble_force(pos, rad, rdot, vel, mg, mv, Re, rho, cell, i, q_prim_vf) result(force)
         $:GPU_ROUTINE(parallelism='[seq]')
-        real(wp), intent(in) :: pos, rad, vel, mg, mv, Re, rho
+        real(wp), intent(in) :: pos, rad, rdot, mg, mv, Re, rho
+        real(wp), intent(in), dimension(3) :: vel
         integer, dimension(3), intent(in) :: cell
         integer, intent(in) :: i
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
 
-        real(wp) :: a
-        real(wp) :: dp, vol, force, v_rel
+        real(wp) :: a, dp, vol, force, Re_b, C_d, v_rel_mag, rho_b
+        real(wp), dimension(3) :: v_rel
         real(wp), dimension(fd_order - 1) :: xi, eta, L
+        integer :: j
 
         if (fd_order == 2) then
             if (i == 1) then
@@ -610,19 +613,43 @@ contains
             dp = L(1)*eta(1) + L(2)*eta(2) + L(3)*eta(3)
         end if
 
-        vol = (4._wp/3._wp)*pi*rad**3._wp
-        force = -1._wp*vol*dp
+        vol = (4._wp/3._wp) * pi * (rad**3._wp)
+        v_rel_mag = 0._wp
 
-        v_rel = vel - f_interpolate_velocity(pos, cell, i, q_prim_vf)
+        do j = 1, num_dims
+            v_rel(j) = vel(j) - f_interpolate_velocity(pos, cell, j, q_prim_vf)
+            v_rel_mag = v_rel_mag + v_rel(j)**2._wp
+        end do
+
+        force = 0._wp
 
         if (lag_params%drag_model == 1) then ! Free slip Stokes drag
-            force = force - (4._wp*pi*rad*v_rel)/Re
+            force = force - (4._wp*pi*rad*v_rel(i))/Re
         else if (lag_params%drag_model == 2) then ! No slip Stokes drag
-            force = force - (6._wp*pi*rad*v_rel)/Re
+            force = force - (6._wp*pi*rad*v_rel(i))/Re
+        else if (lag_params%drag_model == 3) then ! Levich drag
+            force = force - (12._wp*pi*rad*v_rel(i))/Re
+        else if (lag_params%drag_model > 0) then ! Drag coefficient model
+            v_rel_mag = sqrt(v_rel_mag)
+            Re_b = max(1e-3, rho * v_rel_mag * 2._wp * rad * Re)
+            if (lag_params%drag_model == 4) then ! Mei et al. 1994
+                C_d = 16._wp / Re_b * (1._wp + (8._wp / Re_b + 0.5_wp * (1._wp + 3.315_wp * Re_b**(-0.5_wp))) ** -1._wp)
+            end if
+            force = force - 0.5_wp * C_d * rho * pi * rad**2._wp * v_rel(i) * v_rel_mag
+        end if
+
+        if (lag_params%pressure_force) then
+            force = force - vol * dp
+        end if
+
+        if (lag_params%gravity_force) then
+            force = force + (mg + mv) * accel_bf(i)
         end if
 
-        a = force/(mg + mv)
+        if (lag_params%momentum_transfer_force) then
+            force = force  - 4._wp * pi * rho * rad**2._wp * v_rel(i) * rdot
+        end if
 
-    end function f_get_acceleration
+    end function f_get_bubble_force
 
 end module m_bubbles_EL_kernels
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
@@ -816,6 +816,9 @@ contains
         lag_params%nBubs_glb = dflt_int
         lag_params%vel_model = dflt_int
         lag_params%drag_model = dflt_int
+        lag_params%pressure_force = .true.
+        lag_params%gravity_force = .false.
+        lag_params%momentum_transfer_force = .false.
         lag_params%c_d = dflt_real
         lag_params%epsilonb = 1._wp
         lag_params%charwidth = dflt_real
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
@@ -188,7 +188,8 @@ contains
 
         if (bubbles_lagrange) then
             #:for VAR in [ 'heatTransfer_model', 'massTransfer_model', 'pressure_corrector', &
-                & 'write_bubbles', 'write_bubbles_stats', 'vel_model', 'drag_model']
+                & 'write_bubbles', 'write_bubbles_stats', 'vel_model', 'drag_model', &
+                & 'pressure_force', 'gravity_force', 'momentum_transfer_force']
                 call MPI_BCAST(lag_params%${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
             #:endfor
 
diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py
@@ -163,7 +163,7 @@ def analytic(self):
 
     PRE_PROCESS[f"patch_bc({bc_p_id})%radius"] = ParamType.REAL
 
-for p_id in range(1, 10+1):
+for p_id in range(1, 20+1):
     for attribute, ty in [("geometry", ParamType.INT), ("smoothen", ParamType.LOG),
                       ("smooth_patch_id", ParamType.INT), ("hcid", ParamType.INT)]:
         PRE_PROCESS[f"patch_icpp({p_id})%{attribute}"] = ty
@@ -326,7 +326,8 @@ def analytic(self):
 })
 
 for var in [ 'heatTransfer_model', 'massTransfer_model', 'pressure_corrector',
-             'write_bubbles', 'write_bubbles_stats' ]:
+             'write_bubbles', 'write_bubbles_stats', 'pressure_force', 
+             'gravity_force', 'momentum_transfer_force']:
     SIMULATION[f'lag_params%{var}'] = ParamType.LOG
 
 for var in [ 'solver_approach', 'cluster_type', 'smooth_type', 'nBubs_glb',
