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.

The matrix inside-out warm-start updates tray temperatures from the bubble-sum residual $\sum_i K_{i,j} x_{i,j} - 1$ and cached derivatives of $\ln K$ with respect to temperature:

\[\Delta T_j = -\frac{\sum_i K_{i,j} x_{i,j} - 1} {\sum_i x_{i,j} K_{i,j} \frac{\partial \ln K_{i,j}}{\partial T}}\]

The code limits $\Delta T_j$ to ±8 K and enforces bounds of 50–1000 K for numerical stability during the matrix warm-start stage.

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.
MATRIX_INSIDE_OUT solveMatrixInsideOut() Adaptive matrix inside-out mode that bypasses matrix setup for small columns and otherwise solves component-balance tridiagonal systems as a warm start before rigorous inside-out polishing. Avoids warm-start overhead on small columns while preserving a matrix path for larger hydrocarbon fractionators.
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.
NAPHTALI_SANDHOLM solveNaphtaliSandholm() Inside-out warm start followed by guarded simultaneous Newton correction of liquid component flows, tray temperatures, and vapor flows. Best for rigorous residual-driven MESH convergence on well-conditioned hydrocarbon fractionators.
MESH_RESIDUAL solveMeshResidual() Inside-out initialization followed by full MESH residual evaluation. Best for auditing material, equilibrium, summation, energy, specification, and product-draw residuals.
AUTO ColumnSolverFactory.AutoSolver Runs a feasibility pre-screen, initializes a copied candidate, solves a relaxed damped base case, probes built-in candidate solvers on copies, and accepts the first solved non-fallback candidate or best valid fallback. Useful when an agent or workflow should request robust automatic solver selection while still reporting the concrete solver through getLastSolverTypeUsed().

Sequential substitution details

Automatic solver pipeline

AUTO mode is intended for workflows where the caller wants a robust answer and diagnostics rather than a specific numerical method. The pipeline first calls screenSpecificationFeasibility() and adds a FEASIBILITY_SCREEN entry to the automatic solver summary. It then creates a candidate copy, tries shortcut or thermodynamic-profile initialization, runs a relaxed DAMPED_SUBSTITUTION base solve, and probes an ordered candidate list on copies of that warmed base state.

The candidate order depends on column structure:

Candidate probes do not run their own nested damped fallback solves. AUTO accepts the first solved candidate that did not rely on guarded fallback products; otherwise it scores the available results by residuals and falls back to a live damped-substitution solve when necessary. The public solver request remains AUTO, while getLastSolverTypeUsed() reports the concrete accepted solver and getLastAutoSolverSummary() records the candidate trace.

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

Outer tear variables

Some column features are coupled outside the inner tray solver because they change the traffic or pressure profile that the tray sweeps depend on:

solveWithColumnTearVariables() wraps the selected inner solver for these features. It records getLastColumnTearIterationCount(), getLastColumnTearResidual(), isLastColumnTearConverged(), getLastPumparoundRelativeChange(), getLastHydraulicPressureDropPa(), and getLastHydraulicPressureDropResidual(). If a side-draw target is physically impossible, the draw fraction is bounded by available tray traffic and the latest tear diagnostics report non-convergence instead of allowing an impossible product draw.

Specification homotopy

Purity, recovery, and product-flow specifications can be difficult if the outer loop jumps directly from the current product split to a sharp target. setSpecificationHomotopySteps(steps) stages the effective target over the requested number of continuation steps. Each stage solves with a temporary target between the current product value and the final user target; the stored public ColumnSpecification still contains the final target. getLastSpecificationHomotopyStepCount() reports how many staged targets were completed by the latest run.

AUTO mode enables three specification-homotopy stages automatically for adjustable product specifications when the user has not configured a larger or smaller value. Reflux-ratio and duty specifications are not staged because they are direct operating specifications rather than product targets manipulated through condenser or reboiler temperature.

Matrix solver specifics

Naphtali-Sandholm 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, product-draw, and specification residuals. Public diagnostics expose the full norm and group norms through getLastMeshResidualNorm(), getLastMeshMaterialResidualNorm(), getLastMeshEquilibriumResidualNorm(), getLastMeshSummationResidualNorm(), getLastMeshEnergyResidualNorm(), getLastMeshProductDrawResidualNorm(), and getLastMeshSpecificationResidualNorm().

For MATRIX_INSIDE_OUT, diagnostics also expose the adaptive warm-start decision: wasMatrixInsideOutWarmStartUsed(), wasMatrixInsideOutWarmStartBypassed(), getLastMatrixInsideOutIterationCount(), getLastMatrixInsideOutTemperatureResidual(), and getLastMatrixInsideOutSolveTimeSeconds(). These values make it clear whether a run paid the matrix setup cost or intentionally followed the rigorous inside-out path directly.

By default, the MESH residual vector is diagnostic for sequential solver modes. For NAPHTALI_SANDHOLM and MESH_RESIDUAL, the residual gate is effective unless explicitly disabled. When the gate is effective, solved() requires the latest residual norm to be finite and less than getMeshResidualTolerance(), and also requires the product-draw residual to pass getMeshProductDrawResidualTolerance(). For legacy solver modes, calling setEnforceMeshResidualTolerance(true) opts into the same convergence contract; calling it with false explicitly disables the gate even for the residual-oriented solvers.

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\]

Most NeqSim column solvers use a tray-by-tray flash approach: each tray mixes its input streams and performs a pressure-enthalpy (PH) flash. This provides a robust MESH-consistent update for any equation of state available in NeqSim (SRK, CPA, GERG-2008, etc.). The NAPHTALI_SANDHOLM solver adds a simultaneous residual-correction layer on top of that warm start.

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
MATRIX_INSIDE_OUT Adaptive matrix warm start plus rigorous inside-out polish; bypasses matrix setup for small columns Larger hydrocarbon columns where the matrix warm start can help
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
NAPHTALI_SANDHOLM Simultaneous MESH residual Newton correction with guarded acceptance Rigorous residual convergence checks
MESH_RESIDUAL Inside-out initialization with MESH residual diagnostics Residual auditing and diagnostics
AUTO Feasibility-screened automatic candidate selection with copy-based solver probes and damped fallback Agent workflows and uncertain column setups

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 and evaluates the scaled MESH residual vector without running an additional Newton-polishing solve. When the residual or product-draw gate is not satisfied, the current implementation leaves the inside-out state in place and reports the failed gate through solved(), the residual getters, and getConvergenceDiagnostics(). This keeps MESH_RESIDUAL as a diagnostics-oriented entry point; use NAPHTALI_SANDHOLM when a guarded simultaneous residual correction should be attempted.

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 track three scalar residual metrics, and the residual-oriented modes also gate on the MESH diagnostics:

  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, product-draw, 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