Numerics behind the fluid demos

The Rayleigh–Bénard and urban-flow demos on this site share a single pure-JavaScript NS2D solver. The page below documents what the solver does, why it is a good fit for a browser demo, and — more importantly — why it is not appropriate for research-grade simulation.

Governing equations

The shared NS2D code solves the two-dimensional incompressible Navier-Stokes equations. In dimensional form the velocity u = (u, w) obeys

∂_tu + (u · ∇)u = -∇p + ν∇²u + βgTe_z, ∇ · u = 0,

∂_tT + u · ∇T = κ∇²T.

The buoyancy term is active in the Rayleigh–Bénard demo, where the top and bottom plates set the thermal boundary condition. In the urban-flow demo the temperature equation is passive and the same incompressible momentum solver is driven by inflow/outflow boundary conditions and rectangular solid masks.

The symbols are the velocity u = (u, w), pressure divided by a constant reference density p, temperature T, kinematic viscosity ν, thermal diffusivity κ, buoyancy coefficient βg, and the upward unit vector e_z. The operator ∇· is divergence, ∇ is the gradient, and ∇² is the Laplacian. In the demos βg is only used for Rayleigh–Bénard convection; in the urban-flow demo the temperature field is not part of the visible dynamics.

Numerical discretisation

The browser solver follows the classical finite-volume/finite- difference projection-method ideas described in Ferziger and Peric. It uses a staggered MAC-type grid: horizontal velocity u is stored on vertical cell faces, vertical velocity w on horizontal cell faces, and pressure and temperature at cell centres with one ghost-cell layer. This arrangement makes the discrete pressure gradient and velocity divergence naturally compatible, which is the main reason staggered grids are so common for incompressible flow.

Each time step first transports velocity and temperature using semi-Lagrangian back-tracing with bilinear interpolation. That is more dissipative than a high-order research discretisation, but it is robust for an interactive web page. Diffusive terms use the standard second-order central Laplacian on the same grid. The pressure is then found from the projection equation

∇²p = (1 / Δt) ∇ · u^*, uⁿ⁺¹ = u^* - Δt∇p,

where u^* is the provisional velocity after advection, diffusion and body forces. The Poisson equation is solved by successive over-relaxation (SOR), with a deliberately fixed iteration budget so the animation remains responsive. The projection is therefore good enough to make the flow visually incompressible, but it is not converged to the tolerances one would demand from a DNS or LES calculation.

The demos expose this directly as the max div diagnostic. In an exactly incompressible calculation ∇ · u would be zero everywhere after the pressure projection. Here it is reduced, but never exactly zero. Numerically that means the browser model is still allowing tiny local sources and sinks of volume: air is being slightly created in some cells and slightly removed in others. That is acceptable for an interactive visual model, but it is precisely why these demos should not be used for quantitative mass-conservation arguments.

Why this is fine for a teaching demo

Stable advection — visitors can move parameters and obstacles without the simulation immediately becoming unusable.
Cheap enough for a browser — the SOR projection is simple, predictable, and runs interactively on a 200×100 grid.
Directly connected to the equations — the demos use the same pressure-projection structure that appears in many research flow solvers, only in a much smaller and more dissipative form.
Captures the qualitative physics — onset of convection, plumes, wakes, channelling and stagnation — at the right trends in Ra, Re and Pr.

Why this is not for real applications

Several deliberate simplifications keep the solver fast enough for a browser tab. Each one is a reason not to trust quantitative output:

Numerical dissipation from semi-Lagrangian advection. First-order accurate in space; small-scale features bilinearly average toward the mean within a handful of steps. The effective viscosity is set by the grid spacing, not by the molecular Re or Pr you typed into the slider. So the slider values are aspirational at high Re/Ra.
No subgrid-scale model. At resolutions a fluid mechanician would call "turbulent", everything below the grid is gone — there is no Smagorinsky, no dynamic SGS, no resolved boundary layer.
Crude wall treatment. First-order extrapolation on no-slip walls; no wall-resolved turbulence and no log-law.
Approximate pressure projection. The SOR solve is deliberately capped so each frame remains responsive. It reduces divergence enough for the visual model, but it does not make ∇ · u exactly zero. The displayed nonzero max div is the residual local creation/removal of fluid volume left by the approximate projection.
First-order in time. Phase errors accumulate over many steps; the demo cannot honestly produce energy spectra.
Effective Prandtl number is partly numerical. The numerical viscosity blends with the molecular viscosity to set the effective Pr. So the Pr slider on the Rayleigh–Bénard demo controls a mixed parameter, not the pure molecular ratio.

The net effect is that quantitative outputs (Nu, Cd, shedding frequency) trend correctly with parameters but their absolute values can be off by 20–50% from the literature, depending on the regime. Use the demos to see that Nu grows with Ra; use a research code to argue what Nu is.

What we use for actual research

The lab's urban LES code uDALES (second-order finite-difference, MPI-parallel, immersed boundary method, multiple SGS options) handles the full 3D atmospheric flow problem at research scale. For canonical turbulence work, including Rayleigh–Bénard-style convection, jets and plumes, gravity currents, stratified turbulence, fountains, and clouds, we use SPARKLE, the lab's in-house finite-volume DNS solver for incompressible turbulent flow. SPARKLE uses a fourth-order symmetry-preserving spatial discretisation for the convection and diffusion terms — the symmetry-preserving form ensures that energy and other quadratic invariants are conserved by the discrete operators, so the only numerical dissipation in a simulation is the molecular viscosity we ask for — paired with a third-order adaptive Adams–Bashforth time integrator that adjusts its step to the running CFL condition.

References

For the staggered-grid finite-volume/finite-difference discretisation: Ferziger, J. H., Peric, M., & Street, R. L. (2020). Computational Methods for Fluid Dynamics. 4th ed., Springer.
For the pressure-projection formulation: Chorin, A. J. (1968). Numerical solution of the Navier–Stokes equations. Mathematics of Computation, 22, 745–762.
For successive over-relaxation as an iterative elliptic solver: Young, D. M. (1954). Iterative methods for solving partial difference equations of elliptic type. Transactions of the American Mathematical Society, 76, 92–111.
For the Rayleigh–Bénard physics behind the convection demo: Ahlers, G., Grossmann, S., & Lohse, D. (2009). Heat transfer and large-scale dynamics in turbulent Rayleigh–Bénard convection. Reviews of Modern Physics, 81, 503–537.
For the research-grade urban LES code discussed above: Suter, I., Grylls, T., Sützl, B. S., Owens, S. O., Wilson, C. E., & van Reeuwijk, M. (2022). uDALES 1.0: a large-eddy simulation model for urban environments. Geoscientific Model Development, 15, 5309–5335.