Local simulation workbench

Model scope and interpretation

Reactive CO₂ / water / calcite

A separate PHREEQC 3.8.6 extension uses the pinned USGS property dataset at 25 °C. Initially pure water reacts with finite calcite and a fixed-volume, initially dry CO₂ headspace. Water vapour is explicitly included at equilibrium, with the dataset’s CO₂/H₂O PR interaction parameter 0.19. Total gas pressure and CO₂ partial pressure are distinct absolute quantities. Aqueous concentrations are mol/kg of the final solvent water, not mol/L. Dissolution may exhaust the mineral, in which case undersaturation is allowed but supersaturation is not. Carbon/calcium/hydrogen/oxygen inventories, charge, ionic strength, reaction mass action, calcite complementarity, both gas solubility equilibria and the mixture PR equation must close before export.

Pressure-corrected equilibrium constants and solvent mass remain native property outputs. The browser independently reconstructs balances and mass action, checks constants against the pinned 25 °C references, and screens activity coefficients using rounded water-property constants. This is not an independently validated property package. Very low gas inventories can fail the independent gas checks and are rejected. No kinetics, salt/brine model, hydrocarbon liquid, redox network, reactive oil/water VLLS, mineral solid solution, calorics or guaranteed arbitrary-chemistry stability is supplied.

Transient immiscible two-fluid flow

Foundation OpenFOAM 14 incompressibleVoF advances phase fraction and momentum on a true three-dimensional Cartesian mesh with adaptive Courant-limited time steps. Equal-property periodic slab advection has an exact transported cell-average reference. A separate closed water/air tank release exercises gravity, surface tension, no-slip sidewalls and density contrast. All saved times include complete native fields; the browser verifies the Cartesian cells, initial column, water volume, phase bounds, energy ceiling and available analytical checks. Reduced pressure p_rgh is in Pa, unlike the steady workbench’s kinematic pressure.

The advection test uses fluids with identical density/viscosity and is invariant transversely. It is not density-contrast validation. The water/air tank has mesh and time-step checks but no measured-flow validation. These short isothermal laminar cases do not predict turbulent wellbore slugging, compressible trapped-air impacts, contact-angle hysteresis, dissolution, phase change or reacting transport. The reactive and CFD extensions are separate engines, not a coupled reactive-flow simulator.

USGS PHREEQC release and documentation · Gas-phase equations · Foundation VOF formulation

Configurable fluid and pure-crystal equilibrium

The equilibrium model accepts two to twenty components with explicit critical temperatures, critical pressures, acentric factors, feed mole fractions and a complete symmetric binary-interaction matrix. The native fluid solver is Thermo 0.6.0 FlashVLN with classical Peng–Robinson mixing. A separate joint Gibbs minimization allocates up to eight distinct pure-component crystals while re-solving the depleted fluid composition. Up to four candidate liquids may be requested; a requested phase is not necessarily present at equilibrium.

Each crystal needs a matched melting-temperature/fusion-enthalpy reference. The constant heat-capacity difference is liquid minus solid; the constant molar-volume difference is solid minus hypothetical liquid. Zero explicitly neglects that correction. Results include phase mole amounts, compositions, fluid chemical-potential closure, crystal complementarity and multistart tangent-plane stability diagnostics. Negative stability or poor closure prevents an accepted export. Finite starting points cannot prove a global minimum for arbitrary mixtures.

The five-component benchmark contains propane, heptane, water, octadecane and eicosane. NIST fusion values use matched Barbillon reference rows: 301.5 K and 61.5 kJ/mol for octadecane, 309.6 K and 69.0 kJ/mol for eicosane. Its water/hydrocarbon kij of 0.48 and hydrocarbon kij of zero are illustrative assumptions, not calibrated measurements. This is mathematical model verification, not a validated crude-oil prediction. No caloric or validated bulk-density predictions are inferred. Raw potential values use the dimensionless ln(f/P) reference, not fugacity in pascals or absolute chemical potentials.

This is not unrestricted all-chemistry VLLS. Reactions, electrolytes, compound crystals, general solid solutions, hydrate structures, competing polymorphs and fully solid endpoints remain outside this companion. Near-exhausted fluid bounds are rejected rather than clipped into a plausible-looking allocation. Model and property selection still require independent engineering review.

Independent temperature and pressure sweeps

Optional PR sweeps independently solve 3–21 states at fixed feed, varying temperature or absolute pressure while holding the other coordinate fixed. Each point must pass the original equilibrium checks; a failed point rejects the entire result. The chart shows vapor, combined liquid and combined pure-crystal mole amounts, with per-sample compositions and a CSV in K and Pa absolute. Liquid numbering is local to each state. Connected points are not a traced phase envelope, located transition or energy balance.

Companion 1.3 includes PR volume roots so the browser independently recomputes pressure, mixture fugacity with the full binary-interaction matrix, and fusion/Cp/volume crystal potentials. These are numerical audit quantities, not validated density predictions. Both benchmark sweeps are also checked by a separate 60-digit joint phase-allocation calculation; this does not validate the assumed physical properties. Older imports explicitly identify their weaker audit coverage.

Calibrated quaternary wax solid solution

A separate low-pressure wax model uses an ideal liquid and one mixed Wilson solid for n-C11, n-C14, n-C16 and n-C18. Its empirical interaction is Λij = exp[−a(Ci−Cj)²/T], with a = 131.6938 K per carbon-number squared. This is a fitted Wilson parameterization, not the predictive Wilson or UNIQUAC models in the source paper, and is not coupled into the PR vapor/multiple-liquid flash.

The measurements are the 30 quaternary mixtures in Table 1 of Shariatrad and colleagues (2022), at approximately 0.9 bar. Every third mixture was withheld before fitting: 20 training and 10 held-out mixtures. The training WDT average absolute error is 0.837 K (maximum 2.893 K); held-out error is 0.733 K, with a maximum of 1.227 K. The ideal-solid held-out baseline is 5.671 K. Experimental expanded temperature uncertainty is 0.1 K; these model errors are larger, not experimental-uncertainty agreement.

Input temperature is limited to 275–294 K and the four named components. Printed experimental mass fractions that sum to 0.9999 or 1.0001 are explicitly normalized only when preparing the calibration dataset; raw printed numbers are retained. User case fractions must sum to one. Mixed-solid allocations, trace-component closure and phase disappearance are numerically checked. Below-WDT solid amounts have no measured validation. No crude-oil calibration, solid-solid transition, Cp correction, pressure extrapolation or all-solid endpoint is provided.

Optional temperature sweeps independently solve 3–41 equilibrium samples at fixed feed and approximately 0.9 bar. Every sample must satisfy the balance and chemical-potential checks; a failed sample rejects the entire result. The chart separates solid mole amount from solid mass fraction and marks the independently solved WDT. Connecting samples does not add solved states, and the sweep is not a transient cooling, nucleation, hysteresis or deposition model.

Composition support is checked against two convex hulls: the 20 training mixtures and all 30 measured mixtures. The ten held-outs are excluded from the training hull. Outside a hull means composition extrapolation; inside or on its boundary establishes geometric support only, not measured accuracy or validation of solid fractions. This is stricter than checking each component’s minimum and maximum separately.

Native laminar pipe and concentric annular flow

Laminar templates use Foundation OpenFOAM 14 without a turbulence closure. They solve steady, incompressible, isothermal Newtonian flow with stationary no-slip walls and fully developed axial-periodic conditions. The conservative 1–1,000 hydraulic Reynolds-number screen does not predict transition. Eccentric laminar annuli are not supported. There are no k, omega, eddy-viscosity or y-plus outputs.

The pipe limit is G = 8νU/Ro² and Darcy f = 64/Re; the annular friction factor is not generally 64/Re based on hydraulic diameter. Full native fields must agree with circular Poiseuille theory to less than 5% in gradient, maximum velocity error and volume-weighted RMS velocity error. Chordal inner-wall faces can place cells inside the nominal circular solid; their reference velocity is zero, not a negative profile extrapolation, and their volume fraction is reported. Native residual, two-checkpoint stationarity, local continuity and forcing/wall-drag checks remain mandatory. No analytical velocity field is substituted for an OpenFOAM result.

Eight native cases check three transverse resolutions per geometry, a separate pipe axial-resolution holdout and a Re = 250 holdout. Fine-grid analytical errors must be below 2%, final transverse friction change below 2%, and the holdout friction changes below 0.5% after expected Reynolds scaling. This is analytical numerical verification, not experimental validation or a uniform 3-D grid convergence estimate. Complete laminar verification record

Native three-dimensional turbulent mean flow

Foundation OpenFOAM 14 solves rectangular ducts, smooth circular pipes and stationary concentric or eccentric annuli with its kOmegaSST closure. Pipes and annuli use fully developed axial-periodic conditions and native mean-velocity forcing, not an entrance-flow model. Multiple cells in all three directions and no-slip walls distinguish this from a two-dimensional cross-section calculation. The original duct has a uniform inlet and specified turbulence intensity and length scale, with zero outlet kinematic gauge pressure. Circular cases use turbulence values as initialization only. Wall grading and turbulence boundary conditions are part of the generated case, not friction-factor substitutions.

RANS predicts mean flow rather than resolving turbulent eddies. The Reynolds-number screen does not predict transition. Native mesh checks, steady residual convergence, inlet/outlet volume closure and y-plus diagnostics accompany the imported field. Pressure is p/rho in square metres per square second, not pascals. The view contains sampled native cell fields; positions are arithmetic vertex means, not a separately exported native cell-centre field. The complete native case and VTK export remain on the local machine.

New CFD exports also check volume-averaged absolute finite-volume divergence normalized by U/Dh, with a limit of 10⁻⁶. This comes from the native local-continuity log, not just periodic inlet/outlet balance. Old version 1.0 result files remain readable without this newer diagnostic.

The original duct uses the native k-based wall-function y-plus diagnostic. Circular cases use resolved-wall SST conditions: k = 0, zero wall eddy viscosity via nutLowReWallFunction, and omegaWallFunction. Their native y-plus uses wall velocity gradients, and the maximum must be at most 2 before an accepted export. At least 16 radial cells and graded near-wall spacing are required. The browser samples come from single-precision legacy VTK; higher-precision original fields remain in the native case directory.

This release prepares only steady, isothermal, incompressible Newtonian ducts, fully developed periodic circular pipes and stationary annuli. Relative annular eccentricity is center offset divided by radial clearance and is limited to 0–0.6; inner diameter is limited to 90% of outer diameter. Hydraulic diameter is outer minus inner diameter, not outer diameter. Arbitrary meshes, rotating walls, cuttings, transient flow, LES, DNS, compressible flow, reactions and multiphase CFD are not exposed. Mesh refinement and comparison with applicable measurements remain necessary; convergence is not physical validation. SST variants differ, so another implementation is not interchangeable without checking its precise equations and wall treatment.

Circular cases use conformal four-block annuli or a five-block pipe O-grid, with curved boundary edges and no collapsed axis cells. Finite-volume boundary faces are polygonal: their area deficit relative to the circular section is reported, and flux is checked against the actual polygon area. Periodic lengths are 0.2–5 outer diameters. The driving gradient must settle over the last 20 iterations and its integrated force must agree with native wall shear drag within 1%. Nonzero Ux/k/omega initial residuals must be below 10⁻⁶. Near-zero transverse velocity and periodic pressure use physical scales instead of ill-conditioned relative residual norms: transverse speed / bulk speed and pressure span / bulk speed² must each be below 10⁻⁵. U/k/omega/nut fields must change by less than 10⁻⁴ between saved checkpoints at least 20 iterations apart. Equivalent p/rho is reconstructed from native periodic-pressure fluctuations plus the solved macroscopic gradient, referenced to zero at the outlet plane. The equivalent pressure loss is gradient times periodic length, not entrance loss. These checks do not establish mesh independence or experimental accuracy.

The recorded duct benchmark uses 1,728, 13,824 and 46,656 cells. Successive pressure-loss changes are 0.87% and 0.94%, with a non-monotonic trend. This comparison is not a grid-independence claim or an accuracy bound. Recorded three-mesh diagnostics

Ten additional native cases check pipe, concentric and eccentric annular transverse refinement and a separate pipe axial-resolution holdout. Final friction changes are 1.57%, 2.44% and 2.57%; the axial holdout changes friction by less than 0.001%. Fine-pipe friction differs by 1.44% from the original Princeton Superpipe operating point at Re = 31,577, while maximum axial velocity differs by about 4.50%. This is one measured pipe comparison, not measured annular validation or a universal accuracy bound. Complete pipe and annular validation record

File provenance and reproducibility

Results include the engine version, original case and SHA-256 of the exact case-file bytes. The browser independently checks component or boundary-flow balance and diagnostic consistency. Editing an input clears the old result. Hash matching detects a mismatched case but does not authenticate an engine run or validate properties. Recorded benchmarks are labelled, never presented as fresh calculations. Case and result files are portable backups; no permanent browser archive or cloud backup is implied.

Sources and related calculations

• Foundation 14: linear solvers, preconditioners and relaxation
• Thermo 0.6.0 FlashVLN source: candidate liquids and unsupported native solids
• Foundation 14: laminar and Newtonian transport selection
• MIT Advanced Fluid Mechanics: no-slip concentric annular flow problem
• SciPy Qhull: convex hull supporting-plane equations
• SciPy HiGHS: independent nonnegative linear feasibility checks
• Foundation 14 volume-averaged continuity definition
• Foundation 14 residual normalization source
• Foundation 14 mean-velocity forcing constraint
• Foundation 14 resolved-wall viscosity and y-plus implementation
• NIST: undecane fusion references
• NIST: tetradecane fusion references
• NIST: hexadecane fusion references
• Shariatrad et al. (2022): measured quaternary wax disappearance, Table 1
• Thermo Wilson model: excess Gibbs energy and activity coefficients
• Princeton Superpipe: measured smooth-pipe operating point Re 31,577
• Foundation blockMesh: curved edges and conformal multi-block meshes
• Thermo: multiphase flash (solids require the companion's outer minimization)
• NIST: n-octadecane phase-change data
• NIST: n-eicosane phase-change data
• Foundation OpenFOAM 14 installation and native solver
• Foundation OpenFOAM 14: k-based wall-function y-plus definition
• Turbulence Modeling Resource: precise SST model variants