Invicid two-way fluid-structure interaction

examples/2D_mibm_cylinder_in_cross_flow/case.py

@@ -0,0 +1,105 @@
+import json
+import math
+
+Mu = 1.84e-05
+gam_a = 1.4
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            # For these computations, the cylinder is placed at the (0,0,0)
+            # domain origin.
+            # axial direction
+            "x_domain%beg": 0.0e00,
+            "x_domain%end": 6.0e-03,
+            # r direction
+            "y_domain%beg": 0.0e00,
+            "y_domain%end": 6.0e-03,
+            "cyl_coord": "F",
+            "m": 200,
+            "n": 200,
+            "p": 0,
+            "dt": 6.0e-6,
+            "t_step_start": 0,
+            "t_step_stop": 10000,  # 10000,
+            "t_step_save": 100,
+            # Simulation Algorithm Parameters
+            # Only one patches are necessary, the air tube
+            "num_patches": 1,
+            # Use the 5 equation model
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            # One fluids: air
+            "num_fluids": 1,
+            # time step
+            "mpp_lim": "F",
+            # Correct errors when computing speed of sound
+            "mixture_err": "T",
+            # Use TVD RK3 for time marching
+            "time_stepper": 3,
+            # Use WENO5
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "weno_Re_flux": "T",
+            "weno_avg": "T",
+            "avg_state": 2,
+            "mapped_weno": "T",
+            "null_weights": "F",
+            "mp_weno": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            # We use ghost-cell
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            "bc_y%beg": -3,
+            "bc_y%end": -3,
+            # Set IB to True and add 1 patch
+            "ib": "T",
+            "num_ibs": 1,
+            "viscous": "T",
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "E_wrt": "T",
+            "parallel_io": "T",
+            # Patch: Constant Tube filled with air
+            # Specify the cylindrical air tube grid geometry
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 3.0e-03,
+            # Uniform medium density, centroid is at the center of the domain
+            "patch_icpp(1)%y_centroid": 3.0e-03,
+            "patch_icpp(1)%length_x": 6.0e-03,
+            "patch_icpp(1)%length_y": 6.0e-03,
+            # Specify the patch primitive variables
+            "patch_icpp(1)%vel(1)": 0.05e00,
+            "patch_icpp(1)%vel(2)": 0.0e00,
+            "patch_icpp(1)%pres": 1.0e00,
+            "patch_icpp(1)%alpha_rho(1)": 1.0e00,
+            "patch_icpp(1)%alpha(1)": 1.0e00,
+            # Patch: Cylinder Immersed Boundary
+            "patch_ib(1)%geometry": 2,
+            "patch_ib(1)%x_centroid": 1.5e-03,
+            "patch_ib(1)%y_centroid": 4.5e-03,
+            "patch_ib(1)%radius": 0.3e-03,
+            "patch_ib(1)%slip": "F",
+            "patch_ib(1)%moving_ibm": 2,
+            "patch_ib(1)%vel(2)": -0.1,
+            "patch_ib(1)%angles(1)": 0.0,  # x-axis rotation in radians
+            "patch_ib(1)%angles(2)": 0.0,  # y-axis rotation
+            "patch_ib(1)%angles(3)": 0.0,  # z-axis rotation
+            "patch_ib(1)%angular_vel(1)": 0.0,  # x-axis rotational velocity in radians per second
+            "patch_ib(1)%angular_vel(2)": 0.0,  # y-axis rotation
+            "patch_ib(1)%angular_vel(3)": 0.0,  # z-axis rotation
+            "patch_ib(1)%mass": 0.0001,  # z-axis rotation
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50(Not 1.40)
+            "fluid_pp(1)%pi_inf": 0,
+            "fluid_pp(1)%Re(1)": 2500000,
+        }
+    )
+)

src/common/m_derived_types.fpp

@@ -332,6 +332,8 @@ module m_derived_types
 
         !! Patch conditions for moving imersed boundaries
         integer :: moving_ibm ! 0 for no moving, 1 for moving, 2 for moving on forced path
+        real(wp) :: mass, moment ! mass and moment of inertia of object used to compute forces in 2-way coupling
+        real(wp), dimension(1:3) :: force, torque ! vectors for the computed force and torque values applied to an IB
         real(wp), dimension(1:3) :: vel
         real(wp), dimension(1:3) :: step_vel ! velocity array used to store intermediate steps in the time_stepper module
         real(wp), dimension(1:3) :: angular_vel

src/common/m_mpi_common.fpp

@@ -472,6 +472,30 @@ contains
 
     end subroutine s_mpi_allreduce_sum
 
+    !>  This subroutine follows the behavior of the s_mpi_allreduce_sum subroutine
+    !>  with the additional feature that it reduces an array of vectors.
+    impure subroutine s_mpi_allreduce_vectors_sum(var_loc, var_glb, num_vectors, vector_length)
+
+        integer, intent(in) :: num_vectors, vector_length
+        real(wp), dimension(:, :), intent(in) :: var_loc
+        real(wp), dimension(:, :), intent(out) :: var_glb
+
+#ifdef MFC_MPI
+        integer :: ierr !< Generic flag used to identify and report MPI errors
+
+        ! Performing the reduction procedure
+        if (loc(var_loc) == loc(var_glb)) then
+            call MPI_Allreduce(MPI_IN_PLACE, var_glb, num_vectors*vector_length, &
+                              mpi_p, MPI_SUM, MPI_COMM_WORLD, ierr)
+        else
+            call MPI_Allreduce(var_loc, var_glb, num_vectors*vector_length, &
+                              mpi_p, MPI_SUM, MPI_COMM_WORLD, ierr)
+        end if
+
+#endif
+
+    end subroutine s_mpi_allreduce_vectors_sum
+
     !>  The following subroutine takes the input local variable
         !!      from all processors and reduces to the sum of all
         !!      values. The reduced variable is recorded back onto the

src/pre_process/m_global_parameters.fpp

@@ -560,6 +560,8 @@ contains
             patch_ib(i)%vel(:) = 0._wp
             patch_ib(i)%angles(:) = 0._wp
             patch_ib(i)%angular_vel(:) = 0._wp
+            patch_ib(i)%mass = dflt_real
+            patch_ib(i)%moment = dflt_real
 
             ! sets values of a rotation matrix which can be used when calculating rotations
             patch_ib(i)%rotation_matrix = 0._wp

src/simulation/m_ibm.fpp

@@ -101,7 +101,6 @@ contains
             end if
             call s_update_ib_rotation_matrix(i)
         end do
-
         $:GPU_ENTER_DATA(copyin='[patch_ib]')
 
         ! Allocating the patch identities bookkeeping variable
@@ -197,6 +196,21 @@ contains
         type(ghost_point) :: gp
         type(ghost_point) :: innerp
 
+        ! set the Moving IBM interior Pressure Values
+        do patch_id = 1, num_ibs
+            if (patch_ib(patch_id)%moving_ibm == 2) then
+                do j = 0, m
+                    do k = 0, n
+                        do l = 0, p
+                            if (ib_markers%sf(j,k,l) == patch_id) then
+                                q_prim_vf(E_idx)%sf(j, k, l) = 1._wp
+                            end if
+                        end do
+                    end do
+                end do
+            end if
+        end do
+
         if (num_gps > 0) then
             $:GPU_PARALLEL_LOOP(private='[i,physical_loc,dyn_pres,alpha_rho_IP, alpha_IP,pres_IP,vel_IP,vel_g,vel_norm_IP,r_IP, v_IP,pb_IP,mv_IP,nmom_IP,presb_IP,massv_IP,rho, gamma,pi_inf,Re_K,G_K,Gs,gp,innerp,norm,buf, radial_vector, rotation_velocity, j,k,l,q,qv_K,c_IP,nbub,patch_id]')
             do i = 1, num_gps
@@ -244,6 +258,17 @@ contains
                 if (surface_tension) then
                     q_prim_vf(c_idx)%sf(j, k, l) = c_IP
                 end if
+
+                ! set the pressure
+                do q = 1, num_fluids
+                    if (patch_ib(patch_id)%moving_ibm == 0) then
+                        q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
+                    else
+                      ! TODO :: improve for two fluid
+                      q_prim_vf(E_idx)%sf(j, k, l) = pres_IP / (1._wp - 2._wp*abs(levelset%sf(j,k,l,patch_id)*alpha_rho_IP(q)/pres_IP) * dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, levelset_norm%sf(j,k,l,patch_id,:)))
+                    end if
+                end do
+
                 if (model_eqns /= 4) then
                     ! If in simulation, use acc mixture subroutines
                     if (elasticity) then
@@ -969,6 +994,67 @@ contains
 
     end subroutine s_update_mib
 
+    ! compute the surface integrals of the IB via a volume integraion method described in
+    ! "A coupled IBM/Euler-Lagrange framework for simulating shock-induced particle size segregation"
+    ! by Archana Sridhar and Jesse Capecelatro
+    subroutine s_compute_ib_forces(pressure)
+
+        real(wp), dimension(0:m, 0:n, 0:p), intent(in) :: pressure
+
+        integer :: i, j, k, ib_idx
+        real(wp), dimension(num_ibs, 3) :: forces, torques
+        real(wp), dimension(1:3) :: pressure_divergence, radial_vector
+        real(wp) :: cell_volume
+
+        forces = 0._wp
+        torques = 0._wp
+
+        ! TODO :: This is currently only valid inviscid, and needs to be extended to add viscocity
+        $:GPU_PARALLEL_LOOP(private='[ib_idx,radial_vector,pressure_divergence,cell_volume]', copy='[forces,torques]', copyin='[ib_markers]', collapse=3)
+        do i = 0, m
+            do j = 0, n
+                do k = 0, p
+                    ib_idx = ib_markers%sf(i, j, k)
+                    if (ib_idx /= 0) then ! only need to compute the gradient for cells inside a IB
+                        if (patch_ib(ib_idx)%moving_ibm == 2) then
+                            if (num_dims == 3) then
+                                radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid] ! get the vector pointing to the grid cell
+                            else 
+                                radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, patch_ib(ib_idx)%y_centroid, 0._wp] ! get the vector pointing to the grid cell
+                            end if
+                            pressure_divergence(1) = (pressure(i+1, j, k) - pressure(i-1, j, k)) / (2._wp * x_cc(i))
+                            pressure_divergence(2) = (pressure(i, j+1, k) - pressure(i, j-1, k)) / (2._wp * y_cc(j))
+                            cell_volume = x_cc(i) * y_cc(j)
+                            if (num_dims == 3) then
+                                pressure_divergence(3) = (pressure(i, j, k+1) - pressure(i, j, k-1)) / (2._wp * z_cc(k))
+                                cell_volume = cell_volume * z_cc(k)
+                            else
+                                pressure_divergence(3) = 0._wp
+                            end if
+
+                            forces(ib_idx, :) = forces(ib_idx, :) - (pressure_divergence * cell_volume)
+                            torques(ib_idx, :) = torques(ib_idx, :) + (cross_product(radial_vector, pressure_divergence) * cell_volume)
+                        end if
+                    end if
+                end do
+            end do
+        end do
+        $:END_GPU_PARALLEL_LOOP()
+
+        ! reduce the forces across all MPI ranks
+        call s_mpi_allreduce_vectors_sum(forces, forces, num_ibs, 3)
+        call s_mpi_allreduce_vectors_sum(torques, torques, num_ibs, 3)
+
+        ! apply the summed forces
+        do i = 1, num_ibs
+            patch_ib(i)%force(:) = forces(i, :)
+            patch_ib(i)%torque(:) = matmul(patch_ib(i)%rotation_matrix_inverse, torques(i, :)) ! torques must be computed in the local coordinates of the IB
+        end do
+
+        print *, patch_ib(1)%force(1:2)
+
+    end subroutine s_compute_ib_forces
+
     !> Subroutine to deallocate memory reserved for the IBM module
     impure subroutine s_finalize_ibm_module()
 
@@ -978,6 +1064,78 @@ contains
 
     end subroutine s_finalize_ibm_module
 
+    subroutine s_compute_moment_of_inertia(ib_marker, axis)
+
+        real(wp), dimension(3), optional :: axis !< the axis about which we compute the moment. Only required in 3D.
+        integer, intent(in) :: ib_marker
+
+        real (wp) :: moment, distance_to_axis, cell_volume
+        real (wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis
+        integer :: i, j, k, count
+
+        if (p == 0) then
+            axis = [0, 1, 0]
+        else if (sqrt(sum(axis**2)) == 0) then
+            ! if the object is not actually rotating at this time, return a dummy value and exit
+            patch_ib(ib_marker)%moment = 1._wp
+            return
+        else
+            axis = axis / sqrt(sum(axis))
+        end if
+
+        ! if the IB is in 2D or a 3D sphere, we can compute this exactly
+        if (patch_ib(ib_marker)%geometry == 2) then ! circle
+            patch_ib(ib_marker)%moment = 0.5 * patch_ib(ib_marker)%mass * (patch_ib(ib_marker)%radius)**2
+        elseif (patch_ib(ib_marker)%geometry == 3) then ! rectangle
+            patch_ib(ib_marker)%moment = patch_ib(i)%mass * (patch_ib(ib_marker)%length_x**2 + patch_ib(ib_marker)%length_y**2) / 6._wp                
+        elseif (patch_ib(ib_marker)%geometry == 8) then ! sphere
+            patch_ib(ib_marker)%moment = 0.4 * patch_ib(ib_marker)%mass * (patch_ib(ib_marker)%radius)**2
+
+        else ! we do not have an analytic moment of inertia calculation and need to approximate it directly
+            count = 0
+            moment = 0._wp
+            cell_volume = (x_cc(1) - x_cc(0))*(y_cc(1) - y_cc(0)) ! computed without grid stretching. Update in the loop to perform with stretching
+            if (p /= 0) then
+              cell_volume = cell_volume * (z_cc(1) - z_cc(0))
+            end if
+
+
+            $:GPU_PARALLEL_LOOP(private='[position,closest_point_along_axis,vector_to_axis,distance_to_axis]', copy='[moment,count]', copyin='[ib_marker,cell_volume,axis]', collapse=3)
+            do i = 0, m
+                do j = 0, j
+                    do k = 0, p
+                        if (ib_markers%sf(i, j, k) == ib_marker) then
+                            $:GPU_ATOMIC(atomic='update')
+                            count = count + 1 ! increment the count of total cells in the boundary
+
+                            ! get the position in local coordinates so that the axis passes through 0, 0, 0
+                            if (p == 0) then
+                                position = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, 0._wp]
+                            else
+                                position = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid]
+                            end if
+
+                            ! project the position along the axis to find the closest distance to the rotation axis
+                            closest_point_along_axis = axis*dot_product(axis, position)
+                            vector_to_axis = position - closest_point_along_axis
+                            distance_to_axis = sum(vector_to_axis) ! saves the distance to the axis squared
+
+                            ! compute the position component of the moment
+                            $:GPU_ATOMIC(atomic='update')
+                            moment = moment + distance_to_axis
+                        end if
+                    end do
+                end do
+            end do
+            $:END_GPU_PARALLEL_LOOP()
+
+            ! write the final moment assuming the points are all uniform density
+            patch_ib(ib_marker)%moment = moment * patch_ib(ib_marker)%mass / (count * cell_volume)
+            $:GPU_UPDATE(device='[patch_ib(ib_marker)%moment]')
+        end if
+
+    end subroutine s_compute_moment_of_inertia
+
     function cross_product(a, b) result(c)
         implicit none
         real(wp), intent(in) :: a(3), b(3)

src/simulation/m_mpi_proxy.fpp

@@ -203,7 +203,7 @@ contains
 
         do i = 1, num_ibs
             #:for VAR in [ 'radius', 'length_x', 'length_y', 'length_z', &
-                & 'x_centroid', 'y_centroid', 'z_centroid', 'c', 'm', 'p', 't', 'theta', 'slip']
+                & 'x_centroid', 'y_centroid', 'z_centroid', 'c', 'm', 'p', 't', 'theta', 'slip', 'mass']
                 call MPI_BCAST(patch_ib(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
             #:endfor
             #:for VAR in ['vel', 'angular_vel', 'angles']

src/simulation/m_time_steppers.fpp

@@ -611,6 +611,8 @@ contains
             if (ib) then
                 ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms
                 if (moving_immersed_boundary_flag) then
+                    call s_compute_ib_forces(q_prim_vf(E_idx)%sf(0:m, 0:n, 0:p)) ! compute the force and torque on the IB from the fluid
+
                     do i = 1, num_ibs
                         if (s == 1) then
                             patch_ib(i)%step_vel = patch_ib(i)%vel
@@ -621,14 +623,22 @@ contains
                             patch_ib(i)%step_z_centroid = patch_ib(i)%z_centroid
                         end if
 
-                        if (patch_ib(i)%moving_ibm == 1) then
-                            do j = 1, 3
-                                patch_ib(i)%vel(j) = (rk_coef(s, 1)*patch_ib(i)%step_vel(j) + rk_coef(s, 2)*patch_ib(i)%vel(j) + rk_coef(s, 3)*0._wp*dt)/rk_coef(s, 4) ! 0.0 is a placeholder for accelerations
-                                patch_ib(i)%angular_vel(j) = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel(j) + rk_coef(s, 2)*patch_ib(i)%angular_vel(j) + rk_coef(s, 3)*0._wp*dt)/rk_coef(s, 4)
+                        if (patch_ib(i)%moving_ibm > 0) then
+                            patch_ib(i)%vel = (rk_coef(s, 1)*patch_ib(i)%step_vel + rk_coef(s, 2)*patch_ib(i)%vel)/rk_coef(s, 4)
+                            patch_ib(i)%angular_vel = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel + rk_coef(s, 2)*patch_ib(i)%angular_vel)/rk_coef(s, 4)
 
-                                ! Update the angle of the IB
-                                patch_ib(i)%angles(j) = (rk_coef(s, 1)*patch_ib(i)%step_angles(j) + rk_coef(s, 2)*patch_ib(i)%angles(j) + rk_coef(s, 3)*patch_ib(i)%angular_vel(j)*dt)/rk_coef(s, 4)
-                            end do
+                            if (patch_ib(i)%moving_ibm == 2) then ! if we are using two-way coupling, apply force and torque
+                                ! update the velocity from the force value
+                                patch_ib(i)%vel = patch_ib(i)%vel + rk_coef(s, 3)*dt*(patch_ib(i)%force/patch_ib(i)%mass)/rk_coef(s, 4)
+
+                                ! update the angular velocity with the torque value
+                                patch_ib(i)%angular_vel = (patch_ib(i)%angular_vel * patch_ib(i)%moment) + (rk_coef(s, 3)*dt*patch_ib(i)%torque/rk_coef(s, 4)) ! add the torque to the angular momentum
+                                call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel) ! update the moment of inertia to be based on the direction of the angular momentum
+                                patch_ib(i)%angular_vel = patch_ib(i)%angular_vel / patch_ib(i)%moment ! convert back to angular velocity with the new moment of inertia
+                            end if
+
+                            ! Update the angle of the IB
+                            patch_ib(i)%angles = (rk_coef(s, 1)*patch_ib(i)%step_angles + rk_coef(s, 2)*patch_ib(i)%angles + rk_coef(s, 3)*patch_ib(i)%angular_vel*dt)/rk_coef(s, 4)
 
                             ! Update the position of the IB
                             patch_ib(i)%x_centroid = (rk_coef(s, 1)*patch_ib(i)%step_x_centroid + rk_coef(s, 2)*patch_ib(i)%x_centroid + rk_coef(s, 3)*patch_ib(i)%vel(1)*dt)/rk_coef(s, 4)