Why Determyn

Swipe an equation sideways to see all of it.

One vessel

The engine’s reference plant is a synthetic gas train. Its inlet separator, V-101, is a horizontal drum 2 m across and 6 m long. PIC-101 holds it at 28 bar absolute on the gas outlet valve PV-101, and LIC-101 holds the level at 0.8 m. The fluid has twenty components, from nitrogen to m-xylene.

The engine models V-101 as a liquid zone and a vapor zone sharing the volume and the pressure. The state is the inventory: each zone’s moles of each component and its internal energy, 42 numbers. Away from the minimum holdups they obey these balances.

\frac{d n_{i}^{L}}{d t}

= \sum_{k} F_{k}^{L} x_{i, k} - N_{i}, i = 1, \dots, 20

\frac{d n_{i}^{V}}{d t}

= \sum_{k} F_{k}^{V} y_{i, k} + N_{i}

\frac{d U^{L}}{d t}

= \sum_{k} F_{k}^{L} h_{k}^{L} - N h^{V} + Q_{int} - W

\frac{d U^{V}}{d t}

= \sum_{k} F_{k}^{V} h_{k}^{V} + N h^{V} - Q_{int} + W, N = \sum_{i} N_{i}

nᵢᴸ, nᵢⱽ: moles of component i in the liquid and vapor zones, mol; differential states
Uᴸ, Uⱽ: internal energy of each zone, J; differential states
Fₖᴸ, Fₖⱽ: molar flows of the liquid and vapor parts of the stream at nozzle k, mol/s, positive into the drum and free to change sign
x_i,k, y_i,k: mole fractions of component i in those parts; an outflow carries its zone’s composition
hₖᴸ, hₖⱽ: molar enthalpies of those parts, J/mol; each pure component as an ideal gas at 298.15 K has h = 0
Nᵢ, N: interphase transfer of component i and in total, mol/s, positive for evaporation
hⱽ: molar enthalpy of the vapor zone, J/mol, carried by the transferred moles
Q_int: heat from the vapor zone to the liquid zone, W
W: work the liquid zone does on the vapor zone as its volume changes, W

The transfer, heat and work terms leave one zone and enter the other, so mass and energy close exactly. V-101 has no heater, its heat loss is not modeled, and there is no reaction term yet.

Temperature and pressure follow

Temperature and pressure are not states. The inventory and the vessel fix them: each zone’s energy must equal its moles times its molar internal energy, the zone volumes must fill the drum, and the liquid volume sets the level through the filled angle β.

0

= U^{L} - n^{L} u^{L} (T^{L}, P, 𝐱)

0

= U^{V} - n^{V} u^{V} (T^{V}, P, 𝐲)

0

= V - n^{L} v^{L} (T^{L}, P, 𝐱) - n^{V} v^{V} (T^{V}, P, 𝐲)

0

= V_{fill} (β) - n^{L} v^{L} (T^{L}, P, 𝐱)

V_{fill} (β)

= \frac{L_{s} R_{s}^{2}}{2} (β - \sin β), ℓ = R_{s} (1 - \cos \frac{β}{2})

Tᴸ, Tⱽ: zone temperatures, K; algebraic unknowns
P: the common pressure, Pa; algebraic
β: filled angle of the horizontal shell, rad; algebraic. Unlike an arccos form of the level, these formulas are defined for any value a Newton iterate takes, and the liquid heel keeps the solution inside (0, 2π)
nᴸ, nⱽ: total moles of each zone, mol; x, y: their compositions, xᵢ = nᵢᴸ/nᴸ and yᵢ = nᵢⱽ/nⱽ
u, v: molar internal energy, J/mol, and molar volume, m³/mol, from the property package, each on its own zone’s root
V: drum volume, m³; with flat ends (heads are not modeled yet) it is πR_s²L_s, 18.85 m³ for V-101
L_s, R_s: shell length, 6 m, and radius, 1 m
ℓ: liquid level, m

These four rows are the UV problem. No separate flash solves them; they join the plant’s other equations in one Newton iteration. A standalone UV flash is planned.

On stable roots the rows are well posed. At fixed molar volume a phase’s internal energy rises with temperature, because c_v is positive, and its pressure falls as its volume grows. The four rows then have a nonsingular Jacobian near a solution, so the inventory fixes the temperatures, the pressure and the level.

Peng-Robinson

The molar properties come from the Peng-Robinson (PR) equation of state in its 1976 form, with the one-fluid quadratic mixing rule.

P

= \frac{R T}{v - b} - \frac{a (T)}{v^{2} + 2 b v - b^{2}}

a_{i} (T)

= 0.45724 \frac{R^{2} T_{c, i}^{2}}{P_{c, i}} α_{i} (T), b_{i} = 0.07780 \frac{R T_{c, i}}{P_{c, i}}

α_{i} (T)

= {[1 + κ_{i} (1 - \sqrt{T / T_{c, i}})]}^{2}, κ_{i} = 0.37464 + 1.54226 ω_{i} - 0.26992 ω_{i}^{2}

a

= \sum_{i} \sum_{j} x_{i} x_{j} (1 - k_{i j}) \sqrt{a_{i} a_{j}}, b = \sum_{i} x_{i} b_{i}

P, T, v: pressure, Pa; temperature, K; molar volume, m³/mol
R: molar gas constant, 8.314462618 J/(mol·K)
T_c,i, P_c,i, ωᵢ: critical temperature, K, critical pressure, Pa, and acentric factor of component i
aᵢ, bᵢ: attraction parameter, Pa·m⁶/mol², and covolume, m³/mol, of component i
αᵢ, κᵢ: alpha function and its slope parameter
a, b: the phase’s mixture attraction parameter and covolume
xᵢ: mole fractions of the phase
kᵢⱼ: binary interaction parameter, symmetric (kᵢⱼ = kⱼᵢ, kᵢᵢ = 0) and independent of temperature

The 1976 κ applies at every ω; the 1978 revision differs only above ω = 0.49, beyond n-decane’s 0.488. The kᵢⱼ are cited to their sources, and a missing pair fails the build unless the fluid system declares it zero, as the reference fluid does for 107 of its 190 pairs.

The engine solves the cubic in Z and takes the enthalpy as a departure from the ideal gas.

0

= Z^{3} - (1 - B) Z^{2} + (A - 3 B^{2} - 2 B) Z - (A B - B^{2} - B^{3})

Z

= \frac{P v}{R T}, A = \frac{a P}{{(R T)}^{2}}, B = \frac{b P}{R T}

h

= h^{ig} + R T (Z - 1) + \frac{T a^{'} - a}{2 \sqrt{2} b} \ln \frac{Z + (1 + \sqrt{2}) B}{Z + (1 - \sqrt{2}) B}, u = h - P v

h^{ig}

= \sum_{i} x_{i} \int_{T_{ref}}^{T} c_{p, i}^{ig} d T, \frac{c_{p, i}^{ig}}{R} = C_{0, i} + C_{1, i} T + C_{2, i} T^{2} + C_{3, i} T^{3} + C_{4, i} T^{4}

Z: compressibility factor
A, B: dimensionless attraction and covolume parameters
h, u: molar enthalpy and molar internal energy of the phase, J/mol
hⁱᵍ: ideal-gas enthalpy, J/mol, zero for each pure component at T_ref = 298.15 K
a′: da/dT, from analytic derivatives of the alpha functions
c_p,iⁱᵍ: ideal-gas heat capacity of component i, J/(mol·K), a polynomial whose coefficients C₀,ᵢ to C₄,ᵢ come from the component database

Cardano’s formula gives the roots, two Newton steps polish each, and only roots with Z > B count: the smallest for the liquid zone, the largest for the vapor. If the requested root is missing or near the spinodal, the call still returns a finite root, flagged as clamped. That keeps NaN out of the solver, but where a root vanishes the switch to the other one is a jump, so the argument above holds only while no root is clamped.

By default PR runs untranslated: the standard Peneloux shift, fitted for Soave-Redlich-Kwong, takes propane’s liquid from 6.5 to 15 percent too dense. Against reference equations for methane, propane and n-butane at reduced temperatures from 0.5 to 0.95, PR saturated liquid densities range from 12.8 percent high to 7.7 percent low, and vapor pressures fall within 0.9 percent above 0.6.

Between the zones

The zones are not forced into equilibrium, which keeps the equations well posed as a phase appears, disappears or runs low. Mass crosses at a finite rate driven by the fugacity difference, each zone’s coefficients taken at its own temperature.

\ln φ_{i}

= \frac{b_{i}}{b} (Z - 1) - \ln (Z - B) - \frac{A}{2 \sqrt{2} B} (\frac{2 \sum_{j} x_{j} a_{i j}}{a} - \frac{b_{i}}{b}) \ln \frac{Z + (1 + \sqrt{2}) B}{Z + (1 - \sqrt{2}) B}

D_{i}

= x_{i} φ_{i}^{L} (T^{L}, P, 𝐱) - y_{i} φ_{i}^{V} (T^{V}, P, 𝐲) = \frac{f_{i}^{L} - f_{i}^{V}}{P}

N_{i}

= k_{mt} [χ^{L} r (D_{i}) - χ^{V} r (- D_{i})]

Q_{int}

= U A_{int} (T^{V} - T^{L})

φᵢ: fugacity coefficient of component i in a phase, from Peng-Robinson at that phase’s root; in the first line a, b, A, B, Z and the xⱼ are that phase’s, and aᵢⱼ = (1 − kᵢⱼ)√(aᵢaⱼ)
Dᵢ: fugacity difference divided by the pressure, positive when component i tends to evaporate
fᵢᴸ, fᵢⱽ: fugacities of component i in the liquid and vapor zones, Pa
k_mt: interphase rate constant, mol/s; by default 100 mol/(s·m³) times the vessel volume, and tunable
r: smoothed positive part, equal to max(D, 0) outside a band of width 10⁻⁶ around zero and continuously differentiable, with r(D) − r(−D) = D exactly
χᴸ, χⱽ: zone availabilities, 0 at the minimum holdup and 1 from twice it up, a smooth step in between
UA_int: conductance between the zones, W/K

With both zones above twice their minimum holdup, χᴸ = χⱽ = 1 and Nᵢ = k_mt Dᵢ. A large k_mt drives the fugacities together, but equilibrium also needs a large UA_int to equalize the temperatures. A smaller k_mt gives real lags, such as slow re-flashing after a pressure drop.

The crossing moles carry the vapor zone’s molar enthalpy, an approximation to each component’s partial molar enthalpy that keeps total energy exact and puts the latent heat on the liquid. The work W = P dVᴸ/dt would raise the index; dummy derivatives, the temperature and pressure rates carried as algebraic unknowns, keep the index at 1 and the energy exact.

The outlet valve

Links set flows from the pressures at their ends. PV-101 follows the ISA-75.01 gas sizing equation, made smooth through zero flow and through choking.

\dot{m}

= \frac{N_{6} F_{p} C_{v} f (λ)}{3600} Y \sqrt{\frac{ρ_{1}}{1000}} s (Δ p_{eff}; δ) \sqrt{\frac{x_{eff}}{x_{p}}}

x_{p}

= \frac{{| Δ p_{eff} |}_{δ}}{p_{1}}, Y = 1 - \frac{x_{eff}}{3 F_{γ} x_{T}}, F_{γ} = \frac{γ}{1.40}

x_{eff}

= \frac{1}{2} (x_{p} + F_{γ} x_{T} - \sqrt{{(x_{p} - F_{γ} x_{T})}^{2} + ε_{x}^{2}})

s (Δ p; δ)

= \frac{Δ p}{{(Δ p^{2} + δ^{2})}^{1 / 4}}, {| Δ p |}_{δ} = \sqrt{Δ p^{2} + δ^{2}}

ṁ: mass flow, kg/s, of either sign
N₆ = 2.73: the standard’s constant for kg/h, kPa and kg/m³; the factors 3600 and 1000 convert to kg/s and Pa
F_p: piping geometry factor, 1; C_v: flow coefficient, 800 for PV-101
f(λ): inherent characteristic at valve position λ, linear for PV-101
p₁, ρ₁: pressure, Pa, and density, kg/m³, of the upstream end, blended smoothly over 10 Pa around zero pressure difference
Δp_eff: pressure difference across the valve, Pa; it equals p_from − p_to unless a connected holdup is at its minimum
x_p: pressure-drop ratio; x_T: pressure-differential ratio factor at choked flow, 0.7
F_γ: specific-heat-ratio factor; γ: ratio of specific heats, taken as c_p/(c_p − R) from the inlet’s molar heat capacity
Y: expansion factor
s, |·|_δ: regularized signed square root and smooth absolute value of a pressure difference Δp
δ = 1 Pa, ε_x = 10⁻⁴: smoothing widths, both registry parameters

The plain signed root of Δp has an infinite slope at zero, so Newton breaks down exactly where flows stop or reverse. The regularized root s has slope 1/√δ there and stays within a relative δ²/(4Δp²) of the plain root once |Δp| ≫ δ. Newton’s method, its reuse of one Jacobian across steps and the error control all need continuous derivatives, so every max, min, absolute value and switch that would put a kink in a residual is smoothed this way, its width a registry parameter. Links are quasi-steady, apart from the compressor’s duct inertia, so water hammer is outside the model.

One system

The compiler joins V-101, its valves, the pump, the compressor and the exchangers into one implicit differential-algebraic system of about 580 unknowns. The controllers run outside it, once per tick, and their outputs enter as inputs held over the next tick.

𝐅 (t, 𝐲, \dot{𝐲}; 𝜽, 𝐮_{k}) = 𝟎, t \in [t_{k}, t_{k + 1}], 𝐲 = (𝐲_{d}, 𝐲_{a})

\det \frac{\partial 𝐅}{\partial ({\dot{𝐲}}_{d}, 𝐲_{a})} \neq 0

{\dot{𝐲}}_{n} \approx \frac{1}{Δ t_{n}} \sum_{j = 0}^{q} α_{n, j} 𝐲_{n - j}, 1 \leq q \leq 5

𝐉 Δ 𝐲^{(m)} = - 𝐆 (𝐲_{n}^{(m)}), 𝐉 = \frac{\partial 𝐅}{\partial 𝐲} + \frac{α_{n, 0}}{Δ t_{n}} \frac{\partial 𝐅}{\partial \dot{𝐲}}

F: every residual row in the plant, one per unknown
y_d: differential variables, such as component moles, zone energies, liquid-full junction pressures, exchanger wall temperatures, valve positions, machine speeds and the compressor’s mass flow
y_a: algebraic variables, such as temperatures, gas and vessel pressures, fill angles, the flows of pipes, valves and pumps, and the vessel’s dummy derivatives
θ: parameters from the registry, constant between events
uₖ: inputs held over tick k (500 ms by default), such as valve and speed commands from the controllers, boundary pressures, temperatures and compositions, and equipment status
det ∂F/∂(ẏ_d, y_a) ≠ 0: the index-1 condition
Δtₙ, q, α_n,j: step size, s; order, 1 to 5; and coefficients of the backward differentiation formula
G: the residual with ẏ replaced by that formula
J, Δy⁽ᵐ⁾: Newton iteration matrix and the update at iteration m

The system is built to be index 1: the matrix in the determinant is nonsingular except at special states, such as a spinodal, where the integrator shortens the step. Elsewhere the implicit function theorem makes the algebraic variables and derivatives continuously differentiable functions of the differential states, inputs and parameters, on the branch through a consistent initial point. Locally the system is then an ordinary differential equation with a locally Lipschitz right side, which a bare square root would break at zero flow, so exactly one solution passes through that point. That is what deterministic means here.

IDA integrates it with backward differentiation formulas of order 1 to 5 and a modified Newton iteration, at a default relative tolerance of 10⁻⁸. The Jacobian is exact, from forward-mode automatic differentiation, and KLU factors it.

Discontinuities, such as a pump start, are applied only at ticks, followed by a consistent re-initialization. A restored snapshot reproduces a run bit for bit on the same build. Steady states come today from a settle run; a steady-state solver and a cold start from rest are still to come.

Deterministic, not predetermined

Predetermined means fixed in advance, whatever anyone does. Deterministic means the outcome follows from the state and the inputs. Change an input and the future changes with it, which is what makes a simulator useful.

An operator moves the PIC-101 set point, LV-101 fails closed, or the compressor trips. The set-point move reaches the equations through the controller’s output, as a change in u at the following ticks; the failure and the trip change u directly. Training, startup and shutdown, and control testing rest on this. A historian cannot show the result of a move nobody made.

What tuning determines

A stable controlled plant forgets small errors in its initial state. Start V-101’s level 5 cm high, and LIC-101 brings it back. Only an inventory that no controller or throughflow corrects keeps its error.

Parameter errors persist. If the model’s Cv for PV-101 is 10 percent high, PIC-101 still holds 28 bar, with the valve about 0.45 open where the plant’s sits at 0.50. Integral control moves the error from the pressure to the controller’s output, where it never decays. A wrong lag or metal mass biases the transients instead.

So tuning has to determine the parameters from the plant’s own records. It is in development: the estimator is built and under test, with no acceptance results yet. It solves weighted least squares in noise units, with priors, inside the registry’s bounds.

\hat{𝝃}

= \underset{𝝃_{\min} \leq 𝝃 \leq 𝝃_{\max}}{arg min} \sum_{e} \sum_{k} {(\frac{{\hat{y}}_{e, k} (𝜽 (𝝃)) - y_{e, k}}{σ_{e, k}})}^{2} + \sum_{j} {(\frac{ξ_{j} - ξ_{j}^{0}}{σ_{j}^{ξ}})}^{2}

ξ_{j}

\in {\ln θ_{j}, θ_{j}}

θ: physical parameters, such as valve Cv, exchanger UA, actuator lags and k_mt
ξ: fit coordinates; ξⱼ = ln θⱼ for most parameters, so that a step is relative, and ξⱼ = θⱼ for the rest
ŷ_e,k: the engine’s prediction of measurement k in experiment e, a steady average or a step-response trace
y_e,k: the prepared historian value
σ_e,k: noise standard deviation, from the record’s spread in steady periods
ξⱼ⁰, σⱼ^ξ: starting value and prior spread, from a datasheet, a document or an engine default
ξ_min, ξ_max: bounds from the engine’s parameter registry

The fit runs in stages: hydraulics, then heat transfer, then dynamics, each carrying the earlier results as priors. Identifiability is judged locally, from linearized sensitivities and the Fisher information. A parameter the data cannot see or separate stays out, and the record names the combination its partners carry, such as the sum of two series resistances with no pressure measured between them.

The limits

Deterministic here describes the averaged quantities a plant measures. Turbulence and bubble nucleation vary at small scales; they enter the equations through averages such as friction factors and transfer coefficients.

No model is perfect. Instruments drift, some inputs go unmeasured, and parameters change as exchangers foul and valves wear. The fit tests for drift and step changes, and it reports a mismatch the equations cannot represent, such as an undocumented bypass, instead of folding it into a Cv. So fidelity has to be measured on events the fit never saw. The planned fidelity report will do that unit by unit and state its limits.

Physics determines what a plant does next.