Part 2 of 3 · Black Oil PVT from Scratch

PR3 EOS Engine in Python: Implementation and Verification

29 May 2026

Part 1 established the characterised fluid — 12 components with T_c, P_c, ω, BIPs and volume shifts computed from MW and SG using Kesler-Lee correlations, all verified against commercial software exactly.

This part builds the PR3 EOS engine and computes the untuned bubble point pressure. The implementation follows Peneloux (1982) and Peng-Robinson (1976) exactly. All results are verified against the commercial software untuned output of 200.84 bar at 107°C.

Three implementation details are not explicitly documented in the software manual or standard references. Each is identified and verified here.

PR3 Pure-Component Parameters

The Peng-Robinson EOS describes fluid pressure as a function of molar volume:

P = RT / (V − b) − a(T) / [V(V + b) + b(V − b)]

Two parameters per component at reservoir temperature T = 380.15 K:

Attractive parameter a_i:

a_i = Ω_A · (R²T_c² / P_c) · α_i α_i = [1 + κ_i(1 − √T_r,i)]²

Covolume b_i:

b_i = Ω_B · RT_c / P_c

with Ω_A = 0.45724, Ω_B = 0.077796, R = 83.14472 cm³·bar/(mol·K).

PR3 kappa — the third parameter:

Standard PR2 uses a quadratic κ that becomes unreliable for heavy components. PR3 switches to a cubic formula at ω > 0.491:

κ = 0.37464 + 1.54226·ω − 0.26992·ω² (ω ≤ 0.491) κ = 0.379642 + 1.48503·ω − 0.164423·ω² + 0.016666·ω³ (ω > 0.491)

For this fluid the cubic formula activates for C16+ (ω = 1.374). This is the primary reason PR3 is used for heavy oils — PR2 underestimates attraction forces for large acentric factors.

Table 1 — Pure-component PR3 parameters at T = 380.15 K

Component	T_r	κ	α	a_i (cm⁶·bar/mol²)	b_i (cm³/mol)
N₂	3.012	0.436	0.462	6.84×10⁵	24.05
CO₂	1.248	0.708	0.841	3.34×10⁶	26.68
C1	1.994	0.395	0.701	1.75×10⁶	26.78
C2	1.245	0.524	0.882	5.33×10⁶	40.45
C3	1.028	0.603	0.983	1.00×10⁷	56.34
iC4	0.932	0.650	1.046	1.51×10⁷	72.37
nC4	0.894	0.674	1.075	1.62×10⁷	72.44
iC5	0.826	0.711	1.134	2.24×10⁷	87.87
nC5	0.810	0.745	1.155	2.39×10⁷	90.13
C6	0.749	0.812	1.230	3.33×10⁷	109.05
C7-C15	0.593	1.061	1.547	8.78×10⁷	181.44
C16+	0.383	2.153	3.316	1.49×10⁹	925.84

Mixing Rules

Van der Waals mixing rules with BIPs from Part 1:

a_mix = Σᵢ Σⱼ zᵢzⱼ √(aᵢaⱼ) · (1 − k_ij) b_mix = Σᵢ zᵢbᵢ

The Peneloux volume shifts from Part 1 enter as:

c_i = s_i · b_i c_mix = Σᵢ zᵢc_i

At P = 213 bar, T = 380.15 K, the dimensionless mixture parameters are:

A = a_mix·P / (RT)² = 18.353 B_raw = b_mix·P / RT = 1.402 B_eff = (b_mix − c_mix)·P / RT = 0.920

Implementation Detail 1 — The Cubic Solver Uses B_raw

The PR cubic in Z:

Z³ + (B−1)·Z² + (A − 2B − 3B²)·Z + (B³ + B² − A·B) = 0

Per the Peneloux (1982) formulation, the cubic is solved with B_raw, the unshifted covolume parameter. The volume shift is not incorporated into the cubic coefficients.

At P = 213 bar, B_raw = 1.402. This gives a single real root Z_raw = 1.624 — physically correct for a compressed liquid above its bubble point. As pressure decreases toward the bubble point, B_raw drops and two roots appear.

The distinction between B_raw and B_eff matters practically: using B_eff in the cubic gives incorrect root locations and therefore incorrect fugacity coefficients.

Implementation Detail 2 — Reported Z and Volume Are Shifted

The commercial software does not report Z_raw. It reports the Peneloux-corrected quantities:

Z_reported = Z_raw − c_mix·P / RT V_reported = V_raw − c_mix

This is verified by back-calculation from the reported Z = 1.0793 at bubble point:

Z_raw = 1.0793 + (71.6 × 200.8) / (83.14 × 380.15) = 1.534

Which is exactly the root of the PR cubic at P = 200.8 bar with B_raw. The density computed from V_reported matches to 0.3%.

Table 2 — Raw vs shifted quantities at bubble point (200.84 bar)

Quantity	Raw (unshifted)	Shifted	Commercial software
Z_liquid	1.534	1.079	1.0793
V_liquid (cm³/mol)	241.4	169.8	169.800
Density (kg/m³)	516	735	732.7

Implementation Detail 3 — Volume Shift Does Not Enter Fugacity

The PR fugacity coefficient for component i:

ln φᵢ = (bᵢ/b_mix)·(Z−1) − ln(Z−B) − A/(2√2·B) · (daᵢ/a_mix − bᵢ/b_mix) · ln[(Z + (1+√2)·B) / (Z + (1−√2)·B)]

where da_i = 2·Σⱼ zⱼaᵢⱼ.

No volume shift term appears. Z and B here are Z_raw and B_raw respectively. Peneloux (1982) proved that the shift cancels exactly in the fugacity ratio φᵢᴸ/φᵢᵛ — it therefore has no effect on K-values or bubble point pressure. This is why the untuned bubble point can be reproduced exactly without knowing the volume shifts.

Table 3 — Fugacity coefficients at 213 bar

Component	ln(φᵢ)	φᵢ	Volatility
N₂	1.247	3.480	High
C1	0.556	1.743	High — primary driver of bubble point
CO₂	−0.070	0.932	Moderate
C3	−1.284	0.277	Low
C6	−3.402	0.033	Very low
C7-C15	−6.298	0.002	Negligible
C16+	−23.776	~0	Anchors liquid phase

Bubble Point Solver

At bubble point the liquid composition equals the feed (x_i = z_i) and the incipient vapour satisfies:

Σᵢ Kᵢzᵢ = 1 Kᵢ = φᵢᴸ / φᵢᵛ

The solution procedure uses Wilson K-values as initial estimates, successive substitution to convergence, and Brent's method on the outer pressure loop.

Result

The residual scan confirms a sign change between 200 and 210 bar:

P (bar)	Residual
150	+0.149
170	+0.078
190	+0.024
200	+0.002
210	−0.018
230	−0.051

Brent's method converges to:

P_bub = 200.84 bar

Table 4 — K-values: Python vs commercial software

Component	z_i	y_i	K_i (Python)	K_i (Reference)
N₂	0.00410	0.01186	2.892	2.8917
CO₂	0.03799	0.05061	1.332	1.3322
C1	0.39916	0.79733	1.998	1.9975
C2	0.06072	0.06293	1.036	1.0364
C3	0.05449	0.03830	0.703	0.7029
C6	0.02342	0.00496	0.212	0.2116
C7-C15	0.19774	0.00829	0.042	0.0419
C16+	0.15903	~0	~0	4.49×10⁻⁶

All K-values match the commercial software output to 4 significant figures.

Summary of Implementation Details

Three details not explicitly stated in the software manual, all consistent with Peneloux (1982):

Detail	Standard assumption	Correct implementation
B in cubic solver	Use B_eff (shifted)	Use B_raw — per Peneloux (1982)
Z and V reported	Raw PR values	Shifted: subtract c_mix·P/RT
Volume shift in fugacity	Applied	Not applied — cancels in φᴸ/φᵛ

The Untuned Gap

The untuned EOS gives 200.84 bar. The lab measurement is 213.1 bar — a gap of 12.3 bar (6.1%). This is expected for a default characterisation with no regression. In Part 3 we close this gap by adjusting the critical properties of the pseudocomponents to match the lab bubble point.

Data and Code

Full Jupyter notebooks: github.com/eskoantg/PVTi_cross_check

Volve dataset: equinor.com/energy/volve-data-sharing

See Part 1 for data sources and disclaimer.

▶ Nomenclature (click to expand)

Symbol	Description	Units
P	Pressure	bar
T	Temperature	K
R	Universal gas constant (83.14472)	cm³·bar/(mol·K)
V	Molar volume	cm³/mol
Z	Compressibility factor	—
a_i	Pure-component attractive parameter	cm⁶·bar/mol²
b_i	Pure-component covolume	cm³/mol
c_i	Peneloux volume shift	cm³/mol
s_i	Dimensionless shift parameter	—
α_i	Temperature correction function	—
κ_i	PR kappa function	—
ω	Acentric factor	—
T_r	Reduced temperature = T/T_c	—
T_c	Critical temperature	K
P_c	Critical pressure	bar
Ω_A	EOS constant (0.45724)	—
Ω_B	EOS constant (0.077796)	—
a_mix	Mixture attractive parameter	cm⁶·bar/mol²
b_mix	Mixture covolume	cm³/mol
c_mix	Mixture volume shift	cm³/mol
A	Dimensionless mixture attractive parameter	—
B_raw	Dimensionless unshifted covolume	—
B_eff	Dimensionless shifted covolume	—
k_ij	Binary interaction parameter	—
da_i	Partial derivative of a_mix w.r.t. composition	cm⁶·bar/mol²
z_i	Overall mole fraction	—
x_i	Liquid mole fraction	—
y_i	Vapour mole fraction	—
K_i	K-value (equilibrium ratio)	—
φ_i	Fugacity coefficient	—