Coupled Modeling of CO₂ Frosting, Fluid Flow, and Heat Transfer under Cryogenic Conditions Using a Nucleation-Based CFD Framework

1 Introduction

As a stable, clean, and efficient low-carbon energy source, natural gas plays a vital role in the global transition toward green and low-carbon energy systems. In recent years, Pressurized Liquefied Natural Gas (PLNG) technology has garnered widespread attention due to its low-energy consumption advantage [1]. Compared to conventional LNG technology, PLNG liquefies and stores natural gas at higher pressures (1~2 MPa), corresponding to a liquefaction temperature of approximately −120 to −100°C, significantly higher than atmospheric pressure LNG (−162°C). At a liquefaction pressure of 1.7 MPa, the liquefaction temperature is approximately −115°C, requiring only 60% of the power needed by conventional LNG technology [2]. Furthermore, the PLNG process eliminates the need for CO2 and heavy hydrocarbon pretreatment, reducing land requirements and simplifying the liquefaction process flow.

However, due to increased permissible CO2 fraction, CO2 desublimation and frost formation readily occur during cryogenic heat exchange, thereby impairing heat transfer efficiency. To address this, researchers have proposed integrating CO2 cryogenic desublimation separation with the PLNG process. This approach utilizes the cold energy from the liquefaction process to remove CO2 via desublimation. Compared to traditional pretreatment methods, it significantly reduces capital investment and energy consumption [3], simplifies the pretreatment workflow, and achieves efficient energy utilization.

At present, most of the existing research on frosting has focused on the phenomenon of water vapor frosting. Hayashi et al. [4] divided the frosting process into three distinct stages: the growth of ice crystals, the development of the frost layer, and the stage of a fully developed frost layer. Lei et al. [5] experimentally investigated the effect of surface temperature on the frosting behavior of inverted cold plates under natural convection conditions. Zhang et al. [6] experimentally analyzed how icing conditions affect frost distribution and growth characteristics across different fin spacings in finned tube heat exchangers. Concurrently, researchers proposed various phase-change numerical models to simulate the icing process, employing a combined approach of theoretical analysis and numerical simulation to predict water vapor icing phenomena. Lee et al. [7] developed a computational model for frost layer growth that accounted for water vapor permeation effects, and proposed a formula to calculate phase-change mass transfer rates for analyzing the evolution of frost layer porosity and thickness. Byun et al. [8] developed a frost formation model to predict frost aggregation, density stratification, and internal property distribution under cryogenic conditions. Qi et al. [9] proposed a frost formation model under cryogenic conditions to explain the phenomenon of trace water vapor frosting on a plate during nitrogen cooling. Wong et al. [10] proposed a mass-transfer-theory-based CFD model using dimensionless numbers to predict frosting on cold surfaces, which avoids empirical constants tied to specific conditions and controls thickness prediction error within 20%.

Regarding experimental studies on CO2 cryogenic desublimation and frost formation, Song et al. [11,12,13,14] designed a three-stage cooling test bench based on a Stirling refrigerator for separating CO2-containing mixed gases, analyzing the effects of gas flow rate and cold head temperature on CO2 desublimation, system energy consumption, and heat transfer performance; Naletov et al. [15] conducted experiments on N2/CO2 mixtures of varying concentrations, revealing the effects of flow rate, concentration, and temperature on CO2 capture efficiency and product purity; Wang et al. [16] obtained crystal morphologies of CO2 desublimation under different temperatures and supercooling conditions through visualization experiments, these experimental studies elucidate the patterns of CO2 frost formation and typical crystal morphologies, providing intuitive evidence for understanding the frosting mechanism. Cai et al. [17] designed a set of high-pressure and low-temperature visualization experimental equipment. Through conducting synchronous experiments of visualization and heat transfer, they explored the desublimation laws of carbon dioxide under pressure and its impact on heat transfer. Overall, however, research on the frosting mechanism of CO2 remains limited.

In the field of CO2 cryogenic desublimation frosting simulation, researchers have primarily established heat and mass transfer models based on the thermodynamic properties of CO2 mixtures, frosting temperature prediction, and phase equilibrium calculations to analyze desublimation behavior under various working conditions. Jensen et al. [18] calculated CO2 frost points for multicomponent systems using Dalton’s partial pressure law and phase equilibrium theory based on PR or SRK equations, though without experimental comparison. Srikanth et al. [19] proposed a thermodynamic framework based on the Peng-Robinson state equation (PR EoS) and the Lagrangian cubic solver, which is used to accurately predict the CO2 frost point. For process modeling, Naletov et al. [20] developed a quasi-transient mathematical model for CO2 desublimation in thermal power plant flue gases, analyzing the effects of temperature distribution, frost layer thickness, and flow field parameters on capture efficiency. Haddad et al. [21] proposed a two-dimensional transient model incorporating the porous medium properties of frost layers, revealing key drivers of frost thickness and density evolution. Lei et al. [22] employed the lattice Boltzmann method to investigate cryogenic on cylindrical surfaces, elucidating flow and mass transfer characteristics within porous structures. He et al. [23] integrated desublimation phase-change and gas-solid equilibrium models into Computational Fluid Dynamics (CFD) simulations to model CO2 frost layers on cold surfaces. Du et al. [24] conducted experimental studies on the influence of CO2 desublimation on heat transfer performance in a binary gas under high flow rates, analyzing the heat transfer coefficient, frost layer thermal resistance, and desublimation mechanism.

Overall, existing simulations primarily focus on frost layer properties and phase equilibrium calculations, lacking in-depth analysis of frost’s impact on mainstream gas flow and heat transfer.

At present, understanding of CO2 frosting patterns under PLNG working conditions remains limited, particularly lacking in-depth research on the evolution of frost layers on heat transfer surfaces and their impact on natural gas flow and heat transfer. Therefore, this study establishes a two-dimensional CFD numerical model based on classical nucleation theory to investigate the heat and mass transfer mechanisms of CO2 desublimation frosting within cryogenic heat exchange channels. The CO2 frost point under various gas compositions was calculated using the PR state equation combined with MATLAB programming, providing a basis for determining frosting conditions. Simultaneously, by incorporating gas-solid phase transition kinetics and multiphysics transport equations, the evolution of the flow field and frost layer growth characteristics were analyzed, revealing the operational mechanism of CO2 cryogenic desublimation under pressurized liquefaction conditions. The research findings will address gaps in understanding the coupled frost layer-flow-heat transfer mechanism within cryogenic heat exchange channels, providing theoretical support for optimizing cryogenic carbon capture and PLNG processes.

2 Modelling and Methodology

2.1 Physical Modelling

The fundamental mechanism of CO2 frost formation is as follows: when a mixed gas stream flows over a cryogenic wall surface, part of the CO2 directly desublimates on the surface to form frost crystals, while another portion diffuses into the frost layer and deposits, causing the frost layer’s thickness and density to continuously increase. As shown in Fig. 1, the growth process of the frost layer is divided into three main stages.

Figure 1: Diagram of the frost layer growth process.

While actual CO2 frost exhibits 3D morphological features, a 2D channel approximation was adopted for this preliminary study. This simplifies the computational domain to represent the mid-plane of a wide rectangular channel, allowing for high-resolution transient coupling of multiphase heat and mass transfer while minimizing computational costs. This approach is sufficient for capturing the macro-scale growth kinetics and flow blockage effects, though it inherently neglects 3D edge effects and lateral morphological variations.

In practical experimental and industrial cryogenic configurations, such heat exchange channels typically feature a high width-to-height aspect ratio (W/H ≥ 10). Therefore, the 2D simulation effectively represents the mid-plane of a wide rectangular channel, where the lateral wall edge effects on the central core flow are negligible, further justifying this dimensional simplification [25].

To describe this process, a two-dimensional channel model was established, as shown in Fig. 2. The channel height is H = 0.015 m, with an adiabatic inlet section (L1 = 0.2 m) and a constant-temperature section (L2 = 0.4 m). The mixed gas enters from the left inlet, undergoes heat exchange with the cryogenic wall surface, and exits from the right outlet.

Figure 2: Two-dimensional simulation domain.

To investigate the effects of mixed-gas inlet velocity, inlet temperature, CO2 mol fraction, and cryogenic wall surface temperature on CO2 desublimation, numerical simulations were conducted for each of the nine working conditions listed in Table 1.

2.2 Numerical Modelling

Within the heat exchange channel, when the gas cools below the CO2 frost point and the CO2 supersaturation exceeds 1, CO2 undergoes a phase transition from gas to solid. These two phases are treated as independent flows, and their interactions are represented through source terms introduced into the continuity, momentum, energy, species, and volume fraction equations.

To ensure computational accuracy while minimizing computational load and enhancing efficiency, the model employs the following assumptions:

• a.The flow rate of the mixed gas is relatively low, keeping the Mach number well below 0.3, meaning dynamic compressibility is negligible. While a 50 K temperature gradient exists within the domain, the localized density variations are assumed to have a secondary effect compared to the dominant forced convection. Additionally, the mass removed via desublimation is small relative to the total bulk flow. Thus, the mixture is treated as an incompressible Newtonian fluid to optimize computational efficiency, with the understanding that neglecting density variations introduces a minor local deviation.
• b.Deublimation and melting processes within the frost layer are disregarded. This assumption is primarily applicable to the early and intermediate crystal growth phases simulated in this 60-min transient forecast. For long-term modeling extending into the ‘Fully Developed Period’, surface desublimation balances internal deposition, and this assumption would need to be relaxed.
• c.The frost layer is characterized by a porous structure, while the influence of velocity is disregarded.
• d.The frost layer temperature is low, allowing radiation heat transfer to be ignored.

The nucleation and growth of CO2 crystals involve mass transfer of CO2 from the gas phase to the solid phase. Nucleation theory, a classical framework for describing gas condensation nucleation processes, is also employed to calculate mass transfer rates during frost layer formation and growth.

• (1)The nucleation process of CO2 crystals

Gibbs and Thomsion provided the free energy required to form a spherical nucleus with a radius of:

The first term of Eq. (1) is the volume free energy, which describes the difference between the formed spherical crystal nucleus with a radius of and the free energy of the parent phase, which is negative. The second term is surface energy, which is the energy required to maintain a surface with a radius of and a surface energy per unit area, and it is positive.

When ∂ΔG/∂r = 0, ΔG reaches the maximum value, the corresponding r is the critical nucleus radius:

The critical nucleus radius is the dividing point of nucleus stability. When r < rc, the crystal nucleus is unstable and its size keeps decreasing until it disappears. When r > rc, the crystal nucleus is stable and continues to grow.

At this time, the nuclear energy barrier is:

Among them:

In the Classical Nucleation Theory (CNT) framework, the nucleation rate is highly sensitive to the solid-vapor surface tension (σ) and the contact angle (θ). Due to the experimental challenges in measuring cryogenic solid-gas interfaces, the surface tension of CO2 in this model was determined based on established theoretical correlations utilized in recent desublimation studies [26,27].

Furthermore, the heterogeneous nucleation energy barrier is modified by the contact angle between the solid CO2 nucleus and the cryogenic wall. Based on recent molecular dynamics simulations of CO2 wettability on cryogenic metallic surfaces [28], a contact angle of θ = 58° was adopted within the User-Defined Function (UDF). This specific wetting characteristic has been shown to minimize lattice deviation and promote stable CO2 crystal nucleation during the early stages of subcooling and desublimation.

• (2)The growth process of CO2 crystals

The classical homogeneous nucleation rate (the number of critical nuclei formed per unit volume per unit time) is controlled by the collision rate between molecules and critical nuclei, and the calculation formula is:

The growth of crystal nuclei can be regarded as a process of increasing radius. Based on the principle of conservation of energy, the growth rate of crystal nuclei is expressed as:

The formula for the latent heat of CO2 desublimation fitted by Azreg-Ainou based on theoretical methods:

The parameters in the formula are shown in Table 2.

• (3)Governing equation

The Eulerian multiphase flow model is adopted to calculate the interaction between the mixed gas phase and the dry ice phase during the flow process. The model contains two mass conservation equations for each phase, two momentum conservation equations for each phase, two energy conservation equations for each phase, one volume fraction conservation equation for the secondary phase, and one component conservation equation for the secondary components of the primary phase, which are respectively solved.

The mass conservation equations are established separately for the mixed gas phase and the dry ice phase as follows:

In the equation, S m g and S m i are the mass source terms for the mixed gas phase and the dry ice phase, respectively.

The momentum conservation equations are established separately for the mixed gas phase and the dry ice phase as follows:

In the equation, S u g and S u i are the momentum source terms for the mixed gas phase and the dry ice phase, respectively. τ ¯ g and τ ¯ i are the sum of the viscosity stress tensor and the turbulent stress tensor of the mixed gas phase and the dry ice phase. R → i g and R → g i are the momentum transfer between the two phases.

The energy conservation equations are established separately for the mixed gas phase and the dry ice phase as follows:

In the equation, S h g and S h i are the energy source terms for the mixed gas phase and the dry ice phase, respectively. q → g and q → i are the internal heat exchange for the mixed gas phase and the dry ice phase, respectively. Q i g and Q g i are the energy transfer between the two phases.

In the Eulerian multiphase flow model, each phase is characterized by its volume fraction. It is necessary to solve the volume fraction conservation equation for the secondary phase, i.e., the solid CO2 particles. The volume fraction conservation equation for the dry ice phase is as follows:

In the equation, S α i is the volumetric fraction source term resulting from phase change. α g and α i are the volume fractions of the mixed gas phase and the dry ice phase, respectively.

Since the mass transfer from the gas mixture phase to the dry ice phase essentially corresponds to the mass transfer of the CO2 gas component from the gas mixture to the dry ice phase, it is necessary to consider the conservation equation for the CO2 gas component in the gas mixture. The CO2 gas component conservation equation is as follows:

In the equation, S ω , g C O 2 is the source term of the component equation for phase change. J ¯ g C O 2 is the diffusion fluct of CO2. ω g N 2 and ω g C O 2 are the mass fractions of N2 and CO2 in the mixed gas, respectively.

Given the channel height (H = 0.015 m) and low inlet velocities (0.4 to 0.8 m/s), the nominal Reynolds number of the flow remains in the laminar to transitional regime (Re ≈ 500~1000). However, the SST k-ω turbulence model was selected because the growth of the porous frost layer introduces significant local surface roughness and flow disturbances. The calculation formula is as follows:

In the equation, k n is the turbulent kinetic energy of phase n, ω n is the specific dissipation rate, μ t , n is the turbulent viscosity, μ n is the molecular dynamic viscosity, G k , n and G ω , n are the production terms of k and ω induced by velocity gradients, respectively, G b , n and G ω b , n are the buoyancy production terms, F 2 is the blending function, F n is the magnitude of the strain rate or vorticity, and σ k , n and σ ω , n are the corresponding turbulent Prandtl numbers.

When a mixed gas phase comes into contact with dry ice, momentum exchange occurs. The momentum transfer between the two phases is:

The interphase momentum transfer coefficient is:

The drag coefficient of solid particles is calculated using the Schiller and Naumann model:

The Reynolds number of solid particles is:

When a mixed gas phase comes into contact with dry ice, energy exchange occurs. The energy transfer between the two phases is:

The interphase heat transfer coefficient is:

The Nussel number is calculated using the Ranz-Marshall model:

The Prandtl number of the mixed gas is:

• (4)Source item

According to the mass transfer model of the classical nucleation theory, the mass transfer between the gas phase and the solid phase during the nucleation and growth process of CO2 crystals can be expressed as:

In the above formula, the left phase represents the mass transfer caused by newly formed solid CO2 particles, and the right phase represents the mass transfer caused by the growth of solid CO2 particles.

This study uses custom functions for programming and incorporates the entire phase transition process into the control equation through source terms, thereby reflecting the phase transition phenomenon of the frosting process in the calculation. When the temperature is lower than the triple point temperature, the source terms in each control equation can be calculated by the following methods; Otherwise, all source items are set to 0.

The source terms for each control equation are shown in Table 3. m ˙ g → i and m ˙ i → g represent the mass transfer rates between the mixed-phase and dry-ice phases, respectively. h g is the specific enthalpy of the mixed-gas phase. L s u b is the latent heat of sublimation of CO2. The interpretation of the symbols related to source terms in the other governing equations is provided in the previous section.

2.3 Solution Programme

The simulation was conducted using the CFD software Fluent 2021 R2. The Eulerian multiphase flow model and the SST-k-ω turbulence model were selected, and the energy equation and component transport model were opened. The setting of boundary conditions and solution methods is shown in Table 4 and Table 5. For the cryogenic wall surface, the mass transfer of the species is driven volumetrically by the source term in the computational cells adjacent to the wall, rather than applying a fixed saturation concentration boundary condition. A custom function UDF written based on classical nucleation theory is loaded for transient calculation. The time step is set to 0.1 s, and the convergence criterion for all variables is selected to be less than 10−6.

As shown in Fig. 3, these disturbances locally elevate the turbulent kinetic energy and trigger transition near the phase interface. The SST model effectively captures this artificially enhanced diffusivity and heat transfer induced by the porous boundary, which a standard laminar model would underpredict. To prevent the artificial overprediction of turbulent viscosity in the undisturbed mainstream flow, the Low-Reynolds number (Low-Re) correction option within the SST k-ω model was activated during the setup. This specific damping function ensures that the model transitions accurately between the laminar bulk flow and the locally turbulent regions induced by the rough, porous frost surface.

Figure 3: Mixed gas CO₂ mol fraction and turbulent kinetic energy with surface position (t = 30 min).

2.4 Independence Validation

The simulation domain is meshed using a structured grid partitioning method, primarily considering the boundary layers on the bottom surface and the frost layer surface. To ensure grid quality and simulation accuracy, a combination of coarse and fine meshes is employed. Since the frost layer continuously grows over time on the cryogenic bottom surface, the boundary layer in the two-dimensional model is meshed with increased density, as shown in Fig. 4.

Figure 4: Schematic diagram of grid structure.

By adjusting the overall cell size and local encryption degree of the grid, five types of grids, namely 6940, 11,140, 13,240, 17,440, and 21,640, were divided. The grid quality was all above 0.9. Meanwhile, five time steps, namely 0.05 s, 0.1 s, 0.15 s, 0.2 s, and 0.25 s, were selected. Under working condition 1, the independence of the grid and time step was verified by detecting the volume fraction and mass transfer rate of the dry ice phase. The results are shown in Fig. 5. When the number of grids reaches 17,440, the results of both no longer change with the increase of the number of grids. When the time step is less than 0.1 s, the differences in the results of the two under different time steps are very small. Therefore, considering the calculation accuracy and efficiency comprehensively, the number of grids of 17,440 and the time step size of 0.1 s were selected for the simulation calculation. This fine temporal resolution is necessary not only for numerical stability but also to physically capture the highly transient nonlinear kinetics of CO2 initial nucleation and the subsequent rapid micro-structural evolution of the frost layer [29,30].

Figure 5: Independence validation (a) Grid (b) Time step.

2.5 Model Validation

The average frost layer thickness and thermal conductivity were selected as comparison metrics. By contrasting with numerical models and experimental results from previous literature, the accuracy of the numerical model constructed in this study was verified. The same experimental conditions as Song et al. [14] (cryogenic wall surface temperature 153.15 K, mixed gas inlet temperature 243.15 K, mixed gas inlet velocity 0.0328 m/s, CO2 mol fraction 13%) were employed to demonstrate the frost layer growth process.

Fig. 6 shows the comparison between Song’s experimental results and the simulation results of this model. The results indicate that the relative error between this model and the experimental results is less than 15%, which can accurately describe the CO2 desublimation and frosting process.

Figure 6: Comparison between simulated values and experimental values (a) Average frost layer thickness (b) Average frost layer thermal conductivity.

3 Results and Discussion

3.1 Calculation of CO₂ Frost Point in the Mixed Gas

In this study, gas-solid phase equilibrium theory is employed to construct a gas-solid escape equation for calculating the CO2 frost point in the CH4-CO2 binary system.

The gas-solid escape rate equation is expressed as follows:

The gas-phase fugacity and solid-phase fugacity of CO2 can be calculated respectively by the following formulas:

By combining Formulas (33) and (34), the core equation for calculating the CO2 frost point is established:

The Antoine equation, which divides the temperature into three intervals to calculate the saturated vapor pressure of CO2 desublimation, respectively:

The parameters in the formula are shown in Table 6.

The standard form of the PR state equation is expressed as:

The cubic form expressed by the compression factor Z is:

Among them,

The classical van der Waals mixing rule is employed to calculate the parameter a in Eq. (39) and the parameter b in Eq. (40):

Among them,

Combining the PR state equation with the vdW mixing rule, the fugacity coefficient of component i in the mixed gas can be expressed as:

The fugacity coefficient of CO2 gas in the saturated state does not require a mixing rule, and the above formula can be simplified as:

The CO2 frost point of the CH4-CO2 binary system under the corresponding working conditions was calculated by the PR equation of state method, and the relative deviation and absolute deviation were calculated. They were compared and analyzed with the experimental values, and the results are shown in Table 7 and Fig. 7.

Figure 7: Analysis of frost point in CH₄-CO₂ binary system (a) Deviation distribution (b) Summary of Deviation Data.

Within the pressure range of 689.29 to 2757.12 kPa, for CO2 mol fraction of 2%, 4%, 8%, and 10%, the results calculated using the PR equation of state are consistently lower than the experimental values, with average relative deviations below 3%, indicating a relatively high accuracy. Therefore, in the numerical simulation process of this paper, the PR state equation method is adopted to calculate the CO2 frost point and load it into the user-defined function.

Natural gas in actual PLNG systems contains trace amounts of heavier hydrocarbons. These trace components subtly alter the phase envelope curve and may act as heterogeneous nucleation sites, causing a slight deviation in the actual CO2 crystallization point compared to the idealized binary model.

3.2 Flow Field Characteristics

The simulation is conducted with a cryogenic wall surface temperature of 150 K, a mixed-gas inlet temperature of 200 K, an inlet velocity of 0.4 m/s, and a CO2 mol fraction of 10%. The results capture the evolution of the temperature field, velocity field, and related parameters within the heat exchange channel during the frosting process, thereby revealing the flow-field characteristics in the vicinity of the cryogenic wall surface under CO2 frosting conditions.

3.2.1 Temperature Distribution

Fig. 8 shows the temperature distribution inside the heat exchange channel at different times. In the initial stage, the cooling effect of the cryogenic wall surface on the gas is confined to the near-wall region, forming a thin low-temperature boundary layer, while the main flow region remains at a relatively high temperature (approximately 195 K). As time progresses, by t = 30 min, the cooling gradually propagates toward the channel core, and the low-temperature zone near the cryogenic wall surface thickens significantly, developing a relatively stable vertical temperature stratification. At this stage, the gas temperature in the lower region decreases to 155~170 K. In addition, the average gas temperature decreases progressively along the flow direction toward the channel outlet. This is mainly because the incoming mixed gas initially remains warm and exchanges heat insufficiently with the cryogenic wall surface, whereas continuous flow enhances convective heat transfer, allowing more effective cooling downstream.

Figure 8: Contours of temperature distribution.

As frosting proceeds, the process becomes relatively stable at t = 50 min. The temperature field indicates a slight expansion of the high-temperature region in the gas space away from the cryogenic wall surface near the outlet. For instance, the gas temperature adjacent to the upper surface at 0.58 m from the inlet rises from about 185 K at t = 30 min to approximately 190 K at t = 50 min. This slight increase results from the growing frost layer, which increases thermal resistance and weakens the cryogenic wall surface’s cooling capability in the later stage of frosting. Consequently, although the average gas temperature near the outlet remains lower and the frost layer thicker, the reduced heat transfer efficiency allows warmer gas to accumulate near the outlet region.

Overall, the evolution of the temperature distribution reflects the variation in local cooling intensity along the flow path during steady gas flow. It captures the dynamic interaction between the cryogenic wall surface and the mixed gas, revealing the time-dependent evolution of the heat transfer process and indirectly illustrating the spatial distribution, growth rate, and staged characteristics of the frost layer.

3.2.2 Velocity Distribution

As shown in Fig. 9, the velocity distribution cloud map illustrates the temporal evolution of the mixed gas flow during different frosting stages. At t = 10 min, the maximum velocity in the channel reached 0.65 m/s, with a relatively uniform distribution and a concentrated mainstream region. The velocity near the cryogenic wall surface was slightly lower due to viscous effects. At this stage, the frost layer remained thin and exerted little disturbance on the gas flow. The effective cross-sectional area of the channel was largely preserved, and the flow resistance was not yet significant.

Figure 9: Contours of velocity distribution.

At t = 30 min, the velocity in the middle and upper regions of the channel increased, and the location of the maximum velocity gradually shifted toward the outlet region. This occurred because the frost layer continuously developed along the cryogenic wall surface, progressively reducing the available flow space. As a result, the mixed gas was diverted to flow above the growing frost layer. With the ongoing frosting process, the effective flow area decreased, causing the gas velocity to increase correspondingly. Moreover, closer to the channel outlet, heat exchange between the gas and the cryogenic wall surface became more sufficient, leading to higher CO2 supercooling and supersaturation, which promoted thicker frost formation. Consequently, regions with thicker frost layers exhibited higher gas velocities above them, and the high-velocity zones expanded both vertically and laterally.

At t = 50 min, the maximum velocity increased to 0.8 m/s, and the high-velocity region became more continuous and concentrated. This indicates that the frost layer imposed a stronger blocking effect on the flow, forcing the mixed gas through narrower passages and markedly accelerating the flow compared with the initial stage.

Overall, the evolution of the velocity field reveals the feedback effect of frost growth on the flow dynamics: as the frost layer thickens, the effective flow channel narrows, the local gas velocity intensifies, and the overall flow resistance increases, thereby reshaping the velocity distribution within the heat exchange channel.

3.3 Frost Growth Characteristics

In this study, the characteristic parameters of frost layer growth-namely thickness, density, and thermal resistance-were dynamically evaluated using an unsteady-state model. The effects of cryogenic wall surface temperature, mixed-gas inlet velocity, inlet temperature, and CO2 mol fraction on the evolution of these parameters were systematically analyzed.

3.3.1 Frost Layer Thickness

The frost layer thickness serves as an indicator of the extent of frost formation and enables quantitative analysis of its growth behavior. The average frost layer thickness is defined as:

In the formula, Vf represents the total volume of the frost layer; αi is the volume fraction of the dry ice phase; Vj represents the volume of the computing domain unit.

Based on the classical nucleation theory, under isothermal conditions, molecular aggregation must overcome a high energy barrier at low supersaturation, resulting in a low nucleation rate. In contrast, at high supersaturation, the phase-change driving force becomes sufficient for the spontaneous formation of stable nuclei, leading to a significantly increased nucleation rate that eventually stabilizes after reaching its maximum value. As shown in Fig. 10, a broad region of high nucleation rate appears along the entire cryogenic wall surface, indicating that during the initial cooling stage, the gas near the cryogenic wall surface rapidly attains a high degree of supersaturation, thereby triggering continuous and intense nucleation activity.

Figure 10: Contour of nucleation rate.

At the initial moment, the frost layer thickness is zero. As the cryogenic wall surface temperature decreases to a lower level, the temperature of the mixed gas inside the channel gradually drops under the driving force of the temperature difference, with the gas closer to the cryogenic wall surface exhibiting lower temperatures. Consequently, CO2 in the mixed gas progressively enters a supercooled state. With the rapid increase in subcooling, the supersaturation of CO2 near the cryogenic wall surface becomes significantly greater than 1, as shown in Fig. 11. The spatial distribution of supersaturation closely corresponds to that of the nucleation rate, during the initial stage, the entire cryogenic wall surface region exhibits a high level of supersaturation, corresponding to a broad and continuous zone of high nucleation activity.

Subsequently, CO2 rapidly nucleates and grows on the subcooled surface, leading to the formation of a frost layer. As the gas continues to flow, CO2 within the gas mixture continuously nucleates and grows, resulting in a gradual increase in frost layer thickness. With increasing time, the growth rate of the frost layer gradually decreases until the thickness reaches a steady value. This behavior can be explained as follows: during the initial stage of CO2 desublimation, the mixed gas directly contacts the cryogenic wall surface, where heat exchange is highly efficient. The near-wall region exhibits high CO2 supersaturation, promoting rapid nucleation and crystal growth, thereby causing the frost layer to thicken quickly. However, as the frost layer develops, it forms a thermal barrier that hinders heat transfer between the cryogenic wall surface and the gas. Meanwhile, the gas located farther from the cryogenic wall surface experiences weaker subcooling and lower supersaturation, which reduces the local nucleation and growth rates of CO2, thus slowing the overall frost growth.

Figure 11: Contour of supersaturation.

In addition, as the frost layer thickens, the effective flow passage within the heat exchange channel becomes narrower, increasing the gas velocity above the frost layer and shortening its residence time. This shortens the duration for heat and mass transfer between the gas and the cryogenic wall surface, further limiting frost growth. As frosting continues, the thicker the frost layer becomes, the slower its subsequent growth rate. Eventually, the frost layer thickness approaches a quasi-steady state. At this stage, the frost layer is relatively thick, the thermal resistance is high, and the gas velocity is elevated, making further CO2 nucleation and growth on the frost surface negligible. Consequently, the average frost layer thickness remains essentially unchanged.

As shown in Fig. 12, a lower cryogenic wall surface temperature, reduced inlet velocity of the mixed gas, lower inlet temperature, and higher CO2 mol fraction all contribute to an increase in the average frost layer thickness. While the frost layer thickness approaches a quasi-steady state for the 6% and 8% CO2 cases within the 60-min timeframe, the scenario with the highest CO2 fraction (10%) still exhibits a slight upward trajectory. This indicates that the strict ‘Fully Developed Period’-where surface desublimation perfectly balances continuous deposition-requires extended simulation times to be fully attained under high-concentration conditions. Consequently, when comparing the final structural values at exactly t = 60 min, readers should be mindful of this dynamic discrepancy between the transient state (10%) and the quasi-steady states (6% and 8%). The observed disparities in frost thickness and density among these concentrations would likely amplify if the simulation were extended until all cases reached their absolute fully developed periods.

Figure 12: Variation of the average frost layer thickness under the influence of different factors (a) Cryogenic wall surface temperature (b) Mixed gas inlet velocity (c) Mixed gas inlet temperature (d) Mixed gas CO₂ mol fraction.

Specifically, lowering the cryogenic wall surface temperature markedly enhances the heat exchange between the frost layer and the surface, leading to a rapid decrease in the gas temperature within the near-wall region. As a result, both the degree of subcooling and the degree of supersaturation increase significantly, thereby accelerating the desublimation of CO2 and promoting frost layer growth. Reduced cryogenic wall surface temperatures not only strengthen the heat transfer driving force, enabling the gas to reach desublimation conditions more quickly, but also lead to the formation of a denser frost structure, further increasing the overall layer thickness.

Decreasing the inlet velocity of the mixed gas weakens the fluid shear effect on the frost surface, thereby reducing the peeling or erosion of the nascent frost layer and facilitating its continuous growth under stable conditions. A lower velocity also extends the residence time of the gas near the cryogenic wall surface, enhancing heat and mass exchange with the frost layer and providing greater opportunity for CO2 to desublimate at the surface. This mechanism allows the frost layer to achieve greater thickness.

Lowering the inlet temperature of the mixed gas reduces its initial enthalpy, shortens the time required to cool to the CO2 frost point, and enables the frost surface temperature to be maintained at a relatively low level. This condition substantially increases both the CO2 concentration gradient and the degree of supersaturation at the gas-solid interface, thereby accelerating the phase-change mass transfer rate and promoting rapid frost growth. Furthermore, lower inlet temperatures favor the development of frost layers with more stable microstructures and reduced porosity.

Increasing the CO2 mol fraction in the mixed gas elevates the amount of condensable component per unit volume, directly raising the vapor pressure of CO2 and enhancing the mass transfer driving force at both the surface and interior of the frost layer. Elevated CO2 concentrations not only intensify local supersaturation and the driving force for desublimation but also promote the diffusion and deposition of CO2 within frost pores, thereby significantly increasing the average frost layer thickness and deposition density.

3.3.2 Frost Layer Density

The variation in frost layer density is mainly caused by the change in frost layer structure. The frost layer density is defined as:

In the formula, ρi represents the density of dry ice; ρg represents the density of the mixed gas. αi represents the volume fraction of the dry ice phase, and its variation directly reflects the increase in the frost layer thickness and density.

As shown in Fig. 13, the local frost layer density above the cryogenic wall surface was analyzed at frosting times of 20, 40, and 60 min, with the coordinate range x = 0.2~0.6 m. The results indicate that, for a given frosting time, the local density at the downstream end of the cryogenic heat exchange channel increases significantly as the desublimation and heat transfer processes continue.

Figure 13: The local frost layer density with surface position.

Fig. 3 presents the spatial variations of CO2 mol fraction and turbulent kinetic energy in the mixed gas at t = 30 min. Upon entering the heat exchange channel, the heat transfer between the gas and the cryogenic wall surface is initially limited, such that only CO2 in the near-wall region undergoes phase change to form a frost layer. Consequently, changes in CO2 fraction at the front end of the channel are relatively minor. At this stage, the initially deposited frost layer exhibits high porosity and a loose structure, allowing gas to flow through its pores unevenly. This results in higher turbulent kinetic energy and stronger flow disturbances at the channel inlet.

As the gas continues along the cryogenic wall surface, CO2 gradually reaches a subcooled state, promoting rapid desublimation and frost layer growth. With increasing CO2 deposition, the mole fraction in the mixed gas decreases to a minimum of 0.028 at x = 0.56 m. Simultaneously, the frost layer density increases and the structure becomes more compact, leading to a smoother frost surface that reduces flow disturbances. Accordingly, the turbulent kinetic energy at the rear end of the heat exchange channel diminishes, reflecting a more stable flow over the developed frost layer.

The average frost layer density is defined as:

In the formula, mf represents the total mass of the frost layer; mj is the unit mass of the computational domain.

This study quantitatively analyzes the temporal evolution of frost layer density under different frosting durations. The results indicate that the average frost layer density gradually increases over time. During the initial stage (t = 0~10 min), the gas flow is in direct contact with the cryogenic wall surface, resulting in intense heat transfer. The near-wall gas exhibits high subcooling and supersaturation, promoting rapid CO2 nucleation and growth on the cryogenic wall surface. Consequently, the frost layer undergoes rapid formation, with its mass and density increasing sharply.

In the intermediate stage (t = 10~30 min), the frost layer continues to develop and its thickness increases rapidly, corresponding to a period of accelerated growth, although the average density continues to rise, the growth rate diminishes relative to the initial stage. This deceleration is attributed to the combined effects of decreasing CO2 mol fraction and the evolving temperature field. On one hand, ongoing desublimation reduces the CO2 fraction in the mixed gas, lowering its partial pressure and diminishing the supersaturation of the gas-solid phase transition. On the other hand, the increasing frost layer thickness elevates its thermal resistance, gradually raising the surface temperature of the frost layer and reducing the CO2 subcooling, thereby weakening the phase-change driving force.

During the stabilization stage (t = 30~60 min), the frost layer continues to form, but the growth rate of its thickness gradually declines. At this stage, CO2 desublimation occurs within the frost layer itself, leading to gradual densification of the structure and a continued increase in the average density.

As shown in Fig. 14, a lower cryogenic wall surface temperature, a moderately increased inlet velocity of the mixed gas, a lower inlet temperature, and a higher CO2 mol fraction all contribute to increasing the average frost layer density.

Figure 14: Variation of the average frost layer density under the influence of different factors (a) Cryogenic wall surface temperature (b) Mixed gas inlet velocity (c) Mixed gas inlet temperature (d) Mixed gas CO₂ mol fraction.

Specifically, decreasing the cryogenic wall surface temperature markedly enhances the degree of subcooling in the near-wall region [31]. When the cryogenic wall surface temperature decreases from 160 K to 150 K, the average frost density at t = 60 min increases by approximately 78%. A lower cryogenic wall surface temperature expands the region within the frost layer where the local temperature falls below the CO2 frost point, strengthens local supersaturation, and promotes denser CO2 desublimation and deposition per unit volume, thereby effectively increasing the overall frost density.

The influence of inlet velocity on the average frost density exhibits a nonlinear behavior. When the velocity increases from 0.4 m/s to 0.8 m/s, the average frost density at t = 60 min decreases by only 15 kg/m3, indicating a relatively limited variation. During the early stage of frosting, a moderate increase in velocity enhances convective heat and mass transfer between the gas and solid phases, facilitating the formation of a denser frost structure. However, at later stages, excessively high velocities intensify fluid shear effects, which may disturb the deposited frost, increase porosity, and consequently suppress further density growth [32].

To quantitatively elucidate the non-monotonic effect of inlet velocity on the average frost density, the dimensionless mass-transfer Péclet number (Pe) is introduced. It characterizes the ratio of the convective mass transport rate to the diffusive mass transport rate, defined as: Pe=uinHDAB where uin is the inlet velocity, H is the channel height representing the characteristic length, and DAB is the mass diffusion coefficient of CO2 in the binary gas mixture.

As the inlet velocity uin increases from 0.4 m/s to 0.8 m/s, the Pe number increases proportionally, indicating a transition to a heavily convection-dominated regime. At excessively high velocities (high Pe), the intense shear flow rapidly transports the CO2 species downstream, severely stripping the concentration boundary layer above the frost surface. This drastic reduction in local residence time limits the ability of CO2 molecules to diffuse deeply into the existing porous frost structure. Consequently, the rapid convective flow suppresses the internal diffusive densification process-which requires sufficient time for molecules to penetrate the pores and undergo internal deposition. This transport limitation results in a frost layer with higher internal porosity and, paradoxically, a lower average macroscopic density at higher velocities.

Lowering the inlet temperature of the mixed gas significantly favors the formation of a high-density frost layer. When the inlet temperature increases from 200 K to 210 K, the average frost density at t = 60 min decreases by approximately 69%, indicating that higher inlet temperatures strongly inhibit frost densification. A lower inlet temperature allows the gas mixture to reach high subcooling and supersaturation levels more rapidly in the near-wall zone, thereby promoting uniform and accelerated CO2 nucleation and the formation of a compact frost structure with fewer pores, resulting in a marked increase in average frost density.

An increase in the CO2 mol fraction of the mixed gas substantially enhances both the desublimation driving force and the total mass of deposited frost. When the CO2 mol fraction decreases from 8% to 6%, the average frost density at t = 60 min drops to only 36 kg/m3. A higher CO2 fraction not only increases the density of nucleation sites and the desublimation rate but also strengthens diffusion and mass transfer within the frost layer, leading to a higher solid CO2 fraction per unit volume. Consequently, the frost layer becomes more compact, and its average density increases significantly.

3.3.3 Frost Layer Thermal Resistance

Previous studies have demonstrated that the thermal conductivity of frost layers is primarily governed by their density, while the microstructural characteristics of the frost layer also play an important role. This paper employs the thermal conductivity model proposed by Lee et al. [33] to quantitatively analyze variations in the average frost layer thermal conductivity:

The calculation formula for the average frost layer thermal resistance is as follows:

The average frost layer thermal resistance initially increases rapidly over time, reaches a peak, and subsequently decreases gradually. In the early stage of frost formation, the rapid thickening of the frost layer leads to a corresponding fast increase in its average thermal resistance. However, as the frost layer accumulates, the growth rate of its thickness gradually slows due to the increasing thermal resistance, resulting in a deceleration of the average thermal resistance growth.

During the stable growth and deposition stage, CO2 in the mixed gas predominantly undergoes phase change within the pores of the frost layer, thereby increasing its density. The proportion of low-thermal-conductivity pores decreases, enhancing the overall frost layer thermal conductivity. As a result, the average thermal resistance gradually declines and eventually stabilizes.

As shown in Fig. 15, a lower cryogenic wall surface temperature, a higher inlet velocity of the mixed gas, a lower inlet temperature, and a higher CO2 mol fraction all contribute to an increase in the peak value of the average frost layer thermal resistance, while simultaneously shortening the time required to reach the peak. Specifically:

Figure 15: Variation of the average frost layer thermal resistance under the influence of different factors (a) Cryogenic wall surface temperature (b) Mixed gas inlet velocity (c) Mixed gas inlet temperature (d) Mixed gas CO₂ mol fraction.

Lowering the cryogenic wall surface temperature significantly increases both the degree of subcooling and the degree of supersaturation in the near-wall region, thereby accelerating the CO2 desublimation rate. At lower cryogenic wall surface temperatures, the frost layer grows to a greater thickness, resulting in a higher peak thermal resistance. Meanwhile, the accelerated deposition rate at lower temperatures promotes the formation of a denser frost structure with improved thermal conductivity, leading to a more pronounced decline in thermal resistance after the peak is reached.

Increasing the inlet velocity of the mixed gas enhances both the impingement effect of the flow on the frost surface and the overall mass transfer efficiency. Higher flow velocities facilitate CO2 transport toward the cryogenic wall surface, accelerating frost deposition and densification, and causing the thermal resistance to reach its peak earlier. However, under relatively high-velocity conditions, the frost layer tends to form a compact structure at an early growth stage, leaving limited potential for further densification. Consequently, the thermal resistance decreases more slowly after the peak, and the overall variation becomes smoother.

Lowering the inlet temperature of the mixed gas markedly enhances the degree of supersaturation near the cryogenic wall surface, strengthens the nucleation driving force, and accelerates the CO2 desublimation process. As a result, the frost layer grows to a larger thickness within a shorter time, and the thermal resistance peak occurs earlier. Moreover, frost layers formed at lower inlet temperatures exhibit a more uniform and compact structure, resulting in a more distinct decrease in thermal resistance after the peak.

An increase in the CO2 mol fraction directly raises the partial pressure of the condensable component in the mixed gas, significantly enhancing both the nucleation rate and the desublimation driving force. Under high CO2 concentration conditions, the frost layer grows rapidly and reaches a higher thermal resistance peak earlier. During the deposition stage, the accelerated filling of pores within the frost structure under high CO2 concentration leads to a more pronounced decline in thermal resistance after the peak.

4 Conclusion

In this study, a CO2 frosting phase-change model based on the classical nucleation theory was established to investigate the flow and heat transfer processes of CO2 in natural gas and its desublimation-frosting behavior under cryogenic conditions. The kinetic mechanism of CO2 frosting and the associated heat transfer characteristics on the cryogenic wall surface were elucidated. The main conclusions are as follows:

• (1)Based on gas-solid phase equilibrium theory, a gas-solid fugacity equation was developed to calculate the CO2 frost point of natural gas under different working conditions. The calculated results were compared with published experimental data, showing a deviation of less than 3%. Furthermore, a CO2 frosting phase-change model was constructed based on the classical nucleation theory, and its accuracy was verified against available experimental results.
• (2)During the frosting process, the flow characteristics of the mixed gas and the growth behavior of the frost layer exhibit strong coupling and mutual restriction. As frosting progresses, the high-temperature region in the gas phase near the outlet of the heat exchange channel expands slightly due to the increasing thermal resistance of the frost layer in the later stage. In addition, frost growth imposes a more pronounced flow resistance, forcing the mixed gas to pass through narrowed regions, which leads to a significant velocity increase compared with the initial stage.
• (3)The growth characteristics of the frost layer are significantly influenced by the cryogenic wall surface temperature, the inlet velocity of the mixed gas, the inlet temperature, and the CO2 mol fraction. A lower cryogenic wall surface temperature, a lower inlet temperature, a higher CO2 mol fraction, and a moderate inlet velocity of the mixed gas all favor an increase in the driving force and rate of CO2 desublimation, thereby promoting rapid frost formation and densification. This behavior manifests as a larger frost layer thickness and density, as well as a higher and earlier peak thermal resistance.

It should be noted that under the combined effects of various working conditions, the frost layer structure and heat transfer characteristics exhibit complex evolutionary behavior. Therefore, comprehensive optimization is required to achieve a balance between heat transfer efficiency and frost control. Finally, while this preliminary study employs a binary idealization to capture the macroscopic growth kinetics, actual PLNG systems involve complex gas mixtures. Future studies should further investigate the competitive nucleation mechanisms in multicomponent systems, particularly assessing how trace heavier hydrocarbons act as heterogeneous nucleation sites and alter the kinetic growth and internal morphology of the CO2 frost layer.