Petroleum starter - Nodal analysis

MODEL "Petroleum starter - Nodal analysis"
DESCRIPTION "Vogel inflow meets a quadratic outflow surrogate. Pressures are absolute psi, rate is STB/day; this is not a validated multiphase tubing model."
VERSION "1.0.0"
TAGS "petroleum", "starter", "nodal", "scenarios"

strict

# Synthetic saturated solution-gas-drive inflow and illustrative tubing inputs.
input ReservoirPressure_psia = 3000
  VALIDATE BETWEEN 100 AND 20000 MESSAGE "Reservoir pressure must be 100 to 20,000 absolute psi in this example."
input AOF_stbd = 1351.3513513513512
  VALIDATE BETWEEN 1 AND 100000 MESSAGE "AOF must be positive (1 to 100,000 STB/day)."
input StaticLiftPressure_psia = 500
  VALIDATE BETWEEN 0 AND 20000 MESSAGE "Static lift pressure must be 0 to 20,000 absolute psi."
input TubingCoefficient = 0.001
  VALIDATE BETWEEN 0 AND 1 MESSAGE "The quadratic coefficient must be nonnegative, in psi/(STB/day)^2."
input BracketLower_stbd = 0
  VALIDATE BETWEEN 0 AND 100000 MESSAGE "The lower bracket must be nonnegative."
input BracketUpper_stbd = 1300
  VALIDATE BETWEEN 0 AND 100000 MESSAGE "The upper bracket must be nonnegative."

# Enforce the inflow correlation's pressure domain inside the callable.
outflow_pressure = (rate) => StaticLiftPressure_psia + TubingCoefficient * rate^2
inflow_rate = (pressure) => pressure >= 0 AND pressure <= ReservoirPressure_psia THEN AOF_stbd * (1 - 0.2 * pressure / ReservoirPressure_psia - 0.8 * (pressure / ReservoirPressure_psia)^2) ELSE #N/A
balance = (rate) => inflow_rate(outflow_pressure(rate)) - rate
output OperatingRate_stbd IS number = ROOT(balance, BracketLower_stbd, BracketUpper_stbd, 1e-8, 200)
output OperatingPressure_psia IS number = outflow_pressure(OperatingRate_stbd)
output RateResidual_stbd IS number = balance(OperatingRate_stbd)
output LowerResidual_stbd = balance(BracketLower_stbd)
output UpperResidual_stbd = balance(BracketUpper_stbd)
output Assumptions = "Vogel saturated solution-gas-drive screening; quadratic outflow is an illustrative surrogate. Valid inflow pressure is 0 to reservoir pressure. Bracket endpoints outside that domain fail explicitly."
output ScenarioGuide = "Save Baseline. Change TubingCoefficient to 0.0005 and save Lower tubing loss. Restore Baseline or reset inputs; verify the residual and bracket."

# Plot only rates for which the pressure curves are physically in range.
ChartMaxRate_stbd = TubingCoefficient > 0 THEN MIN(AOF_stbd, SQRT(MAX(0, ReservoirPressure_psia - StaticLiftPressure_psia) / TubingCoefficient)) ELSE AOF_stbd
Rates_stbd = MAP(SEQUENCE(14, 1, 0, 1), (index) => ChartMaxRate_stbd * index / 13)
OutflowPressures_psia = MAP(Rates_stbd, (rate) => outflow_pressure(rate))
InflowPressures_psia = MAP(Rates_stbd, (rate) => ReservoirPressure_psia * (SQRT(0.04 + 3.2 * (1 - rate / AOF_stbd)) - 0.2) / 1.6)
Curve = MAP(Rates_stbd, (rate) => {
  rate: rate,
  outflow: outflow_pressure(rate),
  inflow: ReservoirPressure_psia * (SQRT(0.04 + 3.2 * (1 - rate / AOF_stbd)) - 0.2) / 1.6
})
A1 = "Operating oil rate (STB/day)"
B1 = OperatingRate_stbd
A2 = "Operating bottomhole pressure (absolute psi)"
B2 = OperatingPressure_psia
A3 = "Equation residual (STB/day)"
B3 = RateResidual_stbd
A4 = Assumptions

Charts_1!type = "charts"
Charts_1!config = <toml>
  version = 1
  kind = "charts"
  title = "Nodal analysis"

  [layout]
  title = "Synthetic well operating point"
  x_axis = "Oil rate (STB/day)"
  y_axis = "Bottomhole pressure (absolute psi)"

  [[charts]]
  type = "line"
  bind = "Curve"
  x = "rate"
  y = "outflow"
  label = "Quadratic outflow surrogate"

  [[charts]]
  type = "line"
  bind = "Curve"
  x = "rate"
  y = "inflow"
  label = "Vogel inflow"
</toml>

END MODEL