@@ -802,6 +802,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca
 | `lag_mg_wrt`            | Logical | Add the bubble gas mass to the database file |
 | `lag_betaT_wrt`         | Logical | Add the bubble heat flux model coefficient to the database file |
 | `lag_betaC_wrt`         | Logical | Add the bubble mass flux model coefficient to the database file |
+| `lag_voidfrac_wrt`      | Logical | Add the particle volume fraction in the host cell to the database file (`particles_lagrange` only) |
 
 The table lists formatted database output parameters. The parameters define variables that are outputted from simulation and file types and formats of data as well as options for post-processing.
 
@@ -1084,6 +1085,41 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 
 - `kahan_summation` uses Kahan compensated summation when smearing the bubble contributions onto the Eulerian void fraction, reducing the round-off sensitivity of the accumulation to the summation order. It is not compatible with `--mixed` precision builds.
 
+#### 9.3 Euler-Lagrange Solid Particle Model
+
+| Parameter                                 | Type    | Description                                                         |
+| ---:                                      | :---:   | :---                                                                |
+| `particles_lagrange`                      | Logical | Lagrangian solid particle model switch                              |
+| `particle_pp%%rho0ref_particle`           | Real    | Particle material density                                           |
+| `particle_params%%nparticles_glb`         | Integer | Global number of particles                                          |
+| `particle_params%%input_path`             | String  | Path to the particle input file                                     |
+| `particle_params%%solver_approach`        | Integer | 1: One-way coupling, 2: two-way coupling                            |
+| `particle_params%%stationary`             | Logical | Keep the particles fixed in space (default false)                   |
+| `particle_params%%qs_force`                | Integer | Quasi-steady drag: 0 off, 1 Gidaspow, 2 Parmar, 3 Osnes             |
+| `particle_params%%qs_fluct_force`         | Logical | Quasi-steady drag fluctuation force                                 |
+| `particle_params%%pressure_gradient_force`| Logical | Pressure-gradient force                                             |
+| `particle_params%%added_mass_force`       | Integer | Added-mass force: 0 off, 1 on                                       |
+| `particle_params%%mu_ref(i)`              | Real    | Reference viscosity of fluid $i$ for the drag (inviscid cases)      |
+| `particle_params%%suth(i)`                | Real    | Sutherland constant of fluid $i$ (optional)                         |
+| `particle_params%%interpolation_order`    | Integer | Even order (2 to 8) of the fluid-to-particle barycentric interpolation |
+| `particle_params%%epsilonb`               | Real    | Standard deviation scaling for the Gaussian kernel                  |
+| `particle_params%%charwidth`              | Real    | Domain virtual depth (z direction, for 2D simulations)              |
+| `particle_params%%valmaxvoid`             | Real    | Maximum particle volume fraction permitted                          |
+| `particle_params%%write_particles`        | Logical | Write the particle evolution to `D/lag_particle_evol_<rank>.dat`    |
+| `particle_params%%write_void_evol`        | Logical | Write the volume fraction evolution over time                       |
+
+- `particles_lagrange` activates the Euler-Lagrange solid particle model: rigid spherical particles are tracked individually and projected onto the grid with the Gaussian kernel of \cite Maeda18. It requires a 2D or 3D case, `model_eqns = 2` and `fd_order = 2` or `4` (the particle gradient stencil and halo width), and cannot be combined with `bubbles_lagrange`, `igr` or `cyl_coord`. The added-mass force needs two-way coupling (`solver_approach = 2`), which computes the fluid acceleration it uses, and `qs_fluct_force` needs a quasi-steady drag model (`qs_force > 0`). Restarts and particle post-processing output (`lag_db_wrt`) need ``parallel_io = 'T'``, which writes the particle restart file. A restart must use the same number of MPI ranks as the run that wrote the particle restart file (each rank reads the particles that rank wrote, as for Lagrangian bubbles); the simulation aborts otherwise. The restart file keeps each particle's drag-fluctuation state and random seed, so a restarted run continues the same random sequence. Known limitation: the added-mass force uses the fluid acceleration from the previous Runge-Kutta stage, which is not saved, so the first stage after a restart uses zero fluid acceleration. With `added_mass_force = 1` a restarted run therefore differs slightly from an uninterrupted one (about 5e-5 relative in the particle velocities over the next 10 steps of the regression case); without the added-mass force, restarts reproduce the uninterrupted run to round-off. The interpolation stencil must stay within the ghost layers holding valid data, so `interpolation_order` is at most 8, and at most `4 + fd_order` with the pressure-gradient or added-mass force (6 for `fd_order = 2`). Periodic and reflective (symmetry) boundaries and particle collisions are not yet supported.
+
+- `input_path` Path to the particle input file. Each row specifies one particle, with columns `x  y  z  u  v  w  radius`; any extra columns are ignored.
+
+- `qs_force` selects the quasi-steady drag correlation: Gidaspow \cite Gidaspow1994, Parmar et al. \cite Parmar2010 with the volume-fraction correction of \cite Sangani1991, or Osnes et al. \cite Osnes2023, whose Loth et al. \cite Loth2021 single-particle drag uses the continuous reformulation of \cite Daoud2026. In inviscid cases without chemistry, the drag viscosity of fluid $i$ is `mu_ref(i)`, corrected with Sutherland's law when `suth(i)` is given.
+
+- `pressure_gradient_force` and `added_mass_force` use the fluid pressure, density and velocity gradients at the particles, taken from the reconstructed cell-face states. `fd_order` sets the stencil of the gradients in the two-way source terms and the particle halo width.
+
+- `qs_fluct_force` adds the stochastic quasi-steady force fluctuations of \cite Osnes2023 and \cite Lattanzi2022, advanced once per time step as an Ornstein-Uhlenbeck process whose relaxation rate uses the radial distribution function of \cite MaAhmadi1986.
+
+- The Lagrangian post-processing flags (`lag_db_wrt`, `lag_pos_wrt`, `lag_vel_wrt`, ...) also apply to particles, which are written as the lag_particles point mesh in the Silo output. `lag_voidfrac_wrt` adds the particle volume fraction in each particle's host cell.
+
 ### 10. Velocity Field Setup {#sec-velocity-field-setup}
 
 | Parameter              | Type    | Description |
 
@@ -29,6 +29,9 @@
             "m_bubbles_EE",
             "m_bubbles_EL",
             "m_bubbles_EL_kernels",
+            "m_particles_EL",
+            "m_particles_EL_kernels",
+            "m_euler_lagrange",
             "m_qbmm",
             "m_hypoelastic",
             "m_phase_change",
 
@@ -195,6 +195,84 @@ @article{Maeda18
   doi     = {10.1016/j.jcp.2018.05.029}
 }
 
+% --- Lagrangian particle drag ---
+
+@article{Daoud2026,
+  author  = {T. Daoud and T. Jackson and S. Balachandar},
+  title   = {A careful examination of closure models in {Euler--Lagrange} simulations of compressible multiphase flow in a planar shock particle curtain problem},
+  journal = {International Journal of Multiphase Flow},
+  pages   = {105740},
+  year    = {2026},
+  doi     = {10.1016/j.ijmultiphaseflow.2026.105740}
+}
+
+@article{Osnes2023,
+  author  = {A. N. Osnes and M. Vartdal and M. Khalloufi and J. Capecelatro and S. Balachandar},
+  title   = {Comprehensive quasi-steady force correlations for compressible flow through random particle suspensions},
+  journal = {International Journal of Multiphase Flow},
+  volume  = {165},
+  pages   = {104485},
+  year    = {2023},
+  doi     = {10.1016/j.ijmultiphaseflow.2023.104485}
+}
+
+@article{Loth2021,
+  author  = {E. Loth and J. T. Daspit and M. Jeong and T. Nagata and T. Nonomura},
+  title   = {Supersonic and hypersonic drag coefficients for a sphere},
+  journal = {AIAA Journal},
+  volume  = {59},
+  number  = {8},
+  pages   = {3261--3274},
+  year    = {2021},
+  doi     = {10.2514/1.J060153}
+}
+
+@article{Parmar2010,
+  author  = {M. Parmar and A. Haselbacher and S. Balachandar},
+  title   = {Improved drag correlation for spheres and application to shock-tube experiments},
+  journal = {AIAA Journal},
+  volume  = {48},
+  number  = {6},
+  pages   = {1273--1276},
+  year    = {2010}
+}
+
+@article{Sangani1991,
+  author  = {A. S. Sangani and D. Z. Zhang and A. Prosperetti},
+  title   = {The added mass, {Basset}, and viscous drag coefficients in nondilute bubbly liquids undergoing small-amplitude oscillatory motion},
+  journal = {Physics of Fluids A},
+  volume  = {3},
+  number  = {12},
+  pages   = {2955--2970},
+  year    = {1991}
+}
+
+@book{Gidaspow1994,
+  author    = {D. Gidaspow},
+  title     = {Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions},
+  publisher = {Academic Press},
+  year      = {1994}
+}
+
+@article{MaAhmadi1986,
+  author  = {D. Ma and G. Ahmadi},
+  title   = {An equation of state for dense rigid sphere gases},
+  journal = {Journal of Chemical Physics},
+  volume  = {84},
+  number  = {6},
+  pages   = {3449--3450},
+  year    = {1986}
+}
+
+@article{Lattanzi2022,
+  author  = {A. M. Lattanzi and V. Tavanashad and S. Subramaniam and J. Capecelatro},
+  title   = {Stochastic model for the hydrodynamic force in {Euler--Lagrange} simulations of particle-laden flows},
+  journal = {Physical Review Fluids},
+  volume  = {7},
+  pages   = {014301},
+  year    = {2022}
+}
+
 % --- Elasticity ---
 
 @article{Rodriguez19,
 
@@ -0,0 +1,18 @@
+# 2D Shock–Particle Curtain Interaction
+
+A Mach 1.66 planar shock, held by a Dirichlet inflow, hits a 2 mm thick curtain of glass particles (radius 57.5 µm, volume fraction 0.21) in 82.7 kPa air between slip walls. It uses the Euler-Lagrange solid-particle solver (`particles_lagrange`) with two-way coupling, Osnes quasi-steady drag, pressure-gradient force and added mass.
+
+The setup follows the multiphase shock tube experiment of Wagner et al. (2012):
+> J. L. Wagner, S. J. Beresh, S. P. Kearney, W. M. Trott, J. N. Castaneda, B. O. Pruett, and M. R. Baer, "A multiphase shock tube for shock wave interactions with dense particle fields", Experiments in Fluids, vol. 52, no. 6, pp. 1507–1517, 2012. https://doi.org/10.1007/s00348-012-1272-x
+
+`gen_particles.py` writes the particles (seeded, uniform in the curtain band) to `input/particles.dat`; `case.py` calls it when pre_process runs, or run it yourself with `python3 gen_particles.py`. In 2D each particle stands for a slab of depth `charwidth`, so projected particles may overlap.
+
+```shell
+./mfc.sh run examples/2D_particle_curtain/case.py -n 4
+```
+
+## Result
+
+Density-gradient schlieren |∇ρ|/ρ with the particles (white) at t = 500 µs.
+
+<img src="result.png"/>
@@ -0,0 +1,130 @@
+#!/usr/bin/env python3
+# Planar shock (driven by a Dirichlet inflow) hitting a dense particle curtain between slip walls
+# (Euler-Lagrange, two-way coupling). The particles come from gen_particles.py.
+import argparse
+import json
+import math
+import os
+
+import gen_particles
+
+parser = argparse.ArgumentParser(prog="2D_particle_curtain", formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT", help="MFC's toolchain's internal state.")
+args = parser.parse_args()
+
+# Write input/particles.dat only when ./mfc.sh run runs pre_process (validate, test and docs also load this file)
+if args.mfc.get("command") == "run" and "pre_process" in args.mfc.get("targets", []):
+    gen_particles.write_particles(os.path.dirname(os.path.abspath(__file__)))
+
+# Air (ideal gas)
+gamma = 1.4
+R_air = 287.0  # J/(kg K)
+
+
+def normal_shock(M, p1, rho1):
+    """Pressure, density and lab-frame velocity behind a normal shock of Mach number M moving into still gas (p1, rho1)."""
+    p2 = p1 * (1.0 + 2.0 * gamma / (gamma + 1.0) * (M**2 - 1.0))
+    rho2 = rho1 * (gamma + 1.0) * M**2 / ((gamma - 1.0) * M**2 + 2.0)
+    u2 = 2.0 / (gamma + 1.0) * math.sqrt(gamma * p1 / rho1) * (M - 1.0 / M)
+    return p2, rho2, u2
+
+
+# Ambient air (Wagner et al. 2012) and the state behind the incident shock
+M_shock = 1.66
+p_amb, T_amb = 82700.0, 296.4
+rho_amb = p_amb / (R_air * T_amb)
+p_post, rho_post, u_post = normal_shock(M_shock, p_amb, rho_amb)
+
+# Grid: x in [-0.2, 0.3], y in [0, 0.015] (slip walls at y = 0 and y = 0.015)
+xb, xe, yb, ye = -0.2, 0.3, 0.0, 0.015
+Nx, Ny = 2000, 30
+dx = (xe - xb) / Nx
+
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational domain
+            "x_domain%beg": xb,
+            "x_domain%end": xe,
+            "y_domain%beg": yb,
+            "y_domain%end": ye,
+            "m": Nx - 1,
+            "n": Ny - 1,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.4,
+            "n_start": 0,
+            "t_stop": 5.0e-4,
+            "t_save": 5.0e-5,
+            # Simulation algorithm
+            "model_eqns": 2,
+            "num_fluids": 1,
+            "num_patches": 2,
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "mapped_weno": "T",
+            "mp_weno": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -17,
+            "bc_x%end": -8,
+            "bc_y%beg": -15,
+            "bc_y%end": -15,
+            # Output
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "T",
+            "lag_db_wrt": "T",
+            "lag_voidfrac_wrt": "T",
+            # Patch 1: ambient air
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5 * (xb + xe),
+            "patch_icpp(1)%y_centroid": 0.5 * (yb + ye),
+            "patch_icpp(1)%length_x": xe - xb,
+            "patch_icpp(1)%length_y": ye - yb,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": p_amb,
+            "patch_icpp(1)%alpha_rho(1)": rho_amb,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Patch 2: post-shock state on the left (x < -5 mm); the Dirichlet inflow holds it
+            "patch_icpp(2)%geometry": 3,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%x_centroid": -0.1025,
+            "patch_icpp(2)%y_centroid": 0.5 * (yb + ye),
+            "patch_icpp(2)%length_x": 0.195,
+            "patch_icpp(2)%length_y": ye - yb,
+            "patch_icpp(2)%vel(1)": u_post,
+            "patch_icpp(2)%vel(2)": 0.0,
+            "patch_icpp(2)%pres": p_post,
+            "patch_icpp(2)%alpha_rho(1)": rho_post,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            # Fluid: air
+            "fluid_pp(1)%eos": "ideal_gas",
+            "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(1)%cv": 717.5,
+            # Lagrangian particles
+            "particles_lagrange": "T",
+            "fd_order": 2,
+            "particle_pp%rho0ref_particle": 2520.0,
+            "particle_params%input_path": "input/particles.dat",
+            "particle_params%nparticles_glb": gen_particles.n_particles,
+            "particle_params%solver_approach": 2,
+            "particle_params%qs_force": 3,
+            "particle_params%pressure_gradient_force": "T",
+            "particle_params%added_mass_force": 1,
+            "particle_params%mu_ref(1)": 1.716e-5,
+            "particle_params%suth(1)": 110.4,
+            "particle_params%interpolation_order": 2,
+            "particle_params%epsilonb": 1.0,
+            "particle_params%valmaxvoid": 0.9,
+            "particle_params%charwidth": gen_particles.charwidth,
+        },
+        indent=4,
+    )
+)
@@ -0,0 +1,31 @@
+#!/usr/bin/env python3
+# Writes <outdir>/input/particles.dat for the particle curtain (default outdir: this example's directory).
+# case.py calls it when pre_process runs; the test suite calls it for the Example test.
+import math
+import os
+import random
+import sys
+
+# Glass particles, radius 57.5 um, volume fraction 0.21 in the band x in [0, 2 mm] across the channel y in [0, 15 mm].
+# In 2D each particle stands for a slab of depth charwidth, so projected particles may overlap (collisions are not modeled).
+rp = 57.5e-6
+x_curtain = (0.0, 2.0e-3)
+y_channel = (0.0, 0.015)
+vf = 0.21
+charwidth = 2.5e-4  # the grid spacing
+n_particles = round(vf * (x_curtain[1] - x_curtain[0]) * (y_channel[1] - y_channel[0]) * charwidth / (4.0 / 3.0 * math.pi * rp**3))
+
+
+def write_particles(outdir):
+    """Uniform random (seeded) positions in the curtain band, one line per particle: x, y, z, u, v, w, radius."""
+    random.seed(1)
+    os.makedirs(os.path.join(outdir, "input"), exist_ok=True)
+    with open(os.path.join(outdir, "input", "particles.dat"), "w") as f:
+        for _ in range(n_particles):
+            x = random.uniform(*x_curtain)
+            y = random.uniform(y_channel[0] + rp, y_channel[1] - rp)
+            f.write(f"{x:.16e} {y:.16e} 0.0 0.0 0.0 0.0 {rp:.16e}\n")
+
+
+if __name__ == "__main__":
+    write_particles(sys.argv[1] if len(sys.argv) > 1 else os.path.dirname(os.path.abspath(__file__)))
@@ -0,0 +1,15 @@
+# 2D Blast Wave Through a Particle Half-Ring
+
+A Mach 10 cylindrical blast wave, released from a high-pressure driver at the origin, hits a half-ring of 1612 magnesium particles (radius 50 µm) in air. It uses the Euler-Lagrange solid-particle solver (`particles_lagrange`) with two-way coupling, Osnes quasi-steady drag, pressure-gradient force and added mass.
+
+`gen_particles.py` writes the particles (seeded random placement without overlaps, radius 11.05–12.65 mm) to `input/particles.dat`; `case.py` calls it when pre_process runs, or run it yourself with `python3 gen_particles.py`.
+
+```shell
+./mfc.sh run examples/2D_particle_hemisphere/case.py -n 4
+```
+
+## Result
+
+Density-gradient schlieren |∇ρ|/ρ with the particles (white) at t = 50 µs.
+
+<img src="result.png"/>