← The Process Simulation Handbook

Thermodynamic & Physical Property Models

Chapter 1 solved the flowsheet as a graph. But every unit in that graph — every flash, column stage, and reactor outlet — is only as correct as the property model behind it. Choosing that model is the single most consequential decision in a simulation, and the one most likely to be wrong. This chapter is how you choose, and how MaximaLabs helps you get it right.

2.1 The most important choice you make

The same flowsheet, solved with two different property methods, can give two completely different answers — not slightly different, qualitatively different. The textbook case is ethanol–water: it forms an azeotrope at about 89 mol% ethanol, a hard ceiling that ordinary distillation cannot cross. A model that captures the strong non-ideality of the liquid (an activity-coefficient model) predicts that azeotrope; a bare cubic equation of state, tuned for hydrocarbons, does not — so it will happily (and wrongly) report a column reaching 99.9% ethanol overhead. Same column, same specs; the thermodynamics decide whether the answer is physical.

2.2 Equations of state

An equation of state (EoS) relates pressure, volume, and temperature for both phases from one set of parameters. The workhorses are the cubics — Peng-Robinson (PR) and Soave-Redlich-Kwong (SRK):

$P,\;T,\;v$

Pressure, temperature, and molar volume — the three state variables the EoS relates.

$R$

Universal gas constant.

$a,\;b$

Peng-Robinson pure-component energy and co-volume parameters, computed once from each component's critical temperature/pressure and acentric factor.

$\alpha(T)$

A temperature-dependent correction factor on a, fit so the EoS reproduces the pure component's own vapor-pressure curve near its critical point.

The EoS gives fugacity coefficients

in each phase, and phase equilibrium is where they match, which fixes the K-value :

$\hat\varphi_i^{\,V},\;\hat\varphi_i^{\,L}$

Fugacity coefficient of component i in the vapor and liquid mixture — both computed from the same EoS at the mixture's T, P, and composition, so no separate liquid-phase model is needed.

$y_i,\;x_i$

Vapor- and liquid-phase mole fraction of component i.

$K_i$

Equilibrium ratio K_i = y_i/x_i, fixed once the fugacities of both phases match.

Cubics are the right default for non-polar, hydrocarbon, and high-pressure gas systems — natural gas, refinery cuts, compressors. MaximaLabs ships PR and SRK plus more specialized equations of state (PC-SAFT, CPA for associating fluids, GERG-2008 for LNG). You pick the property method from the status bar; the menu spans equations of state, activity models, electrolytes, steam, and polymer models:

2.3 Activity-coefficient models

For a strongly non-ideal liquid — alcohols, water, acids, anything with hydrogen bonding — the right tool is an activity-coefficient model. Here the liquid non-ideality lives in an activity coefficient

, and VLE is the modified Raoult's law:

$y_i,\;x_i$

Vapor- and liquid-phase mole fraction of component i.

$\varphi_i^{\,V}$

Vapor-phase fugacity coefficient of i (≈1 near atmospheric pressure, where this method is normally used).

$P$

System (total) pressure.

$\gamma_i$

Liquid activity coefficient of i — the non-ideality correction the activity model below supplies; 1 is the ideal (Raoult's-law) limit.

$P_i^{\text{sat}}(T)$

Pure-component vapor pressure of i at the mixture temperature T.

The default in MaximaLabs for polar systems is NRTL (Non-Random Two-Liquid), whose activity coefficient for component

is

$\gamma_i$

Activity coefficient of component i being solved for.

$x_i,\;x_j,\;x_k,\;x_m$

Liquid mole fractions — j, k, m range over every component (the sums); i is the one component gamma is evaluated for.

$\tau_{ij}$

Dimensionless NRTL interaction energy for the ordered pair i,j (asymmetric — tau_ij is not tau_ji).

$G_{ij}$

NRTL local-composition weighting factor for the pair i,j — how strongly j clusters around i relative to i itself.

$\alpha_{ij}$

Non-randomness parameter for the pair (a third fitted constant, often held near 0.2-0.3 across many families rather than regressed).

$a_{ij},\;b_{ij}$

The two fitted binary interaction constants — a is temperature-independent, b/T carries the temperature dependence.

$T$

Temperature (K).

The binary interaction parameters

are the data dependency — MaximaLabs carries a DECHEMA/databank set (and you can regress or override your own).

2.4 The rest of the activity-model family

NRTL isn't the only local-composition model, and it isn't predictive — every one of these needs fitted binary parameters. Four more are worth knowing by their actual equations, not just their names, because each makes a different trade between simplicity, physical grounding, and whether it needs data at all.

Van Laar is the oldest of the family — built directly on the van der Waals equation for regular solutions, two parameters, no temperature dependence. For a binary:

$G^{ex}$

Excess Gibbs energy of mixing (J/mol) — the departure from an ideal solution.

$x_1,\;x_2$

Liquid mole fractions of components 1 and 2.

$R,\;T$

Gas constant and temperature (van Laar's constants carry no further T dependence themselves).

$A_{12},\;A_{21}$

The two fitted, temperature-independent van Laar binary constants.

which splits into the two per-component activity coefficients MaximaLabs actually evaluates:

$\gamma_1,\;\gamma_2$

Activity coefficients of components 1 and 2.

$A_{12},\;A_{21}$

Same two fitted constants as above — note the asymmetry: A_12 weights gamma_1 with A_21, and vice versa.

$x_1,\;x_2$

Liquid mole fractions.

Its simplicity is also its limit — no temperature dependence means it can't track VLE over a wide range, but it's cheap to fit and still shows up for quick screening. It's a real, selectable MaximaLabs thermo package (

, regressed to the same cited DECHEMA ethanol-water points as NRTL/UNIQUAC/Wilson) — not just the textbook equation above, though its 2-constant fit is honestly a touch looser (it places the azeotrope at ~91.6 mol% vs. the 89.4% reference, a known trait of having no temperature dependence).

Wilson introduced local composition — the idea that the composition immediately around a molecule differs from the bulk because of energetic preference — which every model after it (NRTL, UNIQUAC) inherits:

$\gamma_1$

Activity coefficient of component 1 (ethanol here); gamma_2 follows the same form with subscripts 1 and 2 swapped.

$x_1,\;x_2$

Liquid mole fractions.

$\Lambda_{12},\;\Lambda_{21}$

Wilson local-composition parameters, Lambda_ij = (v_j/v_i)·exp(-a_ij/T), from the pure-liquid molar volumes v_i and two fitted binary energies a_ij.

Wilson fits VLE well with only two parameters, but structurally cannot predict liquid-liquid splitting — for that you need NRTL or UNIQUAC.

UNIQUAC (UNIversal QUAsi-Chemical) splits the excess Gibbs energy into a combinatorial term (molecule size/shape, from pure-component data alone) and a residual term (the energetic interactions, from two fitted parameters per pair):

$G^{ex}_{comb}$

Combinatorial part of the excess Gibbs energy — accounts for molecule size/shape differences alone, from pure-component data, independent of energetics.

$x_1,\;x_2$

Liquid mole fractions.

$\phi_1,\;\phi_2$

Segment (volume) fractions, from r_i below.

$\theta_1,\;\theta_2$

Surface-area fractions, from q_i below.

$q_1,\;q_2$

Pure-component surface-area parameters.

$z$

Lattice coordination number, fixed at 10 in MaximaLabs' implementation.

$G^{ex}_{res}$

Residual part of the excess Gibbs energy — the energetic contribution, carrying the two fitted binary parameters.

$\tau_{12},\;\tau_{21}$

UNIQUAC binary interaction parameters (below), asymmetric between the two directions.

$\theta_1,\;\theta_2$

Surface-area fractions, same as above.

$q_1,\;q_2$

Pure-component surface-area parameters.

where the segment fraction

and surface-area fraction come from pure-component size () and surface-area () parameters, and carries the two fitted binary energies. The coordination number .

UNIFAC (UNIversal Functional Activity Coefficient) is the predictive answer to "what if I have no binary data at all": it reuses UNIQUAC's exact combinatorial/ residual form, but builds

, , and the interaction terms from the molecule's functional groups rather than the whole molecule —

$r_i,\;q_i$

Molecule i's overall volume and surface-area parameters — the same r_i, q_i UNIQUAC's equations use, now built up rather than looked up.

$\nu_k^{(i)}$

Number of group k that molecule i contains (an integer count from its structure — e.g. ethanol has one -CH3, one -CH2-, one -OH).

$R_k,\;Q_k$

Tabulated volume/area parameters for functional group k, from the UNIFAC group table — not fitted per molecule or per pair.

$\gamma_i^{res}$

Residual (energetic) part of component i's activity coefficient.

$\nu_k^{(i)}$

Count of group k in molecule i, same as above.

$\Gamma_k$

Group k's activity coefficient evaluated in the actual mixture (all groups from every component, pooled).

$\Gamma_k^{(i)}$

Group k's activity coefficient evaluated in a reference solution of pure component i — subtracting this out is what makes the term vanish for a molecule made of only one group type.

where

counts how many of group molecule contains, are tabulated group volume/area parameters, is group 's activity coefficient in the actual mixture, and is the same group's coefficient in a reference mixture of pure . Break ethanol into -CH₃, -CH₂-, and -OH groups, water into H₂O, look up the tabulated group interaction energies, and UNIFAC predicts the γ's — no fitted binary data for this specific pair required. It's the fallback MaximaLabs' Property Estimator (§2.6) reaches for when no measured VLE exists for a new mixture; the trade is accuracy — a regressed NRTL/UNIQUAC fit to real data beats a UNIFAC prediction whenever that data exists.

2.5 Choosing a property method

The classic decision reduces to a short tree: polar / hydrogen-bonding liquid at low-to- moderate pressure → an activity-coefficient model (NRTL/UNIQUAC); non-polar / hydrocarbon / high-pressure gas → a cubic EoS (PR/SRK); water + gas at high pressure → an associating model (CPA); aqueous salts → an electrolyte model. MaximaLabs encodes that judgment: it looks at your component list and, if the active method is a poor fit, surfaces a one-click recommendation. Put a polar ethanol-water mixture on Peng-Robinson and this banner appears over the canvas — the decision tree, made actionable, before a wrong answer is ever computed:

2.6 Seeing it: the phase diagram

The honest test of a property method is whether it reproduces measured phase equilibrium. MaximaLabs plots the binary T-xy / P-xy diagram straight from the active thermo package, so you can look at the VLE before trusting a column. For ethanol-water under NRTL, the vapor and liquid curves pinch together and cross the diagonal at the azeotrope — exactly the feature a cubic EoS misses:

This is also where trust comes from: MaximaLabs pins the ethanol-water NRTL prediction against published data on its validation dashboard (bubble-point AAD < 0.2 K, azeotrope at x ≈ 0.894 / 351.3 K).

2.7 When a component isn't in the databank

Real work runs into molecules the databank has never heard of. Rather than stop, you estimate the missing pure-component properties from structure — group contribution (Joback) for critical constants and normal boiling point, Ambrose-Walton for vapor pressure. MaximaLabs does this in the Thermo tools panel — Analysis & reports ▸ Thermo tools ▸ Property estimator. Enter the molecular groups (or look the component up by CAS number) and it returns the critical constants, normal boiling point, and acentric factor; the estimated component then drops straight into the same flowsheet:

2.8 Try it: watch the thermodynamics change the answer

The lesson of the chapter: the solver from Chapter 1 will faithfully converge whatever thermodynamics you give it — so the property method is yours to get right. Next up, the calculation every one of these models feeds: the flash.

2.9 Exercises