Devlog #85 – Wave 64: Drag Coefficient, Geodesic Dome, Stable Matching, Hologram, Louvain & Grain Structure

New Simulations

🪂 Drag Coefficient — Aerodynamics Across Shape and Reynolds Number

The drag force equation

The drag force on any body moving through a fluid is described by a single dimensionless coefficient C_d that absorbs all the geometric complexity of flow separation, wake formation, and pressure recovery:

F_d = ½ · ρ · v² · C_d · A

where:
  ρ  = fluid density (kg/m³)
  v  = relative velocity (m/s)
  A  = reference area (frontal area for bluff bodies, planform for airfoils)
  C_d = drag coefficient (dimensionless)

Decomposition:
  C_d = C_d,pressure + C_d,friction
  Pressure (form) drag dominates for bluff bodies
  Friction (viscous) drag dominates for streamlined bodies and flat plates

The simulation evaluates F_d in real time as you adjust velocity and fluid density sliders, with the formula displayed numerically. Six body shapes are available, each with its characteristic C_d(Re) curve derived from experimental data.

Shape comparison and C_d values

At intermediate Reynolds numbers (Re ≈ 10⁴–10⁵), the six bodies span nearly two orders of magnitude in drag coefficient:

Body shape         C_d (Re ~ 10⁵)    Notes
─────────────────────────────────────────────────────────
Flat plate         ~1.17             Pure pressure drag, massive wake
Cube               ~1.05             Separation at sharp corners
Cylinder           ~1.0              Wide wake, separation at ~80°
Sphere             ~0.47             Subcritical (laminar BL)
Streamlined body   ~0.04–0.10        Attached flow, thin wake
Airfoil (NACA 0012) ~0.006–0.012     Mostly friction drag at low AoA

The multi-shape overlay panel plots all six C_d(Re) curves simultaneously on a log-log axis, making it easy to see where curves converge, diverge, and where the drag crisis occurs.

The drag crisis

The most striking feature of the drag coefficient simulation is the drag crisis for a sphere. At subcritical Reynolds numbers (Re < 3×10⁵), the boundary layer is laminar and separates early, producing a wide wake and C_d ≈ 0.47. As Re crosses the critical value, the boundary layer transitions to turbulent before it reaches the separation point. Turbulent boundary layers have higher momentum and can remain attached further around the sphere, dramatically narrowing the wake:

Subcritical (Re < 3×10⁵):
  Laminar boundary layer → early separation (~80° from front stagnation)
  Wide wake, large pressure drag → C_d ≈ 0.47

Critical transition (Re ≈ 3×10⁵):
  BL transition to turbulent → delayed separation (~120° from front)
  Narrow wake, small pressure drag → C_d ≈ 0.10

Effect on sports:
  Golf ball dimples: artificially trip BL turbulent at lower Re
  → drag crisis at Re ≈ 4×10⁴ instead of 3×10⁵
  → 2× longer range vs smooth ball at typical golf speeds

The flow-separation visualisation draws the separation point as an animated marker on the body outline and sketches the wake width as a shaded region, giving an intuitive sense of what the C_d number represents physically.

🔺 Geodesic Dome — From Icosahedron to Spherical Triangulation

The icosahedron as starting point

A regular icosahedron has 20 equilateral triangular faces, 30 edges, and 12 vertices, all on a circumscribed sphere. It is the Platonic solid with the most faces, making it the natural starting point for a sphere approximation. Its vertex coordinates are:

Golden ratio φ = (1 + √5) / 2 ≈ 1.6180

12 vertices (normalised to unit sphere):
  (0, ±1, ±φ)  →  normalised: (0, ±1/√(1+φ²), ±φ/√(1+φ²))
  (±1, ±φ, 0)  →  normalised similarly
  (±φ, 0, ±1)  →  normalised similarly

Face count:  F = 20
Edge count:  E = 30
Vertex count: V = 12
Euler: V - E + F = 12 - 30 + 20 = 2  ✓

Frequency-ν subdivision

To subdivide each triangular face at frequency ν, each edge of the original icosahedron is divided into ν equal parts, creating a grid of ν² smaller triangles per face. The new interior vertices are placed by linear interpolation and then projected outward onto the circumscribed sphere:

For each original triangle (A, B, C) at frequency ν:
  For i = 0 to ν:
    For j = 0 to ν - i:
      k = ν - i - j
      Point P = (i·A + j·B + k·C) / ν    // barycentric interpolation
      P_sphere = P / |P|                   // project to unit sphere

Result counts:
  New vertices per face:  (ν+1)(ν+2)/2
  New triangles per face: ν²
  Total vertices: V = 10ν² + 2
  Total faces:    F = 20ν²
  Total edges:    E = 30ν²
  Euler check:    V - E + F = 10ν² + 2 - 30ν² + 20ν² = 2  ✓ for all ν

The live Euler formula display in the simulation updates as you drag the frequency slider, confirming V−E+F = 2 at every step — a satisfying check that the topology is consistent.

Dome mode and structural properties

Dome mode clips the subdivided sphere at the equator, retaining only the lower hemisphere. This halves the vertex count roughly, but the boundary ring requires reinforcing struts. The simulation calculates the maximum strut length variance — geodesic domes built from a single frequency have all struts nearly equal in length, which simplifies construction. At ν = 1 there is only one strut length; by ν = 6 there are still only 9 distinct lengths for the entire dome structure.

💍 Stable Matching — Gale-Shapley and the Deferred Acceptance Algorithm

The stable matching problem

Given N proposers and N acceptors, each with a complete strict preference ranking over the other side, a stable matching is a perfect matching (every person matched exactly once) with no blocking pair: a pair (p, a) where proposer p prefers acceptor a to their current match and acceptor a prefers p to their current match. An unstable matching would "unravel" in practice — the blocking pair would defect.

Formal definition:
  Proposers P = {p₁, …, pN}
  Acceptors A = {a₁, …, aN}
  Each pᵢ has ranking: aσ(i,1) > aσ(i,2) > … > aσ(i,N)
  Each aⱼ has ranking: pτ(j,1) > pτ(j,2) > … > pτ(j,N)

Blocking pair: (pᵢ, aⱼ) such that
  pᵢ prefers aⱼ to μ(pᵢ)   AND
  aⱼ prefers pᵢ to μ⁻¹(aⱼ)
where μ is the current matching.

Gale-Shapley deferred acceptance

The algorithm (Gale and Shapley, 1962) runs in rounds. Each unmatched proposer proposes to their most-preferred acceptor not yet rejected. Each acceptor tentatively holds the best proposal received so far and rejects all others. The key insight is that proposals only go down the proposer's preference list, while acceptors can only trade up — so the algorithm always terminates and produces a stable matching:

Gale-Shapley (proposer-optimal):
  while ∃ unmatched proposer pᵢ with non-empty list:
    aⱼ = pᵢ's top remaining acceptor
    if aⱼ is unmatched:
      tentatively match (pᵢ, aⱼ)
    elif aⱼ prefers pᵢ to current holder pₖ:
      unmatch pₖ (pₖ becomes free again)
      tentatively match (pᵢ, aⱼ)
    else:
      aⱼ rejects pᵢ; pᵢ removes aⱼ from list

Complexity: O(N²) proposals in the worst case
Terminates in at most N² rounds

The simulation's step-by-step mode highlights the current proposer in yellow, the target acceptor in blue, and rejected proposals in red. A proposal-round counter shows how many rounds were needed, which varies between N (best case) and N² (worst case).

Optimality and blocking-pair verification

Gale-Shapley produces the proposer-optimal stable matching: every proposer gets the best partner they could receive in any stable matching. Acceptors receive the worst stable partner. The simulation's block-pair detector exhaustively checks all N² pairs after convergence, confirming zero blocking pairs exist. The 2000-run verification test randomises all preference lists and checks stability on each — confirming correctness even for adversarial preference structures.

🎭 Hologram — Interference Recording and Wavefront Reconstruction

Recording: interference on the holographic plate

Unlike a photograph, which records only intensity, a hologram records both the amplitude and phase of the object wavefront by interfering it with a coherent reference beam. The intensity pattern on the holographic plate is:

Recording geometry:
  E_ref(x) = A_r · exp(i k x sin θ_r)     (plane reference wave)
  E_obj(x) = Σⱼ aⱼ / rⱼ · exp(i k rⱼ)    (spherical waves from object points)

Recorded intensity (hologram transmittance ∝ exposure):
  I(x) = |E_ref + E_obj|²
       = |E_ref|² + |E_obj|² + E_ref* E_obj + E_ref E_obj*
         ───────────────────────────────────────────────────
         DC bias    speckle    signal term   conjugate term

The signal term E_ref* · E_obj encodes the object's phase map.

The simulation renders I(x) as a grey-scale fringe pattern on the holographic plate in real time as you drag object point sources around the scene. Multiple object points produce complex overlapping fringe sets that together encode the 3D object.

Reconstruction: diffraction orders

Illuminating the processed hologram with the original reference beam diffraction-reads the stored wavefront. The transmitted field decomposes into three orders:

Hologram transmittance t(x) ∝ I(x)   (linear recording regime)

Reconstructed field = t(x) · E_ref(x)
  = (|E_ref|² + |E_obj|²) · E_ref    → zero order (DC, straight-through)
  +  E_ref* · E_obj · E_ref           → +1 order: E_obj reconstructed  (virtual image)
  +  E_ref · E_obj* · E_ref           → −1 order: E_obj* conjugate wave (real image)

+1 order: diverging wave appears to come from original object positions → virtual image
−1 order: converging wave focuses to a real image (pseudoscopic, depth-reversed)

The simulation animates the reconstruction by propagating waves from the hologram plane. The +1 order virtual image appears behind the plate in the correct 3D position; the −1 order real image forms in front. A slider adjusts the reconstruction wavelength, showing how wavelength mismatch with the recording wavelength causes magnification and chromatic distortion.

🔵 Louvain — Community Detection and Modularity

Modularity as an objective

The modularity Q of a graph partition measures how much more within-community edges there are compared to a random null model with the same degree sequence:

Q = (1/2m) · Σᵢⱼ [Aᵢⱼ - kᵢkⱼ/(2m)] · δ(cᵢ, cⱼ)

where:
  m     = total number of edges
  Aᵢⱼ  = adjacency matrix entry
  kᵢ   = degree of node i
  cᵢ   = community label of node i
  δ(a,b) = 1 if a = b, 0 otherwise

Q range: [−0.5, 1)
  Q ~ 0: no community structure (random graph)
  Q > 0.3: significant community structure
  Q > 0.7: very strong communities

The two-phase Louvain algorithm

The Louvain algorithm (Blondel et al., 2008) alternates two phases until no further improvement is possible:

Phase 1 — Local Optimisation:
  Initialise: each node in its own community
  For each node i (random order):
    Compute ΔQ for moving i to each neighbour community c:
      ΔQ = [Σ_in + 2k_i,in) / (2m) - ((Σ_tot + kᵢ) / (2m))²]
           − [Σ_in / (2m) − (Σ_tot / (2m))² − (kᵢ / (2m))²]
    where Σ_in  = edges within c
          Σ_tot = total degree of nodes in c
          k_i,in = edges from i to c
    Move i to community with max ΔQ > 0
  Repeat until no move increases Q

Phase 2 — Coarsening (super-node aggregation):
  Build new graph: nodes = communities from Phase 1
  Edge weights = sum of edges between communities
  Self-loops = sum of edges within each community
  Repeat Phase 1 on the coarsened graph

Complexity: O(n log n) in practice on sparse graphs
Terminates when Q stops increasing between full passes

Graph generation and verification

The simulation uses a Stochastic Block Model to generate test graphs with known ground-truth community structure. The SBM places n_communities groups of equal size; within-group edge probability p_in exceeds between-group probability p_out, creating the planted partition. After Louvain converges, the simulation computes the Normalised Mutual Information (NMI) between the detected communities and the ground truth, and displays the final Q value. For well-separated communities (p_in/p_out > 4) Louvain reliably recovers the planted structure with NMI > 0.95.

🔩 Grain Structure — Solidification, Coarsening, and the Hall-Petch Relation

Voronoi nucleation

Solidification begins with nucleation: crystalline seeds form at random positions in the melt. Each grain grows outward until it meets a neighbouring grain, forming a grain boundary. The resulting tessellation is a Voronoi diagram — each grain cell contains all points closer to its nucleus than to any other:

Voronoi cell for nucleus i:
  V(i) = { x ∈ ℝ² : |x − nᵢ| ≤ |x − nⱼ| for all j ≠ i }

Grain boundary between cells i and j:
  The perpendicular bisector of segment nᵢnⱼ
  meets at a triple junction (three grains) with angle ~120°
  (minimises grain boundary energy per unit length)

The simulation nucleates grains uniformly at random and renders the Voronoi tessellation with each grain assigned a random crystallographic orientation (colour-coded by hue).

Potts model annealing and grain growth

Grain growth is driven by curvature: grain boundaries migrate toward their centre of curvature, reducing total boundary length and energy. The Potts model discretises this on a grid: each lattice site carries a spin (grain ID); the Monte Carlo Metropolis step flips boundary sites to the orientation of a neighbouring grain if it reduces total boundary energy:

Potts model energy:
  H = J · Σ_{⟨i,j⟩} (1 − δ(sᵢ, sⱼ))

  J = grain boundary energy per bond
  Sum over nearest-neighbour pairs only
  δ(a,b) = 1 within same grain, 0 across boundary

Metropolis update at site i:
  Choose random neighbour grain orientation s'
  ΔH = H(s') − H(sᵢ)
  Accept if ΔH ≤ 0, or with probability exp(−ΔH / kT)

Grain growth law (parabolic growth):
  ⟨D⟩² − ⟨D₀⟩² = K · t
  ⟨D⟩ ∝ tⁿ  with n ≈ 0.5 for ideal normal grain growth

The Hall-Petch relation

The yield strength of a polycrystalline metal increases as grain size decreases, because grain boundaries act as obstacles to dislocation motion. The Hall-Petch relation quantifies this:

Hall-Petch relation:
  σ_y = σ₀ + k / √d

where:
  σ_y = yield strength (MPa)
  σ₀  = lattice friction stress (yield strength of a single crystal)
  k   = Hall-Petch slope (strengthening coefficient, MPa·m^½)
  d   = mean grain diameter (m)

Example values (steel):
  σ₀ ≈ 70 MPa
  k  ≈ 0.74 MPa·m^½
  At d = 10 μm: σ_y ≈ 70 + 0.74/√(10×10⁻⁶) ≈ 70 + 234 = 304 MPa
  At d = 1  μm: σ_y ≈ 70 + 740 = 810 MPa

The simulation tracks mean grain diameter D(t) as annealing proceeds and plots σ_y vs D^−½ in real time alongside the grain map. The four material presets (pure iron, aluminium alloy, copper, austenitic steel) supply realistic σ₀ and k values so the Hall-Petch plot has physically accurate units.

Up Next

Wave 65 is in development with six more simulations covering capillary action (Jurin's law and meniscus shape), A* pathfinding with three heuristics, parametric curves from de Casteljau Bézier through B-splines to NURBS, Laguerre-Gaussian vortex beams with orbital angular momentum, colloidal stability via DLVO theory, and internet routing protocols including Bellman-Ford and Dijkstra link-state. The library is on track to exceed 600 simulations in June.