Flowsheet simulator — validation

Each check below runs the live calculation engine (mb-engine 0.17.0) and compares it with a reference computed outside the engine. The engine numbers are generated by a build script, not typed in.

1 · Peng–Robinson flash vs an independent implementation

Reference: the open-source thermo 0.6.1 library (Gibbs-stability flash with PR gas and liquid phases), given the same critical constants, acentric factors, binary interaction parameters and H₂ effective constants as the engine. This verifies the φ–φ flash itself; how well Peng–Robinson matches experiment depends on the k_ij values, which can be tuned.

2 · Peng–Robinson flash vs DWSIM

Reference: DWSIM 9.0.5, the open-source process simulator (Peng–Robinson property package, default flash), run through its automation interface on the same six feeds.

Same input data

DWSIM given the engine's critical constants, acentric factors, H₂ effective constants and k_ij — a second, separately written implementation of the same flash.

DWSIM's own data — ammonia loop

The two ammonia-loop cases again, this time with DWSIM's built-in constants and k_ij, as a new DWSIM user would get them. The flash method is the same; these inputs are not:

With matching inputs the two programs agree to about 6.1e-5, so the gap below comes from input data, not from the flash. With DWSIM's own data the separator vapour carries 2.842 mol % NH₃ against the engine's 2.555 mol %, and the H₂ and N₂ dissolved in the liquid differ by a factor of about 1.7 and 3.4. Dissolved N₂ follows the N₂–NH₃ k_ij, which the engine sets to zero; dissolved H₂ follows how H₂ is represented (effective or true critical constants); the NH₃ in the vapour moves with both, in opposite directions. This comparison does not show which set is closer to a real loop; for loop work, check these numbers against plant separator data.

3 · NRTL activity coefficients

γ from the engine's NRTL flash at its converged liquid composition, against thermo's NRTL model at the same composition and temperature (Aspen form τ_ij = a_ij + b_ij/T, α = 0.3; parameters are test inputs, not a recommended data set).

4 · Energy balance with condensation

20 mol % steam in N2/CO2 at 150 °C / 5 bar + N2 at 25 °C → cooler 40 °C → knock-out 4.8 bar. Reference: CoolProp 8.0 (IAPWS-95 water, reference equations of state for N₂ and CO₂) evaluated at the component partial pressures with a Raoult split at the drum.

This case exercises the dew-point rule: wet process gas above its dew point is treated as all vapour, with water evaluated at its partial pressure.

5 · Flowsheet solver

6 · Property methods — flashes vs DWSIM

The rigorous units use a property method chosen for the whole flowsheet: Peng–Robinson, Soave–Redlich–Kwong, NRTL, UNIQUAC, Wilson, UNIFAC, Dortmund UNIFAC or ideal, with pure-component data for 432 components and binary parameters from the ChemSep databank. Each case below is a flash or bubble / dew point run through the engine's own heater and flash units, against DWSIM 9.0.5 with the same property package and parameters (Poynting and Henry corrections off on both sides).

UNIQUAC is not listed: DWSIM's built-in UNIQUAC parameters for methanol–water and acetone pairs differ from the ChemSep set the engine uses; against thermo with identical parameters the engine's activity coefficients agree to machine precision for all five activity models.

7 · Three-phase (VLLE) and liquid–liquid flash

The three-phase separator, and every flash when the method's valid phases include a second liquid, tests the two-phase answer for another liquid (Michelsen tangent-plane distance), converges vapour, light liquid and heavy liquid with a three-phase Rachford–Rice, and keeps the state with the lowest Gibbs energy. References: thermo's FlashVLN (one Peng–Robinson gas and two Peng–Robinson liquids, same constants and k_ij) for water–hydrocarbon systems, and isoactivity equations solved separately with thermo's activity models for partially miscible liquids.

8 · Rigorous distillation columns vs DWSIM

Equilibrium-stage columns solved by simultaneous correction (Naphtali–Sandholm Newton on the full MESH equations), against DWSIM's rigorous column on the same feed, stages, pressure and specifications (reflux ratio and bottoms rate).

9 · Gibbs reactor vs NIST-JANAF

Equilibrium by Gibbs-energy minimisation (the Gordon–McBride element-potential iteration of NASA CEA, so species whose element is nearly exhausted still converge), with no reactions written, against equilibrium constants from the NIST-JANAF formation Gibbs energies: water-gas shift at 1000 K and ammonia synthesis — stoichiometric at 700 K / 150 bar and with hydrogen nearly used up at 400 K / 100 bar.

10 · Phase envelopes vs DWSIM

Bubble and dew curves traced by Michelsen's continuation method (Newton on ln K, ln T and ln P with the most sensitive variable specified, stepping through the critical point), with the critical point located on the curve and the cricondenbar / cricondentherm refined. Reference: DWSIM 9.0.5 with the same critical constants and k_ij — its critical point from the Heidemann–Khalil method, its cricondentherm as the hottest point of its own dew curve (sampled every 0.5 bar). Points along the traced curves were also confirmed with independent bubble / dew calculations in thermo to within 2·10⁻⁴ K.

Property table

A temperature / pressure sweep tabulates phase fractions, enthalpy, entropy, Cp and the density, viscosity and thermal conductivity of each phase. Reference for the vapour density and compressibility: thermo's Peng–Robinson gas root with the same constants and k_ij.

11 · Binary T–x–y diagrams

Bubble and dew curves from the flowsheet property method, azeotropes located by bisection on y₁ − x₁, and — for partially miscible pairs — the three-phase line found from the three-phase flash, with the stable diagram drawn across the liquid–liquid gap. Reference: modified Raoult's law with thermo's NRTL / UNIFAC activity coefficients and the same vapour-pressure equations, solved independently with scipy (isoactivity liquid–liquid equilibrium plus the bubble condition for the heteroazeotropes).

12 · Kinetic plug-flow and CSTR reactors vs exact solutions

Power-law kinetics in the chosen phase of the property method. The plug-flow reactor integrates molar flows and enthalpy flow (Dormand–Prince RK5(4), switching to the L-stable Rosenbrock method RODAS3 when the kinetics are stiff); the CSTR solves the extents and temperature by Newton from several starts, reports every steady state it finds, and takes very fast kinetics to the steady state by pseudo-transient continuation. References: Levenspiel's closed forms for first-order reactions with volume change, scipy LSODA and root-bracketing on the same ideal-gas balances (adiabatic, cooled and fixed-duty reactors, including all three steady states of an exothermic CSTR), the closed-form solution of a linear system with a stiffness ratio of 3·10⁶, and the NIST-JANAF water-gas-shift equilibrium for reactors whose reverse rate comes from the equilibrium constant.

With very fast reversible kinetics the plug-flow and CSTR outlets also reproduce the engine's Gibbs reactor (adiabatic, same species) to 10⁻⁷ in mole fraction for steam reforming + shift, butene isomerisation and ammonia synthesis — two independent routes to the same equilibrium.

13 · Catalytic (LHHW / Hougen–Watson) rate laws

Rates of the form k(T)·[Π a^order − reverse]/(Σ K_t(T) Π a^e)^m, with a thermodynamically consistent driving force when the reverse rate comes from the equilibrium constant (the databank's or your own ln K correlation). References: the closed-form Langmuir–Hinshelwood plug-flow and stirred-tank solutions, and the published Xu & Froment (1989) steam-reforming model — three reversible reactions with negative hydrogen orders and a squared adsorption term — integrated independently with scipy, isothermal and adiabatic.

14 · Vapour association of carboxylic acids

With an activity-coefficient method, acids dimerise in the vapour (chemical theory: ideal gas of monomers, dimers and cross dimers). The measured vapour pressure and heat of vaporisation refer to this associated vapour, so the liquid enthalpy built from them is now consistent with the formation enthalpies — without association it is about 28 kJ/mol too high for acetic acid, which made liquid heats of reaction and adiabatic temperatures wrong. References: NIST liquid enthalpies of formation, the DIPPR vapour-pressure equation, and an independently written chemical-theory solution for vapour compressibilities.

15 · CO₂ absorption in aqueous MEA

The aqueous-MEA method solves the carbamate (2 MEA + CO₂ ⇌ MEAH⁺ + MEACOO⁻) and bicarbonate (MEA + CO₂ + H₂O ⇌ MEAH⁺ + HCO₃⁻) equilibria of the published Akula et al. model as implemented in the open-source IDAES process-modelling framework; CO₂ follows Henry's law on the free molecular CO₂ (N₂O analogy), water and MEA Raoult's law on the true species. Enthalpies come from the temperature dependence of the same fugacities, so the heat of absorption is consistent with the solubility. References: the IDAES property-package test values, an independent scipy solution of the published equations, and the IDAES heat of absorption.

Against measured vapour pressures

Total pressure over 30 wt % MEA (7.017 mol/kg water) measured by Arshad et al. (J. Chem. Eng. Data 59, 2014, 764; NIST ThermoML), against the engine's bubble pressure. At 120 °C the model is within 2 % up to a loading of 0.32 mol/mol and then over-predicts more and more (+19 % at 0.40, +37 % at 0.44); at 80 °C it is within 10 % to 0.44, and at 40 °C within 20 % at 0.5. Above about 0.55 it over-predicts the CO₂ pressure about twofold — a known limit of ideal-solution chemical-equilibrium models. In a regenerator this makes CO₂ look slightly easier to strip at the rich end, so reboiler duties come out somewhat low.

Mean absolute deviation 8.2 % for loadings up to 0.45 (10 points) and 112 % above (13 points). The test vectors below hold four of these points with tolerances set at the model's accuracy.

16 · Hot potassium carbonate (Benfield-type solutions)

A Pitzer electrolyte model (Harvie–Møller–Weare form with unsymmetrical mixing terms) with the parameters and equilibrium constants of the USGS PHREEQC pitzer.dat database: species K⁺, CO₃²⁻, HCO₃⁻, CO₂(aq), OH⁻ and H⁺, the water activity from the osmotic coefficient, and the water consumed by bicarbonate formation. Reference: PHREEQC itself, adding K₂CO₃ and CO₂ to 1 kg of water.

Full grid against PHREEQC

20, 30 and 40 wt % K₂CO₃ at conversions to bicarbonate of 0.05, 0.3, 0.6 and 0.9. Above 100 °C PHREEQC evaluates the solution at the water vapour pressure instead of 1 atm, which accounts for the small difference there.

The pitzer.dat KHCO₃ parameters are fitted to about 3 mol/kg. In strong, highly converted solutions they make the model thermodynamically unstable (a spurious liquid–liquid split) close to where KHCO₃ would in fact crystallise; the engine maps that onset over 0–250 °C and, past 95 % of it, holds the properties at the boundary and extrapolates smoothly (the stream report flags such states). PHREEQC evaluates those states without a stability check, so they are listed separately and not graded — all are cold or 40 wt % solutions converted 60–90 % to bicarbonate, beyond normal hot-carbonate practice:

Against handbook data

Water activity, heat capacity and density of potassium carbonate solutions from the Armand Products potassium carbonate handbook. The density correlation is fitted to the handbook's 15.6 °C table, so its temperature dependence is the test here; the heat capacity follows from the Pitzer parameters' temperature dependence and is not fitted.

Phase-boundary edge cases

Gas and condensate around a carbonate absorber and regenerator: wet CO₂ at its dew point with the salt listed but absent, the absorber overhead carrying a numerical trace of K₂CO₃, and the same gas condensing water. Reference: IAPWS-IF97 saturation of water at its partial pressure.

At high concentration the Pitzer water activity falls below the handbook's 20 °C vapour pressures (−8 % at 46 wt %) and the heat capacity of 40 wt % solution is 10 % low; between 20 and 30 wt % — the Benfield range — both agree within 1–4 %. The Pitzer parameters carry fitted temperature terms to about 90–100 °C and are extrapolated above that. Promoters (DEA, glycine, vanadium) change absorption rates, not the equilibrium, and are represented through tray efficiencies.

17 · Murphree tray efficiency

A tray with Murphree vapour efficiency E leaves its vapour at y = y_in + E(K x − y_in) instead of in equilibrium; the efficiency can be set for every tray or per component. Reference: the exact solution of the tray equations for a dilute absorber with a straight equilibrium line (numpy), which equals the Lewis–Kremser closed form. The regenerator case checks the condenser-temperature specification used for amine and carbonate strippers.

Scope

These checks verify the implementation against independent references. Agreement with another simulator shows the equations are solved correctly for the given parameters; it does not make a property method suitable for every system — Peng–Robinson without fitted k_ij over-predicts the solubility of water in hydrocarbon liquids, ChemSep NRTL / UNIQUAC parameters are fitted to vapour–liquid data, and the urea units remain empirical correlations. Use the results for preliminary engineering and confirm them against plant data before design use. Kinetic reactors are only as good as the rate expressions entered, and the acid association model covers dimerisation only (no acid–water cross association; binary parameters fitted without association may need refitting). The CO₂ packages cover CO₂ only (no H₂S, NH₃ or MDEA / activated MDEA chemistry), and tray efficiencies stand in for rate-based mass transfer. The inputs for every case above are in the engine test set, cases MB-089 to MB-195.