Add packed bed thermal conductivity benchmark

Cargo.toml

@@ -126,6 +126,10 @@ path = "examples/example_template/main.rs"
 name = "cone_hopper"
 path = "examples/cone_hopper/main.rs"
 
+[[example]]
+name = "bench_thermal_bed"
+path = "examples/bench_thermal_bed/main.rs"
+
 [profile.release]
 opt-level = 3
 lto = "fat"

examples/bench_thermal_bed/README.md

@@ -0,0 +1,66 @@
+# Packed Bed Thermal Conductivity Benchmark
+
+## Physics
+
+This benchmark validates DEM contact-based heat conduction through a packed bed of monodisperse spheres. Two flat walls at different temperatures (hot bottom at 400 K, cold top at 300 K) drive heat flow through particle-particle and particle-wall contacts in vacuum (no fluid phase).
+
+The heat transfer model uses Hertzian contact mechanics:
+
+```
+Q = k · 2a · (T_j - T_i)
+```
+
+where `a = sqrt(R_eff · δ)` is the contact radius, `k` is the thermal conductivity, and `δ` is the overlap.
+
+## Validation
+
+The benchmark checks against the Zehner-Bauer-Schlunder (ZBS) framework for effective thermal conductivity in packed beds:
+
+1. **Linear temperature profile** — At steady state with constant effective conductivity, Fourier's law predicts a linear temperature distribution (R² > 0.8).
+2. **Steady-state convergence** — Wall heat flux variation < 10% in the last 20% of the simulation.
+3. **Reasonable k_eff** — The effective thermal conductivity ratio k_eff/k_s falls in the physically expected range (0.001 to 1.0) for contact-only conduction in vacuum.
+4. **Temperature bounds** — All particle temperatures stay between T_cold and T_hot.
+5. **Average temperature** — The bed average temperature is near the midpoint (350 K).
+
+## Setup
+
+- **Particles**: 200 monodisperse steel spheres (R = 1 mm, ρ = 7800 kg/m³)
+- **Domain**: 20 × 20 mm periodic in x,y; walls at z = 0 (bottom) and z = 22 mm (top)
+- **Material**: k = 50 W/(m·K), cₚ = 500 J/(kg·K), E = 10 MPa (soft for fast packing)
+- **Stage 1**: FIRE energy minimization to create a dense random packing under gravity
+- **Stage 2**: Thermal conduction with wall temperatures T_hot = 400 K, T_cold = 300 K
+
+## Running
+
+```bash
+# Run simulation (release mode recommended, ~2-3 minutes)
+cargo run --release --example bench_thermal_bed -- examples/bench_thermal_bed/config.toml
+
+# Validate results
+cd examples/bench_thermal_bed
+python3 validate.py
+
+# Generate plots
+python3 plot.py
+```
+
+## Expected Output
+
+- `data/ThermalProfile.csv` — Per-atom z-position and temperature at steady state
+- `data/WallHeatFlux.txt` — Time series of bottom wall heat flux
+- `temperature_profile.png` — Steady-state T(z) with linear fit
+- `wall_heat_flux.png` — Heat flux convergence to steady state
+- `avg_temperature.png` — Average temperature evolution
+
+## Plots
+
+![Temperature Profile](temperature_profile.png)
+![Wall Heat Flux](wall_heat_flux.png)
+![Average Temperature](avg_temperature.png)
+
+## References
+
+- Zehner, P. & Schlunder, E.U. (1970). *Chemie Ingenieur Technik*, 42(14), 933-941.
+- Bauer, R. & Schlunder, E.U. (1978). *Int. Chem. Eng.*, 18(2), 189-204.
+- Batchelor, G.K. & O'Brien, R.W. (1977). *Proc. R. Soc. Lond. A*, 355, 313-333.
+- Vargas, W.L. & McCarthy, J.J. (2001). *AIChE Journal*, 47(5), 1052-1059.

examples/bench_thermal_bed/config.toml

@@ -0,0 +1,107 @@
+# ============================================================================
+# Packed Bed Thermal Conductivity Benchmark
+# ============================================================================
+# Validates DEM heat conduction through a packed bed against expected physics.
+#
+# Stage 1 (settling): Particles fall under gravity onto the bottom wall.
+# Stage 2 (thermal):  Top wall moved to bed surface, temperatures set,
+#                      heat conducts to steady state.
+# ============================================================================
+
+[comm]
+processors_x = 1
+processors_y = 1
+processors_z = 1
+
+# Domain: periodic in x,y; fixed walls in z
+[domain]
+x_low  = 0.0
+x_high = 0.012         # 12 mm [m]
+y_low  = 0.0
+y_high = 0.012         # 12 mm [m]
+z_low  = -0.001        # Below bottom wall [m]
+z_high = 0.040         # Well above settling region [m]
+periodic_x = true
+periodic_y = true
+periodic_z = false
+
+[neighbor]
+skin_fraction = 1.1    # Neighbor list skin multiplier [dimensionless]
+bin_size = 0.005       # Bin width [m]
+
+# Gravity for settling
+[gravity]
+gz = -9.81             # [m/s²]
+
+# Material properties
+[[dem.materials]]
+name = "steel"
+youngs_mod = 2e8       # Young's modulus [Pa]
+poisson_ratio = 0.3    # Poisson's ratio [dimensionless]
+restitution = 0.2      # Low COR for fast energy dissipation
+friction = 0.5         # Sliding friction [dimensionless]
+rolling_friction = 0.1 # Rolling friction [dimensionless]
+
+# 100 spheres inserted in a tall region — they settle under gravity
+[[particles.insert]]
+material = "steel"
+count = 100
+radius = 0.001         # Particle radius [m] (d = 2 mm)
+density = 7800.0       # Steel density [kg/m³]
+region = { type = "block", min = [0.001, 0.001, 0.002], max = [0.011, 0.011, 0.014] }
+
+# Thermal properties — high k and low cp for fast equilibration in benchmark
+[thermal]
+conductivity = 500.0      # Thermal conductivity [W/(m·K)] — enhanced for fast SS
+specific_heat = 50.0      # Specific heat capacity [J/(kg·K)] — reduced for speed
+initial_temperature = 350.0  # Initial temperature [K]
+
+# ── Walls ──────────────────────────────────────────────────────────────────
+# Bottom wall at z=0 (normal up)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 0.0
+normal_y = 0.0
+normal_z = 1.0
+material = "steel"
+name = "bottom"
+
+# Top wall — placed very high during settling, moved down at thermal stage start
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.035
+normal_x = 0.0
+normal_y = 0.0
+normal_z = -1.0
+material = "steel"
+name = "top"
+
+# Output
+[output]
+dir = "examples/bench_thermal_bed"
+
+[thermo]
+columns = ["step", "atoms", "ke", "walltime"]
+
+[dump]
+interval = 200000      # Write dump every 200k steps
+format = "text"
+
+# ── Stage 1: Gravity settling ─────────────────────────────────────────────
+[[run]]
+name = "settling"
+steps = 500000         # 0.25 s at dt=5e-7 — enough for 100 particles to settle
+dt = 5e-7              # Timestep [s]
+thermo = 50000
+
+# ── Stage 2: Thermal conduction ───────────────────────────────────────────
+# Top wall moved to bed surface, temperatures set, gravity disabled
+[[run]]
+name = "thermal"
+steps = 4000000        # 2.0 s for thermal equilibration
+dt = 5e-7              # Timestep [s]
+thermo = 20000         # Output every 0.01 s
+gravity.gz = 0.0       # No gravity during thermal stage