Numerical Methods¶

Overview¶

The MSHX solver relies on several numerical techniques to handle the nonlinear, multi-variable nature of multi-stream heat exchange with phase change. This section documents the key algorithms and their convergence properties.

Q-Bisection Solver¶

The primary solution strategy is bisection on total heat duty \(Q\). This approach was chosen over temperature-based iteration because:

\(Q\) varies monotonically with the UA or MITA constraint
\(Q\) is smooth even when temperature has discontinuities (phase change)
Bisection is globally convergent — it never diverges

Algorithm¶

Q_lo ← 0
Q_hi ← Q_max  (minimum of hot-side and cold-side maximum heat)

for iter = 1 to 300:
    Q_trial ← (Q_lo + Q_hi) / 2

    Distribute Q_trial among streams
    PH flash → compute outlet temperatures
    Build composite curves
    Calculate UA_calc and MITA_calc

    if |UA_calc - UA_spec| / UA_spec < 0.001:  → converged
    if UA_calc < UA_spec:  Q_lo ← Q_trial  (need more heat)
    else:                  Q_hi ← Q_trial  (need less heat)

    if |Q_hi - Q_lo| / Q_hi < 1e-5:  → converged (bounds collapsed)

Convergence Rate¶

Bisection has linear convergence with a guaranteed reduction factor of 2 per iteration:

\[ |Q^{(n)} - Q^*| \le \frac{Q_{\text{max}}}{2^n} \]

After 300 iterations, the relative precision is:

\[ \frac{Q_{\text{max}}}{2^{300}} \approx 10^{-90} \]

In practice, convergence to 0.1% typically occurs within 20–40 iterations.

PH Flash Integration¶

At each bisection iteration, pressure-enthalpy (PH) flashes are performed using DWSIM's built-in thermodynamic engine:

\[ T = f_{\text{PH}}(P, H) \]

This is the key enabler for handling phase change: the PH flash automatically determines the equilibrium temperature and phase fractions for any enthalpy value, including the two-phase region.

Flash Types Used

Flash	When Used	Purpose
PT (Pressure-Temperature)	Known outlet T	Calculate outlet enthalpy
PH (Pressure-Enthalpy)	Known outlet H	Calculate outlet temperature

The PropertyPackage (e.g., Raoult's Law, Peng-Robinson, NRTL) handles the equilibrium calculation. The MSHX solver is property-package agnostic — it works with any package supported by DWSIM.

Heat Distribution for N-Stream Mode¶

For 3+ stream configurations, the total heat duty must be distributed among individual streams. The MSHX uses proportional capacity distribution:

\[ Q_i = Q_{\text{total}} \cdot \frac{Q_{i,\text{max}}}{\sum_j Q_{j,\text{max}}} \]

where \(Q_{i,\text{max}}\) is the maximum heat that stream \(i\) could exchange if brought to the extreme temperature of the opposite side.

This approach ensures:

Each stream receives a share proportional to its thermal capacity
No stream exceeds its maximum possible heat exchange
The distribution remains physically consistent at every bisection step

Interpolation¶

Linear interpolation is used to evaluate composite curve temperatures at arbitrary heat duty fractions:

\[ T(Q) = T_k + \frac{Q - Q_k}{Q_{k+1} - Q_k} \cdot (T_{k+1} - T_k) \]

Phase Change Regions¶

During phase change, the T-Q curve has a flat region where \(T_{k+1} \approx T_k\) but \(Q_{k+1} \neq Q_k\). The interpolation handles this correctly:

For T→Q interpolation: returns the midpoint Q when \(\Delta T < 0.0001\) K
For Q→T interpolation: standard linear interpolation (temperature is constant)

Numerical Safeguards¶

Safeguard	Implementation
Temperature crossover check	If \(\Delta T \le 0.001\) K at any LMTD point, reduce \(Q_{\text{hi}}\)
NaN/Infinity protection	Skip zone if LMTD is NaN, Infinity, or ≤ 0
Maximum iterations	Hard limit of 300 iterations prevents infinite loops
Q bounds collapse	Exit if \(\\|Q_{\text{hi}} - Q_{\text{lo}}\\| / Q_{\text{hi}} < 10^{-5}\)
Flash failure fallback	Linear T interpolation if PH flash fails at a segment point
Energy balance check	Post-solve validation: net imbalance must be < 1% of max stream Q