Generate a lift table with NeqSim

A lift table (Flow’s VFPPROD) tells the simulator which bottom-hole pressure a well needs to lift a given rate to the surface at a given tubing-head pressure (THP). You can type one in by hand (ps.lift.ProductionTable, see examples/lift_table.py), or let NeqSim compute it from the well’s tubing and the fluid it flows. ps.lift.generate does the second: for every point of the axes you give, NeqSim marches the stream up the tubing with a multiphase-flow correlation and PetroScript solves for the bottom-hole pressure that delivers that THP. What comes back is an ordinary ps.lift.ProductionTable.

You need five things:

A model the table is for: an ordinary model that checks clean, in METRIC or FIELD units. Its well gives the table its depth (ps.Well(..., bhp_reference_depth=...); a well that states none needs bhp_reference_depth= on generate), and its units are the table’s.
The stream, a ps.fluid.Compositional: the components, the equation of state and a named composition. It only describes what flows up the tubing; your model can stay black-oil.
The tubing, a ps.lift.Tubing: its stations, bore, roughness and temperature.
The axes: rates, THPs, water fractions and gas fractions, each in the model’s units.
A NeqSim runner, named by PETROSCRIPT_NEQSIM (Install and engines).

1. The model

Any model works, loaded or built. This one is a small dead-oil field with one producer; the producer’s bhp_reference_depth is the depth the table’s pressures refer to.

import petroscript as ps

m = ps.Model(
    name="LIFT_NEQSIM",
    units="metric",
    start="2025-01-01",
    grid=ps.grid.Cartesian(shape=(10, 10, 3), dx=300.0, dy=300.0, dz=20.0, top=2000.0),
    rock=ps.Rock(porosity=0.2, permeability=300.0, vertical_ratio=0.1,
                 compressibility=5.0e-5, reference_pressure=230.0),
    fluid=ps.fluid.OilWater(
        oil=ps.pvt.DeadOil(pressure=[50.0, 230.0, 400.0], oil_fvf=[1.12, 1.10, 1.08],
                           oil_viscosity=[1.8, 2.0, 2.3]),
        water=ps.pvt.Water(reference_pressure=230.0, fvf=1.01, compressibility=4.5e-5,
                           viscosity=0.4, viscosibility=0.0),
        oil_surface_density=850.0,
        water_surface_density=1020.0,
    ),
    relperm=[ps.relperm.OilWater(sw=[0.2, 0.4, 0.6, 0.8], krw=[0.0, 0.05, 0.2, 0.5],
                                 krow=[0.9, 0.4, 0.1, 0.0], pcow=[0.0, 0.0, 0.0, 0.0])],
    initial=ps.initial.Equilibrium(datum_depth=2000.0, datum_pressure=230.0,
                                   oil_water_contact=2100.0),
    wells=[ps.Well(name="PROD", i=5, j=5, bhp_reference_depth=2000.0)],
    schedule=ps.Schedule(
        end="2028-01-01",
        report="monthly",
        events=[ps.Event("2025-01-01",
                         ps.Perforate(well="PROD", k=(1, 3), diameter=0.2),
                         ps.Produce(well="PROD", liquid_rate=400.0, bhp_limit=50.0))],
    ),
)

2. The stream

The components’ constants are in the model’s units: molecular weight g/mol (lb/lb-mol in FIELD), critical temperature K (degR), critical pressure barsa (psia), critical volume m3/kg-mol (ft3/lb-mol). They come from the fluid’s characterisation (a PVT report, or published tables for pure components: methane is 16.04 g/mol, 190.6 K, 46.0 bar, 0.0986 m3/kg-mol and an acentric factor of 0.013; n-decane 142.29, 617.7, 21.1, 0.6 and 0.49). The composition’s mole fractions sum to 1. Water is a component too: a water-cut axis above 0 needs H2O in the composition, or NeqSim has no water to flow.

oil = ps.fluid.Compositional(
    equation_of_state="peng_robinson",             # or "soave_redlich_kwong"
    components=ps.pvt.Components(
        name=["C1", "C2", "C3", "NC4", "C7P", "C15P", "H2O"],
        molecular_weight=[16.043, 30.07, 44.097, 58.123, 100.0, 210.0, 18.015],
        critical_temperature=[190.56, 305.32, 369.83, 425.12, 540.0, 700.0, 647.3],
        critical_pressure=[45.99, 48.72, 42.48, 37.96, 27.4, 15.0, 220.89],
        critical_volume=[0.0986, 0.1455, 0.2, 0.255, 0.4, 0.8, 0.056],
        acentric_factor=[0.011, 0.099, 0.152, 0.2, 0.35, 0.65, 0.344],
    ),
    binary_interaction_coefficients={("C1", "C7P"): 0.03},
    compositions=[ps.pvt.Composition(
        name="RES",
        mole_fraction={"C1": 0.4, "C2": 0.05, "C3": 0.05, "NC4": 0.04, "C7P": 0.18,
                       "C15P": 0.18, "H2O": 0.1},
    )],
    surface_temperature=15.0,                      # degC: surface volumes are counted here
    surface_pressure=1.01325,                      # barsa
)

This fluid is passed to generate only. Making it the model’s own fluid would need the reservoir temperature (temperature=), and OPM Flow runs no equation-of-state model.

3. The tubing

Stations go down the string, the first at the wellhead; each segment between two stations has one bore and one roughness. Depths are in the model’s depth reference and units.

tubing = ps.lift.Tubing(
    measured_depth=[0.0, 2000.0],                  # m along the string
    true_vertical_depth=[0.0, 2000.0],             # m
    inner_diameter=[0.124],                        # m, a list: one per segment
    roughness=[1.5e-5],                            # m, a list: one per segment
    wellhead_temperature=40.0,                     # degC
    temperature_gradient=0.03,                     # degC warmer per m down
)

4. Generate the table

Each axis rises and is in the model’s units. rate_kind says which rate the rate axis is ("oil", "liquid", "gas" or "water"); the water fraction is a water cut (a fraction) unless water_fraction_kind says otherwise, and the gas fraction a gas-oil ratio (sm3/sm3 in METRIC, Mscf/stb in FIELD) unless gas_fraction_kind says otherwise. Left out, the THP axis is sized from the model: six values from the lowest THP it states for the well up to its initial pressure.

table = ps.lift.generate(
    tubing,                                        # the pipe the stream flows up
    model=m, well="PROD",
    fluid=oil, composition="RES",
    correlation="beggs_brill",                     # or "hagedorn_brown", "gray", "mukherjee_brill"
    number=1,                                      # the table's number in the deck
    rate=[50.0, 200.0, 400.0, 800.0], rate_kind="liquid",   # sm3/day of liquid
    thp=[10.0, 20.0, 40.0, 80.0],                  # barsa
    water_fraction=[0.0, 0.5],                     # water cut
    gas_fraction=[0.0],                            # gas-oil ratio, sm3/sm3
    folder="lift-neqsim",                          # NeqSim's request, result and log go here
)
print(table.bhp.shape)      # (4, 2, 1, 1, 4): thp, water fraction, gas fraction, lift, rate
print(table.frame())        # one row per point, in long form, with its bhp
print(table.source)         # NeqSim's version, the correlation, the composition, the datum, ...

Choosing the correlation:

beggs_brill takes any inclination and keeps water beside oil and gas. It leaves the kinetic-energy (momentum) term out, so it refuses a point where that term would matter: the acceleration rule, Ek above 0.05, met where the gas moves fast (a high gas fraction, a high rate, a low THP).
hagedorn_brown and gray take only tubing within 15 degrees of vertical, and drop the water where gas, oil and water flow together. They keep the momentum term where liquid flows, but leave it out where gas flows alone, so the acceleration rule applies to them there too.
mukherjee_brill leaves the momentum term out, like beggs_brill, and drops that water too.

Measured on NeqSim 3.25.0, 2,500 m of vertical 0.1 m tubing: a 50/50 methane and n-decane oil at 200 sm3/sm3 and THP 20 bar, where beggs_brill refused 2 of 18 points by the acceleration rule, filled all 18 with hagedorn_brown; pure methane at THP 10 bar was refused by all three. For oil with gas, hagedorn_brown can fill what beggs_brill refuses; for gas alone, change the axes.

5. Use it

The table is like any other: the model holds it at lift_tables/<number>, and a control names the table itself. Here the producer keeps its rate until its THP falls to 70 barsa:

m = m.set({
    "lift_tables/1": table,
    "schedule/2025-01-01/PROD/control": ps.Produce(
        well="PROD", liquid_rate=400.0, bhp_limit=50.0, thp_limit=70.0, lift_table=table),
})
deck = ps.compile(m)                               # the deck writes VFPPROD 1 from the table
assert "VFPPROD" in deck.files["LIFT_NEQSIM.DATA"]

From here the model is checked, compiled and run like any other (deck.run()). examples/lift_neqsim.py is this page as one script, run in OPM Flow, printing the date the THP limit takes over.

When NeqSim refuses points

A table needs every entry, so generate raises ps.RunError when NeqSim refuses any point, and lists each refused point with its axis values and the reason:

the fluid makes no aqueous phase at standard conditions: water-cut values above 0 with no H2O in the composition.
the tubing cannot carry the point: the acceleration rule: ... (choked) (the flow reaches sonic speed at the wellhead: too high a rate for the bore at that THP), or the well cannot lift this point (even a bottom-hole pressure 1000 bar above the THP does not deliver it). With the THP axis stated, NeqSim knows nothing of the reservoir: it searches the bottom-hole pressure between the THP and THP + 1000 bar. Lower the top rate, widen the bore, or raise the THP.
no convergence ... or the march ... fails: the correlation’s pressure march jumps at this point (a flow-regime boundary, a dew point), so no bottom-hole pressure in the search matches the THP. Move the axis value a little, or try another correlation.
the acceleration rule: Ek ... > 0.05 where beggs_brill leaves the momentum term out: the gas moves too fast at that point for a correlation without the momentum term (a high gas fraction, a high rate or a low THP). Lower the gas-fraction or rate axis, or raise the lowest THP. For oil with gas in vertical tubing, hagedorn_brown may fill the point; for gas alone it refuses it too.
the three-phase rule: hagedorn_brown, gray or mukherjee_brill with gas, oil and water flowing together. Use beggs_brill, within its acceleration rule.

Other refusals come before NeqSim runs, each a plain ValueError: a well the model does not hold, a well with no bhp_reference_depth (state it on the well, or bhp_reference_depth= on generate), no THP axis and nothing to size one from (state thp=), and a LAB model.