Solid particle implementation for E-L solver

Author: sbryngelson

docs/documentation/case.md

@@ -266,6 +266,16 @@ Setup: Only requires specifying `init_dir` and filename pattern via `zeros_defau
 Implementation: All variables and file handling are managed in `src/common/include/ExtrusionHardcodedIC.fpp` with no manual grid configuration needed.
 Usage: Ideal for initializing simulations from lower-dimensional solutions, enabling users to add perturbations or modifications to the base extruded fields for flow instability studies.
 
+The following parameters support hardcoded initial conditions that read interface data from files:
+
+| Parameter        | Type    | Description                                               |
+| ---:             | :---:   | :---                                                      |
+| `interface_file` | String  | Path to interface geometry data file                      |
+| `normFac`        | Real    | Interface normalization factor                            |
+| `normMag`        | Real    | Interface normal magnitude                                |
+| `g0_ic`          | Real    | Initial gas volume fraction for interfacial IC            |
+| `p0_ic`          | Real    | Initial pressure for interfacial IC                       |
+
 #### Parameter Descriptions
 
 - `num_patches` defines the total number of patches defined in the domain.
@@ -788,7 +798,7 @@ Details of the transducer acoustic source model can be found in \cite Maeda17.
 | ---:               | :----:  |          :---                                  |
 | `bubbles_euler`    | Logical | Ensemble-averaged bubble modeling |
 | `bubbles_lagrange` | Logical | Volume-averaged bubble modeling |
-| `bubble_model`     | Integer | [1] Gilmore; [2] Keller--Miksis; [3] Rayleigh-Plesset |
+| `bubble_model`     | Integer | [0] Particle; [1] Gilmore; [2] Keller--Miksis; [3] Rayleigh-Plesset |
 | `Ca`               | Real    | Cavitation number |
 | `Web`              | Real    | Weber number |
 | `Re_inv`           | Real    | Inverse Reynolds number |
@@ -892,6 +902,13 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 | `epsilonb`            | Real    | Standard deviation scaling for the gaussian function      |
 | `charwidth`           | Real    | Domain virtual depth (z direction, for 2D simulations)    |
 | `valmaxvoid`          | Real    | Maximum void fraction permitted                           |
+| `drag_model`          | Integer | Drag model for bubble dynamics                            |
+| `vel_model`           | Integer | Velocity model for bubble interface                       |
+| `charNz`              | Integer | Characteristic size parameter                             |
+| `input_path`          | String  | Path to bubble input file (default: `./input`)            |
+| `pressure_force`      | Logical | Enable pressure gradient force                            |
+| `gravity_force`       | Logical | Enable gravitational force                                |
+| `write_void_evol`     | Logical | Write void fraction evolution data                        |
 
 - `nBubs_glb` Total number of bubbles. Their initial conditions need to be specified in the ./input/lag_bubbles.dat file. See the example cases for additional information.
 
@@ -905,6 +922,39 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 
 - `massTransfer_model` Activates the mass transfer model at the bubble's interface based on (\cite Preston07).
 
+#### 9.3 Lagrangian Solid Particle Model
+
+| Parameter                        | Type    | Description                                                     |
+| ---:                             | :---:   | :---                                                            |
+| `particles_lagrange`             | Logical | Lagrangian solid particle model switch                          |
+| `nParticles_glb`                 | Integer | Global number of particles                                      |
+| `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) |
+| `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             |
+| `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                                 |
+| `particle_pp%%rho0ref_particle`  | Real    | Reference particle material density                             |
+| `particle_pp%%cp_particle`       | Real    | Particle specific heat capacity                                 |
+
+- `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.
+
+- `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.
+
 ### 10. Velocity Field Setup {#sec-velocity-field-setup}
 
 | Parameter              | Type    | Description |

docs/documentation/equations.md

@@ -514,6 +514,46 @@ with \f$\sigma = \varepsilon_b \max(\Delta x^{1/3}_\text{cell},\;R_\text{bubble}
 
 Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order adaptive Runge-Kutta time integration.
 
+### 6.3 Euler-Lagrange Solid Particles (`particles_lagrange = .true.`)
+
+**Source:** `src/simulation/m_particles_EL.fpp`, `src/simulation/m_particles_EL_kernels.fpp`
+
+The Euler-Lagrange particle solver tracks non-deformable solid particles in a compressible carrier flow.
+The volume-averaged carrier flow equations use the same source term framework as the bubble model (Section 6.2),
+with the particle volume fraction \f$\alpha_p\f$ replacing the bubble void fraction \f$\alpha\f$.
+
+**Particle volume fraction via Gaussian kernel:**
+
+\f[\alpha_p(\mathbf{x}) = \sum_n V_{p,n}\,\delta_\sigma(\mathbf{x} - \mathbf{x}_n)\f]
+
+where \f$V_{p,n} = \frac{4}{3}\pi R_n^3\f$ is the volume of particle \f$n\f$ and \f$\delta_\sigma\f$ is the Gaussian regularization kernel from Section 6.2.
+
+**Particle equation of motion:**
+
+\f[m_p \frac{d\mathbf{u}_p}{dt} = \mathbf{F}_\text{drag} + \mathbf{F}_\text{pressure} + \mathbf{F}_\text{AM} + \mathbf{F}_\text{gravity} + \mathbf{F}_\text{collision}\f]
+
+**Quasi-steady drag force:**
+
+\f[\mathbf{F}_\text{drag} = \beta\,(\mathbf{u}_f - \mathbf{u}_p)\f]
+
+where \f$\beta\f$ is a drag coefficient computed from one of several correlations selected via `qs_drag_model`:
+- Parmar et al. (2010): Re and Ma corrections with Sangani et al. (1991) volume fraction correction
+- Modified Parmar: Re and Ma corrections with Osnes et al. (2023) volume fraction correction
+- Osnes et al. (2023): comprehensive correlation for compressible flow through random particle suspensions
+- Gidaspow (1994): Ergun/Wen-Yu correlation for dense suspensions
+
+**Pressure gradient force** (`pressure_force = .true.`):
+
+\f[\mathbf{F}_\text{pressure} = -V_p\,\nabla p\f]
+
+**Added mass force** (`added_mass_model > 0`):
+
+\f[\mathbf{F}_\text{AM} = C_\text{AM}\,\rho_f\,V_p\left(\frac{D\mathbf{u}_f}{Dt} - \frac{d\mathbf{u}_p}{dt}\right)\f]
+
+**Collision force** (`collision_force = .true.`):
+
+Soft-sphere DEM model with Hertzian contact stiffness and viscous damping.
+
 ---
 
 ## 7. Fluid-Structure Interaction

docs/module_categories.json

@@ -19,6 +19,8 @@
             "m_bubbles_EE",
             "m_bubbles_EL",
             "m_bubbles_EL_kernels",
+            "m_particles_EL",
+            "m_particles_EL_kernels",
             "m_qbmm",
             "m_hyperelastic",
             "m_hypoelastic",