Numerical Methods¶
Overview¶
The PHE solver employs an iterative property evaluation loop to determine the outlet temperatures, heat duty, and pressure drops from the specified plate geometry and inlet conditions. The core challenge is that the fluid transport properties (\(\rho\), \(\mu\), \(c_p\), \(k\)) depend on temperature, but the outlet temperatures are unknown at the start of the calculation. The solver resolves this circular dependency through successive substitution on the average bulk temperature of each stream, using DWSIM's thermodynamic engine for property evaluation at each iteration.
Iterative Property Evaluation¶
Average Temperature Approach¶
Fluid properties for each stream are evaluated at the arithmetic mean of the inlet and outlet bulk temperatures:
At the first iteration, the outlet temperatures are unknown. The solver uses an initial guess to start the process:
This assumes approximately 30% of the maximum temperature span as an initial estimate, which provides a reasonable starting point for most operating conditions.
Successive Substitution¶
At each iteration \(k\):
- Compute average temperatures: \(T_{h,\text{avg}}^{(k)}\) and \(T_{c,\text{avg}}^{(k)}\)
- Evaluate fluid properties at \(T_{h,\text{avg}}^{(k)}\) and \(T_{c,\text{avg}}^{(k)}\) via PH flash
- Compute \(\text{Re}\), \(\text{Pr}\), \(\text{Nu}\), \(h\) for each stream using the Martin (1996) correlation
- Assemble \(U\) and compute \(UA\)
- Apply the e-NTU method to obtain \(T_{h,\text{out}}^{(k+1)}\) and \(T_{c,\text{out}}^{(k+1)}\)
- Check convergence: \(\left|T_{\text{out}}^{(k+1)} - T_{\text{out}}^{(k)}\right| < \text{tol}\)
Algorithm Flowchart¶
The following diagram illustrates the complete solution procedure:
flowchart TD
A[Read Inputs:<br/>Geometry, Inlet T, Flow Rates,<br/>Fouling, Flow Arrangement] --> B[Initial Guess:<br/>Outlet Temperatures]
B --> C[Compute Average<br/>Temperatures]
C --> D[Evaluate Fluid Properties<br/>at T_avg via PH Flash]
D --> E[Compute Re, Pr<br/>for Each Stream]
E --> F[Martin 1996 Correlation:<br/>Compute Nu, h for Each Stream]
F --> G[Assemble Overall U<br/>and Compute UA]
G --> H[e-NTU Method:<br/>Compute ε, Q, T_out]
H --> I{Converged?<br/>ΔT < tolerance}
I -->|No| J[Update Average<br/>Temperatures]
J --> C
I -->|Yes| K[Compute Pressure Drop:<br/>Kumar 1984 Correlation]
K --> L[Apply Heat Leak<br/>50/50 Split]
L --> M[Output Results:<br/>Q, T_out, U, UA, ε,<br/>ΔP_hot, ΔP_cold]
style A fill:#e3f2fd,stroke:#1565c0,color:#000
style M fill:#e8f5e9,stroke:#2e7d32,color:#000
style I fill:#fff3e0,stroke:#e65100,color:#000PH Flash Integration¶
At each iteration, the solver calls DWSIM's built-in thermodynamic engine to evaluate fluid properties via pressure-enthalpy (PH) flash or pressure-temperature (PT) flash calculations:
| Flash Type | When Used | Purpose |
|---|---|---|
| PT flash | At average temperature | Evaluate \(\rho\), \(\mu\), \(c_p\), \(k\) at known T and P |
| PH flash | At computed outlet | Determine outlet temperature from enthalpy balance |
The solver is property-package agnostic — it works with any thermodynamic package supported by DWSIM (Peng-Robinson, NRTL, UNIQUAC, Raoult's Law, etc.). The property package handles the equilibrium calculation and returns the required transport properties.
Property Evaluation
For single-phase applications, the PT flash is the primary operation: given the average temperature and stream pressure, DWSIM returns all transport properties needed for the heat transfer and pressure drop correlations. The PH flash is used for the final outlet state determination where the enthalpy is known from the energy balance.
Heat Leak Implementation¶
The optional heat leak \(Q_{\text{leak}}\) represents thermal energy exchanged between the PHE and its surroundings. The implementation splits the heat leak equally between the two streams:
The adjusted energy balances become:
Sign Convention
A positive \(Q_{\text{leak}}\) means heat flows into the exchanger from the surroundings (net heat gain). A negative value means heat flows out to the surroundings (net heat loss). The 50/50 split distributes the effect symmetrically between the two streams.
Convergence Criteria¶
The solver uses the following convergence test at each iteration:
where \(\varepsilon_T\) is the temperature convergence tolerance.
| Parameter | Default Value | Description |
|---|---|---|
| Temperature tolerance \(\varepsilon_T\) | 0.01 K | Maximum change in outlet temperature between successive iterations |
| Maximum iterations | 100 | Hard limit to prevent infinite loops |
| Energy balance check | 1% | Post-solve validation: net imbalance must be less than 1% of heat duty |
In practice, convergence is typically achieved within 5 to 15 iterations for well-conditioned problems with moderate temperature differences.
Numerical Safeguards¶
The solver includes several protective measures to handle edge cases and prevent numerical failures:
| Safeguard | Implementation | Purpose |
|---|---|---|
| Temperature crossover detection | If \(T_{h,\text{out}} < T_{c,\text{in}}\) or \(T_{c,\text{out}} > T_{h,\text{in}}\), flag infeasible | Prevents thermodynamically impossible results |
| Maximum iteration limit | Hard cap of 100 iterations | Prevents infinite loops on non-converging cases |
| Property bound checking | Verify \(\mu > 0\), \(k > 0\), \(\rho > 0\), \(c_p > 0\) after each flash | Catches flash failures or invalid property returns |
| Reynolds number floor | Enforce \(\text{Re} \ge 0.1\) | Prevents division by zero in friction factor correlation |
| NTU ceiling | Cap NTU at 1000 | Prevents numerical overflow in the exponential terms of the e-NTU formulas |
| Effectiveness bounds | Clamp \(0 \le \varepsilon \le 1\) | Ensures physically meaningful heat transfer |
| Energy balance validation | Post-solve check: \(\|Q_h - Q_c\| / Q_{\max} < 0.01\) | Confirms consistency of the converged solution |
| Damping factor | Under-relaxation \(\alpha = 0.7\) on outlet temperature updates if oscillation detected | Stabilises convergence for sensitive cases |