Exercise 15.2
Gas Property Module - Z, Bg, Density, Viscosity & Compressibility Table
Package the chapter's property functions into a single gas_properties(P, T, gamma_g) function that returns a dict of {Z, Bg, density, viscosity, cg} in one call, then add a gas_property_table(P_range, T, gamma_g) wrapper that returns a formatted DataFrame. Validate against published data for a 0.7 gravity gas at 250°F.
---
You'll assemble the chapter's gas-PVT toolkit into a single property table. The six verified correlations are embedded for you; do not modify them or re-derive the physics. The Sutton pseudo-criticals, the Dranchuk–Abou-Kassem Z-factor, gas_fvf, gas_density, the Lee–Gonzalez–Eakin viscosity, and the numerical gas_compressibility are all in scope already.
Your job is the one missing piece: the table builder.
Write gas_property_table(P_range, T_F, gamma_g) that returns a pandas.DataFrame with one row per pressure in P_range and these columns, in this exact spelling and order:
'P (psia)', 'Z', 'Bg (rcf/scf)', 'rho (lb/ft3)', 'mu (cp)', 'cg (1/psi)'For each pressure P:
- Get
Ppc, Tpc = pseudo_critical_properties(gamma_g). Tpr = (T_F + 459.67) / Tpc.Z = z_factor_DAK(P / Ppc, Tpr).Bg = gas_fvf(P, T_F, Z)→ column'Bg (rcf/scf)'rho = gas_density(P, T_F, Z, gamma_g)→ column'rho (lb/ft3)'mu = gas_viscosity_lee(P, T_F, Z, gamma_g)→ column'mu (cp)'cg = gas_compressibility(P, T_F, Z, gamma_g)→ column'cg (1/psi)'
Store the raw P in 'P (psia)' and the same Z you computed in 'Z'. Keep the numbers as floats (no rounding/formatting; the grader reads the raw values).
The embedded constants are the chapter's published validation case (a 0.70-gravity gas at 250 °F):
GAMMA_G = 0.70
T_F = 250.0
P_RANGE = [2000.0, 4000.0, 6000.0] # psiaThen build the table and expose these exact names the tests read:
props = gas_property_table(P_RANGE, T_F, GAMMA_G) # the DataFrame
z_2000 = float(props.loc[0, 'Z'])
bg_2000 = float(props.loc[0, 'Bg (rcf/scf)'])
rho_4000 = float(props.loc[1, 'rho (lb/ft3)'])> Think about it: as you walk up the pressure range, Bg should fall > (gas compresses into less reservoir volume per scf) while rho rises. The > ideal-gas skeleton Bg = 0.02827 · Z · T_R / P means Bg · P / Z is the same > constant on every row, a quick sanity check that your Z is threaded through > correctly. Why does cg shrink toward roughly 1/P as pressure climbs?
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
import pandas as pd
from scipy.optimize import brentq
# ── Verified gas-PVT correlations from Chapter 15 (do not edit) ──────────
def pseudo_critical_properties(gas_sg):
"""
Estimate pseudo-critical properties from gas specific gravity.
Uses the Sutton (1985) correlations.
Parameters
----------
gas_sg : float
Gas specific gravity (air = 1.0). Typical range: 0.55 to 0.90.
Returns
-------
tuple
(Ppc in psia, Tpc in deg R)
"""
Ppc = 756.8 - 131.0 * gas_sg - 3.6 * gas_sg**2
Tpc = 169.2 + 349.5 * gas_sg - 74.0 * gas_sg**2
return Ppc, Tpc
def z_factor_DAK(Ppr, Tpr):
"""
Calculate gas Z-factor using the Dranchuk-Abou-Kassem (1975)
equation of state.
This correlation fits the Standing-Katz chart with 11 empirical
coefficients. Valid for 1.0 <= Tpr <= 3.0 and 0.2 <= Ppr <= 15.
Parameters
----------
Ppr : float
Pseudo-reduced pressure (dimensionless).
Tpr : float
Pseudo-reduced temperature (dimensionless).
Returns
-------
float
Compressibility factor Z.
"""
# DAK (1975) coefficients
A1, A2, A3 = 0.3265, -1.0700, -0.5339
A4, A5, A6 = 0.01569, -0.05165, 0.5475
A7, A8 = -0.7361, 0.1844
A9, A10, A11 = 0.1056, 0.6134, 0.7210
def residual(rho_r):
"""DAK residual function. Z = 0.27 * Ppr / (rho_r * Tpr)"""
Z = 0.27 * Ppr / (rho_r * Tpr)
F = (
A1 + A2/Tpr + A3/Tpr**3 + A4/Tpr**4 + A5/Tpr**5
) * rho_r + (
A6 + A7/Tpr + A8/Tpr**2
) * rho_r**2 - A9 * (
A7/Tpr + A8/Tpr**2
) * rho_r**5 + A10 * (
1 + A11*rho_r**2
) * (rho_r**2 / Tpr**3) * np.exp(-A11*rho_r**2) + 1 - Z
return F
# Solve for reduced density
rho_r = brentq(residual, 0.001, 5.0)
Z = 0.27 * Ppr / (rho_r * Tpr)
return Z
def gas_fvf(P_psia, T_F, Z):
"""Gas formation volume factor (rcf/scf -- reservoir ft^3 per scf).
The 0.02827 constant gives reservoir CUBIC FEET per scf. For reservoir
bbl/scf, divide by 5.615 (i.e. use 0.005035).
"""
T_R = T_F + 459.67
return 0.02827 * Z * T_R / P_psia
def gas_density(P_psia, T_F, Z, gas_sg):
"""Gas density (lb/ft3)."""
T_R = T_F + 459.67
return 2.7 * gas_sg * P_psia / (Z * T_R)
def gas_viscosity_lee(P_psia, T_F, Z, gas_sg):
"""
Gas viscosity using Lee-Gonzalez-Eakin (1966) correlation.
Returns viscosity in centipoise (cp).
"""
T_R = T_F + 459.67
Mg = 28.97 * gas_sg # molecular weight
rho_g = gas_density(P_psia, T_F, Z, gas_sg) # lb/ft3
rho_gcc = rho_g / 62.428 # convert to g/cc
K = (9.4 + 0.02 * Mg) * T_R**1.5 / (209 + 19*Mg + T_R)
X = 3.5 + 986/T_R + 0.01*Mg
Y = 2.4 - 0.2*X
return 1e-4 * K * np.exp(X * rho_gcc**Y)
def gas_compressibility(P_psia, T_F, Z, gas_sg, dP=10):
"""
Gas isothermal compressibility cg (1/psi).
Calculated numerically from Z(P).
"""
Ppc, Tpc = pseudo_critical_properties(gas_sg)
T_R = T_F + 459.67
Tpr = T_R / Tpc
Ppr_plus = (P_psia + dP) / Ppc
Ppr_minus = (P_psia - dP) / Ppc
try:
Z_plus = z_factor_DAK(Ppr_plus, Tpr)
Z_minus = z_factor_DAK(Ppr_minus, Tpr)
dZ_dP = (Z_plus - Z_minus) / (2 * dP)
cg = 1/P_psia - 1/Z * dZ_dP
except (ValueError, RuntimeError):
cg = 1/P_psia # ideal gas approximation as fallback
return cg
# ── Validation constants: 0.70-gravity gas at 250 deg F (do not edit) ────
GAMMA_G = 0.70 # gas specific gravity (air = 1.0)
T_F = 250.0 # reservoir temperature, deg F
P_RANGE = [2000.0, 4000.0, 6000.0] # pressures, psia
def gas_property_table(P_range, T_F, gamma_g):
"""Build a gas-property table, one row per pressure.
Columns (exact spelling/order):
'P (psia)', 'Z', 'Bg (rcf/scf)', 'rho (lb/ft3)', 'mu (cp)', 'cg (1/psi)'
"""
Ppc, Tpc = pseudo_critical_properties(gamma_g)
T_R = T_F + 459.67
Tpr = T_R / Tpc
rows = []
for P in P_range:
Z = z_factor_DAK(P / Ppc, Tpr)
rows.append({
"P (psia)": P,
"Z": Z,
"Bg (rcf/scf)": gas_fvf(P, T_F, Z),
"rho (lb/ft3)": gas_density(P, T_F, Z, gamma_g),
"mu (cp)": gas_viscosity_lee(P, T_F, Z, gamma_g),
"cg (1/psi)": gas_compressibility(P, T_F, Z, gamma_g),
})
return pd.DataFrame(rows)
props = gas_property_table(P_RANGE, T_F, GAMMA_G)
z_2000 = float(props.loc[0, "Z"])
bg_2000 = float(props.loc[0, "Bg (rcf/scf)"])
rho_4000 = float(props.loc[1, "rho (lb/ft3)"])
print(f"GAS PROPERTY TABLE - gamma_g = {GAMMA_G}, T = {T_F} deg F")
print(props.to_string(index=False))
print()
print("Z @ 2000 psia :", z_2000)
print("Bg @ 2000 psia :", bg_2000, "rcf/scf")
print("rho @ 4000 psia :", rho_4000, "lb/ft3")
lockCopying code is a Full Access feature.