Gas Properties
Z-factor (Dranchuk-Abou-Kassem)
DAK eleven-parameter EOSZ = 1 + (A1 + A2/Tpr + A3/Tpr³ + A4/Tpr⁴ + A5/Tpr⁵)·ρr
+ (A6 + A7/Tpr + A8/Tpr²)·ρr²
− A9·(A7/Tpr + A8/Tpr²)·ρr⁵
+ A10·(1 + A11·ρr²)·(ρr²/Tpr³)·exp(−A11·ρr²)
where Tpr = T / Tpc, Ppr = P / Ppc, ρr = 0.27·Ppr/(Z·Tpr)
Pseudo-critical properties (Sutton, dry gas)
Tpc = 169.2 + 349.5·γg − 74.0·γg²
Ppc = 756.8 − 131.0·γg − 3.6·γg²
Viscosity (Lee-Gonzalez-Eakin)
μg = K · exp(X · ρg^Y)
K = (9.4 + 0.02·M)·T^1.5 / (209 + 19·M + T)
X = 3.5 + 986/T + 0.01·M
Y = 2.4 − 0.2·X
Pseudo-pressure (real gas)
m(p) = 2·∫[from p_ref to p] (p' / (μ(p')·Z(p'))) dp'
Integrated by Simpson's rule over 200 intervals, then stored as a lookup table.
Reservoir — Material Balance
Volumetric form (dry gas)
G·(Bg(p) − Bgi) = Gp·Bg(p) + We − Wp·Bw
For volumetric reservoirs (no aquifer): We = 0, Wp = 0. Solved bisectionally for p given Gp.
p/Z linearization
p/Z = (pi/Zi)·(1 − Gp/G)
Used for history matching: linear regression of observed (Gp, p/Z) yields G as the intercept.
Gas formation volume factor
Bg = 0.02827 · Z·T / p [rcf/scf]
Well — Inflow Performance
Darcy productivity index (pseudo-pressure)
q = J · [m(p̄) − m(pwf)]
J = k·h / (1422·T·[ln(re/rw) − 0.75 + S]) [Mscf/d / psi²·cp⁻¹]
Houpeurt deliverability (back-pressure)
q = C · (p̄² − pwf²)^n
C and n fitted from multi-rate test log-log regression. Equivalent J extracted at average drawdown for use in coupled forecast.
Surface Network
Weymouth pipe (horizontal gas)
p1² − p2² = (25 · γg · T̄ · Z̄ · L · f / d⁵) · q²
f = 0.032 / d^(1/3) [Weymouth friction]
q in MMscf/d
p in psia
L in ft
d in inches
T̄ in °R
Cullender-Smith VLP (gas tubing)
pwf² = pwh²·S + B·q²
S = exp(2·γg·H / (53.34·T̄·Z̄)) static head factor
B = 25·γg·T̄·Z̄·f·L_tubing / d_tubing⁵ friction
Network solver
Forward (flow): Kahn topological sort, sources → sinks
flow_downstream += Σ flow_upstream
Backward (press): reverse topo, sinks → sources
p_upstream = √(p_downstream² + K·q²)
Coupled Iteration
Nodal intersection (per well, per timestep)
repeat (damped fixed-point):
pwf = VLP(q, pwh) ← tubing
q' = J · [m(p̄) − m(pwf)] ← inflow
q = 0.5·(q + q') ← damping
until |q' − q| < tolerance
Outer network iteration
repeat (up to 6 outer iterations):
{pwh_i} = solveNetwork(model, {q_i})
for each well i:
q_i = NodalIntersection(reservoir, well_i, pwh_i)
update model average pressure via MBE
until max rate change < 1% of mean rate
All constants assume field units. The factor 25 in Weymouth/Cullender-Smith assumes q in MMscf/d; the implementation uses q in Mscf/d and divides q² by 10⁶ to compensate.