NeqSim traces the PT phase envelope using Michelsen’s continuation method (Michelsen & Mollerup, Thermodynamic Models: Fundamentals & Computational Aspects, 2nd ed., 2007). This document describes the mathematical formulation, the Newton-Raphson solver, critical point detection, quality line tracing, and the adaptive step-size control implemented in PTPhaseEnvelopeMichelsen and SysNewtonRhapsonPhaseEnvelope.

1. Problem Formulation

At every point on the phase boundary, the following $n_c + 2$ equations hold simultaneously:

\[\ln K_i + \ln \hat{\varphi}_i^{\text{vap}}(T, P, \mathbf{y}) - \ln \hat{\varphi}_i^{\text{liq}}(T, P, \mathbf{x}) = 0, \quad i = 1, \dots, n_c\] \[\sum_i y_i - \sum_i x_i = 0\] \[u_{s} - S = 0 \qquad \text{(specification equation)}\]

where:

The unknowns are collected into a vector $\mathbf{u}$:

\[\mathbf{u} = [\ln K_1, \ln K_2, \dots, \ln K_{n_c}, \ln T, \ln P]^T\]

Using logarithmic variables ensures positivity of $K$, $T$, and $P$ naturally.

Phase Compositions

The liquid and vapor mole fractions are computed from the Rachford-Rice equation at a given vapor mole fraction $\beta$:

\[x_i = \frac{z_i}{1 - \beta + \beta K_i}, \qquad y_i = \frac{K_i z_i}{1 - \beta + \beta K_i}\]

For the dew point curve, $\beta \to 1$ (i.e., phaseFraction $\approx 1 - 10^{-10}$). For the bubble point curve, $\beta \to 0$ (i.e., phaseFraction $\approx 10^{-10}$). For quality lines, $0 < \beta < 1$.

2. Jacobian Matrix

The Jacobian $\mathbf{J}$ of the system has $(n_c + 2) \times (n_c + 2)$ entries:

Rows 1 to $n_c$ (fugacity equations):

\[J_{ij} = \delta_{ij} + \sum_k \left( \frac{\partial \ln \hat{\varphi}_i^{\text{vap}}}{\partial y_k} \frac{\partial y_k}{\partial \ln K_j} - \frac{\partial \ln \hat{\varphi}_i^{\text{liq}}}{\partial x_k} \frac{\partial x_k}{\partial \ln K_j} \right)\] \[J_{i, n_c+1} = T \left( \frac{\partial \ln \hat{\varphi}_i^{\text{vap}}}{\partial T} - \frac{\partial \ln \hat{\varphi}_i^{\text{liq}}}{\partial T} \right)\] \[J_{i, n_c+2} = P \left( \frac{\partial \ln \hat{\varphi}_i^{\text{vap}}}{\partial P} - \frac{\partial \ln \hat{\varphi}_i^{\text{liq}}}{\partial P} \right)\]

Row $n_c + 1$ (material balance):

\[J_{n_c+1, j} = \frac{\partial y_j}{\partial \ln K_j} - \frac{\partial x_j}{\partial \ln K_j}\]

Row $n_c + 2$ (specification):

\[J_{n_c+2, s} = 1, \quad J_{n_c+2, j \ne s} = 0\]

3. Continuation Method

3.1 Sensitivity Vector

At each converged point, the sensitivity vector $d\mathbf{u}/ds$ is computed by solving:

\[\mathbf{J} \cdot \frac{d\mathbf{u}}{ds} = \mathbf{e}_{n_c+2}\]

where $\mathbf{e}_{n_c+2}$ is a unit vector with 1 in the last row (specification equation).

3.2 Specification Variable Selection

Michelsen (2007) recommends choosing the specification variable as the component of $d\mathbf{u}/ds$ with the largest absolute value:

\[s = \arg\max_j \left| \frac{du_j}{ds} \right|\]

This keeps the Jacobian well-conditioned near turning points (cricondenbar, cricondentherm) and the critical region where the standard choices ($\ln P$ or $\ln T$) would lead to singular or near-singular systems.

In NeqSim, the first 5 points use $\ln P$ as the specification variable (simple pressure stepping from the starting low pressure), then switch to the sensitivity-based selection.

3.3 Polynomial Extrapolation

The next point estimate uses cubic polynomial extrapolation from the last 4 converged points. Given stored solutions $\mathbf{u}_1, \mathbf{u}_2, \mathbf{u}_3, \mathbf{u}_4$ at specification values $s_1, s_2, s_3, s_4$, the cubic polynomial coefficients are found by solving:

\[\begin{bmatrix} 1 & s_1 & s_1^2 & s_1^3 \\ 1 & s_2 & s_2^2 & s_2^3 \\ 1 & s_3 & s_3^2 & s_3^3 \\ 1 & s_4 & s_4^2 & s_4^3 \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \end{bmatrix} = \begin{bmatrix} u_j^{(1)} \\ u_j^{(2)} \\ u_j^{(3)} \\ u_j^{(4)} \end{bmatrix}\]

The next specification value is $s_5 = s_4 + \Delta s \cdot (d\mathbf{u}/ds)_s$, and the initial estimate is:

\[u_j^{(5)} = c_0 + c_1 s_5 + c_2 s_5^2 + c_3 s_5^3\]

3.4 Adaptive Step Size

The arc-length step $\Delta s$ is adjusted adaptively based on Newton convergence:

\[\Delta s_{\text{new}} = \Delta s_{\text{old}} \cdot \sqrt{\frac{n_{\text{target}}}{n_{\text{actual}}}}\]

where $n_{\text{target}} = 4$ (target Newton iterations) and $n_{\text{actual}}$ is the actual iteration count. The scale factor is clamped to $[0.25, 2.0]$ to avoid excessively large or small steps.

Additionally, absolute step limits enforce:

\[|\Delta \ln T| \le \ln\left(1 + \frac{\Delta T_{\max}}{T}\right), \qquad |\Delta \ln P| \le \ln\left(1 + \frac{\Delta P_{\max}}{P}\right)\]

with $\Delta T_{\max} = 10$ K and $\Delta P_{\max} = 10$ bar.

4. Newton-Raphson Solver

Each point on the envelope is converged using Newton-Raphson iteration with damped steps and iterative backtracking:

4.1 Damped Newton Step

At each iteration $k$:

\[\mathbf{u}^{(k+1)} = \mathbf{u}^{(k)} - \lambda \cdot \mathbf{J}^{-1} \mathbf{f}\]

The damping factor is $\lambda = \min(1, 1/\max(1, |\Delta \mathbf{u}|_\infty))$, which limits each logarithmic correction to one before re-evaluating the residual. It prevents a single Newton update from making an unbounded multiplicative change in temperature, pressure or K-values.

4.2 Convergence Criteria

Convergence requires both norms to be small (dual criterion):

\[\|\mathbf{J}^{-1} \mathbf{f}\|_2 < 10^{-5} \quad \text{AND} \quad \|\mathbf{f}_{\mathrm{updated}}\|_2 < 10^{-8}\]

The dual check prevents false convergence where the correction is small but the residual is not.

4.3 Backtracking

If Newton iteration does not converge within 50 iterations, the solver backtracks to the last converged point and halves the step size:

\[\Delta s_{\text{bt}} = \Delta s \cdot 0.5^k, \quad k = 1, 2, \dots, 15\]

Retries use the tangent predictor until four history points are available, then use the cubic predictor. A singular cubic history falls back to the current tangent, never coefficients from a different variable. Non-finite residuals, trivial multicomponent K=1 solutions and stalled specification steps cannot be accepted. Only converged states update the solution history. Corrected points may not jump more than three configured temperature or pressure steps from the previous accepted point. A pressure or point-limit exit throws and leaves extrema unavailable; the cutoff is not a computed cricondenbar.

After 15 failed backtrack attempts, an IllegalStateException signals failure to PTPhaseEnvelopeMichelsen, which then either terminates the current tracing pass or starts a second pass from the opposite end.

5. Two-Pass Envelope Tracing

The envelope is traced in up to two passes:

Pass 1: Start from the specified end (dew or bubble) at low pressure but trace towards and through the critical point to the other branch.

Pass 2 (if needed): If Pass 1 fails before completing (solver divergence near the critical region), the algorithm:

  1. Saves all cricondenbar/cricondentherm data from Pass 1
  2. Flips bubblePointFirst and phaseFraction
  3. Restarts from the opposite end

The second pass can recover the other branch when continuation fails near the critical point, but does not guarantee that the two branches meet. If neither pass yields any converged point, tracing throws IllegalStateException and returns no invented extrema. The physical branch labels are independent of which side was traced first. The restart is not stopped by comparing its temperature with the last temperature of the failed pass. A dew branch can be warmer than the entire low-pressure bubble segment without overlapping it, as in the SAFT-VR Mie methane/propane regression. Each pass follows its own continuation and pressure termination criteria.

6. Critical Point Detection

6.1 Criterion

For the first multicomponent crossing, the tracer monitors the lightest and heaviest components: $K_{\mathrm{light}} < 1.05$ and $K_{\mathrm{heavy}} > 0.95$. After that crossing, another critical point requires both ratios to approach unity from either direction: $|\ln K| < \ln(1.05)$. The ratios must first diverge, at least 50 points must separate the detections, and pressure must exceed 5 bara. A light-component ratio far below unity and a heavy-component ratio far above unity are not evidence of another critical point. These heuristics do not replace a stability analysis. For a pure component, $K=1$ holds along the entire saturation curve and cannot identify criticality. The tracer therefore skips this K-based detection and rejects candidate saturation points whose two volume roots coincide (relative volume separation at most 1e-6).

6.2 Newton Refinement

After detection, the critical point estimate from polynomial extrapolation is refined by Newton-Raphson iteration. The specification equation (last row) is replaced with the critical constraint:

\[f_{n_c+2} = \sum_{i=1}^{n_c} (\ln K_i)^2 = 0\]

The corresponding Jacobian row becomes:

\[J_{n_c+2, j} = \begin{cases} 2 \ln K_j & j \le n_c \\ 0 & j = n_c+1 \text{ or } n_c+2 \end{cases}\]

This gives a well-conditioned system at the critical point where the standard specification equation becomes degenerate (the sensitivity of any single variable vanishes as all K-values converge).

Up to 10 Newton iterations are performed with convergence tolerance $|\mathbf{f}|_2 < 10^{-8}$. Refinement evaluates trial states on a private clone. Only the estimated critical temperature and pressure are transferred back; the accepted continuation state, phase properties and material-balance sums are preserved.

7. Envelope Closure Check

After tracing completes, the envelope is considered closed if both branches contain at least 3 finite points (NaN separators do not count):

\[\text{closed} = (n_{\text{dew}} \ge 3) \quad \text{AND} \quad (n_{\text{bub}} \ge 3)\]

This legacy check is available via isEnvelopeClosed(). It also returns false after a pressure or point-limit exit. The point-count heuristic alone does not prove that independently restarted segments meet or that the full boundary is stable. Inspect the structured segments before treating a partial trace as a complete dew-bubble-critical loop.

8. Quality Lines

Quality lines are curves of constant molar vapor fraction $\beta$ inside the two-phase region. They connect the bubble curve ($\beta = 0$) to the dew curve ($\beta = 1$).

8.1 Tracing Algorithm

Each quality line is traced using the same continuation method as the phase boundary, but with fixed $\beta$:

  1. Initialization: Wilson K-value correlation estimates the starting temperature at low pressure for the specified $\beta$. A flash calculation (bubble-point if $\beta < 0.5$, dew-point if $\beta \ge 0.5$) refines the K-values.

  2. Continuation: The Newton-Raphson solver traces the curve in $P$-$T$ space using the same Jacobian, but the Rachford-Rice equation uses the specified $\beta$ instead of $\beta \to 0$ or $\beta \to 1$.

  3. Exit criteria: Tracing stops when any of:

    • Pressure exceeds the maximum envelope pressure
    • Pressure drops below minimum after passing a peak (the curve turns around)
    • The solver fails to converge
    • 5000 points are traced

8.2 Volume and Mass Fractions

At each point on a quality line, the volume fraction $\beta_V$ and mass fraction $\beta_W$ are computed from the molar vapor fraction $\beta$:

\[\beta_V = \frac{\beta \cdot V_m^{\text{vap}}}{\beta \cdot V_m^{\text{vap}} + (1-\beta) \cdot V_m^{\text{liq}}}\] \[\beta_W = \frac{\beta \cdot M_w^{\text{vap}}}{\beta \cdot M_w^{\text{vap}} + (1-\beta) \cdot M_w^{\text{liq}}}\]

where $V_m = M_w / \rho$ is the molar volume and $M_w$ is the phase molar mass.

These three fractions differ because the vapor and liquid phases have different densities and molar masses. For example, at $\beta = 0.5$ (50% of moles are vapor):

8.3 Data Access

Quality line data is accessed via the get() method with string keys:

Key Pattern Returns Example
"qualityT_X" Temperatures (K) along quality line $\beta = X$ "qualityT_0.5"
"qualityP_X" Pressures (bara) along quality line $\beta = X$ "qualityP_0.5"
"qualityVolFrac_X" Volume fractions along the line "qualityVolFrac_0.5"
"qualityMassFrac_X" Mass fractions along the line "qualityMassFrac_0.5"

9. Starting Temperature from Wilson K-Values

The initial temperature estimate at a given pressure $P$ and vapor fraction $\beta$ is found by solving:

\[g(\beta, T) = \sum_i \frac{z_i (K_i^W - 1)}{1 + \beta(K_i^W - 1)} = 0\]

where the Wilson K-values are:

\[K_i^W = \frac{P_{c,i}}{P} \exp\left[5.373(1 + \omega_i)\left(1 - \frac{T_{c,i}}{T}\right)\right]\]

Newton’s method is used: starting from $T_0 \approx T_{c,\text{ref}} \cdot 5.373(1+\omega_{\text{ref}}) / [5.373(1+\omega_{\text{ref}}) - \ln(P/P_{c,\text{ref}})]$ where the reference component is the lightest (for bubble) or heaviest (for dew).

10. Cricondenbar and Cricondentherm Tracking

These extrema are tracked incrementally during tracing:

Both include the phase compositions $\mathbf{x}$ and $\mathbf{y}$ at the extremum, accessible via "cricondenbarX", "cricondenbarY", "cricondenthermX", "cricondenthermY".

11. Direct Cricondenbar and Cricondentherm Refinement

After the phase envelope has been traced, the discrete maximum-pressure and maximum-temperature estimates can be refined to machine precision using the Michelsen simultaneous Newton method (Michelsen, 1980; Michelsen & Mollerup, 2007, Chapter 12).

11.1 Mathematical Formulation

At the cricondenbar (maximum pressure on the envelope), the temperature is stationary: $dT/dP = 0$. At the cricondentherm (maximum temperature), the pressure is stationary: $dP/dT = 0$. These extremum conditions can be expressed as additional equations in a simultaneous Newton system.

The unknowns are the same $(n_c + 2)$ variables as the envelope tracer:

\[\mathbf{u} = [\ln K_1, \ln K_2, \dots, \ln K_{n_c}, \ln T, \ln P]^T\]

The equations are:

\[g_i = \ln K_i + \ln \hat{\varphi}_i^{\text{vap}} - \ln \hat{\varphi}_i^{\text{liq}} = 0, \quad i = 1, \dots, n_c\] \[g_{n_c+1} = \sum_i \frac{z_i(K_i - 1)}{1 + \beta(K_i - 1)} = 0 \quad \text{(Rachford-Rice)}\] \[g_{n_c+2} = S = 0 \quad \text{(extremum specification)}\]

where the specification $S$ differs for each extremum type.

11.2 Cricondenbar Specification ($S_P = 0$)

At the cricondenbar, if we were to use $\ln P$ as the specification variable in the envelope equations, the resulting $(n_c+1) \times (n_c+1)$ system in $(\ln K, \ln T)$ would be singular — there is no unique temperature increment for a given pressure increment at the turning point. This singularity condition is:

\[S_P = \sum_{i=1}^{n_c} s_i \cdot \frac{\partial g_i}{\partial \ln T} = 0\]

where $s_i = (y_i - x_i) / |\mathbf{y} - \mathbf{x}|$ is a normalized direction vector and:

\[\frac{\partial g_i}{\partial \ln T} = T \left( \frac{\partial \ln \hat{\varphi}_i^{\text{vap}}}{\partial T} - \frac{\partial \ln \hat{\varphi}_i^{\text{liq}}}{\partial T} \right)\]

11.3 Cricondentherm Specification ($S_T = 0$)

Symmetrically, at the cricondentherm the specification equation uses the pressure derivative:

\[S_T = \sum_{i=1}^{n_c} s_i \cdot \frac{\partial g_i}{\partial \ln P} = 0\]

where:

\[\frac{\partial g_i}{\partial \ln P} = P \left( \frac{\partial \ln \hat{\varphi}_i^{\text{vap}}}{\partial P} - \frac{\partial \ln \hat{\varphi}_i^{\text{liq}}}{\partial P} \right)\]

11.4 Jacobian Structure

The $(n_c+2) \times (n_c+2)$ Jacobian is:

\[\mathbf{J} = \begin{bmatrix} \frac{\partial g_1}{\partial \ln K_1} & \cdots & \frac{\partial g_1}{\partial \ln K_{n_c}} & \frac{\partial g_1}{\partial \ln T} & \frac{\partial g_1}{\partial \ln P} \\ \vdots & \ddots & \vdots & \vdots & \vdots \\ \frac{\partial g_{n_c}}{\partial \ln K_1} & \cdots & \frac{\partial g_{n_c}}{\partial \ln K_{n_c}} & \frac{\partial g_{n_c}}{\partial \ln T} & \frac{\partial g_{n_c}}{\partial \ln P} \\ \frac{\partial g_{\text{RR}}}{\partial \ln K_1} & \cdots & \frac{\partial g_{\text{RR}}}{\partial \ln K_{n_c}} & 0 & 0 \\ \frac{\partial S}{\partial \ln K_1} & \cdots & \frac{\partial S}{\partial \ln K_{n_c}} & \frac{\partial S}{\partial \ln T} & \frac{\partial S}{\partial \ln P} \end{bmatrix}\]

The first $(n_c+1)$ rows use fully analytical derivatives from system.init(3) (fugacity coefficient derivatives getdfugdt(), getdfugdp(), getdfugdx(j)). The last row (extremum specification) uses numerical perturbation ($\epsilon = 10^{-6}$) because the mixed second derivatives $\partial^2 g_i / (\partial \ln T \, \partial \ln K_j)$ involve complex chain-rule expressions.

The composition derivatives of the K-value equations use the Rachford-Rice differentiation:

\[\frac{\partial y_j}{\partial \ln K_j} = \frac{z_j K_j (1 - \beta)}{(1 - \beta + \beta K_j)^2}, \qquad \frac{\partial x_j}{\partial \ln K_j} = \frac{-z_j \beta K_j}{(1 - \beta + \beta K_j)^2}\]

11.5 Convergence Properties

The simultaneous Newton system gives quadratic convergence, typically converging in 3-8 iterations from the discrete envelope estimate. This is a major improvement over the previous alternating Newton approach which had linear convergence and could require hundreds of iterations.

Key numerical features:

11.6 Usage

The refinement is invoked automatically when calling ThermodynamicOperations.calcCricoP() or ThermodynamicOperations.calcCricoT() after a phase envelope has been traced. It requires an initial estimate from the envelope and the corresponding phase compositions.

// After phase envelope calculation
ThermodynamicOperations ops = new ThermodynamicOperations(fluid);
ops.calcPTphaseEnvelope();

// Refine cricondenbar (internally calls CricondenBarFlash)
ops.calcCricoP(dewT, dewP, bubT, bubP, cricondenbar, cricondenbarX, cricondenbarY);

// Refine cricondentherm (internally calls CricondenThermFlash)
ops.calcCricoT(dewT, dewP, bubT, bubP, cricondentherm, cricondenthermX, cricondenthermY);

12. Implementation Classes

Class Responsibility
PTPhaseEnvelopeMichelsen Orchestration: two-pass tracing, point storage, critical/cricondenbar/cricondentherm tracking, quality line management, data access via get()
SysNewtonRhapsonPhaseEnvelope Newton-Raphson solver: Jacobian assembly, sensitivity vector, specification variable selection, polynomial extrapolation, step-size control, critical point refinement
CricondenBarFlash Direct cricondenbar refinement: Michelsen simultaneous $(n_c+2)$ Newton system with $S_P = 0$ specification
CricondenThermFlash Direct cricondentherm refinement: Michelsen simultaneous $(n_c+2)$ Newton system with $S_T = 0$ specification
ThermodynamicOperations Public API: calcPTphaseEnvelope(), calcCricoP(), calcCricoT(), and calcPTphaseEnvelopeWithQualityLines()

References

  1. Michelsen, M.L. & Mollerup, J.M., Thermodynamic Models: Fundamentals & Computational Aspects, 2nd ed., Tie-Line Publications, 2007.
  2. Michelsen, M.L., “Calculation of phase envelopes and critical points for multicomponent mixtures”, Fluid Phase Equilibria, 4(1-2), 1-10, 1980.
  3. Wilson, G.M., “A modified Redlich-Kwong equation of state, application to general physical data calculations”, Paper No. 15C, AIChE 65th National Meeting, Cleveland, Ohio, 1968.
  4. Hoteit, H., Segtovich, I., Freistühler, A., “An efficient and robust algorithm for the calculation of gas-liquid critical point of multicomponent fluid mixtures”, Fluid Phase Equilibria, 241(1-2), 186-195, 2006.
  5. Nikolaidis, I.K., Experiment-based assessment of the accuracy of several EoS in predicting the phase envelope of natural gas, Paper presented at SPE Annual Technical Conference and Exhibition, 2016.

See Also