A Small Whole Refinery#
This notebook strings the difflow_refinery units together into a small refinery, from the
crude assay to the finished-product pools and the hydrogen header:
+-> offgas + unstabilised naphtha -> gas plant -+-> fuel gas
| +-> LPG
| +-> stabilised naphtha
| |
crude -> CDU ------+ naphtha hydrotreater (NHT)
| |
| naphtha fractionator
| | |
| light naphtha heavy naphtha -> reformer -+-> reformate
| | +-> net H2 --+
| +--------> gasoline pool <---------+ |
| v
+-> kerosene + diesel + AGO -> distillate hydrotreater (DHT) -> fractionator -+-> jet hydrogen
| +-> ULSD header
+-> atmospheric residue -> residue desulfurizer (RDS) ------------------------> fuel oil (+ H2 plant)
Every box is a unit of the plugin, solved in its own way:
Unit |
Class |
Model |
How far it is validated |
|---|---|---|---|
Crude unit (CDU) |
|
furnace + equation-oriented MESH column, side strippers, pumparounds |
against an independent Pyomo/IPOPT column and IDAES property packages |
Gas plant |
|
Peng-Robinson stage columns, O’Connell tray efficiencies |
debutanizer against IDAES |
Naphtha and distillate hydrotreaters |
|
trickle beds (LHHW HDS by sulfur class, HDN, aromatics saturation), PR separator, H2 recycle, steam stripper, charge heater |
balances and gradients only; rate constants illustrative |
Product fractionators |
|
smooth TBP split (a sigmoid of 8 K width), not a column |
balances only |
Catalytic reformer |
|
three adiabatic beds, equilibrium from Gibbs energies, PR separator, H2 recycle, trace sulfur |
balances and gradients only; rate constants illustrative |
Residue desulfurizer |
|
trickle beds (HDS by class, HDM, CCR reduction, conversion), once-through treat gas, ideal product split |
balances and gradients only; rate constants illustrative |
Hydrogen header |
|
header balance with a swing hydrogen plant; makeup purity fed back into the hydrotreaters |
balances only |
Product pools |
|
Ethyl RT-70 octane, RVP index, mass-basis sulfur, cetane index from composition; flash, freeze, smoke, viscosity and straight-run octane estimated |
rules against published worked examples; the estimates mostly unverified |
How to read the numbers. The balances, the column equations and the equilibrium
thermochemistry are physics. The kinetic constants of the hydrotreaters, the reformer and
the residue desulfurizer are illustrative: they give the right trends and exact gradients,
but the yields, the product sulfur and the hydrogen numbers are the shape of the answer, not
a prediction for any catalyst. So are the property estimates of section 8, which are marked
unverified in difflow_refinery.properties. The assay is the made-up test crude used
throughout the plugin’s examples (examples/35, 38), given 1.8 wt% sulfur, 1500 wppm
nitrogen, 5 wt% CCR and (new here, for the residue desulfurizer) 40 wppm Ni+V.
What is still assumed (each says so where it is used): the 2 mol % H2S in the crude unit’s offgas; the fractionators’ sigmoid width; the residue desulfurizer’s once-through treat gas; the iso/normal splits of the reformer feed; the hydrogen plant as an unlimited swing source. The last section lists what the earlier version of this notebook had to write as stand-ins, and which issue replaced each.
The run is sized for under twenty minutes on four cores (17.6 min when this copy was executed): a coarse cut grid (30 C cuts below 400 C), short gas-plant columns, single-bed naphtha hydrotreater. Most of the time is compilation: each unit compiles once, and the cross-unit gradient of section 10 once more.
import dataclasses
import time
import warnings
import jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp
import numpy as np
from IPython.display import Markdown, display
from tabulate import tabulate
import difflow_refinery as dr
import difflow_refinery.hydrogen as h2
from difflow_refinery import Assay, BlendCharacterization, BlendComponent, BlendPool, column as cc
from difflow_refinery.composition import CompositionRangeWarning
from difflow_refinery.gasplant import (
AmineTreater, AmineTreaterParams, GasCompressor, GasCompressorParams, GasPlantColumn,
absorber_deethanizer, debutanizer, fuel_gas, gas_plant_feed, lpg_quality)
from difflow_refinery.hydroprocessing.reactor import ReactorOptions
from difflow_refinery.hydrotreating import NAPHTHA_HDT_PARAMS, Hydrotreater, HydrotreaterParams
from difflow_refinery.plant import Chain, Stage, central_difference
from difflow_refinery.reforming import CatalyticReformer, NaphthaFeed, ReformerParams
from difflow_refinery.reforming.thermo import std_liquid_volume
from difflow_refinery.residue import ResidueDesulfurizer, fuel_oil_blend
T_START = time.time()
C = 273.15
BBL_D = 86400.0 / cc.BARREL # m3/s -> bbl/d
MW_S = 32.065e-3 # kg/mol
def table(rows, headers, floatfmt=".4g"):
display(Markdown(tabulate(rows, headers, tablefmt="github", floatfmt=floatfmt)))
def tick(label, t0):
print(f"{label}: {time.time() - t0:.0f} s")
1. Crude and crude unit#
The crude and the column are those of examples/38_refinery_gas_plant.ipynb: 95 000 bbl/d,
a 30-stage atmospheric column with kerosene, diesel and AGO side strippers, two pumparounds,
a 5 % overflash closing the furnace, and a partial condenser held at 40 C so that the
column makes an offgas for the gas plant. The cut grid is coarser than the default (30 C
cuts to 400 C instead of 20 C), which gives 17 pseudo-components and shortens every
compile downstream.
characterize(..., composition=True) adds the hydrocarbon types, hydrogen and the sulfur
and nitrogen classes per pseudo-component that the hydrotreaters, the reformer and the
residue desulfurizer need. The crude unit characterizes the same assay on the same cut
points itself, so its product streams are on the same component names. The Ni+V only
changes the residue (the metals sit in the heaviest cuts).
PCT = [0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100]
T_C = [-10, 60, 95, 150, 205, 260, 315, 370, 430, 500, 600, 680, 850]
assay = Assay(PCT, [t + C for t in T_C], sg=0.86,
light_ends={"ethane": 0.05, "propane": 0.5, "isobutane": 0.3,
"n_butane": 1.0, "isopentane": 0.8, "n_pentane": 1.5},
sulfur_wt=1.8, nitrogen_wppm=1500.0, ccr_wt=5.0, nickel_vanadium_wppm=40.0)
cut_points = dr.default_cut_points(assay, ((673.15, 30.0), (873.15, 60.0), (np.inf, 150.0)))
with warnings.catch_warnings():
warnings.simplefilter("ignore", CompositionRangeWarning)
char = dr.characterize(assay, cut_points, composition=True)
BPD = 95_000.0
Vf = BPD * cc.BARREL / 86400.0
cdu_params = cc.CrudeColumnParams(
n_stages=30, feed_stage=27, P_top=1.5e5, P_bottom=1.9e5, P_condenser=1.3e5,
bottom_steam=150.0, steam_T=C + 260.0, condenser="partial",
side_products=(cc.SideProduct("kero", 9, 4, steam=40.0),
cc.SideProduct("diesel", 16, 4, steam=40.0),
cc.SideProduct("ago", 22, 3, steam=20.0)),
pumparounds=(cc.Pumparound("pa1", 12, 10), cc.Pumparound("pa2", 19, 17)),
specs=(cc.product_rate("naphtha", 0.20 * Vf), cc.product_rate("kero", 0.11 * Vf),
cc.product_rate("diesel", 0.17 * Vf), cc.product_rate("ago", 0.05 * Vf),
cc.pumparound_duty("pa1", 15e6), cc.pumparound_delta_t("pa1", 60.0),
cc.pumparound_duty("pa2", 20e6), cc.pumparound_delta_t("pa2", 60.0),
cc.overflash(0.05), cc.stage_temperature(0, C + 40.0)),
)
t0 = time.time()
cdu_unit = dr.CrudeUnit(assay, cdu_params, cut_points=cut_points)
assert tuple(cdu_unit.crude.names) == tuple(char.names)
cdu = cdu_unit.solve(BPD, T=C + 240.0, P=6e5)
tick(f"CDU converged: {bool(cdu.converged)}", t0)
print(cdu.table())
P = {name: {k: v for k, v in s.items()} for name, s in cdu.products.items()}
CDU converged: True: 47 s
product bbl/d vol % wt % API TBP5 C TBP50 C TBP95 C
offgas 136 0.1 0.1 122.3 -89 -10 89
naphtha 19000 20.0 16.7 65.4 -5 96 162
kero 10450 11.0 10.2 46.7 126 183 236
diesel 16150 17.0 16.5 37.6 195 260 326
ago 4750 5.0 5.1 30.6 271 333 383
residue 44514 46.9 51.4 18.4 322 476 759
coil outlet 313.9 C, 69.7 mol% vaporised; furnace 42.8 MW absorbed, 50.3 MW fired
condenser 31.9 MW, pumparounds 15.0, 20.0 MW
2. Gas plant#
The offgas and the unstabilised naphtha go to the saturated gas plant as in
examples/38_refinery_gas_plant.ipynb: a two-stage compressor to 14.5 bar, an amine
treater on the compressed vapour, an absorber-deethanizer with the whole naphtha as lean oil
(overhead: fuel gas), and a debutanizer (overhead: LPG; bottoms: stabilised naphtha). The
columns are shorter than in example 38 (12 and 16 trays) to save compile time.
gas_plant_feed (#326) puts the crude-unit streams on the gas plant’s component table:
it keeps the cuts that carry more than 0.1 % of the naphtha, folds the rest into the
heaviest kept cut (moles conserved; the mass change is reported), and drops the water the
accumulator decants.
Still assumed: the offgas H2S. The crude unit makes no H2S (the assay’s sulfur stays on
the cuts), so the offgas is given 2 mol % H2S (h2s={"offgas": 0.02}) for the amine to
remove. evolved_h2s(products, char, fraction) would put a sulfur balance behind it, but
the fraction of a crude’s sulfur that evolves in the furnace is crude-specific and no value
is sourced, so it would be an assumption too.
LIGHT = ["hydrogen_sulfide", "ethane", "propane", "isobutane", "n_butane", "isopentane", "n_pentane"]
gp = gas_plant_feed(P, char, LIGHT, min_fraction=1e-3, h2s={"offgas": 0.02}, # H2S ASSUMED
T=C + 40.0, P=1.3e5)
print(gp.summary())
comps = gp.components
gp_flows = gp.flows
def gp_mass(s): # kg/s
return float(gp.mass(s))
offgas, naphtha, H2S = gp["offgas"], gp["naphtha"], float(gp.h2s_mol["offgas"])
t0 = time.time()
comp_gas, condensate, ci = GasCompressor(GasCompressorParams(comps, outlet_P=14.5e5, n_stages=2))(offgas)
sweet, acid, ai = AmineTreater(AmineTreaterParams())(comp_gas)
gp_feed = {k: sweet[k] + condensate[k] for k in sweet if k.startswith("F_")}
gp_feed.update(T=C + 40.0, P=14.5e5)
deeth = GasPlantColumn(absorber_deethanizer(comps, n_trays=12, feed_tray=4, c2_in_bottoms=0.002))
fgas, deeth_btms, di = deeth(dict(naphtha, P=14.5e5), gp_feed)
tick(f"deethanizer converged: {bool(di['converged'])}", t0)
t0 = time.time()
debut = GasPlantColumn(debutanizer(comps, n_trays=16, feed_tray=8, c5_in_lpg=0.01, c4_in_naphtha=0.01))
lpg, stab, bi = debut(deeth_btms)
tick(f"debutanizer converged: {bool(bi['converged'])}", t0)
od, ob = di["outputs"], bi["outputs"]
q = lpg_quality(gp_flows(lpg), comps, grade="commercial_propane")
fg = fuel_gas(gp_flows(fgas), comps)
table([("fuel gas (kg/s, LHV MJ/kg)", f"{gp_mass(fgas):.3f}, {float(fg['lhv_mass']) / 1e6:.1f}"),
("LPG (kg/s); propane / butanes+ (mol frac)",
f"{gp_mass(lpg):.2f}; {float(q['values']['propane']):.2f} / {float(q['values']['butanes_plus']):.2f}"),
("stabilised naphtha (kg/s), RVP (kPa)", f"{gp_mass(stab):.2f}, {float(ob['bottoms.rvp']) / 1e3:.1f}"),
("compressor (MW)", f"{float(ci['power']) / 1e6:.3f}"),
("deethanizer / debutanizer reboiler (MW)",
f"{float(od['reboiler.duty']) / 1e6:.2f} / {float(ob['reboiler.duty']) / 1e6:.2f}")],
["gas plant", "value"])
into = gp_flows(offgas) + gp_flows(naphtha)
out = sum(gp_flows(s) for s in (fgas, acid, lpg, stab))
print(f"largest component imbalance across the gas plant: {float(jnp.max(jnp.abs(out - into))):.1e} mol/s")
components: hydrogen_sulfide, ethane, propane, isobutane, n_butane, isopentane, n_pentane, pc01, pc02, pc03, pc04
cuts kept (naphtha > 0.001): pc01..pc04; folded into pc04: pc05, pc06, pc07, pc08, pc09, pc10, pc11, pc12, pc13, pc14, pc15, pc16, pc17
folded: 4.9e-05 of the naphtha (molar)
offgas: folded 2.212e-07 mol/s (3.309e-08 kg/s); mass change +7.83e-07 kg/s; dropped 0.00279 kg/s; H2S added 0.05143 mol/s
naphtha: folded 1.342e-02 mol/s (2.007e-03 kg/s); mass change -2.67e-04 kg/s; dropped 0 kg/s
deethanizer converged: True: 54 s
debutanizer converged: True: 40 s
| gas plant | value | |——————————————-|——————-| | fuel gas (kg/s, LHV MJ/kg) | 0.021, 46.8 | | LPG (kg/s); propane / butanes+ (mol frac) | 1.65; 0.32 / 0.66 | | stabilised naphtha (kg/s), RVP (kPa) | 23.56, 32.3 | | compressor (MW) | 0.019 | | deethanizer / debutanizer reboiler (MW) | 8.41 / 11.70 |
largest component imbalance across the gas plant: 2.8e-13 mol/s
3. Naphtha hydrotreater#
The stabilised naphtha goes to a Hydrotreater set up as a naphtha unit at a naphtha
unit’s conditions: one bed, 320 C inlet, 30 bar, LHSV 4 h-1, 100 Nm3/m3 of treat gas, the
feed entering the charge heater at 250 C after a feed/effluent exchanger, and a cooler
stripper (150 C feed; the default 230 C does not converge on a naphtha, and says so). The
kinetics are NAPHTHA_HDT_PARAMS (#332), an illustrative naphtha set: the diesel
constants with the light sulfur classes and HDN faster and aromatics saturation slower.
Two changes of #332 are why this works at these conditions: the reactor’s phase model in an all-vapour bed (the earlier version of this notebook saw the aromatics equilibrium run backwards at 30 bar, which was the phase model, not the chemistry), and the charge heater, whose fired duty now goes in the energy table.
The unit’s beds are integrated with ReactorOptions(adjoint="forward")
(diffrax.ForwardMode): the same values, but forward-mode derivatives instead of reverse,
which is what lets section 10 differentiate it together with the reformer in one
jax.jacfwd (see difflow_refinery.plant).
nht_feed = {k: float(v) for k, v in stab.items() if k.startswith("F_") and k[2:] in char.names}
nht_feed.update(T=C + 40.0, P=3e5)
nht_params = HydrotreaterParams(T_in=(C + 320.0,), quench=None, bed_fractions=(1.0,), P=30e5, lhsv=4.0,
h2_oil=100.0, stripper_feed_T=C + 150.0, heater_inlet_T=C + 250.0,
kinetics=NAPHTHA_HDT_PARAMS, reactor=ReactorOptions(adjoint="forward"))
t0 = time.time()
nht = Hydrotreater(char, nht_feed, nht_params)
nht_res = nht.solve(nht_feed)
tick(f"naphtha hydrotreater converged: {bool(nht_res.converged)}", t0)
print(nht_res.table())
print(f"charge heater {float(nht_res.outputs['heater.fired_duty']) / 1e6:.1f} MW fired, "
f"feed/effluent exchanger {float(nht_res.outputs['feed_effluent.duty']) / 1e6:.1f} MW")
print("closures:", {k: f"{float(v):.0e}" for k, v in nht_res.balances.items()})
naphtha hydrotreater converged: True: 86 s
WABT 321.4 C, bed dT 2.1 K
product: S 0.0 wppm, N 0.3 wppm, SG 0.7491, aromatics 19.8 vol%, cetane index 27.8
yields (mass): product 82.35 %, wild naphtha 17.10 %, gas 0.566 %
H2: chemical 1.7 Nm3/m3 (10 scf/bbl), makeup 11.0 Nm3/m3, purity 97.3 %
recycle compressor 50.3 kW, purge 6.913 mol/s
charge heater 5.3 MW fired, feed/effluent exchanger 19.1 MW
closures: {'C': '1e-16', 'H': '1e-16', 'N': '4e-16', 'S': '3e-16', 'mass': '1e-16'}
At 320 C the sulfur is removed (0.0002 wppm; with the illustrative first-order classes, a figure far below 1 wppm only says “removed”) and the nitrogen is down to 0.3 wppm, both under a reformer’s 0.5 wppm. The bed rises 2.1 K on 1.7 Nm3/m3 of chemical hydrogen, the gas yield is positive (0.57 %), and with the feed/effluent exchanger doing 19.1 MW the charge heater fires 5.3 MW.
The naphtha split#
The light naphtha (C5s and the first cut, holding the benzene precursors) goes to gasoline
and only the heavy naphtha is reformed. res.fractionate(...) (#328) splits the stripper
bottoms and the wild naphtha together at 85 C TBP on the hydrotreater’s treated product
grid, with a smooth sigmoid (8 K wide, illustrative) rather than a column, and sends the
H2, H2S, NH3 and C1-C2 dissolved in the wild naphtha to an off-gas, so the products are
blendable and the mass closes. A rigorous alternative is a gas-plant splitter on
gasplant.hydroprocessed_feed(...) (#327; one more column compile).
NAPHTHA_CUT = C + 85.0
def naphtha_split(res, cut=NAPHTHA_CUT):
return res.fractionate(cut_points=(cut,), products=("light_naphtha", "heavy_naphtha"),
feeds=("product", "wild_naphtha"))
fr = naphtha_split(nht_res)
print(fr.table())
light_naphtha: 8.2534 kg/s (35.21 wt%)
heavy_naphtha: 15.1792 kg/s (64.75 wt%)
off_gas: 0.0110 kg/s (0.05 wt%)
mass closure 1.5e-16
4. Catalytic reformer#
NaphthaFeed.from_hydrotreater(fr, "heavy_naphtha") (#327) maps the heavy naphtha onto the
reformer’s P/N/A lumps by carbon number from the treated grid: its molar masses (so the
feed’s mass is the fractionator’s to round-off), its hydrocarbon types after saturation,
and its sulfur, which the reformer carries as a trace element to the reformate and, as
H2S, the net gas and the fuel gas (#330; the conversion to H2S is illustrative). The
iso/normal and MCP/cyclohexane splits are from_characterization’s illustrative defaults.
ref_feed = NaphthaFeed.from_hydrotreater(fr, "heavy_naphtha")
g = ref_feed.group_fractions("volume")
print(f"reformer feed {float(ref_feed.mass_flow):.3f} kg/s (fractionator heavy naphtha "
f"{float(fr.rates['heavy_naphtha']):.3f}): P/N/A {100 * float(g['paraffins']):.0f}/"
f"{100 * float(g['naphthenes']):.0f}/{100 * float(g['aromatics']):.0f} vol%, "
f"N+2A {float(ref_feed.n_plus_2a()):.0f}, S {float(ref_feed.sulfur_wppm):.3g} wppm")
t0 = time.time()
reformer = CatalyticReformer(ReformerParams()) # WAIT 500 C, 12 bar, H2/HC 5, LHSV 1.5
ref = reformer.solve(ref_feed)
tick(f"reformer converged: {bool(ref.converged)} in {ref.iterations} recycle iterations", t0)
print(ref.summary())
print("closures:", {k: f"{float(v):.0e}" for k, v in ref.balances().items()})
reformer feed 15.179 kg/s (fractionator heavy naphtha 15.179): P/N/A 56/23/21 vol%, N+2A 64, S 0.000219 wppm
reformer converged: True in 24 recycle iterations: 44 s
WAIT 500.0 C, WABT 474.9 C, P_sep 12.0 bar, H2/HC 5.00
reactor dT (K): -62.6, -57.5, -42.7
reformate: 81.6 vol% (81.5 C5+), RON 97.9 (linear 90.1), MON 87.5, aromatics 69.7 vol%, benzene 3.08 vol%, RVP 2.3 psi
net H2 2.51 wt% of feed at 87.7 mol%; LPG 0.420 kg/s, fuel gas 0.044 kg/s
fired 53.87 MW, compressor 0.640 MW, coke 0.78 kg/h, cycle 272 d
closures: {'mass': '-6e-13', 'carbon': '-5e-13', 'hydrogen': '-7e-13', 'energy': '3e-13', 'sulfur': '1e-16', 'adiabatic_1': '-1e-16', 'adiabatic_2': '0e+00', 'adiabatic_3': '3e-15'}
The reformer’s fired duty is large for its size because the model has no feed-effluent exchanger: the charge heater does all the work of heating naphtha and recycle gas from the separator to the first reactor (see the reformer docs).
5. Distillate hydrotreater#
Kerosene, diesel and AGO are combined and hydrotreated together: two beds with a 15 %
quench, 340 C first-bed inlet, 60 bar, LHSV 1 h-1, 300 Nm3/m3, and the default (diesel)
kinetics. The crude-unit products (stripping water and light-end traces included) go in as
they are. The charge heater has no feed/effluent exchanger here (heater_inlet_T=None), so
its duty is the whole heat-up from 200 C.
The product is fractionated into jet and diesel at 240 C TBP by the same smooth split
(#328); the stripper’s wild naphtha goes to the gasoline pool on its blend grid
(product_stream("wild_naphtha", gases=False), #333), and the gas dissolved in it to fuel
gas.
dht_feed = {}
for name in ("kero", "diesel", "ago"):
for k, v in P[name].items():
if k.startswith("F_"):
dht_feed[k] = dht_feed.get(k, 0.0) + float(v)
dht_feed.update(T=C + 200.0, P=3e5)
dht_params = HydrotreaterParams(T_in=(C + 340.0,), P=60e5, lhsv=1.0, h2_oil=300.0)
t0 = time.time()
dht = Hydrotreater(char, dht_feed, dht_params)
dht_res = dht.solve(dht_feed)
tick(f"distillate hydrotreater converged: {bool(dht_res.converged)}", t0)
print(dht_res.table())
print(f"charge heater {float(dht_res.outputs['heater.fired_duty']) / 1e6:.1f} MW fired")
print("closures:", {k: f"{float(v):.0e}" for k, v in dht_res.balances.items()})
JET_CUT = C + 240.0
def jet_diesel(res):
return res.fractionate(cut_points=(JET_CUT,), products=("jet", "diesel"))
fd_dht = jet_diesel(dht_res)
print(fd_dht.table())
distillate hydrotreater converged: True: 132 s
WABT 350.6 C, bed dT 11.5, 8.0 K
product: S 8.9 wppm, N 97.3 wppm, SG 0.8215, aromatics 18.6 vol%, cetane index 51.3
yields (mass): product 96.71 %, wild naphtha 2.87 %, gas 0.742 %
H2: chemical 29.1 Nm3/m3 (173 scf/bbl), makeup 48.4 Nm3/m3, purity 96.2 %
recycle compressor 124.6 kW, purge 35.428 mol/s
charge heater 35.8 MW fired
closures: {'C': '1e-16', 'H': '3e-16', 'N': '6e-16', 'S': '5e-16', 'mass': '1e-16'}
jet: 21.3452 kg/s (46.26 wt%)
diesel: 24.7974 kg/s (53.74 wt%)
off_gas: 0.0000 kg/s (0.00 wt%)
mass closure 1.5e-16
6. Residue desulfurizer#
The atmospheric residue carries most of the crude’s sulfur (3 wt%), against a 0.5 wt% fuel
oil spec, and cutter blending cannot fix that: the lever rule would need about 85 % diesel
in the “fuel oil” (residue.cutter_fraction_for_sulfur). The ResidueDesulfurizer (#331)
is the conversion unit: three beds at 150 bar, LHSV 0.25 h-1, HDS by sulfur class with a
refractory residue class, HDM, CCR reduction and some conversion of the 538 C+ cuts, with a
once-through treat gas (90 % H2, 1000 Nm3/m3) and an ideal product split at 350 C into a
distillate and the desulfurized residue. Its constants are illustrative, set to land in
the ARDS ranges usually quoted (unverified).
t0 = time.time()
rds = ResidueDesulfurizer(char, P["residue"])
rds_res = rds.solve(P["residue"])
tick(f"residue desulfurizer finished: {bool(rds_res.converged)}", t0)
print(rds_res.table())
print("closures:", {k: f"{float(v):.0e}" for k, v in rds_res.balances.items()})
print(f"residue mass the RDS layout does not carry: {100 * rds.dropped_mass_fraction:.2g} %")
residue desulfurizer finished: True: 19 s
WABT 391.1 C, total bed dT 64.2 K
feed: S 3.06 wt%, Ni+V 78 wppm, CCR 9.7 wt%
desulfurized residue: S 0.332 wt%, Ni+V 18.2 wppm, CCR 5.42 wt%, SG 0.9294, yield 85.4 %
HDS 90.0 %, HDM 79.8 %, CCR reduction 51.4 %, conversion 11.8 %
distillate 12.6 % (S 1747 wppm), H2 chemical 109 Nm3/m3 (1.04 wt%)
closures: {'C': '2e-16', 'H': '2e-16', 'N': '1e-16', 'NiV': '1e-16', 'S': '6e-16', 'mass': '0e+00'}
residue mass the RDS layout does not carry: 0 %
The residue goes from 3.06 wt% S to 0.332 wt% in the desulfurized residue (90.0 % HDS), with 80 % of its Ni+V and half of its CCR removed. The 12.6 % that converts to distillate keeps 1750 wppm S: its fragments inherit their parent’s sulfur, so a refinery would send it to a distillate hydrotreater, and here it goes to the fuel oil, where it still fits. The price is 109 Nm3/m3 of chemical hydrogen (1.04 wt% of the feed).
7. The hydrogen header#
The reformer’s net gas, at 12 bar and the purity it reports, goes to a 10 bar header
(HydrogenNetwork, #329). The hydrotreaters draw their makeup from it, and so does the
residue desulfurizer, which is by far the largest consumer; a hydrogen plant (SMR + PSA,
99.9 mol %, an illustrative swing source with no capacity limit here) fills the deficit.
Every consumer on a header gets its gas at the header’s purity, and the hydrotreaters’
makeup demand depends on that purity (a leaner makeup carries methane that the purge has to
take out, and hydrogen with it), so close_hydrotreater_loop re-solves each hydrotreater
with the header’s composition until the purity stops moving (each re-solve of a compiled
unit takes a second or two).
The residue desulfurizer is on the header as a consumer whose demand is its chemical hydrogen consumption. Its model has a once-through treat gas of fixed composition, so the purity it receives does not feed back into it, and the hydrogen a real RDS loop loses with its purge and in solution is not counted: its demand is a lower bound.
The reformer was solved on the naphtha hydrotreater’s product at the default makeup (97 % H2); the loop changes that makeup, so the treated naphtha is fractionated again and the reformer re-solved from its converged recycle, and the header is balanced once more on the new net gas. The last line says how far that outer pass moved the header purity.
MOL_PER_NM3 = 101325.0 / (8.314462618 * 273.15)
RDS_P = 150e5
header = h2.Header("main", P=10e5, swing=(h2.H2Plant(capacity=float("inf")),))
net0 = h2.HydrogenNetwork(
producers=[h2.Producer.from_reformer(ref)],
consumers=[h2.Consumer.from_hydrotreater("nht", nht_res, nht_params, min_purity=0.85),
h2.Consumer.from_hydrotreater("dht", dht_res, dht_params, min_purity=0.85),
h2.Consumer("rds", rds_res.outputs["h2.chemical"], P=RDS_P)],
headers=[header])
t0 = time.time()
loop = h2.close_hydrotreater_loop(net0, {"nht": (nht, nht_feed, nht_params), "dht": (dht, dht_feed, dht_params)},
response_step=0.01)
tick(f"hydrogen loop closed: {loop.converged} in {loop.passes} passes", t0)
S_before = {"nht": float(nht_res.outputs["product.S_wppm"]), "dht": float(dht_res.outputs["product.S_wppm"])}
nht_res, dht_res = loop.units["nht"], loop.units["dht"]
nht_params, dht_params = loop.params["nht"], loop.params["dht"]
for k, r in (("nht", nht_res), ("dht", dht_res)):
print(f" {k}: product S {S_before[k]:.3g} -> {float(r.outputs['product.S_wppm']):.3g} wppm at the header's makeup")
t0 = time.time()
fr = naphtha_split(nht_res)
ref_feed = NaphthaFeed.from_hydrotreater(fr, "heavy_naphtha")
ref = reformer.solve(ref_feed, tear_initial=ref.tear)
fd_dht = jet_diesel(dht_res)
net = loop.network.replace(producers=[h2.Producer.from_reformer(ref)])
H = net.solve()
tick(f"reformer re-solved: {bool(ref.converged)}", t0)
print(H.table())
print(f"header purity moved by {abs(float(H.outputs['main.purity']) - float(loop.header.outputs['main.purity'])):.1e} "
"(mole fraction) when the reformer saw the loop's naphtha")
print("feasible:", H.feasible)
o = H.outputs
rows = [("reformer net gas", float(o["reformer.h2"]), f"{100 * float(o['reformer.purity']):.1f}", ""),
("hydrogen plant (swing)", float(o["h2_plant.h2"]), "99.9", "fills the deficit")]
for c, name in (("nht", "naphtha HDT"), ("dht", "distillate HDT"), ("rds", "residue desulfurizer")):
note = ""
if c in loop.units:
u = loop.units[c].outputs
note = (f"chemical {float(u['h2.chemical']):.1f}, dissolved {float(u['h2.dissolved']):.1f}, "
f"purge {float(u['h2.purge']):.1f}")
else:
note = "chemical consumption only (lower bound)"
rows.append((f"{name} makeup", -float(o[f"{c}.makeup_h2"]), f"{100 * float(o[f'{c}.purity']):.1f}", note))
rows.append(("purge to fuel gas", -float(o["main.h2_surplus"]), "", f"{float(o['main.h2_surplus']) / MOL_PER_NM3 * 3.6:.2f} kNm3/h"))
table(rows, ["hydrogen", "mol/s H2", "purity (mol %)", "where it goes (mol/s)"], floatfmt=".1f")
print("header closures:", {k: f"{float(v):.0e}" for k, v in H.balances.items()})
hydrogen loop closed: True in 2 passes: 14 s
nht: product S 0.000203 -> 0.000199 wppm at the header's makeup
dht: product S 8.86 -> 8.62 wppm at the header's makeup
reformer re-solved: True: 5 s
header main: purity 95.24 mol%, H2 surplus 0.00 mol/s (0.00 kNm3/h) to fuel
nht: makeup 16.58 mol/s, H2 15.79 mol/s at 95.24 mol%, margin +10.24 mol%
dht: makeup 131.90 mol/s, H2 125.62 mol/s at 95.24 mol%, margin +10.24 mol%
rds: makeup 418.48 mol/s, H2 398.56 mol/s at 95.24 mol%
fuel gas 0.00 mol/s (0.000 kg/s, H2 0.00 mol/s); export 0.00 mol/s
closures: total 0e+00, hydrogen 0e+00, mass 0e+00, makeup_h2 0e+00
header purity moved by 1.7e-07 (mole fraction) when the reformer saw the loop's naphtha
feasible: {'main.balanced': True, 'main.pressure': True, 'nht.purity': True, 'dht.purity': True}
| hydrogen | mol/s H2 | purity (mol %) | where it goes (mol/s) | |—————————–|————|——————|——————————————-| | reformer net gas | 189.2 | 87.7 | | | hydrogen plant (swing) | 350.8 | 99.9 | fills the deficit | | naphtha HDT makeup | -15.8 | 95.2 | chemical 2.4, dissolved 6.7, purge 6.7 | | distillate HDT makeup | -125.6 | 95.2 | chemical 75.5, dissolved 16.1, purge 34.0 | | residue desulfurizer makeup | -398.6 | 95.2 | chemical consumption only (lower bound) | | purge to fuel gas | -0.0 | | 0.00 kNm3/h |
header closures: {'total': '0e+00', 'hydrogen': '0e+00', 'mass': '0e+00', 'makeup_h2': '0e+00'}
The residue desulfurizer takes 399 of the 540 mol/s of hydrogen the refinery consumes. The reformer’s net gas (189 mol/s at 87.7 %) covers about a third; the hydrogen plant supplies the other 350 mol/s (28 kNm3/h). Mixing the two puts the header at 95.2 mol %, inside both hydrotreaters’ 85 % spec. Re-solved at that purity, the hydrotreaters change little (the distillate unit’s product goes from 8.9 to 8.6 wppm S). Most of a hydrotreater’s makeup is not chemistry: the purge and the hydrogen dissolved in the separator liquid are more than four fifths of the naphtha unit’s and 40 % of the distillate unit’s.
8. Product pools#
Each pool uses BlendPool’s default specs (illustrative: US regular summer gasoline,
Jet A, ULSD S15, and the very-low-sulfur fuel oil of residue.VLSFO_SPECS). The recipe is
“everything the refinery makes of that product goes into its pool”, so the margins say
whether this refinery’s own streams make spec without purchased blendstocks.
What is estimated rather than given (#330, difflow_refinery.properties): the straight-run
octane of the light naphtha and of the distillate hydrotreater’s wild naphtha (pure-compound
octanes of model compounds blended by Ethyl RT-70 from the P/N/A composition), the jet’s
flash, freeze and smoke points, the diesel’s flash point, and the fuel oil’s viscosity. All
but the freeze point’s melting-point data are marked unverified there; a measured value
would be passed as an override and win. The reformate carries its own RON/MON (the
reformer’s) and its sulfur.
One pool takes one mode, so the gasoline pool is in property mode (the reformate is a property-mode component) and the straight-run components’ composition-derived and estimated properties are copied into property-mode components.
GASOLINE_KEYS = ("SG", "RVP_psi", "S_ppm", "RON", "MON", "aromatics_vol", "naphthenes_vol",
"paraffins_vol", "olefins_vol")
def property_mode(c, keys=GASOLINE_KEYS):
return BlendComponent.from_properties(c.name, **{k: c.properties[k] for k in keys})
reformate = ref.blend_component("reformate")
V_ref = std_liquid_volume(ref.flows("reformate"))
lsr = BlendComponent.from_stream("light naphtha", fr.products["light_naphtha"], fr.char)
dht_char = dht_res.product_char
dht_wild = dht_res.product_stream("wild_naphtha", gases=False)
wn = BlendComponent.from_stream("DHT wild naphtha", dht_wild, dht_char)
gasoline = BlendPool("gasoline")([property_mode(reformate), property_mode(lsr), property_mode(wn)],
[V_ref, lsr.available_volume, wn.available_volume], basis="volume_flow")
table([(c.name, f"{float(c.properties['RON']):.1f}", f"{float(c.properties['MON']):.1f}",
f"{float(c.properties['S_ppm']):.2g}", src)
for c, src in ((reformate, "reformer"), (lsr, "estimated"), (wn, "estimated"))],
["gasoline component", "RON", "MON", "S (ppm)", "octane"])
jet_c = BlendComponent.from_stream("jet", fd_dht.products["jet"], dht_char)
jet = BlendPool("jet")([jet_c], [1.0], basis="split")
diesel_c = BlendComponent.from_stream("diesel", fd_dht.products["diesel"], dht_char)
ulsd = BlendPool("ulsd")([diesel_c], [1.0], basis="split")
fuel_oil = fuel_oil_blend([rds_res.blend_component("residue"), rds_res.blend_component("distillate")],
[rds_res.volume("residue"), rds_res.volume("distillate")])
pools = (("gasoline", gasoline), ("jet", jet), ("ULSD", ulsd), ("fuel oil", fuel_oil))
for name, r in pools:
print(f"\n{name}: {float(r.volume) * BBL_D:,.0f} bbl/d")
table([(s, f"{float(r.properties[s.split(' ')[0]]):.4g}",
f"{float(m):+.4g}", "on spec" if float(m) >= 0 else "OFF SPEC") for s, m in r.margins.items()],
["spec", "value", "margin", ""])
| gasoline component | RON | MON | S (ppm) | octane | |———————-|——-|——-|———–|———–| | reformate | 97.9 | 87.5 | 1.2e-05 | reformer | | light naphtha | 65.7 | 64.3 | 9.4e-05 | estimated | | DHT wild naphtha | 28.8 | 28.6 | 4 | estimated |
gasoline: 16,429 bbl/d
| spec | value | margin | | |————–|———|———-|———-| | RON >= 91 | 83.66 | -7.341 | OFF SPEC | | MON >= 82 | 77.19 | -4.808 | OFF SPEC | | RVP_psi <= 9 | 6.076 | 2.924 | on spec | | S_ppm <= 10 | 0.2497 | 9.75 | on spec |
jet: 14,568 bbl/d
| spec | value | margin | | |—————–|———|———-|———| | S_ppm <= 3000 | 3.937 | 2996 | on spec | | freeze_C <= -40 | -64.79 | 24.79 | on spec | | smoke_mm >= 18 | 23.95 | 5.954 | on spec | | flash_C >= 38 | 51.15 | 13.15 | on spec |
ULSD: 15,983 bbl/d
| spec | value | margin | | |——————–|———|———-|———| | S_ppm <= 15 | 12.65 | 2.355 | on spec | | cetane_index >= 40 | 57.14 | 17.14 | on spec | | flash_C >= 52 | 105 | 52.99 | on spec | | T90_d86_C <= 338 | 331.2 | 6.846 | on spec |
fuel oil: 44,766 bbl/d
| spec | value | margin | | |———————-|———|————|———| | S_ppm <= 5000 | 3120 | 1880 | on spec | | viscosity_cSt <= 380 | 68.73 | 311.3 | on spec | | SG <= 0.991 | 0.92 | 0.07096 | on spec | | CCR_wt <= 18 | 4.81 | 13.19 | on spec |
What the pools say about this refinery:
Gasoline is short of octane: RON 83.7 against 91, MON 77.2 against 82. The reformate is at RON 97.9, but the straight-run components pull the pool down, and their octanes are now estimated rather than assumed: the light naphtha at RON 65.7 / MON 64.3 (the first version of this notebook assumed 70/68), and the distillate hydrotreater’s wild naphtha at RON 28.8 (it assumed 60; a heavy, paraffinic straight-run naphtha is a poor gasoline component, and the estimate is marked unverified). A hydroskimming refinery fixes that with a C5/C6 isomerization unit on the light naphtha (
IsomerizationUnit,examples/39_refinery_isomerization.ipynb), by sending the wild naphtha to the naphtha hydrotreater and the reformer rather than to the pool, with a lower naphtha cut point (section 10 gives the derivative), a higher reformer severity, or bought octane. Sulfur is far inside the spec (0.25 ppm, almost all of it the wild naphtha’s 4.0 ppm).Jet and ULSD make spec. The jet’s flash (51 C), freeze (-65 C) and smoke point (24 mm) are estimates now, not assumptions restated, so the margins say something; but the estimates are unverified correlations, and a -65 C freeze point is what the n-paraffin solubility model gives for this cut, not a measured value. ULSD meets 15 ppm with 2 ppm to spare, cetane index 57.
The fuel oil makes spec: 0.31 wt% sulfur against 0.5, 69 cSt at 50 C (estimated), SG 0.920, CCR 4.8 wt%. That is the residue desulfurizer’s doing, and it is the refinery’s largest hydrogen consumer by far (section 7) and its largest source of H2S (section 9).
9. The refinery in one table#
Product rates, the qualities the specs are written on, and the energy the units report. “Fuel gas” collects the gas plant’s deethanizer overhead, the reformer’s fuel gas, the hydrotreaters’ purge and stripper off-gas, the gas dissolved in their wild naphthas, the residue desulfurizer’s gas, and the hydrogen header’s purge (all as reported, untreated). “LPG” adds the reformer’s LPG to the gas plant’s.
Then two checks that the connections lose nothing: the refinery’s overall mass balance over every terminal stream of the process units (water and steam included; each unit closes its own to round-off, so this tests the connections), and where the crude’s sulfur ends up. The hydrogen header is balanced on its own (section 7): in the unit balance the hydrotreaters’ makeup gas and the residue desulfurizer’s treat gas are inputs, and the reformer’s net gas an output.
def flows_mass(res, unit, *names):
return sum(float(res.streams[n].mass(unit.layout)) for n in names)
def char_mass(stream, bchar):
return float(sum(stream.get(f"F_{n}", 0.0) * m for n, m in zip(bchar.names, np.asarray(bchar.mw)))) / 1000.0
ro = ref.outputs()
dissolved = float(dht_res.outputs["naphtha.dissolved_gas_rate"])
fuel_kg_s = (gp_mass(fgas) + float(ro["fuel_gas.kg_s"])
+ flows_mass(nht_res, nht, "purge", "off_gas") + flows_mass(dht_res, dht, "purge", "off_gas")
+ float(fr.rates["off_gas"]) + dissolved + float(rds_res.streams["gas"].mass(rds.layout))
+ float(o["fuel_gas.kg_s"]))
lpg_kg_s = gp_mass(lpg) + float(ro["lpg.kg_s"])
t_h = lambda kg_s: kg_s * 3.6
products = [
("fuel gas", "", t_h(fuel_kg_s), "includes the RDS gas (H2S, NH3, its once-through H2)"),
("LPG", "", t_h(lpg_kg_s), f"gas-plant LPG vapour pressure margin {float(q['margins']['vapor_pressure']) / 1e3:.0f} kPa"),
("gasoline", float(gasoline.volume) * BBL_D, t_h(float(gasoline.mass)),
f"RON {float(gasoline.properties['RON']):.1f}, MON {float(gasoline.properties['MON']):.1f}, "
f"RVP {float(gasoline.properties['RVP_psi']):.1f} psi, S {float(gasoline.properties['S_ppm']):.2g} ppm"),
("jet", float(jet.volume) * BBL_D, t_h(float(jet.mass)),
f"S {float(jet.properties['S_ppm']):.2g} ppm, flash {float(jet.properties['flash_C']):.0f} C, "
f"freeze {float(jet.properties['freeze_C']):.0f} C, smoke {float(jet.properties['smoke_mm']):.0f} mm"),
("ULSD", float(ulsd.volume) * BBL_D, t_h(float(ulsd.mass)),
f"S {float(ulsd.properties['S_ppm']):.2g} ppm, cetane index {float(ulsd.properties['cetane_index']):.1f}, "
f"flash {float(ulsd.properties['flash_C']):.0f} C, T90 {float(ulsd.properties['T90_d86_C']):.0f} C"),
("fuel oil", float(fuel_oil.volume) * BBL_D, t_h(float(fuel_oil.mass)),
f"S {float(fuel_oil.properties['S_ppm']) / 1e4:.2f} wt%, {float(fuel_oil.properties['viscosity_cSt']):.0f} cSt "
f"at 50 C, SG {float(fuel_oil.properties['SG']):.3f}, CCR {float(fuel_oil.properties['CCR_wt']):.1f} wt%"),
("H2 from the hydrogen plant", "", t_h(float(o["h2_plant.h2"]) * 2.01588e-3), "an input, not a product"),
]
table(products, ["product", "bbl/d", "t/h", "key qualities"], floatfmt=",.1f")
# Overall mass balance over the terminal streams of the units (water and steam included).
W = 18.01528e-3
crude_feed = cdu_unit.feed(BPD, T=C + 240.0, P=6e5)
crude_kg_s = sum(float(crude_feed[f"F_{n}"]) * mw for n, mw in zip(char.names, np.asarray(char.component_MW))) / 1000.0
cdu_steam = (cdu_params.bottom_steam + sum(sp.steam for sp in cdu_params.side_products)) * W
cdu_water = sum(float(P[n].get("F_water", 0.0)) for n in ("water", "offgas", "naphtha", "residue")) * W
fo_char = BlendCharacterization.from_characterization(char)
HDT_IN, HDT_OUT = ("makeup", "steam"), ("off_gas", "purge", "acid_gas", "hps_water", "sour_water")
mass_in = {"crude": crude_kg_s, "CDU stripping steam": cdu_steam, "assumed offgas H2S": H2S * 34.081e-3,
"HDT makeup gas and steam": flows_mass(nht_res, nht, *HDT_IN) + flows_mass(dht_res, dht, *HDT_IN),
"RDS treat gas": float(rds_res.streams["treat_gas"].mass(rds.layout))}
mass_out = {"gas plant fuel gas, acid gas, LPG": gp_mass(fgas) + gp_mass(acid) + gp_mass(lpg),
"light naphtha + NHT fractionator off-gas": float(fr.rates["light_naphtha"] + fr.rates["off_gas"]),
"reformate, LPG, fuel gas, net gas": sum(float(ro[k]) for k in ("reformate.kg_s", "lpg.kg_s",
"fuel_gas.kg_s", "net_gas.kg_s")),
"jet + diesel": float(fd_dht.rates["jet"] + fd_dht.rates["diesel"]),
"DHT wild naphtha and its dissolved gas": float(dht_res.stream_mass(dht_res.product_stream("wild_naphtha"))),
"HDT gases and waters": flows_mass(nht_res, nht, *HDT_OUT) + flows_mass(dht_res, dht, *HDT_OUT),
"RDS gas, distillate, residue": sum(float(rds_res.streams[n].mass(rds.layout))
for n in ("gas", "distillate", "residue")),
"CDU water (decanted and in products)": cdu_water}
m_in, m_out = sum(mass_in.values()), sum(mass_out.values())
table([(k, t_h(v)) for k, v in mass_in.items()] + [("**in**", t_h(m_in))]
+ [(k, t_h(v)) for k, v in mass_out.items()] + [("**out**", t_h(m_out))],
["refinery mass balance", "t/h"], floatfmt=",.3f")
# where the remainder comes from: the two places a stream changes component table
handover = {
"gas_plant_feed: folded cuts and light-end molar masses (reported mass_change)":
sum(float(v) for v in gp.mass_change.values()),
"gas plant -> NHT: the stabilised naphtha on the NHT layout's molar masses, its H2S trace dropped":
float(nht_res.streams["feed"].mass(nht.layout)) - gp_mass(stab),
"CDU kero/diesel/AGO -> DHT: on the DHT layout's molar masses":
float(dht_res.streams["feed"].mass(dht.layout))
- sum(dht_feed.get(f"F_{n}", 0.0) * mw for n, mw in zip(char.names, np.asarray(char.component_MW))) / 1000.0
- dht_feed.get("F_water", 0.0) * W,
"CDU residue -> RDS: components the RDS layout does not carry":
float(rds_res.streams["feed"].mass(rds.layout)) - char_mass(P["residue"], fo_char),
}
print(f"imbalance {(m_out - m_in) / m_in:+.1e} of the inputs ({t_h(m_out - m_in):+.4f} t/h), of which:")
for k, v in handover.items():
print(f" {k}: {t_h(v):+.4f} t/h")
print(f" unexplained: {t_h(m_out - m_in - sum(handover.values())):+.1e} t/h")
| product | bbl/d | t/h | key qualities | |—————————-|———-|——-|——————————————————| | fuel gas | | 52.7 | includes the RDS gas (H2S, NH3, its once-through H2) | | LPG | | 7.5 | gas-plant LPG vapour pressure margin 812 kPa | | gasoline | 16,429.0 | 82.9 | RON 83.7, MON 77.2, RVP 6.1 psi, S 0.25 ppm | | jet | 14,568.3 | 76.8 | S 3.9 ppm, flash 51 C, freeze -65 C, smoke 24 mm | | ULSD | 15,982.8 | 89.3 | S 13 ppm, cetane index 57.1, flash 105 C, T90 331 C | | fuel oil | 44,766.0 | 272.6 | S 0.31 wt%, 69 cSt at 50 C, SG 0.920, CCR 4.8 wt% | | H2 from the hydrogen plant | | 2.5 | an input, not a product |
| refinery mass balance | t/h | |——————————————|———| | crude | 540.687 | | CDU stripping steam | 16.214 | | assumed offgas H2S | 0.006 | | HDT makeup gas and steam | 4.538 | | RDS treat gas | 44.974 | | in | 606.419 | | gas plant fuel gas, acid gas, LPG | 6.029 | | light naphtha + NHT fractionator off-gas | 29.822 | | reformate, LPG, fuel gas, net gas | 54.645 | | jet + diesel | 166.100 | | DHT wild naphtha and its dissolved gas | 5.242 | | HDT gases and waters | 5.316 | | RDS gas, distillate, residue | 323.051 | | CDU water (decanted and in products) | 16.214 | | out | 606.418 |
imbalance -1.8e-06 of the inputs (-0.0011 t/h), of which:
gas_plant_feed: folded cuts and light-end molar masses (reported mass_change): -0.0010 t/h
gas plant -> NHT: the stabilised naphtha on the NHT layout's molar masses, its H2S trace dropped: -0.0001 t/h
CDU kero/diesel/AGO -> DHT: on the DHT layout's molar masses: -0.0000 t/h
CDU residue -> RDS: components the RDS layout does not carry: -0.0000 t/h
unexplained: -1.9e-07 t/h
# Where the crude's sulfur goes (the assumed offgas H2S is left out: it is not the assay's).
crude_S = crude_kg_s * 0.018
def hdt_S(r):
o_ = r.outputs
feed = float(o_["feed.S_wppm"] * o_["feed.rate"]) * 1e-6
liquid = float(o_["product.S_wppm"] * o_["product.rate"] + o_["naphtha.S_wppm"] * o_["naphtha.rate"]) * 1e-6
return feed, liquid
nS_feed, nS_liq = hdt_S(nht_res)
dS_feed, dS_liq = hdt_S(dht_res)
fo_S = float(fuel_oil.properties["S_ppm"]) * 1e-6 * float(fuel_oil.mass)
rds_h2s_S = float(rds_res.outputs["h2s.make"]) * MW_S
rS = ref.sulfur()
s_rows = [("fuel oil (desulfurized residue + RDS distillate)", fo_S),
("H2S from the residue desulfurizer", rds_h2s_S),
("naphtha HDT liquids: light naphtha (gasoline)",
float(lsr.properties["S_ppm"]) * 1e-6 * float(fr.rates["light_naphtha"])),
("naphtha HDT liquids: reformate", float(rS["S.reformate"])),
("naphtha HDT liquids: H2S in the reformer's net gas and fuel gas", float(rS["S.net_gas"] + rS["S.fuel_gas"])),
("distillate HDT liquids (jet, ULSD, wild naphtha)", dS_liq),
("H2S from both hydrotreaters (to sulfur recovery)", nS_feed - nS_liq + dS_feed - dS_liq),
("on the cuts in the gas plant's fuel gas and LPG (not tracked)",
crude_S - (fo_S + rds_h2s_S) - nS_feed - dS_feed)]
table([(k, 3.6 * v, 100 * v / crude_S) for k, v in s_rows],
[f"crude sulfur, {3.6 * crude_S:.2f} t/h", "t/h", "% of crude S"], floatfmt=".4g")
print(f"naphtha HDT liquid sulfur against what the split and the reformer account for: "
f"{3.6 * (nS_liq - s_rows[2][1] - s_rows[3][1] - s_rows[4][1]):+.1e} t/h")
energy = [
("CDU furnace", float(cdu.column.furnace_fired_duty) / 1e6, "fired"),
("naphtha HDT charge heater", float(nht_res.outputs["heater.fired_duty"]) / 1e6, "fired (after a feed/effluent exchanger)"),
("distillate HDT charge heater", float(dht_res.outputs["heater.fired_duty"]) / 1e6, "fired (no exchanger modelled)"),
("reformer heaters", float(ro["heaters.fired_MW"]), "fired"),
("gas-plant reboilers", float(od["reboiler.duty"] + ob["reboiler.duty"]) / 1e6, "absorbed (not fired)"),
("compressors (gas plant, reformer, 2 x HDT)",
(float(ci["power"]) + float(nht_res.outputs["compressor.power"]) + float(dht_res.outputs["compressor.power"])) / 1e6
+ float(ro["compressor.power_MW"]), "shaft"),
]
table(energy, ["energy", "MW", ""], floatfmt=".1f")
print("not in the table: the residue desulfurizer's heater (its bed-1 inlet is a spec), makeup "
"compressors, the hydrogen plant")
| crude sulfur, 9.73 t/h | t/h | % of crude S | |—————————————————————–|———–|—————-| | fuel oil (desulfurized residue + RDS distillate) | 0.8505 | 8.739 | | H2S from the residue desulfurizer | 7.656 | 78.67 | | naphtha HDT liquids: light naphtha (gasoline) | 2.802e-09 | 2.879e-08 | | naphtha HDT liquids: reformate | 5.871e-10 | 6.032e-09 | | naphtha HDT liquids: H2S in the reformer’s net gas and fuel gas | 1.115e-08 | 1.146e-07 | | distillate HDT liquids (jet, ULSD, wild naphtha) | 0.001452 | 0.01492 | | H2S from both hydrotreaters (to sulfur recovery) | 1.224 | 12.58 | | on the cuts in the gas plant’s fuel gas and LPG (not tracked) | 1.385e-05 | 0.0001423 |
naphtha HDT liquid sulfur against what the split and the reformer account for: +3.0e-24 t/h
| energy | MW | | |——————————————–|——|—————————————–| | CDU furnace | 50.3 | fired | | naphtha HDT charge heater | 5.3 | fired (after a feed/effluent exchanger) | | distillate HDT charge heater | 35.9 | fired (no exchanger modelled) | | reformer heaters | 53.9 | fired | | gas-plant reboilers | 20.1 | absorbed (not fired) | | compressors (gas plant, reformer, 2 x HDT) | 0.8 | shaft |
not in the table: the residue desulfurizer's heater (its bed-1 inlet is a spec), makeup compressors, the hydrogen plant
The overall balance closes to 1.8e-6 of the inputs. All of that except 3e-10 is the change
gas_plant_feed reports when it folds the heavy cut traces into the heaviest kept cut (the
moles are kept, at a lighter cut’s molar mass). The first version of this notebook closed to
1.7e-4, through its bridges. The sulfur table closes by construction (its last line is the
remainder, sulfur on the cuts that leave in the gas plant’s fuel gas and LPG). The naphtha
hydrotreater’s liquid sulfur is accounted for by the split and the reformer to round-off.
Nine-tenths of the crude’s sulfur leaves as H2S, almost all of it from the residue
desulfurizer, for a sulfur plant this refinery does not model.
10. A gradient across units#
difflow_refinery.plant (#334) composes library units into one differentiable function.
The units do not all differentiate the same way:
from difflow_refinery.plant import ad_mode_table
display(Markdown(ad_mode_table()))
| unit | object | modes | why | other mode | |—|—|—|—|—| | crude unit | unit.CrudeUnit / column.CrudeColumn | fwd + rev | EO MESH Newton in lax.while_loop on stop-gradient inputs, then one Newton step with the converged Jacobian (implicit-function derivative) | – | | vacuum unit | vacuum.VacuumColumn / vacuum.StageColumn | fwd + rev | stage-network Newton, implicit step reusing the converged Jacobian | – | | gas plant columns | gasplant.GasPlantColumn | fwd + rev | the vacuum StageColumn machinery on CubicThermo: three Newton passes, then the implicit step | – | | gas compressor, amine treater | gasplant.GasCompressor, gasplant.AmineTreater | fwd + rev | closed-form stages and removal fractions (pure JAX) | – | | hydrotreater | hydrotreating.Hydrotreater | rev | beds integrated by diffrax with RecursiveCheckpointAdjoint (a custom_vjp, reverse only); the recycle tear’s Newton takes its Jacobian from a forward-adjoint copy and ends with one implicit step on the reverse-mode residual | HydrotreaterParams(reactor=ReactorOptions(adjoint=’forward’)) integrates the beds with diffrax.ForwardMode: then forward mode only (jax.jacfwd), same values | | residue desulfurizer | residue.ResidueDesulfurizer | rev | trickle beds with RecursiveCheckpointAdjoint (reverse only), quench mixing by implicit Newton | RDSParams(reactor=ReactorOptions(adjoint=’forward’)) should give forward mode only, by the hydrotreater’s construction (not tested) | | hydrocracker | hydrocracking.Hydrocracker | rev | checkpointed bed adjoints, and the UCO recycle’s Anderson fixed point with a GMRES adjoint (a custom_vjp); forward mode is not available | – | | catalytic reformer | reforming.CatalyticReformer | fwd | beds integrated with diffrax.ForwardMode inside the Flowsheet’s traced recycle path (optimistix fixed point, implicit differentiation, which needs JVPs of the loop) | none for the whole unit: CatalyticReformer(adjoint=’reverse’) makes a stand-alone reactor reverse-capable, but the recycle’s implicit derivative still needs JVPs | | FCC | fcc.FCCUnit | fwd + rev | riser by constant-step Tsit5 with diffrax.DirectAdjoint (both modes), riser-regenerator heat balance by optimistix Newton with implicit adjoint | – | | isomerization | isomerization.IsomerizationUnit | fwd | with a DIH the recycle is converged by the Flowsheet’s Python Anderson loop and differentiated at its solution by a jax.custom_jvp | – | | alkylation | alkylation.AlkylationUnit | fwd | DIB-overhead tear: the Flowsheet’s traced optimistix fixed point; the planning block differentiates it with jax.jacfwd | – | | preheat train | preheat.PreheatedCrudeUnit | fwd + rev | outer Newton in lax.while_loop with an implicit step (preheat._newton.implicit_step) | – | | blend pool | blending.BlendPool / BlendComponent.from_stream | fwd + rev | closed-form blending rules and property estimates (pure JAX) | – | | adapters | gasplant.gas_plant_feed, hydrotreating.HydrotreaterResult.fractionate, reforming.NaphthaFeed.from_hydrotreater, residue.fuel_oil_blend | fwd + rev | sums, sigmoid splits and lumping (pure JAX); gas_plant_feed’s cut SELECTION is static (pass cuts= under a transform) | – | | hydrogen network | hydrogen.HydrogenNetwork | fwd + rev | header balances, swing clips and the PSA purity target in pure JAX (the clips are kinks); close_hydrotreater_loop is a concrete Python iteration that returns a linear purity response | – |
The chain differentiated here runs from the naphtha hydrotreater to the gasoline pool, three stages:
nht: the hydrotreater at bed inlet temperature
T(with the header’s makeup), the naphtha fractionator at cut pointT_cut,NaphthaFeed.from_hydrotreaterand the light naphtha’s estimated properties;reformer: the reformer on that feed (warm-started from its converged recycle);
pool: the gasoline pool, with the distillate hydrotreater’s wild naphtha held fixed.
The reformer is forward mode only. The naphtha hydrotreater’s default bed adjoint is
reverse only, and a chain of the two could not be traced end to end in either mode: then
Chain.jacobian falls back to the chain rule by unit Jacobians (method="chain"). This
notebook’s hydrotreater was built with ReactorOptions(adjoint="forward"), so every stage is
forward-capable and the whole chain is one jax.jacfwd: two inputs, two tangents. The two unit stages are
jit=True: the reformer’s recycle is a Python loop that, under a transform, would otherwise be traced
and compiled again on every call.
nht_p = nht_params
wild_fixed = property_mode(wn)
V_wild = wn.available_volume
KEEP = GASOLINE_KEYS
tear = ref.tear
def nht_stage(x):
r = nht.solve(nht_feed, params=dataclasses.replace(nht_p, T_in=(x[0],)), warn=False)
f = naphtha_split(r, x[1])
c = BlendComponent.from_stream("light naphtha", f.products["light_naphtha"], f.char)
return {"feed": NaphthaFeed.from_hydrotreater(f, "heavy_naphtha"),
"light": {k: c.properties[k] for k in KEEP}, "V_light": c.available_volume}
def reformer_stage(d):
r = reformer.solve(d["feed"], tear_initial=tear, tol=1e-11, max_iter=400, on_nonconvergence="ignore")
c = r.blend_component("reformate")
return {"reformate": {k: c.properties[k] for k in KEEP}, "V_ref": std_liquid_volume(r.flows("reformate")),
"light": d["light"], "V_light": d["V_light"], "h2": r.outputs()["h2.net_mol_s"]}
def pool_stage(d):
comps_ = [BlendComponent.from_properties("reformate", **d["reformate"]),
BlendComponent.from_properties("light naphtha", **d["light"]), wild_fixed]
gpool = BlendPool("gasoline")(comps_, [d["V_ref"], d["V_light"], V_wild], basis="volume_flow")
return jnp.stack([gpool.properties["RON"], gpool.properties["S_ppm"], gpool.volume * BBL_D,
d["reformate"]["RON"], d["h2"]])
OUT_NAMES = ["gasoline RON", "gasoline S (ppm)", "gasoline (bbl/d)", "reformate RON", "reformer net H2 (mol/s)"]
chain = Chain(Stage("nht", nht_stage, modes="fwd", jit=True),
Stage("reformer", reformer_stage, modes="fwd", jit=True),
Stage("pool", pool_stage))
x0 = jnp.asarray([float(nht_p.T_in[0]), NAPHTHA_CUT])
t0 = time.time()
# central differences of the same three functions, unjitted: four solves of the compiled units
fd = central_difference(lambda x: pool_stage(reformer_stage(nht_stage(x))), x0, [0.25, 0.25])
tick("central differences (4 chain evaluations)", t0)
# Every unit above is compiled and held in memory, and the derivative compiles the hydrotreater
# and the reformer again (their JVPs, each stage as one jitted function). Dropping the compiled
# solves first keeps the peak inside this 4-core, 15 GB machine.
jax.clear_caches()
t0 = time.time()
J = chain.jacobian(x0)
tick(f"jacobian by {J.method} (trace, compile and run)", t0)
t0 = time.time()
J2 = chain.jacobian(x0)
tick("the same jacobian again (compiled)", t0)
table([(n, float(J.jacobian[i, 0]), fd[i, 0], float(J.jacobian[i, 1]), fd[i, 1])
for i, n in enumerate(OUT_NAMES)],
["output", "d/dT_NHT (AD)", "FD, h = 0.25 K", "d/dT_cut (AD)", "FD, h = 0.25 K"], floatfmt=".5g")
central differences (4 chain evaluations): 38 s
jacobian by fwd (trace, compile and run): 435 s
the same jacobian again (compiled): 29 s
| output | d/dT_NHT (AD) | FD, h = 0.25 K | d/dT_cut (AD) | FD, h = 0.25 K | |————————-|—————–|——————|—————–|——————| | gasoline RON | -0.0011095 | -0.0011096 | -0.27534 | -0.27535 | | gasoline S (ppm) | -2.9661e-05 | -2.9798e-05 | -0.00024999 | -0.00025 | | gasoline (bbl/d) | 0.0038774 | 0.0038707 | 22.875 | 22.876 | | reformate RON | -0.0016442 | -0.0016444 | 0.0069574 | 0.0069629 | | reformer net H2 (mol/s) | 0.0059628 | 0.0059633 | -1.6928 | -1.6928 |
The derivatives agree with the central differences to 0.1 % or better. The exception is the gasoline sulfur with respect to the hydrotreater temperature, at 0.5 %: the sulfur falls exponentially with temperature, and that is the truncation error of a 0.25 K difference.
What they say:
The hydrotreater’s bed temperature barely matters downstream of it at these conditions: the sulfur is already gone, and a hotter bed saturates a little more aromatic, which costs the reformate 0.0016 RON per K.
The naphtha cut point is the lever on the octane deficit. One kelvin more cut point sends more light naphtha straight to the pool: gasoline RON falls by 0.275, the pool gains 23 bbl/d (the light naphtha does not lose the reformer’s 18 % volume), and the reformer makes 1.7 mol/s less hydrogen. Read linearly, a 10 K lower cut point would gain about 2.8 RON for 230 bbl/d of gasoline and 17 mol/s of hydrogen. It would not close the 7-point deficit, which is the isomerization unit’s job.
The cost: 435 s for the first Jacobian (trace, compile and run: mostly the reformer’s JVP), 29 s
for a repeat, against 38 s for the four central-difference evaluations of units that were already
compiled. For one gradient, finite differences are cheaper on this machine. What the AD Jacobian
buys is exactness: no step to choose, as the sulfur row shows. A planner that relinearizes a
compiled chain at every iteration also pays the compile only once, and then 29 s per Jacobian.
The section’s earlier cell had to free the compiled solves (jax.clear_caches()) to fit the
derivative into this machine’s memory: see “Chaining units” in the refinery docs for the
measurements, and for the mixed-mode pattern (method="chain") when a unit is reverse-only.
print(f"whole notebook: {(time.time() - T_START) / 60:.1f} min")
whole notebook: 16.5 min
11. Gaps this refinery ran into, and what closed them#
The first version of this notebook had to write a bridge, marked STAND-IN, wherever the API had no connection, and it listed nine gaps. Each became an issue (#326 to #334), and this version uses what they added:
# |
Gap in the first version |
Issue |
Now |
|---|---|---|---|
1 |
No adapter from crude-unit products to a |
#326 |
|
2 |
Hydrotreater product grid not readable by a splitter or the reformer; |
#327 |
|
3 |
No jet/diesel split after a distillate hydrotreater |
#328 |
|
4 |
No hydrogen network |
#329 |
|
5 |
Straight-run octane, flash, freeze and smoke points, viscosity assumed; no sulfur through the reformer |
#330 |
|
6 |
No residue upgrading: fuel oil at 3 wt% S |
#331 |
|
7 |
No hydrotreater charge heater; negative |
#332 |
charge heater duty in the energy table; |
8 |
A NaN instead of |
#332, #333 |
a failed flash reports itself; dissolved gases carried; the balance now closes to 1.8e-6, all but 3e-10 of it the reported fold of |
9 |
No plant-level way to differentiate across units |
#334 |
|
What is still assumed or simplified, and why:
Illustrative constants: every kinetic constant of the hydrotreaters, the reformer and the residue desulfurizer, the reformer’s sulfur conversion and iso/normal splits, and the fractionators’ width. The trends and the gradients are exact for the model; the levels are not predictions for any catalyst or column.
The offgas H2S (above).
The residue desulfurizer has a once-through treat gas and no charge heater, and its hydrogen demand on the header is its chemical consumption only (no purge or solution losses): a lower bound on the largest consumer.
The hydrogen plant is an unlimited swing source; makeup compressors are not modelled, nor is a header purge (the swing makes up the balance exactly).
The gasoline pool is off spec. That is a finding about this refinery, not a gap: it needs an isomerization unit, which the plugin has (
examples/39), and which this notebook leaves out to keep its run time.