Exercise 8.4
Z-Factor Table via Dranchuk-Abou-Kassem
The OML 58 associated gas (gamma_g = 0.65) feeds the field's gas-lift and sales-gas accounting. Both jobs need the gas compressibility factor Z, the correction that turns the ideal-gas law into something that actually balances at reservoir and surface conditions. Z is never a single number: it lives on a P-T surface, and the engineer who books gas volumes needs to know where on that surface the gas behaves least ideally.
You'll build a Z-factor lookup table with the Dranchuk-Abou-Kassem (DAK) correlation, then find the worst-case (most non-ideal) corner of it.
The DAK machinery is already provided for you: pseudocritical_properties, z_factor_DAK, and z_factor(gamma_g, P_psia, T_F). Do not re-derive them. Your job is to drive them across a grid and reason about the result.
Write two functions:
z_table(pressures, temps_f): return adictkeyed by(P, T)tuples mapping
to the DAK Z at gas gravity GAMMA_G = 0.65. Build it over every pressure in pressures and every temperature in temps_f. Use (int(P), int(T)) for the keys so they look up cleanly.
max_deviation(table): scan the table and return the single
(P, T, Z) tuple whose Z is furthest from 1.0 (largest |Z - 1|). This is the corner where the gas is most non-ideal, where assuming Z = 1 would cost you the most barrels-equivalent.
Then evaluate over the field's screening grid:
PRESSURES = range(500, 6001, 500)(500, 1000, ... 6000 psia)TEMPS_F = [150, 200, 250](degrees Fahrenheit)
Store the built table in z_grid and the worst corner in worst as a (P, T, Z) tuple.
Physics to expect: along any temperature row, Z starts near 1 at low pressure, dips below 1 through the mid-pressure range as the gas compresses more than ideal, then climbs back above 1 at high pressure as repulsion dominates. The coldest row (150 F) dips the deepest (colder gas is more non-ideal), so the worst corner lands on the 150 F row.
Stuck? Reveal hints one at a time — they progress from nudge to near-solution.
visibilityReveal reference solutionexpand_more
Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.
import numpy as np
# OML 58 associated gas + the screening grid the gas-accounting team uses.
GAMMA_G = 0.65
PRESSURES = range(500, 6001, 500) # 500, 1000, ... 6000 psia
TEMPS_F = [150, 200, 250] # degrees Fahrenheit
# ── Provided DAK machinery - do NOT modify or re-derive these ────────────
def pseudocritical_properties(gamma_g):
Ppc = 756.8 - 131.0 * gamma_g - 3.6 * gamma_g ** 2
Tpc = 169.2 + 349.5 * gamma_g - 74.0 * gamma_g ** 2
return Ppc, Tpc
def z_factor_DAK(Ppr, Tpr, tol=1e-8, max_iter=100):
A1, A2, A3, A4 = 0.3265, -1.0700, -0.5339, 0.01569
A5, A6, A7, A8 = -0.05165, 0.5475, -0.7361, 0.1844
A9, A10, A11 = 0.1056, 0.6134, 0.7210
Z = 1.0
for _ in range(max_iter):
rho_r = 0.27 * Ppr / (Z * Tpr)
c1 = A1 + A2 / Tpr + A3 / Tpr ** 3 + A4 / Tpr ** 4 + A5 / Tpr ** 5
c2 = A6 + A7 / Tpr + A8 / Tpr ** 2
c3 = A9 * (A7 / Tpr + A8 / Tpr ** 2)
c4 = A10 * (1 + A11 * rho_r ** 2) * (rho_r ** 2 / Tpr ** 3) * np.exp(-A11 * rho_r ** 2)
f = Z - 1 - c1 * rho_r - c2 * rho_r ** 2 + c3 * rho_r ** 5 - c4
d_rho_dZ = -rho_r / Z
df = (1 - c1 * d_rho_dZ - 2 * c2 * rho_r * d_rho_dZ + 5 * c3 * rho_r ** 4 * d_rho_dZ
- A10 / Tpr ** 3 * d_rho_dZ * rho_r
* (2 * rho_r * (1 + A11 * rho_r ** 2) - 2 * A11 * rho_r ** 3) * np.exp(-A11 * rho_r ** 2))
Z_new = Z - f / df
if abs(Z_new - Z) < tol:
return Z_new
Z = Z_new
return Z
def z_factor(gamma_g, P_psia, T_F):
Ppc, Tpc = pseudocritical_properties(gamma_g)
return z_factor_DAK(P_psia / Ppc, (T_F + 460.0) / Tpc)
# ── Your work ─────────────────────────────────────────────────────────────
def z_table(pressures, temps_f):
"""dict (int P, int T) -> DAK Z at GAMMA_G over the P x T grid."""
table = {}
for P in pressures:
for T in temps_f:
table[(int(P), int(T))] = z_factor(GAMMA_G, P, T)
return table
def max_deviation(table):
"""Return the (P, T, Z) entry with the largest |Z - 1|."""
(P, T), Z = max(table.items(), key=lambda kv: abs(kv[1] - 1.0))
return (P, T, Z)
z_grid = z_table(PRESSURES, TEMPS_F)
worst = max_deviation(z_grid)
print("grid points:", len(z_grid))
print("most non-ideal corner (P, T, Z):", worst)
lockCopying code is a Full Access feature.