Ib collisions (#1348)

Co-authored-by: Spencer Bryngelson <shb@gatech.edu>

Author: danieljvickers

docs/documentation/case.md

@@ -331,6 +331,10 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | `moving_ibm`         | Integer | Sets the method used for IB movement. |
 | `vel(i)`             | Real    | Initial velocity of the moving IB in the i-th direction. |
 | `angular_vel(i)`     | Real    | Initial angular velocity of the moving IB in the i-th direction. |
+| `coefficient_of_restitution`     | Real    | A number 0 to 1 describing how elastic IB collisions are |
+| `collision_model`     | Integer    | Integer to select the collision model being used for IB collisions. |
+| `collision_time`     | Real    | Amount of simulation time used to resolve collisions |
+| `ib_coefficient_of_friction`     | Real    | Coefficient of friction used in IB collisions |
 
 These parameters should be prepended with `patch_ib(j)%` where $j$ is the patch index.
 
@@ -361,6 +365,14 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
 - `angular_vel(i)` is the initial angular velocity of the IB about the x, y, z axes for i=1, 2, 3 in radians per second. When `moving_ibm` equals 2, this rotation rate is just the starting rate of the object, which will then change due to external torques. If `moving_ibm` equals 1, then this is constant if it is a number, or can be described analytically with an expression. 
 
+- `coefficient_of_restitution` is a number from 0 (exclusive) to 1 (inclusive) describing how elastic IB collisions are. 0 is for perfectly inellastic collisions while 1 is for perfectly ellastic collisions.
+
+- `collision_model` is an integer to select the collision model being used for IB collisions. Using 0 disables collisions and collisiono checking. 1 enables the soft-sphere collision model, where all IBs must be circles or sphere and those IBs can collide with each other as well as walls.
+
+- `collision_time` is approximately the amount of simulation time used to resolve collisions. This is handled by modifying the spring gonstant used to apply collision forces.
+
+- `ib_coefficient_of_friction` is the coefficient of friction used in IB collisions.
+
 ### 5. Fluid Material's {#sec-fluid-materials}
 
 | Parameter | Type   | Description                                    |

docs/module_categories.json

@@ -26,7 +26,8 @@
             "m_chemistry",
             "m_acoustic_src",
             "m_body_forces",
-            "m_pressure_relaxation"
+            "m_pressure_relaxation",
+            "m_collisions"
         ]
     },
     {

examples/3D_mibm_sphere_head_on_collision/case.py

@@ -0,0 +1,129 @@
+import json
+import math
+
+# This case file is based on verification cases run in
+# "Analytical solution for the problem of frictional-elastic collisions
+# of spherical particles using the linear model"
+# by Francesco Paolo Di Maio and Alberto Di Renzo
+
+# This can be used to check collision rebound angles to recreate
+# figure 4 of that  paper
+
+Mu = 1.84e-05
+gam_a = 1.4
+
+# lead-up-properties
+velocity = 3.9
+dt = 5.0e-6
+collision_time = 20.0 * dt
+
+# parerticle properties
+radius = 5.0
+collision_angle_degrees = 30.0
+collision_angle_radians = collision_angle_degrees * math.pi / 180.0
+domain_size = 4.0 * radius
+lead_distance = 0.2 * radius
+
+# simulation runs long enough to collide and travel about lead distance away again
+simulation_time = 2.0 * (lead_distance / velocity) + collision_time
+num_time_steps = int(simulation_time / dt)
+num_saves = 10
+t_step_save = int(num_time_steps / num_saves)
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": -0.5 * domain_size,
+            "x_domain%end": 0.5 * domain_size,
+            "y_domain%beg": 0.0,
+            "y_domain%end": domain_size,
+            "z_domain%beg": -0.5 * domain_size,
+            "z_domain%end": 0.5 * domain_size,
+            "cyl_coord": "F",
+            "m": 60,
+            "n": 60,
+            "p": 60,
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": num_time_steps,
+            "t_step_save": t_step_save,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            # Use the 5 equation model
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            # One fluids: air
+            "num_fluids": 1,
+            "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_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": -15,
+            "bc_y%end": -15,
+            "bc_z%beg": -3,
+            "bc_z%end": -3,
+            # Set IB to True and add 1 patch
+            "ib": "T",
+            "num_ibs": 1,
+            # 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": 9,
+            "patch_icpp(1)%x_centroid": 0.0,
+            "patch_icpp(1)%y_centroid": 0.5 * domain_size,
+            "patch_icpp(1)%z_centroid": 0.0,
+            "patch_icpp(1)%length_x": domain_size,
+            "patch_icpp(1)%length_y": domain_size,
+            "patch_icpp(1)%length_z": domain_size,
+            # Specify the patch primitive variables
+            "patch_icpp(1)%vel(1)": 0.0e00,
+            "patch_icpp(1)%vel(2)": 0.0e00,
+            "patch_icpp(1)%vel(3)": 0.0e00,
+            "patch_icpp(1)%pres": 1.0e00,
+            "patch_icpp(1)%alpha_rho(1)": 1.0e00,
+            "patch_icpp(1)%alpha(1)": 1.0e00,
+            # Patch: Sphere Immersed Boundary
+            "patch_ib(1)%geometry": 8,
+            "patch_ib(1)%x_centroid": -1.0 * lead_distance * math.sin(collision_angle_radians),  # get a lead up distance to the collision
+            "patch_ib(1)%y_centroid": radius + lead_distance * math.sin(collision_angle_radians),
+            "patch_ib(1)%z_centroid": 0.0,
+            "patch_ib(1)%radius": radius,
+            "patch_ib(1)%slip": "F",
+            "patch_ib(1)%mass": 1.0e6,  # arbitrarily high mass to ignore fluid
+            "patch_ib(1)%vel(1)": velocity * math.sin(collision_angle_radians),
+            "patch_ib(1)%vel(2)": -velocity * math.cos(collision_angle_radians),
+            "patch_ib(1)%moving_ibm": 2,
+            # Collisions
+            "collision_model": 1,  # soft-sphere collision model
+            "ib_coefficient_of_friction": 0.092,
+            "collision_time": collision_time,
+            "coefficient_of_restitution": 0.98,  # almost perfectly elastic
+            # Fluids Physical Parameters
+            "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50(Not 1.40)
+            "fluid_pp(1)%pi_inf": 0,
+        }
+    )
+)

src/common/m_constants.fpp

@@ -23,7 +23,8 @@ 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 = 1000              !< Maximum number of IC patches
+    integer, parameter  :: num_patches_max = 10                !< Maximum number of IC patches
+    integer, parameter  :: num_ib_patches_max = 50000          !< Maximum number of immersed boundary patches (patch_ib)
     integer, parameter  :: num_bc_patches_max = 10             !< Maximum number of boundary condition patches
     integer, parameter  :: max_2d_fourier_modes = 10           !< Max Fourier mode index for 2D modal patch (geometry 13)
     integer, parameter  :: max_sph_harm_degree = 5             !< Max degree L for 3D spherical harmonic patch (geometry 14)

src/common/m_helper.fpp

@@ -18,7 +18,7 @@ module m_helper
     public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, &
         & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, &
         & f_create_bbox, s_print_2D_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, &
-        & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data
+        & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product
 
 contains
 
@@ -297,7 +297,22 @@ contains
 
     end function f_cross
 
-    !> Swap two real numbers.
+    !> @brief Computes the cross product c = a x b of two 3D vectors.
+    subroutine s_cross_product(a, b, c)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+        real(wp), intent(in)  :: a(3), b(3)
+        real(wp), intent(out) :: c(3)
+
+        c(1) = a(2)*b(3) - a(3)*b(2)
+        c(2) = a(3)*b(1) - a(1)*b(3)
+        c(3) = a(1)*b(2) - a(2)*b(1)
+
+    end subroutine s_cross_product
+
+    !> This procedure swaps two real numbers.
+    !! @param lhs Left-hand side.
+    !! @param rhs Right-hand side.
     elemental subroutine s_swap(lhs, rhs)
 
         real(wp), intent(inout) :: lhs, rhs

src/post_process/m_data_input.f90

@@ -108,12 +108,13 @@ impure subroutine s_read_ib_data_files(file_loc_base, t_step)
         integer, intent(in), optional               :: t_step
         character(LEN=len_trim(file_loc_base) + 20) :: file_loc
         logical                                     :: file_exist
-        integer                                     :: ifile, ierr, data_size, var_MOK
+        integer                                     :: ifile, ierr, data_size
 
 #ifdef MFC_MPI
         integer, dimension(MPI_STATUS_SIZE) :: status
         integer(KIND=MPI_OFFSET_KIND)       :: disp
-        integer                             :: m_MOK, n_MOK, p_MOK, MOK, WP_MOK, save_index
+        integer(KIND=MPI_OFFSET_KIND)       :: m_MOK, n_MOK, p_MOK, MOK, WP_MOK, var_MOK
+        integer                             :: save_index
 #endif
 
         if (.not. ib) return

src/pre_process/m_global_parameters.fpp

@@ -153,11 +153,11 @@ module m_global_parameters
 
     !> @name Immersed Boundaries
     !> @{
-    logical                                               :: ib        !< Turn immersed boundaries on
-    integer                                               :: num_ibs   !< Number of immersed boundaries
-    integer                                               :: Np
-    type(ib_patch_parameters), dimension(num_patches_max) :: patch_ib  !< Immersed boundary patch parameters
-    type(vec3_dt), allocatable, dimension(:)              :: airfoil_grid_u, airfoil_grid_l
+    logical                                                  :: ib        !< Turn immersed boundaries on
+    integer                                                  :: num_ibs   !< Number of immersed boundaries
+    integer                                                  :: Np
+    type(ib_patch_parameters), dimension(num_ib_patches_max) :: patch_ib  !< Immersed boundary patch parameters
+    type(vec3_dt), allocatable, dimension(:)                 :: airfoil_grid_u, airfoil_grid_l
     !> @}
 
     !> @name Non-polytropic bubble gas compression
@@ -435,7 +435,7 @@ contains
         ib = .false.
         num_ibs = dflt_int
 
-        do i = 1, num_patches_max
+        do i = 1, num_ib_patches_max
             patch_ib(i)%geometry = dflt_int
             patch_ib(i)%x_centroid = dflt_real
             patch_ib(i)%y_centroid = dflt_real

src/pre_process/m_mpi_proxy.fpp

@@ -112,8 +112,9 @@ contains
             if (chemistry) then
                 call MPI_BCAST(patch_icpp(i)%Y, size(patch_icpp(i)%Y), mpi_p, 0, MPI_COMM_WORLD, ierr)
             end if
-            ! Broadcast IB variables: patch_ib is indexed 1:num_patches_max, not 1:num_bc_patches_max, so these must live in the
-            ! num_patches_max loop.
+        end do
+
+        do i = 1, num_ibs
             #:for VAR in ['vel', 'angular_vel', 'angles']
                 call MPI_BCAST(patch_ib(i)%${VAR}$, size(patch_ib(i)%${VAR}$), mpi_p, 0, MPI_COMM_WORLD, ierr)
             #:endfor