BoilingSim

BoilingSim is a physics-based boiling engine. It solves, on a graphics processing unit, the problem of what happens when heat is driven through a solid vessel into a liquid until that liquid boils: how the heat spreads, how the fluid circulates, when and where bubbles form, and what that thermal history does to anything suspended in the liquid or immersed in it.

The engine is built around one idea: heating, boiling and material properties are a single coupled system, not three problems solved in sequence. Heat entering the vessel base is conducted through the wall and into the liquid as a conjugate problem, so the solid and the fluid are advanced together rather than the wall being pinned at an assumed temperature. Once the wall superheat crosses the onset of nucleate boiling, bubbles nucleate, grow, detach and rise as tracked particles, and the latent heat they carry away and the buoyancy they impart are scattered back onto the grid. That feedback changes the temperature field, which changes where the next bubbles form. The properties of the vessel material set the pace of the whole loop.

Materials are inputs, not assumptions

A materials database supplies density, specific heat and thermal conductivity for 304 stainless steel, cast iron, aluminum and copper, alongside temperature-resolved properties for water at saturation. Swapping the vessel material is a one-line change to a scenario file. The engine has been validated across a 25-fold span in thermal conductivity, from steel at 16 W/m·K to copper at 401 W/m·K, and it reproduces a result the physics demands: bubble departure diameter comes out material-independent at 2.93 mm, because departure depends on contact angle, surface tension and density difference, which are fluid properties. The metal on the other side of the wall changes how fast the wall gets hot, not what the bubbles do once it is.

The demonstration: household cooking and nutrient retention

To show the engine end to end, the reference case is a carrot boiling in a pot on a stove. It was chosen because it exercises every part of the physics at once and because published retention data exist to test against. Inside the carrot, represented as a tetrahedral mesh, nutrient loss is modelled as reaction, diffusion and leaching together: Arrhenius kinetics degrade the nutrient at the local temperature, diffusion moves it through the tissue, a Sherwood correlation carries it across the surface into the water, and a conservative advection scheme then transports it with the flow the fluid solver is already computing.

That structure is what makes the result meaningful rather than fitted. Every molecule is accounted for in one of four buckets, retained, leached, degraded or precipitated, and the four sum to 100.00 percent at every output step. Beta-carotene, which is effectively insoluble in water, comes out almost entirely degraded with negligible leaching. Vitamin C, which is highly water-soluble, loses a fifth of its mass to the water. The engine produces both behaviours from the same equations with different material parameters, and can track the two solutes concurrently through a single thermal and fluid field.

Why it matters

Most models of cooking treat the pot as a well-mixed bath at a single temperature. BoilingSim resolves the spatial temperature field and the bubble dynamics that drive it, which is what allows nutrient loss to be predicted from first principles rather than fitted after the fact.

Nothing in the solver is specific to carrots or to stoves. The same coupled engine applies wherever heat, flow, phase change and a chemical or biological change interact: blanching, pasteurization, sterilization, fermentation, drying, and industrial vessel design. Cooking is the demonstration because it is the hardest available test with public data, not the limit of the scope.

Try it: demonstrations on request

BoilingSim runs as a live interactive session rather than as software you install. On request the lab provisions the full stack, solver, streaming server and dashboard, on a cloud GPU and shares a session URL for a limited window. During that window you drive the simulation yourself from the browser: raise or lower the stove heat flux, swap the pot between stainless steel, aluminum and copper, resize the carrot, switch the tracked nutrient, and watch the 3D scene and the instrumentation respond while the run continues.

To arrange a demonstration, contact Dr. Ebenezer Miezah Kwofie at ebenezer.kwofie@mcgill.ca.

Status and team

BoilingSim was developed by Etornam Celestine Tsyawo and Emmanuel Amankrah Kwofie at the SASEL Lab, McGill University, under the direction of Dr. Ebenezer Miezah Kwofie. Manuscripts covering the heat transfer model, the nutrient retention coupling, and the software itself are in preparation. The code is held under an all-rights-reserved licence as a lab research asset.

The coupled solve

Liquid motion is solved with the Navier-Stokes equations under the Boussinesq approximation on a staggered marker-and-cell grid, using semi-Lagrangian advection, an implicit diffusion step and a pressure projection to enforce incompressibility. Heat is solved as a conjugate problem across the solid-fluid boundary, so the vessel wall and the liquid form one thermal system. An evaporative enthalpy sink at the free surface lets vapour leave the domain, which is what pins the bulk liquid near saturation instead of letting a sealed domain drift several degrees above it.

Boiling is resolved per bubble rather than through a bulk correlation. Once local wall superheat crosses the onset of nucleate boiling, sites nucleate and each bubble is carried as a Lagrangian particle through Mikic-Rohsenow growth, Plesset-Zwick condensation, Fritz departure and terminal-slip rise. The latent heat each bubble absorbs is deposited back into the grid as an energy sink and its buoyancy as a momentum source, both by trilinear scatter, closing the loop between the bubble field and the temperature and velocity fields that created it.

Solute transport

Nutrient loss is a reaction-diffusion-leaching problem rather than a decay curve. Arrhenius kinetics act on the solute inside the immersed body and, separately, on whatever has already dissolved into the liquid. Diffusion redistributes solute within the tissue under a zero-flux surface condition. A Sherwood correlation sets the surface mass-transfer rate, sampled against free-stream velocity a few cells off the surface rather than at the no-slip face. A conservative finite-volume upwind scheme then advects the dissolved solute with the velocity field. A solubility cap prevents the liquid being pushed past saturation, and any mass the cap clips is routed to an explicit precipitation bucket rather than being discarded.

A partition coefficient controls how strongly the solute prefers the tissue over the liquid, and it is the single parameter that separates a fat-soluble pigment from a water-soluble vitamin. Two solute slots run concurrently against the same thermal, fluid and bubble field, so a pair of nutrients with opposite transport behaviour can be validated simultaneously rather than in separate runs.

Validation

Boiling is checked against the Rohsenow nucleate boiling correlation across all three pot materials, landing at 0.92 to 1.04 times the correlation at a stove flux of 80 kW/m², inside the 15 to 30 percent scatter that the pool-boiling literature reports for Rohsenow itself. Mean bubble departure diameter of 2.93 mm, drawn from tens of thousands of sampled departure events per run, sits mid-band in the published 1.5 to 4.0 mm range for saturated water at atmospheric pressure. The conjugate heat transfer solver is validated independently of the boiling model by the conductive drop across the pot base, which reads 26.7 K measured against 26.7 K analytic for steel.

Nutrient retention is validated against published cooking experiments. Beta-carotene retention after 600 s in a 25 mm carrot is 88.72 percent against an 84 percent experimental reference, inside the 80 to 90 percent target band, with the predicted curve crossing the reference value within 20 s of the reference cook time. Vitamin C, a much more soluble solute, reaches 65.80 percent against a 64 percent reference in its own 40 to 70 percent band, and loses 21 percent of its mass to the water where beta-carotene loses effectively none. Running both solutes together in one pot reproduces each single-solute result to within 0.28 percentage points.

A physical instrumented boiling station is being assembled in the lab, with multi-channel temperature sensing, Python-based data streaming and storage, and automated control of the heating level on an electric stove, so the engine can be tuned against measurements taken under the same conditions it simulates.

Engineering for production

The whole pipeline runs as NVIDIA Warp kernels written in Python and compiled just in time for the GPU, which stays the canonical reference implementation. Where profiling identified the pressure Poisson solve as the dominant cost, a replacement kernel was written by hand in CUDA, using shared-memory tiling and coalesced access, and reached from Python through a Rust extension built with PyO3 and maturin. On an RTX 6000 Ada at 2 mm grid spacing that path cuts pressure projection time by 37.8 percent and end-to-end wall time by 20.8 percent for a steel pot, with similar figures for aluminum and copper.

The discipline around those optimizations is the part that makes the engine trustworthy. Every accelerated path sits behind an environment flag and a parity test that holds it to the reference implementation, at a relative tolerance of 1e-5 for a single step and better than 1e-4 across a full 200-iteration projection. Three additional scatter kernels and a full per-bubble kernel were ported, validated and benchmarked, then left switched off by default because the measurements showed the Python-to-Rust call overhead outweighed the kernel-level saving at realistic bubble counts. A Jacobi-preconditioned conjugate gradient pressure solver was built and proved correct, delivering roughly 4.5 times better divergence reduction than the Jacobi loop but running about 8 times slower per projection on this geometry, so it too ships switched off with the conditions under which it would become worthwhile written down. The negative results are recorded next to the successes, which is what lets a reader know exactly where the engine has been checked and where it has not.

Scenarios are declared in YAML and validated through Pydantic models before they reach the solver, so a run can be reconfigured, reviewed and reproduced without touching solver code. The test suite stands at 218 passing tests across the Python, Rust and CUDA layers, covering solver parity, conservation of scattered quantities, sign conventions and multi-step integration. A sign-convention gate in that suite caught a real error in the conjugate gradient derivation before it reached a result.

Deployment and live control

Results stream from the solver over a WebSocket server written in Rust with Tokio and Axum into a React dashboard built on React Three Fiber, which renders the stove, the pot, the liquid, the bubble field and the immersed body in 3D as the run progresses. The link is bidirectional: heat flux, vessel material, body size and the solute being tracked can all be changed from the browser mid-run, and a full scenario can be staged from a configuration page and rebuilt on the fly.

The solver, the streaming server and the dashboard each ship as a Docker image and come up together under Docker Compose on any host with an NVIDIA GPU, behind a GPU precheck that fails loudly rather than falling back to the CPU in silence. A RunPod configuration runs the same stack on rented cloud GPUs on demand, and a tunnel path exposes a stable URL from a workstation. Every published benchmark figure has its exact command, scenario file and expected headline number recorded, so results can be regenerated rather than taken on trust.