Distillation column algorithm

This document describes the mathematical model and solver implementations that power the DistillationColumn class in NeqSim. The implementation is centered on src/main/java/neqsim/process/equipment/distillation/DistillationColumn.java and the package-local solver and residual helper classes in the same package.

Governing equations

Each ideal-equilibrium tray satisfies the familiar MESH relationships:

  1. Total mass balance (tray j)

    $ V_{j-1} + L_{j+1} + F_j = V_j + L_j $

  2. Component balances

    $ V_{j-1} y_{i,j-1} + L_{j+1} x_{i,j+1} + F_j z_{i,j} = V_j y_{i,j} + L_j x_{i,j} $

  3. Phase equilibrium (K-values)

    $ y_{i,j} = K_{i,j} x_{i,j}, \qquad K_{i,j} = \frac{\hat f_{i,j}^{\text{vap}}}{\hat f_{i,j}^{\text{liq}}} $

  4. Energy balance

    $ V_{j-1} h_{j-1}^{V} + L_{j+1} h_{j+1}^{L} + F_j h_j^{F} + Q_j = V_j h_j^{V} + L_j h_j^{L} $

NeqSim evaluates fugacity-based K-values and molar enthalpies through the active SystemInterface. The matrix solver also uses linearized component balances in tridiagonal form:

$ A_j l_{i,j-1} + B_j l_{i,j} + C_j l_{i,j+1} = D_{i,j} $

with stripping factors $S_j = K_{i,j} V_j / L_j$ embedded in the diagonal terms.

Temperature updates rely on the log-Newton step derived from $\sum_i y_{i,j}=1$:

$ \Delta T_j = -\frac{\ln(\sum_i K_{i,j} x_{i,j}) R T_j^2}{h_j^{V} - h_j^{L}} $

The code limits $\Delta T_j$ to ±5 K and enforces bounds of 50–1000 K for numerical stability.

Column preparation

  1. Feed assignment – feeds are attached with addFeedStream; unassigned feeds are auto-placed near matching tray temperatures.
  2. Temperature seeding – init() runs the lowest feed tray, extrapolates temperatures towards condenser and reboiler, and links neighbouring trays with vapour/liquid streams.
  3. Pressure profile – prepareColumnForSolve() imposes a linear pressure drop between the configured bottom and top pressures (or inferred tray values when unspecified).

Solver implementations

Solver Class/Method Strategy Notes
DIRECT_SUBSTITUTION solveSequential() Classic two-sweep sequential substitution (liquids down, vapours up) with adaptive relaxation on temperatures and streams. Converges robustly for well-behaved systems; default choice.
DAMPED_SUBSTITUTION runDamped() Same equations as direct substitution but starts with a user-defined fixed relaxation factor before enabling adaptation. Useful for stiff columns where the default step overshoots.
INSIDE_OUT solveInsideOut() Quadrat-structure inside-out method: streams are relaxed against previous iterates while tray properties update using enthalpy-driven temperature corrections. Balances mass/energy less frequently to reduce cost and supports a polishing phase for tight tolerances.
WEGSTEIN solveWegstein() Wegstein acceleration on the sequential temperature map after direct-substitution warm-up. Speeds up well-conditioned fixed-point problems.
SUM_RATES solveSumRates() Flow-corrected tearing method that adjusts relaxation from tray sum-rate behavior. Useful for absorber and stripper style columns.
NEWTON solveNewton() Finite-difference Newton correction on tray temperatures with line search. A tray-temperature accelerator, not a full simultaneous MESH Newton solver.
MESH_RESIDUAL solveMeshResidual() Inside-out initialization, MESH residual evaluation, and Newton polishing if the residual monitor detects a poor residual state. Best for auditing MESH residuals and preparing cases for future rigorous residual-driven solvers.

Sequential substitution details

Complete usage example

import neqsim.process.equipment.distillation.DistillationColumn;
import neqsim.process.equipment.stream.Stream;
import neqsim.thermo.system.SystemInterface;
import neqsim.thermo.system.SystemSrkEos;

// 1. Create feed fluid (natural gas condensate at 216 K, 30 bar)
SystemInterface feed = new SystemSrkEos(216.0, 30.0);
feed.addComponent("methane", 0.5);
feed.addComponent("ethane", 0.2);
feed.addComponent("propane", 0.15);
feed.addComponent("n-butane", 0.05);
feed.addComponent("n-pentane", 0.05);
feed.addComponent("n-hexane", 0.03);
feed.addComponent("n-heptane", 0.02);
feed.setMixingRule("classic");

// 2. Wrap in a stream and set flow rate
Stream feedStream = new Stream("feed", feed);
feedStream.setFlowRate(100.0, "kg/hr");
feedStream.run();

// 3. Build a 5-tray deethanizer (with reboiler, no condenser)
DistillationColumn column = new DistillationColumn("Deethanizer", 5, true, false);
column.addFeedStream(feedStream, 5);             // feed on tray 5 (top)
column.getReboiler().setOutTemperature(378.15);   // 105 °C
column.setTopPressure(30.0);                      // bar
column.setBottomPressure(32.0);

// 4. Select solver and run
column.setSolverType(DistillationColumn.SolverType.INSIDE_OUT);
column.setMaxNumberOfIterations(50);
column.run();

// 5. Read results
System.out.println("Converged: " + column.solved());
System.out.println("Gas product:    " + column.getGasOutStream().getFlowRate("kg/hr") + " kg/hr");
System.out.println("Liquid product: " + column.getLiquidOutStream().getFlowRate("kg/hr") + " kg/hr");
System.out.println("Iterations:     " + column.getLastIterationCount());
System.out.println("Solve time:     " + column.getLastSolveTimeSeconds() + " s");

Algorithm overview

Matrix solver specifics

Result handling

Once any solver converges, the top gas outlet (gasOutStream) and bottom liquid outlet (liquidOutStream) are cloned from the respective trays. Mass, energy, and iteration statistics are exposed through getters such as getLastIterationCount(), getLastMassResidual(), and getLastEnergyResidual().

Every run() also records a scaled MESH residual vector for the final column state. The vector is assembled by ColumnMeshResidualEvaluator from a ColumnMeshState snapshot and groups equations as material, equilibrium, summation, energy, and specification residuals. Public diagnostics expose the full norm and group norms through getLastMeshResidualNorm(), getLastMeshMaterialResidualNorm(), getLastMeshEquilibriumResidualNorm(), getLastMeshSummationResidualNorm(), getLastMeshEnergyResidualNorm(), and getLastMeshSpecificationResidualNorm().

By default, the MESH residual vector is diagnostic only. setEnforceMeshResidualTolerance(true) makes solved() require the latest residual norm to be finite and less than getMeshResidualTolerance(), which is configured with setMeshResidualTolerance(double).

MESH equations

Each equilibrium stage $j$ in a column with $N$ stages and $n_c$ components satisfies:

Material balance (M):

\[M_{i,j} = L_{j-1} x_{i,j-1} + V_{j+1} y_{i,j+1} + F_j z_{i,j} - (L_j + S_j^L) x_{i,j} - (V_j + S_j^V) y_{i,j} = 0\]

Equilibrium (E):

\[E_{i,j} = K_{i,j}\, x_{i,j} - y_{i,j} = 0\]

Summation (S):

\[S_j = \sum_{i=1}^{n_c} y_{i,j} - \sum_{i=1}^{n_c} x_{i,j} = 0\]

Enthalpy balance (H):

\[H_j = L_{j-1} h_{j-1}^L + V_{j+1} h_{j+1}^V + F_j h_j^F - L_j h_j^L - V_j h_j^V - Q_j = 0\]

Rather than solving these equations algebraically, NeqSim uses a tray-by-tray flash approach: each tray mixes its input streams and performs a pressure–enthalpy (PH) flash. This automatically satisfies M, E, S, and H for any equation of state available in NeqSim (SRK, CPA, GERG-2008, etc.).

Available solver types

Solver Description Best for
DIRECT_SUBSTITUTION Classic tray-by-tray without damping (default) General use
DAMPED_SUBSTITUTION Adaptive relaxation controller Difficult polar/CPA systems
INSIDE_OUT Three-sweep IO with stripping factor correction and K-value tracking Multi-feed, general-purpose, debugging
WEGSTEIN Wegstein acceleration of successive substitution Fast convergence on well-posed problems
SUM_RATES Flow-corrected tearing method Absorbers and strippers
NEWTON Newton-Raphson tray-temperature correction accelerator Difficult temperature convergence cases
MESH_RESIDUAL Inside-out initialization with MESH residual diagnostics and Newton polishing Residual auditing and future rigorous solver preparation

Solver mathematics

DAMPED_SUBSTITUTION applies a relaxation factor $\alpha$ to stream and temperature updates:

\[T_j^{k+1} = T_j^k + \alpha (T_j^{flash} - T_j^k)\]

The factor adapts automatically: if the combined residual grows, $\alpha \leftarrow \max(\alpha_{min},\; 0.5\alpha)$; if it shrinks, $\alpha \leftarrow \min(\alpha_{max},\; 1.2\alpha)$.

WEGSTEIN accelerates successive substitution using the estimated slope of the fixed-point map:

\[T_j^{k+1} = (1 - q_j)\, g(T_j^k) + q_j\, T_j^k, \quad q_j = \frac{s_j}{s_j - 1}, \quad s_j = \frac{g(T_j^k) - g(T_j^{k-1})}{T_j^k - T_j^{k-1}}\]

with $q$ bounded to $[-2, 0]$ (more conservative than the classical $[-5, 0]$) to prevent divergence on multi-component systems. A warm-up phase of direct substitution iterations is run first to establish stable iterates.

SUM_RATES adjusts temperatures using flow correction factors:

\[\alpha_{eff} = \alpha \cdot \theta, \quad \theta = \frac{1}{\bar{r}}, \quad \bar{r} = \frac{1}{N}\sum_{j=1}^{N} \frac{L_j^{out} + V_j^{out}}{F_j^{in}}\]

NEWTON treats all $N$ tray temperatures as simultaneous variables for a temperature-correction accelerator. The residual is $f_i(\mathbf{T}) = T_i^{sweep} - T_i$ and the Jacobian is computed by finite-difference perturbation ($\epsilon = 0.1$ K):

\[J_{ij} \approx \frac{f_i(\mathbf{T} + \epsilon \mathbf{e}_j) - f_i(\mathbf{T})}{\epsilon}\]

A line search ($\lambda = 1, 0.5, 0.25, 0.125$) controls step size, and 2–3 warm-up direct substitution iterations establish the convergence basin.

MESH_RESIDUAL starts from INSIDE_OUT, evaluates the scaled MESH residual vector, and runs NEWTON polishing when the residual vector is non-finite or above the monitor threshold. It does not change the public meaning of NEWTON; instead, it provides a solver entry point that is explicitly organized around residual diagnostics.

Example: selecting a solver

column.setSolverType(DistillationColumn.SolverType.WEGSTEIN);
column.run();

Convergence diagnostics

After solving, the following metrics are available:

column.getLastIterationCount();        // number of iterations
column.getLastTemperatureResidual();   // avg temperature change (K)
column.getLastMassResidual();          // relative mass balance error
column.getLastEnergyResidual();        // relative enthalpy balance error
column.getLastMeshResidualNorm();      // full scaled MESH residual infinity norm
column.getLastSolveTimeSeconds();      // wall-clock time
column.getConvergenceHistory();        // per-iteration [temp, mass, energy] (IO adds K-value residual as 4th element)

Convergence criteria

All solvers use three residual metrics:

  1. Temperature: $\varepsilon_T = \frac{1}{N}\sum_{j=1}^{N} T_j^{new} - T_j^{old} $ (K)
  2. Mass balance: $\varepsilon_M = \max!\Big(\max_j \frac{ M_j^{in} - M_j^{out} }{M_j^{in}},\; \frac{\sum M_j^{in}-M_j^{out} }{\sum M_j^{in}}\Big)$
  3. Energy balance: analogous per-tray and column-wide relative error
  4. Optional MESH residual: infinity norm of the scaled material, equilibrium, summation, energy, and specification residual vector

Tolerances scale with column complexity:

\[\tau = \tau_{base} \cdot \min\!\Big(2.5,\; \max\big(1 + 0.06(N_{stages}-3),\; 1 + 0.25(N_{feeds}-1)\big)\Big)\]

Plotting convergence history

// After column.run()
java.util.List<double[]> history = column.getConvergenceHistory();
for (int k = 0; k < history.size(); k++) {
    double[] h = history.get(k);
    System.out.printf("iter %3d  tempErr=%.2e  massErr=%.2e  energyErr=%.2e%n",
        k + 1, h[0], h[1], h[2]);
}

Inside-Out solver details

The INSIDE_OUT solver performs three sweep phases per iteration (liquid down, vapor up, polish liquid down) with:

\[T_j^{k+1} = T_j^k + \alpha \cdot \beta_j \cdot (T_j^{flash} - T_j^k), \quad \beta_j = 1 + 0.05\,\text{clamp}\!\Big(\ln\frac{V_j}{L_j},\, -1,\, 1\Big)\] \[\varepsilon_K = \max_{i,j} \frac{|K_{i,j}^k - K_{i,j}^{k-1}|}{K_{i,j}^{k-1}}\] \[\Delta T = -\frac{\Sigma - 1}{\sum_i \frac{-b_{i,j}}{T^2} K_i x_i}\]

The convergence history for IO includes a 4th element: [tempErr, massErr, energyErr, kValueResidual].

Configuring the inner loop

column.setSolverType(DistillationColumn.SolverType.INSIDE_OUT);
column.setInnerLoopSteps(3);  // default: 3 cheap iterations between outer flash updates
column.setInnerLoopSteps(0);  // disable simplified model (all iterations use rigorous flash)
column.setInnerLoopSteps(5);  // more inner steps = fewer flashes but less correction per step

Murphree tray efficiency

To model non-ideal trays, set the column-wide Murphree efficiency:

column.setMurphreeEfficiency(0.75); // 75% efficiency

The correction is applied after each tray flash using the standard post-correction formula:

\[y_{i,j}^{out} = y_{i,j-1} + E_{MV} \cdot (y_{i,j}^{eq} - y_{i,j-1})\]

where $y_{i,j}^{eq}$ is the equilibrium composition from the flash, $y_{i,j-1}$ is the inlet vapor from the tray below, and $E_{MV}$ is the Murphree tray efficiency. The reboiler and condenser (when present) are always treated as equilibrium stages.

Newton-Raphson solver details

The NEWTON solver treats all tray temperatures as simultaneous variables and computes a Jacobian by finite-difference perturbation. Each iteration requires $N+1$ column sweeps (1 base + $N$ perturbations). The resulting linear system is solved by Gaussian elimination with partial pivoting and a backtracking line search controls step size.

column.setSolverType(DistillationColumn.SolverType.NEWTON);
column.run();

Suggestions for future improvement