Simulation of Condensation Flow and Heat Transfer in Spiral Tube Heat Exchangers for Non-Azeotropic Hydrocarbon Mixtures

Summary

Here, the article presents a protocol for numerically simulating condensation heat transfer and flow characteristics of non-azeotropic hydrocarbon mixtures in spiral-wound heat exchangers. This method evaluates operational and rolling conditions to predict heat transfer coefficients and pressure drops.

Abstract

As a core component of the natural gas liquefaction process, spiral wound heat exchangers play a crucial role in LNG production. To comprehensively understand the condensation flow and heat transfer characteristics of non-azeotropic hydrocarbon mixtures inside spiral tubes. The high-precision numerical investigation in this study is guided by a well-defined workflow, which involves geometry creation via a professional modeling tool, mesh generation via a dedicated meshing software, numerical computation within a commercial solver (incorporating real-time convergence monitoring), and quantitative post-processing. This integrated approach ensures the high fidelity of the resulting numerical model. The maximum deviations from classical experimental data(Neeraas's experimental data) remain below 15% for the heat transfer coefficient and below 10% for the frictional pressure drop gradient. The simulation results reveal that varying the rolling periods and amplitudes yields similar oscillatory trends in the heat transfer process, exhibiting both enhancement and degradation effects. Specifically, the rolling period alters heat transfer performance by ±20%, whereas the rolling amplitude affects it by ±10%.

Introduction

Natural gas, as a relatively clean fossil fuel, emits significantly less carbon dioxide and other pollutants upon combustion than coal and petroleum. In the global transition toward renewable energy systems, natural gas is often regarded as a "bridge fuel" due to its ability to maintain the stability and reliability of the energy supply¹. Through liquefaction, gaseous natural gas is cooled into a cryogenic liquid (LNG), reducing its volume by a factor of approximately 600, which greatly facilitates transportation and storage². The spiral-wound heat exchanger (SWHE) is a central component in the natural gas liquefaction process. This type of heat exchanger consists of a series of spiral tubes fixed within a cylindrical shell, wound layer by layer in opposite directions around a central mandrel, with spacers separating the layers to ensure sufficient clearance for heat exchange. Due to its spiral configuration, the SWHE provides a large heat transfer surface area within a small footprint³. This compact design makes it highly suitable for integration into large facilities, particularly offshore floating production platforms where space is strictly limited. In the mixed-refrigerant liquefaction processes commonly used in LNG production, non-azeotropic hydrocarbons flow upward inside the tubes, while the shell-side cold fluid flows downward in a countercurrent manner through the gaps within the tube bundle. Under these conditions, the core process on the tube side is the condensation of non-azeotropic hydrocarbons within the spiral tubes, involving complex gas-liquid two-phase flow^4,5.

To accurately predict the flow and heat transfer characteristics of condensation inside the tubes, extensive research has been conducted. For single-component alkanes, Fries et al.⁶ measured the condensation heat transfer characteristics of propane in horizontal tubes, finding that the pressure drop increased as the tube diameter and saturation pressure decreased. They also noted that gravity caused the heat transfer coefficient at the bottom of the tube to be lower than at the top. Zhuang et al.^7,8 studied the condensation of methane and ethane in horizontal tubes, showing that the heat transfer coefficient and frictional pressure drop increased with flow rate and vapor quality. A previous study⁹ analyzed the condensation process of propane in microchannels, confirming that the heat transfer and pressure drop trends were similar to those in conventional channels. For mixed refrigerants, Smit et al.¹⁰ investigated the condensation of R22/R142b mixtures in horizontal tubes, finding that at low mass fluxes, increasing the R142b mass fraction significantly reduced the heat transfer coefficient. Berrada et al.¹¹ studied an R134a/R23 mixture and found that the temperature glide had little effect on heat transfer under different component ratios. Neeraas conducted experiments on ethane/propane mixtures in spiral tubes, noting that the mixing effect significantly impacted the calculation of the condensation heat transfer coefficient¹². In numerical simulations, Li et al.¹³ simulated the ethane/propane condensation process, showing that the heat transfer coefficient and frictional pressure drop decreased with increasing saturation pressure. Qiu et al.¹⁴ introduced the vapor-liquid entrainment effect when simulating propane condensation in spiral tubes; their results showed that considering this effect reduced the deviation between simulation results and experimental data to less than 25%.

Despite extensive research on various working fluids and channel configurations, a substantial gap persists between the pure or binary fluids typically studied in the literature and the multi-component mixtures utilized in industrial LNG production. Specifically, numerical investigations targeting non-azeotropic hydrocarbon mixtures with three or more components within complex flow channels remain remarkably limited¹⁵. Additionally, regarding the unique application of offshore LNG platforms, there is still a lack of a comprehensive understanding of how equipment motion caused by the marine environment alters condensation flow and heat transfer behavior. To address these research gaps, this study combines computational fluid dynamics (CFD) simulations with existing experimental data to develop a detailed three-dimensional, two-phase condensation flow model. Compare and analyze the simulated heat transfer coefficient and frictional pressure drop based on classic experimental data. Based on this model, the paper focuses on simulating the condensation process of representative gas-field compositions in spiral tubes, particularly under rolling conditions, to investigate the underlying mechanisms of complex motion on multi-component condensation heat transfer. This study provides a reliable theoretical basis and engineering guidance for the design and optimization of efficient heat exchangers in offshore natural gas liquefaction processes.

Protocol

Since this work focuses on the local heat transfer and pressure drop characteristics during condensation inside the helical tube, a reduced domain can be employed once the flow is fully developed, enabling an accurate representation of local flow and thermal behavior. For validation against experimental data, a three-section helical tube model is constructed based on the physical model proposed by Neeraas¹², with a tube diameter of 14 mm, a helix angle of 10°, and a coil diameter of 2 m. The model consists of three regions: a fully developed section (0.6 m), a test section (0.2 m), and a pressure stabilization section (0.2 m). The fully developed section ensures that the flow is sufficiently developed before entering the region of interest. The test section is used for comparison with experimental data and for detailed analysis of local flow and heat transfer characteristics. The pressure stabilization section is designed to maintain outlet pressure stability and to prevent backflow, thereby avoiding interference with the results obtained in the test section. The specific modeling software is sourced from the Table of Materials.

1. Physical model and mesh

1. Open modeling software. In the bottom status bar, select Sketch Mode, and click the Z–X plane to enter the sketch environment.
2. In the top toolbar, select the Circle tool. Draw a circle at the origin with a diameter of 14, then press Enter. Click Return to 3D Mode in the top toolbar. The sketch circle will then be converted into a surface.
3. Select the generated circular surface. Click the Move tool (shortcut: M) in the top toolbar. A triad manipulator (three-axis handle) will appear on the surface. Drag the yellow sphere at the center of the manipulator to the global origin (0, 0, 0), which will serve as the reference for rotation and translation.
4. Drag the red arrow along the X-axis, enter 1000 mm, and press Enter. Rotation: Click the rotation ring around the X-axis (blue or green arc), enter 10°, and press Enter.
5. Click the Pull tool (shortcut: P) in the top toolbar and select the circular surface. In the left panel, choose the Revolve option. Then select the Z-axis of the global coordinate system as the axis of rotation. Enable the Helix option in the left panel. Create Volume 1: In the input box or left panel, enter a height of 138.87 mm and an angle of 45.16°, then press Enter. The first fluid domain is generated.
6. Create Volume 2: Select the new end face of Volume 1. Repeat the helical pull operation, again using the Z-axis as the rotation axis. Enter a height of 69.44 and an angle of 22.58°.
7. Create Volume 3: Select the new end face of Volume 2. Create Volume 3 with the same parameters as Volume 2 using the same method.
8. Click the Workbench tab in the top menu bar. Click the Share button. The software will automatically highlight the two faces that intersect among the three volumes. Click the Complete (checkmark) button on the right.
9. Click the Groups tab in the left panel. Select the initial circular face of the first volume, then click Create Named Selection and define it as the inlet (in).
10. Select the final circular face of the third volume, press Ctrl + G to create a group, and define it as the outlet (out).
11. Select the outer cylindrical surfaces of the three volumes and define them as wall boundaries: wall1, wall2, and wall3.
12. In the structure tree on the left, hold Ctrl and select the three solids. Press Ctrl + G to create a group and rename it as fluid.
13. Connect the generated geometry to the Mesh module and double-click to open the Meshing software. In the left tree, click Mesh. In the Details panel at the bottom left, expand Sizing and set the Element Size to 3.
14. Right-click Mesh in the tree, Insert, Sizing. Select the inlet surface (in) as the geometry and click Apply. Set the Element Size to 0.6.
15. Right-click Mesh, Insert, Inflation. Geometry: Select all three fluid domains and click Apply. Boundary: Select the outer wall surfaces defined as wall, then click Apply. Change the option to First Layer Thickness. First Layer Height: 0.01mm. Maximum Layers: 15. Growth Rate: 1.25.
16. Right-click Mesh, Insert, Method. Select the three fluid domains and click Apply. In the Method dropdown menu, choose Sweep. Under this Selection, select Manual Source. Choose the inlet surface (in) as the source face and click Apply.
17. Right-click Mesh in the tree and select Generate Mesh. In this study, the mesh quality was strictly controlled. The minimum Orthogonal Quality of the generated mesh is above 0.90.

2. Simulation software operation

1. Open the solving software. Go to the File tab, and under Read, select Mesh. Then go to Scale Mesh and set Mesh Was Created In to mm.
2. In the Solver settings, select the Pressure-Based solver, choose Absolute for velocity formulation, and enable the Transient option for time.
   NOTE: By superimposing the oscillation equation onto a stationary reference case and implementing it via a user-defined function, the moving coordinate framework can represent the oscillating condition.
3. Click User-Defined, then select Functions. In the Interpreted UDFs section, load the compiled oscillation file.
   NOTE: The resulting motion is expressed as shown in Equation (1). A static-mesh methodology was adopted and implemented using a moving coordinate framework. The fundamental physics of sloshing relies on the relative motion of the fluid with respect to the container boundary. The sloshing excitation is represented as equivalent dynamic acceleration source terms in the momentum equations, enabling a complete reproduction of the dynamic fluid forces on a stationary mesh.
         (1)
   In the equation, Tc represents the swaying period, and X denotes the displacement generated by swaying.
4. Set the Gravitational acceleration in the Y-direction to −9.81 m/s². Under Models, enable Energy and turn on the Energy Equation.
5. Under Models, enable Viscous and select the Reynolds Stress Model (7 equations). In the Reynolds stress model settings, choose Linear Pressure-Strain. For Near-Wall Treatment, select Scalable Wall Functions.
6. In Phases, set Phase-1 (Primary Phase) as gas and Phase-2 (Secondary Phase) as liquid. Under Global Options, enable Surface Tension Force Modeling, and select the Continuum Surface Force model.
   NOTE: An equivalent pseudo-fluid approach based on temperature- and pressure-dependent thermophysical properties was adopted, which is a widely accepted methodology in CFD studies of multi-component mixtures. Given a fixed initial mixture composition, state-dependent thermophysical properties—including density, dynamic viscosity, thermal conductivity, specific heat capacity, and saturation characteristics—were calculated and generated using the NIST REFPROP database across the entire operating temperature and pressure ranges. In the present study, the mixture maintains a homogeneous macro-composition throughout the simulation. Utilizing NIST-derived variable properties accurately captures the nonlinear thermophysical characteristics of the multi-component fluid while avoiding unnecessary computational overhead.
7. Taking an ethane–propane mixture as an example, at a vapor quality of 0.56 and a pressure of 3.2 MPa, define the liquid phase properties in Materials as follows:
   1. Density: 393.06 kg/m³
   2. Specific heat capacity (Cp): 3866.4 J/(kg·K)
   3. Thermal conductivity: 0.078798 W/(m·K)
   4. Viscosity: 5.4796 × 10⁻⁵ Pa·s
   5. Molecular weight: 37.115 kg/kmol
   6. Standard state enthalpy: 0
   7. Reference temperature: 321 K
8. In Materials, define the gas phase properties as follows:
   1. Density: 67.49 kg/m³
   2. Specific heat capacity (Cp): 3488.7 J/(kg·K)
   3. Thermal conductivity: 0.03035 W/(m·K)
   4. Viscosity: 1.129 × 10⁻⁵ Pa·s
   5. Molecular weight: 34.756 kg/kmol
   6. Standard state enthalpy: 0
   7. Reference temperature: 321 K
9. Set the inlet boundary condition as Mass-Flow Inlet(mass flux 300 kg/(m²·s)), the outlet as Pressure Outlet(0 Mpa), and the wall boundary condition as Heat Flux(-10340W/m²).
10. Under Methods, select the PISO algorithm for solution methods. For Volume Fraction, choose Geo-Reconstruct.
    NOTE: Although the volume-of-fluid (VOF) method is widely accepted for tracking macroscale free-surface topological evolution in sloshing and thermal phase-change processes, inherent limitations persist in interface-capturing precision and in capturing microscale interfacial fluctuations. The VOF formulation fundamentally depends on discrete cell phase volume fractions. The Geo-Reconstruct scheme employed herein significantly alleviates numerical diffusion; nevertheless, resolving subgrid microdroplets, spray formation, or microinterfacial structures remains strictly limited by local grid refinement. For the macro-scale sloshing dynamics, bulk thermal convection, and phase-change mass transfer laws prioritized in this investigation, the current VOF framework with approximately 1.42 million grid elements achieves an optimal balance between topological accuracy and computational cost.
11. In Monitors, set up monitoring for:
    1. Pressure at the inlet and outlet of the test section.
    2. Temperature at the inlet and outlet.
    3. Wall temperature.
    4. Volume fraction at the inlet and outlet.
       NOTE: The convergence criterion for the energy residual is set to 1 × 10⁻⁸, while those for the remaining parameters are set to 1 × 10⁻⁴. Crucial global variables, including the area-weighted average temperature and the total pressure drop across the test section, were dynamically monitored. The computation was continued until these variables showed no further fluctuations, ensuring that the flow field had reached a fully developed and stable state.
12. Initialize method selection standard initialization, calculate from all zones. After initialization, in the Run Calculation panel, set: Time Step Size: 1 × 10⁻⁴ s, and Number of Time Steps: 1 × 10⁶.

3. Post-processing and data export configuration

1. In the Calculation Activities panel, click Autosave (Every Flow Time) to open the Autosave window. In the Autosave settings, set Save Data File Every [s] to 0.01 and specify Flow Time as the save interval type. For Save Associated Case Files Type, select Only if Modified, and then click OK.
2. Open the Contours window from the Results panel. In the Contours settings, activate the Filled, Node Values, Boundary Values, Global Range, and Auto Range options.
3. Select Phases as the contour type and Volume Fraction as the variable, then specify phase-1 as the target phase. Finally, click Save/Display to visualize the contour distribution.
   NOTE: The heat transfer coefficient is calculated as the wall heat flux divided by the temperature driving force, obtained from the temperature difference between the inlet and outlet of the test section. Under sloshing conditions, the time-averaged heat transfer coefficient is adopted. The pressure drop is determined by monitoring the difference between the inlet and outlet pressures, and the frictional pressure drop gradient is subsequently calculated as the ratio of this pressure drop to the length of the tube segment.
4. Import the obtained data into Excel, such as inlet and outlet temperature and pressure values.
5. Obtain the temperature difference and pressure difference between the inlet and outlet according to the calculation method in section 3.3.

Results

Using the validated numerical model, the actual liquefaction process was simulated to systematically investigate variations in the heat transfer coefficient and frictional pressure drop across different operating parameters, thereby providing a theoretical basis for the design and optimization of heat exchangers. The main conclusions are summarized as follows: for the condensation of a pure fluid, heat transfer is primarily confined to the liquid film adjacent to the tube wall, where the gas-liquid interface temperature ...

Discussion

Crucial for safeguarding simulation reliability, the three-section configuration establishes fully developed flow conditions upstream of the test section and suppresses backflow at the outlet, thereby optimizing the accuracy of predicted outcomes. This viewpoint has also been reflected in previous studies on horizontal pipes¹. During the grid generation process in this study, the first-layer mesh height, the number of boundary layers, and the minimum orthogonal quality requirement are equally crit...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Fluent	ANSYS	2020r1	simulation software
SpaceClaim	ANSYS	2020r1	modeling software