Solid Particle Implementation for E-L Solver

.gitignore

@@ -114,3 +114,4 @@ cce_*/
 cce_*.log
 run_cce_*.sh
 .ffmt_cache/
+**/.ffmt_cache/

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.
@@ -349,7 +359,7 @@ Definitions for currently implemented immersed boundary patch types are listed i
 
 - `c`, `t`, `p`, and `m` specify the parameters for a NACA airfoil.
 `m` is the maximum camber, `p` is the location of maximum camber, `c` is the coord length, and `t` is the thickness.
-Additional details on this specification can be found in [NACA airfoil](https://en.wikipedia.org/wiki/NACA_airfoil).
+Additional details on this specification can be found in [The Naca Airfoil Series](https://web.stanford.edu/~cantwell/AA200_Course_Material/The%20NACA%20airfoil%20series.pdf)
 
 - `slip` applies a slip boundary to the surface of the patch if true and a no-slip boundary condition to the surface if false.
 
@@ -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,47 @@ 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: 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] 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 |

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",

examples/2D_moving_lag_bubs/case.py

@@ -0,0 +1,202 @@
+#!/usr/bin/env python3
+import argparse
+import json
+import math
+
+import numpy as np
+
+# Domain extents
+xb = -2
+xe = 2
+yb = -2
+ye = 2
+cw = 0.05  # Characteristic width for 3D bubble cloud
+
+# Reference values for nondimensionalization
+L0 = 1e-3  # length - m (min bubble radius)
+rho0 = 950  # density - kg/m3
+c0 = 1449.0  # speed of sound - m/s
+p0 = rho0 * c0 * c0  # pressure - Pa
+T0 = 298  # temperature - K
+
+# Host properties (water)
+gamma_host = 6.12  # Specific heat ratio
+pi_inf_host = 3.43e8  # Stiffness - Pa
+mu_host = 0.001
+c_host = 1449.0  # speed of sound - m/s
+rho_host = 950  # density kg/m3
+T_host = 298  # temperature K
+
+# Lagrangian bubbles' properties
+R_uni = 8314  # Universal gas constant - J/kmol/K
+MW_g = 28.0  # Molar weight of the gas - kg/kmol
+MW_v = 18.0  # Molar weight of the vapor - kg/kmol
+gamma_g = 1.4  # Specific heat ratio of the gas
+gamma_v = 1.333  # Specific heat ratio of the vapor
+pv = 2350  # Vapor pressure of the host - Pa
+cp_g = 1.0e3  # Specific heat of the gas - J/kg/K
+cp_v = 2.1e3  # Specific heat of the vapor - J/kg/K
+k_g = 0.025  # Thermal conductivity of the gas - W/m/K
+k_v = 0.02  # Thermal conductivity of the vapor - W/m/K
+diffVapor = 2.5e-5  # Diffusivity coefficient of the vapor - m2/s
+sigBubble = 0.020  # Surface tension of the bubble - N/m
+mu_g = 1.48e-5
+rho_g = 1  # density kg/m3
+
+# Nondimmensionalization of domain size
+xb = xb / L0
+xe = xe / L0
+yb = yb / L0
+ye = ye / L0
+cw = cw / L0
+
+# patm = 1.0e05  # Atmospheric pressure - Pa
+patm = 1e5
+g0 = 9.81 / (c0 * c0 / L0)
+
+# Patch prim vars
+rho_host = rho_host / rho0
+rho_g = rho_g / rho0
+pres = patm / p0
+
+# Timing
+tend = 0.004 * c0 / L0
+tsave = tend / 50  # save time - sec
+
+c_host = math.sqrt(gamma_host * (pres + pi_inf_host / p0) / rho_host)
+dx = (xe - xb) / (199 + 1)
+dt = 0.05 * dx / c_host
+
+t_step_stop = int(tend / dt)
+t_step_save = int(t_step_stop / 100)
+
+eps = 1e-8
+
+# Configuration case dictionary
+data = {
+    # Logistics
+    "run_time_info": "T",
+    # Computational Domain
+    "x_domain%beg": xb,
+    "x_domain%end": xe,
+    "y_domain%beg": yb,
+    "y_domain%end": ye,
+    "m": 199,
+    "n": 199,
+    "p": 0,
+    "cyl_coord": "F",
+    "dt": 1,
+    "t_step_start": 0,
+    "t_step_stop": 6000,
+    "t_step_save": 60,
+    "adap_dt": "T",
+    "adap_dt_max_iters": 1000,
+    # Simulation Algorithm
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "mixture_err": "F",
+    "mpp_lim": "T",
+    "time_stepper": 3,
+    "weno_order": 5,
+    "mapped_weno": "T",
+    "mp_weno": "T",
+    "avg_state": 2,
+    "weno_eps": 1e-16,
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "bc_x%beg": -1,
+    "bc_x%end": -1,
+    "bc_y%beg": -1,
+    "bc_y%end": -1,
+    "num_patches": 2,
+    "num_fluids": 2,
+    # Database Structure Parameters
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "T",
+    "lag_db_wrt": "T",
+    "lag_pres_wrt": "T",
+    # Fluid Parameters Host
+    "fluid_pp(1)%gamma": 1.0 / (gamma_host - 1.0),
+    "fluid_pp(1)%pi_inf": gamma_host * (pi_inf_host / p0) / (gamma_host - 1.0),
+    # Fluid Parameters Gas
+    "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0),
+    "fluid_pp(2)%pi_inf": 0.0e00,
+    # Bubble parameters
+    "bub_pp%R0ref": 1.0,
+    "bub_pp%p0ref": 1.0,
+    "bub_pp%rho0ref": 1.0,
+    "bub_pp%T0ref": 1.0,
+    "bub_pp%ss": sigBubble / (rho0 * L0 * c0 * c0),
+    "bub_pp%pv": pv / p0,
+    "bub_pp%vd": diffVapor / (L0 * c0),
+    "bub_pp%mu_l": mu_host / (rho0 * L0 * c0),
+    "bub_pp%gam_v": gamma_v,
+    "bub_pp%gam_g": gamma_g,
+    "bub_pp%M_v": MW_v,
+    "bub_pp%M_g": MW_g,
+    "bub_pp%k_v": k_v * (T0 / (L0 * rho0 * c0 * c0 * c0)),
+    "bub_pp%k_g": k_g * (T0 / (L0 * rho0 * c0 * c0 * c0)),
+    "bub_pp%cp_v": cp_v * (T0 / (c0 * c0)),
+    "bub_pp%cp_g": cp_g * (T0 / (c0 * c0)),
+    "bub_pp%R_v": (R_uni / MW_v) * (T0 / (c0 * c0)),
+    "bub_pp%R_g": (R_uni / MW_g) * (T0 / (c0 * c0)),
+    # Viscosity
+    "viscous": "T",
+    "fluid_pp(1)%Re(1)": 1.0 / (mu_host / (rho0 * c0 * L0)),
+    "fluid_pp(2)%Re(1)": 1.0 / (mu_g / (rho0 * c0 * L0)),
+    # Patch for background flow
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": (xb + xe) / 2,
+    "patch_icpp(1)%y_centroid": (yb + ye) / 2,
+    "patch_icpp(1)%length_x": (xe - xb),
+    "patch_icpp(1)%length_y": (ye - yb),
+    "patch_icpp(1)%vel(1)": 0,
+    "patch_icpp(1)%vel(2)": 0,
+    "patch_icpp(1)%pres": pres,
+    "patch_icpp(1)%alpha_rho(1)": (1 - eps) * rho_host,
+    "patch_icpp(1)%alpha_rho(2)": eps * rho_g,
+    "patch_icpp(1)%alpha(1)": 1 - eps,
+    "patch_icpp(1)%alpha(2)": eps,
+    # High pressure
+    "patch_icpp(2)%geometry": 2,
+    "patch_icpp(2)%alter_patch(1)": "T",
+    "patch_icpp(2)%smoothen": "T",
+    "patch_icpp(2)%smooth_patch_id": 1,
+    "patch_icpp(2)%smooth_coeff": 0.5,
+    "patch_icpp(2)%x_centroid": 0,
+    "patch_icpp(2)%y_centroid": 0,
+    "patch_icpp(2)%radius": 0.2 / L0,
+    "patch_icpp(2)%vel(1)": 0,
+    "patch_icpp(2)%vel(2)": 0,
+    "patch_icpp(2)%pres": 10 * pres,
+    "patch_icpp(2)%alpha_rho(1)": (1 - eps) * rho_host,
+    "patch_icpp(2)%alpha_rho(2)": eps * rho_g,
+    "patch_icpp(2)%alpha(1)": 1 - eps,
+    "patch_icpp(2)%alpha(2)": eps,
+    # Lagrangian Bubbles
+    "bubbles_lagrange": "T",
+    "fd_order": 4,
+    "bubble_model": 2,  # (0) Particle (2) KM (3) RP
+    "thermal": 3,
+    "polytropic": "F",
+    "lag_params%nBubs_glb": 5000,  # Max number of bubbles
+    "lag_params%vel_model": 2,
+    "lag_params%drag_model": 1,
+    "lag_params%solver_approach": 2,
+    "lag_params%cluster_type": 2,
+    "lag_params%pressure_corrector": "F",
+    "lag_params%smooth_type": 1,
+    "lag_params%epsilonb": 1.0,
+    "lag_params%valmaxvoid": 0.9,
+    "lag_params%write_bubbles": "T",
+    "lag_params%write_bubbles_stats": "T",
+    "lag_params%write_void_evol": "T",
+    "lag_params%charwidth": cw,
+    "lag_params%charNz": 10,
+}
+
+mods = {}
+
+print(json.dumps({**data, **mods}, indent=4))

examples/3D_lagrange_shbubcollapse/case.py

@@ -155,6 +155,7 @@
             "lag_params%valmaxvoid": 0.9,
             "lag_params%write_bubbles": "F",
             "lag_params%write_bubbles_stats": "F",
+            "lag_params%write_void_evol": "T",
             # Bubble parameters
             "bub_pp%R0ref": 1.0,
             "bub_pp%p0ref": 1.0,