How to model flash drum vapor-liquid equilibrium with Peng-Robinson
Give the drum a temperature and a pressure and the solver runs an isothermal-isobaric flash: Wilson correlation for the first K-value estimates, Rachford-Rice for the vapor fraction, successive substitution on the K-values, then a Newton polish. Peng-Robinson supplies the fugacity coefficients on both sides, so K = phi_liquid / phi_vapor.
A flash drum is one equilibrium stage. Fix any two intensive variables — here temperature and pressure — and the phase split follows; you do not specify recoveries, you specify the conditions and read the split off the thermodynamics.
Peng-Robinson is the usual choice for hydrocarbons and light gases because it is an equation of state: it needs only critical constants and an acentric factor, so it works above the critical point of a component and needs no fitted binary to function. That is also its weakness — see the comparison below.
Show the governing equations
What the solver actually does
- Wilson correlation for a first K-value estimate from Tc, Pc and omega.
- Rachford-Rice for psi at those K-values.
- Successive substitution: recompute phi from the new phase compositions, repeat.
- Newton polish once the iteration is close.
The step most flash implementations get wrong is deciding there is no split at all. Far from the two-phase region successive substitution collapses to K = 1 for every component, which leaves Rachford-Rice indeterminate — any psi satisfies it, and a bracket midpoint of 0.5 leaks out as a phantom half-and-half result. MaximaLabs tests the Rachford-Rice endpoints on the seed instead: if the sum of z*K is at most 1 no vapour can form, and if the sum of z/K is at most 1 no liquid can. A near-critical mixture that fools that test is settled by comparing the Gibbs energy of the two compressibility roots.
Choosing the package
Measured against the NIST TRC ThermoML archive (378 binary pairs, 13,823 measured tie-lines), scoring deviation in vapour mole fraction over the 60 richest pairs:
| Package | Pairs it can evaluate | Median AAD in y | Mean AAD in y |
|---|---|---|---|
| NRTL | 17% | 0.031 | 0.038 |
| Peng-Robinson | 97% | 0.063 | 0.100 |
| SRK | 98% | 0.061 | 1.204 |
Read it as a trade, not a ranking. NRTL is roughly twice as accurate but only on the sixth of pairs it has fitted binaries for; the cubics answer almost everything at about twice the deviation. SRK's mean is worth staring at: 1.204 against a median of 0.061 means a handful of pairs fail catastrophically while the typical one is fine — a package that is usually good and occasionally very wrong, which is exactly why the median is the number the validation report gates on.
Practical rule: hydrocarbons, light gases and anything supercritical, use Peng-Robinson. Polar mixtures, alcohols, water, anything with an azeotrope, use NRTL — and above about 10 bar, where the low-pressure gamma-phi formulation stops being valid, use a Ge mixing rule (pr-mhv1) so the fitted binaries drive the cubic instead.
Numbers are read live from the validation report at /api/validation and move when the corpus or the packages do. Component attribution in the source archive is inferred rather than read, and that caveat travels with the figures.