Updates to E-L Solver (Models/documentation/bugs)

1.
Fixes volume fraction and source term smearing. Previously, these two were combined in a convoluted way. Volume fraction field needs to be computed/communicated at the start of the timestep, source term contributions need to be computed/communicated after computing the fluid force on the particle. These are now split in a clean way.

2.
The filling in the buffer cells has the update that uses sum and replace algorithm (implemented by Ben W.). This was further modified to take in an array as input with the indexes of the variables in q_particles/q_beta that should be updated by the algorithm. This was done so that the volume fraction update could be split from the source term contribution update.

3.
The collision force parameters are now defined in the inputs in the particle physical properties. Includes:

particle_pp%ksp_col, particle_pp%nu_col, particle_pp%E_col, particle_pp%cor_col. Need to be set if collisions is turned on.

The collision forces are not communicated anymore. Instead, each local particle is looped through and checked for overlap with neighbors. Then, only the collision force on this local particle of interest is added to this local particle. This avoids communication of forces.

4.
The sutherland viscosity for air is hardcoded into the force subroutine if "viscous" is turned off. This is a temporary bandaid. Use lag_params%mu_ref = 1.716E-5

5.
Implements quasi-steady drag fluctuations force (logical input : lag_params%qs_fluct_force)
The subroutine is within the kernels file. This uses the random number generator in src/common/m_model.fpp. src/common/m_model.fpp was modified to make the random number generator subroutine public.

6.
The documentation for the qs_drag_model was corrected. New documentation is added for the collision force inputs, quasi-steady fluctuation force, and fluid (air) reference viscosity).

Author: jaguilar

docs/documentation/case.md

@@ -931,30 +931,38 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 | `solver_approach`                | Integer | 1: One-way coupling, 2: Two-way coupling                       |
 | `smooth_type`                    | Integer | Smoothing function. 1: Gaussian, 2: Delta 3x3                  |
 | `stokes_drag`                    | Integer | Stokes drag model flag                                          |
-| `qs_drag_model`                  | Integer | Quasi-steady drag model (0: off, 1: Parmar, 2: Modified Parmar, 3: Osnes, 4: Gidaspow) |
+| `qs_drag_model`                  | Integer | Quasi-steady drag model (0: off, 1: Parmar, 2: Osnes, 3: Modified Parmar, 4: Gidaspow) |
 | `added_mass_model`               | Integer | Added mass model (0: off, >0: active)                           |
 | `interpolation_order`            | Integer | Polynomial order for barycentric field interpolation            |
 | `collision_force`                | Logical | Enable soft-sphere DEM particle-particle collisions             |
+| `qs_fluct_force`                 | Logical | Enable quasi-steady drag force fluctuation contribution         |
 | `pressure_force`                 | Logical | Enable pressure gradient force on particles                     |
 | `gravity_force`                  | Logical | Enable gravitational force on particles                         |
 | `write_void_evol`                | Logical | Write void fraction evolution data                              |
 | `epsilonb`                       | Real    | Standard deviation scaling for the Gaussian kernel              |
 | `valmaxvoid`                     | Real    | Maximum void fraction permitted                                 |
+| `mu_ref`                         | Real    | Fluid reference dynamic viscosity                               |
 | `particle_pp%%rho0ref_particle`  | Real    | Reference particle material density                             |
 | `particle_pp%%cp_particle`       | Real    | Particle specific heat capacity                                 |
+| `particle_pp%%ksp_col`           | Real    | Spring stiffness multiplier for collisions                      |
+| `particle_pp%%nu_col`            | Real    | Poisson's ratio used for collisions                             |
+| `particle_pp%%E_col`             | Real    | Young's modulus [Pa] used for collisions                        |
+| `particle_pp%%cor_col`           | Real    | Coefficient of restitution of particles                         |  
 
 - `particles_lagrange` activates the Euler-Lagrange solid particle solver. Particle initial conditions are read from `./input/lag_particles.dat`. The solver tracks non-deformable spherical particles in a compressible carrier flow using volume-averaged source terms (\cite Maeda18).
 
 - `nParticles_glb` specifies the total number of particles across all MPI ranks. Their initial positions, velocities, and radii must be specified in the input file.
 
 - `solver_approach` specifies the coupling method: [1] one-way coupling where particles are advected by the flow but do not influence it, [2] two-way coupling where particle forces are projected back onto the Eulerian grid as source terms.
 
-- `qs_drag_model` selects the quasi-steady drag correlation: [1] Parmar et al. (2010) with Sangani volume fraction correction, [2] Modified Parmar with Osnes et al. (2023) volume fraction correction, [3] Osnes et al. (2023) full correlation with Loth et al. (2021) rarefied regime, [4] Gidaspow (1994) correlation for dense particle suspensions.
+- `qs_drag_model` selects the quasi-steady drag correlation: [1] Parmar et al. (2010) with Sangani volume fraction correction, [2] Osnes et al. (2023) full correlation with Loth et al. (2021) rarefied regime, [3] Modified Parmar with Osnes et al. (2023) volume fraction correction, [4] Gidaspow (1994) correlation for dense particle suspensions.
 
 - `collision_force` activates soft-sphere DEM collisions using a spring-dashpot contact model with Hertzian stiffness. Collision forces between particles on different MPI ranks are communicated via non-blocking point-to-point messaging.
 
 - `interpolation_order` sets the order of the barycentric Lagrange polynomial used to interpolate Eulerian field quantities (pressure, velocity, density) to particle positions. Must be even; the interpolation stencil uses `N/2` points in each direction.
 
+- `mu_ref` is the fluids reference dynamic viscosity at 273.15 K. Used for particle drag if "viscous" is turned off. Only the air sutherland model is currently implemented. If "viscous" is enabled, the true fluid's viscosity is used.
+
 ### 10. Velocity Field Setup {#sec-velocity-field-setup}
 
 | Parameter              | Type    | Description |