2D Urban Pollutant Dispersion (simpleFoam → scalarTransportFoam)

A simplified 2D case for passive-pollutant (NOx-like tracer) dispersion around three buildings in a left→right atmospheric wind. Derived from the layout of $FOAM_TUTORIALS/basic/scalarTransportFoam/pitzDaily: the backward-facing-step geometry is discarded, the case layout + scalarTransportFoam solver are kept, and the scalar T is reinterpreted as pollutant concentration.

Built and verified against OpenFOAM v2512 (ESI/OpenCFD).

Coordinate convention: x streamwise (left→right), y vertical (buildings rise in +y), z thin spanwise direction (2D, empty). All lengths in metres.

Two-stage workflow

scalarTransportFoam advects T in a frozen, prescribed U/phi — it does not solve momentum. So the flow is produced first by a separate steady RANS run:

1.  cd urbanFlow      && ./Allrun     # blockMesh → checkMesh → simpleFoam (k-ε, ABL inlet)
2.  cd urbanDispersion && ./Allrun     # copies mesh + frozen U/phi/nut, runs scalarTransportFoam

urbanDispersion/Allrun copies the converged mesh and the latest-time U, phi (used as-is, not re-derived) and nut from urbanFlow/. Two separate case directories keep the flow solve and scalar solve cleanly decoupled.

./Allclean in each directory resets it (urbanDispersion/Allclean keeps 0/T + system/ and drops the copied mesh/fields).

Run with Snakemake (alternative driver)

The repo-root Snakefile reproduces both Allruns as a single dependency-aware workflow (makeBlockMeshDict → blockMesh → checkMesh → simpleFoam → copy fields → scalarTransportFoam) and adds three clean jobs. Editing a boundary condition or scheme makes Snakemake re-run only the affected solve. Activate both environments first, from the repo root:

conda activate base                  # provides the snakemake binary
source $WM_PROJECT_DIR/etc/bashrc    # OpenFOAM v2512 (or the `of2512` alias)

snakemake -n                 # dry run: print the plan, change nothing
snakemake --cores 8          # build the whole pipeline; solvers run on 8 cores
snakemake clean              # reset BOTH cases   (= both Allcleans)
snakemake clean_urbanFlow            # reset the flow case only
snakemake clean_urbanDispersion      # reset the dispersion case only

Parallel CFD. --cores N drives the decomposition: only the two solver rules run in parallel, decomposing each case with scotch into N subdomains (numberOfSubdomains in system/decomposeParDict is set to N at run time), running mpirun -np N <solver> -parallel, then reconstructPar. Every other job (mesh, checkMesh, field copy, clean) stays serial. With --cores 1 the solvers run serially with no decomposition.

Mesh refinement (Richardson / GCI). --config refinement_level=L with L ∈ {1,2,3} scales the blockMesh cell counts by √2^(L-1) (gradings unchanged) → a self-similar triplet with constant ratio r=√2: R1=45 520, R2=90 986, R3=182 080 cells (cell size halves R1→R3; effective r≈1.414). Default is 1. Measured solver wall times on 8 cores: R1 simpleFoam ~30 s, R2 ~90 s, R3 ~5.4 min (so the full R3 stage incl. scalarTransport is ~6–7 min — the equal-ratio √2 triplet is kept deliberately, exceeding a 5-min target). Example convergence study:

for L in 1 2 3; do
  snakemake clean
  snakemake --cores 8 --config refinement_level=$L
  postProcess -case urbanDispersion -func airPolution -latestTime  # record the QoI
done

Solver completion is tracked by the urbanFlow.foam / urbanDispersion.foam markers (the solvers' time-directory names are not known in advance). The original Allrun/Allclean scripts are unchanged and still work (serially).

Geometry & mesh

Domain 500 × 300 × 1 m (one cell in z). Four buildings on the floor, 20 m wide, heights 20 / 16 / 14 / 12 m, separated by three 20 m canyons. Domain height = 15× tallest building (blockage 6.7 %); wake = 13 H = 260 m downstream (the 4th building sits in what was the wake).

The mesh is generated programmatically from a column/row occupancy grid by urbanFlow/tools/makeBlockMeshDict.py (re-run automatically by urbanFlow/Allrun if python3 is present). It emits 120 vertices, 35 fluid blocks (10 building-interior blocks omitted), ≈ 45 500 cells, and derives every patch face automatically:

patch	type	role
inlet	patch	west face, x = 0 (ABL profile)
outlet	patch	east face, x = 500
top	symmetryPlane	y = 300 free-slip lid
ground	wall	floor outside the canyons (rough)
emission1 / emission2 / emission3	wall	the three canyon floors (rough + independent source)
buildings	wall	building windward/leeward/roof faces (smooth)
frontAndBack	empty	spanwise faces (2D)

The y-stations are 0, 12, 14, 16, 20, 300 so the four building heights coincide exactly with mesh row-tops (no partial rows). checkMesh is clean (orthogonal, max skewness ≈ 0, total fluid volume 148 760 m³ = 150 000 − 1 240 building volume). Ground-cell centre height ≈ 0.25 m > z0 = 0.1 m, as required for the rough-wall function.

Atmospheric inlet (log-law)

Shared parameters live in urbanFlow/0/include/ABLConditions and are included verbatim by 0/U, 0/k, 0/epsilon:

Uref 5 m/s @ Zref 10 m,  z0 0.1 m,  d 0,  flowDir (1 0 0),  zDir (0 1 0)

⚠️ zDir is the vertical = +y here, despite the historical name.

Derived (and verified by sampling): u* ≈ 0.44 m/s, k ≈ 0.65 m²/s², U(10 m) = 5.0, U(20 m) = 5.7, U(300 m) = 8.7 m/s. Turbulence: standard kEpsilon (RAS). Ground + canyon floors use the rough wall function atmNutkWallFunction (z0 = 0.1, the v2512 name for the former nutkAtmRoughWallFunction); buildings use smooth nutkWallFunction.

Emission sources

Each canyon floor (emission1/2/3) is a fixedGradient on T: at a no-slip wall pollutant enters only by diffusion, so gradient = F / DT (DT = 1 m²/s). The three canyons are independent so each rate is tuned separately — current values span the busy-street NOx range:

patch	canyon	gradient	note
emission1	B1↔B2	1×10⁻⁸	baseline busy street
emission2	B2↔B3	2×10⁻⁸	heavier traffic
emission3	B3↔B4	5×10⁻⁹	lighter traffic

The scalar is passive and linear, so T ∝ F — rescale afterwards for a specific pollutant. DT is a crude constant stand-in for νt/Sct.

Solver settings of note

• Stage 1: SIMPLEC (consistent yes), residualControl 1e-4, bounded Gauss linearUpwind for U. Converged in 1497 iterations.
• Stage 2: ddt steadyState, div(phi,T) bounded Gauss upwind (unconditionally bounded and monotonically convergent so the global balance closes; higher-order limitedLinear/linearUpwind are noted in fvSchemes but their nonlinear limiter stalls the steady residual in a ~1e-5 limit-cycle). T under-relaxed 0.7. Converged in 1696 iterations.
• Both load libs (atmosphericModels) — the copied 0/U carries the atmBoundaryLayerInletVelocity patch field, whose constructor lives there.

Validation (all pass)

check	result
checkMesh	clean; 35 blocks / 120 verts / 45 520 cells; ground cell-centre 0.25 m > z0
Stage-1 convergence	residualControl met @ 1497 it; continuity err ~1e-11
ABL homogeneity	inlet profile matches analytic log-law; preserved aloft to <1 % over the fetch (near-ground deceleration toward the buildings is physical 2D blockage)
canyon / wake recirculation	reversed flow in all three canyons (Ux to −1.4 m/s) and behind B4
Stage-2 boundedness	T ∈ [0, 2.42×10⁻⁷]; no negatives (max at the emission2 floor)
global balance	emission source 7.0×10⁻⁷ ≈ outlet outflux 6.991×10⁻⁷ (0.13 %)
plume	trapped in canyons, advected and diluted downstream

Post-processing helpers: urbanFlow/system/{sampleLines,canyonU}, urbanDispersion/system/{outletFlux,inletFlux,plume,airPolution}. airPolution reports the area-averaged breathing-level concentration (areaAverage of T on a surface 1.5 m above the floor patches emission1/2/3 + ground, via a patchInternalField sampled surface) — it is the objective of the Bayesian optimization below. Run e.g. postProcess -case urbanFlow -func sampleLines -latestTime, or for the balance postProcess -case urbanDispersion -fields '(T phi)' -func outletFlux -latestTime.

Bayesian optimization of the emission split (optimization/)

optimization/ finds how to split a fixed total traffic emission across the three street patches (emission1/2/3) to minimize the worst-case breathing-level pollution over a set of wind speeds (robust / minimax), using BoTorch (sequential, q=1) + Snakemake + OpenFOAM:

• optimization/config.yaml — winds [2,4,6,8] m/s, fixed refinement_level, total emission Φ, aggregation: max, BO budget.
• optimization/Snakefile — flow (one simpleFoam per wind, cached), disperse (per candidate × wind: frozen flow + per-patch gradients → scalarTransportFoam → airPolution), aggregate (max over winds → objective.json). Flow caching makes each candidate cheap.
• optimization/optimize.py — Sobol hot-start then SingleTaskGP + LogExpectedImprovement; the fixed-sum constraint is handled by mapping the box [ε,1]^N to traffic fractions w = x/Σx (the ε floor keeps Σx>0 so the total Φ is conserved — the origin would give zero emission, a spurious zero-pollution optimum). The fitted GP is saved to gp_model.pt after every iteration; a re-launch loads it and skips Sobol. Results logged incrementally to bayesian-optimization.csv.
• optimization/visualize.py — GP posterior mean + std over the traffic-fraction simplex, the evaluated points, and the convergence curve → bo_model.png.

Exact precompute fast-path (default). Because the transport is linear in the per-patch flux (frozen flow, fixedGradient sources), one solve per (patch, wind) with a unit flux gives the full response matrix R[U][j] (response_matrix.json, N·M = 12 solves, rule response_matrix). Any candidate is then A_U = Σ_j R[U][j]·F_j — a dot product, no CFD in the BO loop. optimize.py builds R once and evaluates analytically (BO_EVAL=cfd falls back to per-candidate Snakemake). The full 10+30 campaign runs in ~7 s (vs ~23 min CFD-in-loop); R@F matches the CFD objective to 0.3 %.

conda activate base; source $WM_PROJECT_DIR/etc/bashrc
cd optimization
python optimize.py          # builds R (12 solves) if needed, then BO in seconds
python visualize.py         # bo_model.png

max (worst-case over winds) is used so the objective is robust. Here R shows the downstream canyon emission3 has the lowest response at every wind, so the robust optimum is w ≈ [0,0,1] (all traffic on emission3). A genuine interior trade-off would require the per-patch ranking to flip across winds.

Limitations

Passive, neutral tracer only (no buoyancy/stability). Constant DT under-mixes — for fidelity build DT = νt/Sct from the copied nut (needs a minor solver tweak). With a symmetryPlane top the ABL can still mildly drift; impose the ABL value at the top for strict homogeneity. Blockage 6.7 % (above the <3 % ideal). See the original instruction for the full discussion.