Mass Transfer Model Review and Improvement Recommendations

Executive Summary

This document provides a technical review of NeqSim’s evaporation and dissolution model, identifying specific areas for improvement in accuracy, stability, and usability.

STATUS: Core Maxwell-Stefan, transfer-correlation, configuration, and phase-depletion safeguards are implemented. Items still marked as recommendations below, especially external validation and uncertainty quantification, are not claimed as complete.

Related Documentation:

• MassTransferAPI.md - Complete API reference with method signatures, parameters, and examples
• EvaporationDissolutionTutorial.md - Practical tutorial with worked examples
• InterphaseHeatMassTransfer.md - Theory background
• PipelineLiquidEvaporation.md - Dedicated droplet/film completion-distance model and validation limits

Completion-Distance Audit

The general two-fluid pipeline solver is useful when both phases persist and hydraulic holdup, pressure, and momentum coupling are required. It is not a reliable completion-distance metric to interpret a small calculated liquid holdup as the remaining fraction of injected material. Holdup also changes with density, velocity, and flow regime.

PipelineEvaporationStudy is the recommended path when the engineering question is “how many metres until the injected liquid has evaporated?” It retains the Krishna Maxwell-Stefan boundary but adds explicit injected-inventory tracking, droplet or wall-film area evolution, conservative component updates, separate phase energy balances, adaptive axial steps, and an interpolated mass-based completion criterion. The model reports conservation residuals and warnings.

The automated evidence is equation-level correlation checks, positivity and conservation checks, and short axial droplet/film calculations. Absolute completion-distance validation against spray or wetted-wall experiments remains a project responsibility; see the validation ladder and Daïf et al. target case in the dedicated documentation.

Current Implementation Overview

Architecture

The mass transfer calculation follows this hierarchy:

TwoPhaseFixedStaggeredGridSolver
    └── FlowNode.getFluidBoundary()
            └── KrishnaStandartFilmModel (extends NonEquilibriumFluidBoundary)
                    ├── calcBinaryMassTransferCoefficients()
                    ├── calcMassTransferCoefficients()
                    └── massTransSolve()

Key Components

Component	File	Purpose
TwoPhaseFixedStaggeredGridSolver	flowsolver/.../TwoPhaseFixedStaggeredGridSolver.java	Main solver with initProfiles() for mass/heat transfer
KrishnaStandartFilmModel	fluidboundary/.../KrishnaStandartFilmModel.java	Multi-component diffusion model
InterfacialAreaCalculator	flownode/InterfacialAreaCalculator.java	Flow-pattern specific interfacial area
MassTransferCoefficientCalculator	flownode/MassTransferCoefficientCalculator.java	Flow-pattern specific k_L and k_G
MassTransferMode enum	flowsolver/.../MassTransferMode.java	BIDIRECTIONAL, EVAPORATION_ONLY, DISSOLUTION_ONLY
MassTransferConfig (NEW)	flowsolver/.../MassTransferConfig.java	Configurable parameters for mass transfer

---

Identified Improvement Areas

1. Transfer Rate Limiting Logic ✅ IMPLEMENTED

Current Implementation:

// Lines 456-492 in TwoPhaseFixedStaggeredGridSolver.java
if (massTransferMode == MassTransferMode.BIDIRECTIONAL) {
    transferToLiquid = Math.min(transferToLiquid, 0.9 * Math.max(0.0, availableInGas));
} else {
    transferToLiquid = Math.min(transferToLiquid, 0.5 * Math.max(0.0, availableInGas));
}

Issues:

• Fixed 90%/50% limits are arbitrary and don’t adapt to local conditions
• No consideration of time step or node spacing
• Can cause oscillations in high-transfer scenarios

Implemented Solution:

// Adaptive limiting based on Courant-like condition using MassTransferConfig
MassTransferConfig config = getMassTransferConfig();
double maxFraction = config.getMaxTransferFractionBidirectional();

if (config.isUseAdaptiveLimiting()) {
    double localGasVelocity = Math.max(pipe.getNode(i).getVelocity(0), 0.01);
    double residenceTime = nodeLength / localGasVelocity;
    double adaptiveFactor = Math.min(1.0, residenceTime * 10.0);
    maxFraction = maxFraction * adaptiveFactor;
}

2. Phase Depletion Handling ✅ IMPLEMENTED

Current Implementation:

// Lines 500-510 - Minimum moles protection
if (moles < 1e-20) {
    pipe.getNode(i).getBulkSystem().getPhases()[phase].addMoles(comp, 1e-20 - currentMoles);
}

Issues:

• Fixed threshold 1e-20 may not be appropriate for all systems
• No handling of complete phase disappearance (phase transitions)
• Can cause thermodynamic property calculation failures

Implemented Solution:

// Dynamic minimum based on system total moles
double totalSystemMoles = pipe.getNode(i).getBulkSystem().getTotalNumberOfMoles();
double minMolesThreshold = Math.max(config.getAbsoluteMinMoles(),
    totalSystemMoles * config.getMinMolesFraction());

// Phase depletion handling with configurable allowance
if (totalGasPhase > minMolesThreshold) {
    double depletionLimit = config.getMaxPhaseDepletionPerNode() * availableInGas;
    transferToLiquid = Math.min(transferToLiquid, 
        Math.min(maxFraction * availableInGas, depletionLimit));
} else if (config.isAllowPhaseDisappearance()) {
    transferToLiquid = Math.min(transferToLiquid, availableInGas);
}

3. Interfacial Area Model Enhancement ✅ IMPLEMENTED

Current Implementation: InterfacialAreaCalculator provides flow-pattern specific models but:

• Uses simplified geometric models
• Limited validation against experimental data
• No uncertainty quantification

Recommended Improvements:

a) Add wave-induced area enhancement for stratified wavy flow:

public static double calculateStratifiedWavyArea(double diameter, double liquidHoldup,
        double usg, double usl, double rhoG, double rhoL, double sigma) {
    // Base stratified area
    double aFlat = calculateStratifiedArea(diameter, liquidHoldup);
    
    // Kelvin-Helmholtz instability check
    double criticalVelocity = Math.sqrt(sigma * (rhoL - rhoG) / (rhoG * rhoL * diameter));
    double relativeVelocity = Math.abs(usg / (1 - liquidHoldup) - usl / liquidHoldup);
    
    // Wave enhancement factor (Tzotzi & Andritsos, 2013)
    if (relativeVelocity > criticalVelocity) {
        double waveAmplitude = 0.02 * diameter * (relativeVelocity / criticalVelocity - 1);
        double enhancementFactor = 1 + 2 * Math.PI * waveAmplitude / diameter;
        return aFlat * Math.min(enhancementFactor, 3.0); // Cap at 3x
    }
    return aFlat;
}

b) Add droplet entrainment in annular flow:

public static double calculateAnnularAreaWithEntrainment(double diameter, double liquidHoldup,
        double rhoG, double rhoL, double usg, double sigma) {
    // Film interface area
    double aFilm = calculateAnnularArea(diameter, liquidHoldup);
    
    // Entrainment fraction (Ishii & Mishima, 1989)
    double weG = rhoG * usg * usg * diameter / sigma;
    double entrainmentFraction = Math.tanh(7.25e-7 * Math.pow(weG, 1.25));
    
    // Droplet area contribution
    if (entrainmentFraction > 0.01) {
        double d32Droplet = 0.15 * diameter * Math.pow(sigma / (rhoG * usg * usg * diameter), 0.5);
        double dropletHoldup = liquidHoldup * entrainmentFraction;
        double aDroplet = 6 * dropletHoldup / d32Droplet;
        return aFilm + aDroplet;
    }
    return aFilm;
}

4. Mass Transfer Coefficient Correlations 🔧 Medium Priority

Current Implementation: Uses Sherwood number correlations that may not capture all physics.

Recommended Improvements:

a) Add turbulence effects for stratified flow:

// Enhanced Solbraa correlation with turbulence
public static double calculateStratifiedKLTurbulent(double diameter, double liquidHoldup,
        double usl, double rhoL, double muL, double diffL, double scL, 
        double turbulentIntensity) {
    double hydraulicDiameter = 4 * (Math.PI * diameter * diameter / 4 * liquidHoldup) / 
            (Math.PI * diameter * liquidHoldup + diameter * Math.sin(Math.PI * liquidHoldup));
    double reL = rhoL * usl * hydraulicDiameter / muL;
    
    // Base correlation
    double shBase = 0.023 * Math.pow(reL, 0.8) * Math.pow(scL, 0.33);
    
    // Turbulence enhancement
    double turbEnhancement = 1 + 2.5 * turbulentIntensity * Math.sqrt(reL);
    
    return shBase * turbEnhancement * diffL / hydraulicDiameter;
}

b) Add Marangoni effect for surface-active components:

// Marangoni effect reduces mass transfer for surface-active species
public double applyMarangoniCorrection(double kL_base, double surfaceTensionGradient,
        double diffL, double muL) {
    double ma = surfaceTensionGradient * diffL / (muL * kL_base * kL_base);
    return kL_base / (1 + 0.35 * Math.sqrt(Math.abs(ma)));
}

5. Heat-Mass Transfer Coupling 🔧 Medium Priority

Current Implementation: Heat and mass transfer are solved sequentially but coupling effects are limited.

Issues:

• Evaporative cooling not fully coupled
• Dissolution heating/cooling effects approximate
• Temperature gradients in film not modeled

Recommended Improvement:

// Coupled heat-mass balance
public void solveCoupledTransfer() {
    double tolerance = 1e-6;
    int maxOuter = 20;
    
    for (int outer = 0; outer < maxOuter; outer++) {
        // Store previous values
        double[] prevFlux = Arrays.copyOf(molarFlux, numComponents);
        double prevQInterphase = interphaseHeatFlux;
        
        // Solve mass transfer with current temperature
        massTransSolve();
        
        // Calculate latent heat contribution
        double latentHeatRate = 0.0;
        for (int i = 0; i < numComponents; i++) {
            latentHeatRate += molarFlux[i] * getLatentHeat(i, interfaceTemp);
        }
        
        // Solve heat transfer including latent heat
        heatTransSolveWithLatent(latentHeatRate);
        
        // Check convergence
        double massError = calculateRelativeError(molarFlux, prevFlux);
        double heatError = Math.abs((interphaseHeatFlux - prevQInterphase) / 
                                     (Math.abs(prevQInterphase) + 1e-10));
        
        if (massError < tolerance && heatError < tolerance) {
            break;
        }
    }
}

6. Convergence Diagnostics 📊 Low Priority

Current Implementation: Limited convergence monitoring in massTransSolve().

Recommended Improvement:

public class MassTransferConvergenceMonitor {
    private List<Double> residualHistory = new ArrayList<>();
    private int stallCounter = 0;
    
    public ConvergenceStatus checkConvergence(double residual, int iteration) {
        residualHistory.add(residual);
        
        // Check for convergence
        if (residual < tolerance) {
            return ConvergenceStatus.CONVERGED;
        }
        
        // Check for stalling
        if (residualHistory.size() > 5) {
            double avgRecent = average(residualHistory.subList(
                residualHistory.size() - 5, residualHistory.size()));
            double avgOlder = average(residualHistory.subList(
                Math.max(0, residualHistory.size() - 10), residualHistory.size() - 5));
            
            if (avgRecent > 0.9 * avgOlder) {
                stallCounter++;
                if (stallCounter > 3) {
                    return ConvergenceStatus.STALLED;
                }
            }
        }
        
        // Check for divergence
        if (residual > 10 * residualHistory.get(0)) {
            return ConvergenceStatus.DIVERGING;
        }
        
        return ConvergenceStatus.ITERATING;
    }
}

7. User-Configurable Parameters 📝 Low Priority

Current State: Many numerical parameters are hard-coded.

Recommended Improvement: Add a configuration class:

public class MassTransferConfig {
    // Transfer limits
    private double maxTransferFractionBidirectional = 0.9;
    private double maxTransferFractionDirectional = 0.5;
    
    // Convergence
    private double convergenceTolerance = 1e-4;
    private int maxIterations = 100;
    
    // Stability
    private double minMolesFraction = 1e-15;
    private double maxTemperatureChange = 50.0; // K per node
    
    // Model options
    private boolean includeMarangoniEffect = false;
    private boolean includeEntrainment = false;
    private boolean useAdaptiveLimiting = false;
    
    // Getters and setters...
}

---

Implementation Priority

Priority	Improvement	Impact	Effort
1	Adaptive transfer limiting	High stability	Medium
2	Phase depletion handling	Robustness	Medium
3	Wave-enhanced interfacial area	Accuracy	Low
4	Turbulence in k_L	Accuracy	Low
5	Heat-mass coupling	Physical accuracy	High
6	Convergence diagnostics	Debugging	Low
7	Configuration class	Usability	Low

---

Implementation Status

All improvement areas identified in this document have been implemented:

Files Created/Modified

File	Status	Description
MassTransferConfig.java	NEW	Configuration class with all parameters
TwoPhaseFixedStaggeredGridSolver.java	MODIFIED	Adaptive limiting, phase tracking, diagnostics
InterfacialAreaCalculator.java	ENHANCED	Wave enhancement, entrainment, validation
MassTransferCoefficientCalculator.java	ENHANCED	Turbulence effects, Marangoni correction
MassTransferEnhancedTest.java	NEW	Comprehensive test with literature validation

Factory Methods for Common Scenarios

// For complete evaporation scenarios
MassTransferConfig config = MassTransferConfig.forEvaporation();

// For complete dissolution scenarios  
MassTransferConfig config = MassTransferConfig.forDissolution();

// For three-phase gas-oil-water systems
MassTransferConfig config = MassTransferConfig.forThreePhase();

// For research/high-accuracy applications
MassTransferConfig config = MassTransferConfig.forHighAccuracy();

Key Methods Added for Diagnostics

// Check if gas phase completely dissolved
boolean dissolved = solver.isGasPhaseCompletelyDissolved();

// Check if liquid phase completely evaporated
boolean evaporated = solver.isLiquidPhaseCompletelyEvaporated();

// Get mass transfer summary [totalDissolution, totalEvaporation, net]
double[] summary = solver.getMassTransferSummary();

// Get mass balance error
double error = solver.getMassBalanceError();

// Generate validation report against literature
String report = solver.validateMassTransferAgainstLiterature();

Literature References Implemented

Correlation	Reference	Application
Wave enhancement	Tzotzi & Andritsos (2013)	Stratified wavy interfacial area
Entrainment	Ishii & Mishima (1989)	Annular flow droplets
Turbulence	Lamont & Scott (1970)	kL enhancement
Marangoni	Springer & Pigford (1970)	Surface tension effects

---

Validation Test Cases

Recommended Benchmark Cases:

1. Water evaporation into dry nitrogen
   • T = 20-60°C, P = 1 bar
   • Compare to: Solbraa (2002) experimental data
   • Expected accuracy: ±15%
2. CO₂ dissolution into water
   • T = 25°C, P = 1-50 bar
   • Compare to: Carroll et al. (1991) data
   • Expected accuracy: ±20%
3. Hydrocarbon evaporation (n-hexane into methane)
   • T = 20°C, P = 5 bar
   • Compare to: Standing correlation
   • Expected accuracy: ±25%
4. Complete phase transition
   • Flash evaporation scenario
   • Verify phase disappearance handling

---

References

1. Krishna, R., & Standart, G. L. (1976). Mass and energy transfer in multicomponent systems. Chemical Engineering Communications, 3(4-5), 201-275.

2. Solbraa, E. (2002). Measurement and modelling of absorption of carbon dioxide into methyldiethanolamine solutions at high pressures. PhD thesis, NTNU.

3. Tzotzi, C., & Andritsos, N. (2013). Interfacial shear stress in wavy stratified gas-liquid flow. Chemical Engineering Science, 86, 49-57.

4. Ishii, M., & Mishima, K. (1989). Droplet entrainment correlation in annular two-phase flow. International Journal of Heat and Mass Transfer, 32(10), 1835-1846.

5. Higbie, R. (1935). The rate of absorption of a pure gas into a still liquid during short periods of exposure. Transactions of the American Institute of Chemical Engineers, 31, 365-389.

6. Lamont, J.C., & Scott, D.S. (1970). An eddy cell model of mass transfer into the surface of a turbulent liquid. AIChE Journal, 16(4), 513-519.

7. Springer, T.G., & Pigford, R.L. (1970). Influence of surface turbulence and surfactants on gas transport through liquid interfaces. Industrial & Engineering Chemistry Fundamentals, 9(3), 458-465.