CFD Simulation of Frost on Horizontal Cold Surfaces

Summary

Here, we present a numerical model for simulating frost formation on horizontal cold surfaces using a Eulerian multiphase model with the Lee phase change approach. The model dynamically updates frost volume fraction to capture density variation and is validated against experimental data on frost thickness, density, and distribution.

Abstract

Frost formation is a common phenomenon observed across many different fields, including refrigeration, construction, and natural gas processing. However, due to its complex nature, developing an accurate and reliable numerical model remains a significant challenge. Despite previous efforts to tackle this issue, current models still have certain limitations. This paper introduces a modified numerical model for frost formation, developed based on the fundamental mechanisms underlying frost formation. The model utilizes an Eulerian multiphase flow approach coupled with the Lee phase change model. In addition, the approach for determining the maximum frost volume fraction is updated, enabling the model to consider density variations during the frosting process. The model is rigorously validated by comparing it with experimental data on thickness, density, and distribution from various studies. The results indicate that the mean absolute relative deviation (MARD) for frost thickness is 8.97%, while the MARD for density is 16.06%. Furthermore, the frost morphology predicted by the model closely matches the experimental observations reported in the reference.

Introduction

Frost formation is a common phenomenon observed across many different fields. The accumulation of frost on heat exchanger surfaces significantly impairs heat transfer efficiency¹, obstructs fluid flow, and disrupts the overall performance of heat exchangers², ultimately hindering their normal operation³. Therefore, understanding the mechanisms and behavior of frost formation is critical for addressing this problem in refrigeration systems⁴. In recent decades, a substantial body of research has been dedicated to investigating the causes and characteristics of frost formation in these systems.

Experimental studies have demonstrated that the formation of frost is influenced by various factors, including air temperature, humidity, and the temperature of the cold surface^5,6,7,8,9,10. Numerous experimental findings have indicated that lower incoming air temperatures tend to result in thicker frost layers⁶, while higher humidity levels contribute to the formation of denser frost layers⁷. Song et al. studied frost formation on horizontal surfaces and found that cyclic temperature variations of the cold surface can cause melting at the frost layer interface, which significantly influences the frost formation rate, frost layer thickness, and dynamic frost density⁸. Other research has examined both the morphology and distribution of frost. Jeong et al. observed in their experiments that frost initially forms near the inlet, leading to the occurrence of a phenomenon known as frost hill⁹. Noorshams et al. investigated frost formation on a horizontal circular tube surface and found that the frost layers on the front and rear surfaces of the cylinder were thicker than those on the top surfaces¹⁰. In addition, several studies^{11,12,13,14,15,16} have developed models to predict the one-dimensional frost layer thickness, utilizing experimental frost formation patterns along with both theoretical and empirical approaches. Jones and Parker developed a predictive model for frost thickness based on molecular diffusion theory¹¹. The discrepancy between their model and experimental data remained below 30% over a 3 h period. With the progression of computational technology, an increasing number of researchers have turned to Computational Fluid Dynamics (CFD) for simulating frost formation. Unlike traditional one-dimensional models, CFD simulations provide significant benefits, particularly in visualizing frost thickness distribution and temperature profiles. Cui et al. performed CFD simulations of frost formation based on nucleation theory¹². Their predictions of frost thickness showed a deviation of less than 13% from the experimental data provided by Lenic et al.¹³. In parallel, CFD studies on condensation in structured tubes have shown that geometric features like dimples¹⁴ or helical pitches¹⁵ enhance local heat and mass transfer. Recently, You et al.¹⁶developed a dynamic mesh-based CFD model that characterizes the frost layer as a growing porous medium and incorporates vapor diffusion directly, achieving a relative deviation of less than 5% while maintaining low computational cost. These findings underscore CFD's potential in resolving complex phase change phenomena, offering valuable insights for frost formation modeling.

In conclusion, a significant number of studies^{5,6,7,8,9,10,11,12,13} have explored frost formation on cold surfaces, contributing to an evolving understanding of frost formation patterns under various parameters. Although multiple numerical models have been developed, incorporating various dimensions and mechanisms, they frequently lack comprehensive validation. Most studies^{6,7,8,9,10,11,12,13} primarily validate models using frost thickness, limiting their broader applicability¹⁶. To overcome these limitations, this paper introduces a numerical model that integrates the Lee phase change model and the Eulerian multiphase flow model, with an emphasis on the core mechanisms underlying frost formation. Moreover, a novel method is introduced to calculate the upper limit of the frost volume fraction, taking into account time-dependent changes in frost density, thus addressing the shortcomings of previous models. The proposed model's accuracy and reliability are assessed from various angles, including frost thickness, density fluctuations, and frost formation patterns across different experimental conditions. This extensive validation offers a robust theoretical framework for more accurately predicting frost behavior in real-world applications.

Protocol

1. Physical model and mesh

1. Open SpaceClaim, select the Sketch tab, and choose the Rectangle option under the Create function.
2. Create a 2D geometric model on the XOY plane with a length of 500 mm along the x-axis and a width of 15 mm along the y-axis.
3. Open ICEM, go to the File tab, select Geometry tab, then choose the Open Geometry and import the 2D geometry model.
4. Select the Boundary of the 2D Model, open the Create Part function under the Parts tab, and assign names to the different boundaries.
5. Select the Blocking tab, then click Create Block. In the Create Block window, check the Inherit Part Name option. Then go to the Create Block tab and select Initialize Block.
6. In the Initialize Blocks tab, select 2D Surface Blocking as the type and check Initialize with Settings.
7. Under the Choose tab, select Form Surface for the Method option. Then, under Surfaces, click the Select Surface(s) icon and choose the 2D Plane in the graphics window.
8. Under the Surface Blocking tab, select Mostly Mapped for the Method option, All Quad for the Free Face Mesh Type, and ICEM CFD Quad for the Free Face Mesh Method.
9. Under the Merge Blocks Across Curves tab, select All for the Method, and set Ignore Size to 0.0.
10. In the Pre-Mesh Params tab, under the Meshing Parameters, click Edge Params. Under Edge, click the Select Edge(s) icon and choose the Edge in the X-direction in the graphical interface.
11. Under the Meshing Parameters, set Nodes to 500, Spacing 1 to 1e+10, Ratio 1 to 2, Spacing 2 to 1e+10, Ratio 2 to 2, and Max Space to 1e+10. Then check the Copy Parameters option.
12. Under the Meshing Parameters option, set Mesh Law to BiGeometric and click Apply. In the Pre-Mesh Params tab, under the Meshing Parameters option, click Edge Params.
13. For the Edge option under Meshing Parameters, click the Select Edge(s) icon and choose the edge in the y-direction in the graphical interface.
14. Under the Meshing Parameters, set Nodes to 150, Spacing 1 to 1e+10, Ratio 1 to 2, Spacing 2 to 1e+10, Ratio 2 to 2, and Max Space to 1e+10. Then check the Copy Parameters option.
15. Under the Meshing Parameters option, set the Mesh Law to BiGeometric and click Apply. Under the Blocking option, activate the Pre-Mesh function and click Yes when prompted.
16. Right-click on the Pre-Mesh, then select Convert to Unstructured Mesh from the context menu.
17. Go to the Output Mesh tab and click Solver Setup. In the Solver Setup, select ANSYS Fluent as the solver and click Apply.
18. Under the Output tab, click Write Input. In New Window 1, click Save; in new window 2, click Yes, then Save; in new window 3, click Open; and in new window 4, click Done.

2. Frost formation simulation software operation

1. Open Ansys Fluent. Go to the File tab, then select the Mesh option under the Read. Go to the Scale Mesh, and set mm under Mesh Was Created In. The procedure presented here is based on the experimental conditions defined in Case 1 of Table 1.
2. In the Solver settings, choose Pressure-Based Type, Absolute Velocity Formulation, and Transient Time.
3. Set Gravitational Acceleration in the Y-direction to -9.81. Under Models, click Energy and enable the Energy Equation.
4. Under Models, enable Viscous and select k-epsilon (2 eqn). Under the k-epsilon Model section, select Standard.
5. For Near-Wall Treatment, choose Standard Wall Functions. Under the Turbulence Multiphase Model, choose Mixture.
6. Under Model Constants option, set Cmu, C1-Epsilon, C2-Epsilon, TKE Prandtl Number, TDR Prandtl Number, Dispersion Prandtl Number, Energy Prandtl Number, Wall Prandtl Number, and Turbulent Schmidt Number as 0.09, 1.44, 1.92, 1, 1.3, 0.75, 0.85, 0.85, and 0.7, respectively.
7. Under User-Defined Functions, set Turbulent Viscosity mixture, phase-1, phase-2 to none. Under Models, click Species and enable Species Transport.
8. In the Species Transport Model window, check Diffusion Energy Source and select phase-1.
9. Under Models, enable Multiphase, and choose the Eulerian model. In the Multiphase model window, select Volume Fraction Parameters Formulation Implicit, and set the Number of Eulerian Phases to 2.
10. In the Multiphase model window, select Phases tab, select phase-1 - Primary Phase, set Name as phase-1, and set Phase Material to mixture-template.
11. In the Multiphase model window, select Phases tab, select phase-2 - Secondary Phase, set Name as phase-2, set Phase Material to ice, and enable Granular.
12. In the Phases setup window of the Phases tab, select Phase Property as Granular Temperature Model, set Diameter to 0.0001.
13. In the Granular Properties window of the Phases tab, set Granular Viscosity to 1e-05, set Granular Bulk Viscosity to 30, set Solids Pressure to lun-et-al, set Granular Temperature to algebraic, set Frictional Viscosity to none, set Packing Limit to user-defined, set Radial Distribution to lun-et-al, set Elasticity Modulus to derived.
14. In the Phases Interaction tab of the Multiphase Model window, select the Forces tag. Then select phase-1, phase-2, set the coefficient to wen-yu.
15. In the Phases Interaction tab of the Multiphase Model window, select the Forces tag. Then select phase-2, phase-2, set Restitution Coefficient to 0.9.
16. In the Phases Interaction tab of the Multiphase Model window, select the Interfacial Area tab. Then choose the ia-symmetric option.
17. Under the Boundary Conditions, click Inlet and select the phase-1 tab, set the Phase option to phase-1.
18. In the Velocity Inlet window of the phase-1 tab, select the Momentum tab, set the Velocity Specification Method to Magnitude, Normal to Boundary, set the Reference Frame to Absolute, and set the Velocity Magnitude to 0.6.
19. In the Velocity Inlet window of phase-1, select Thermal tag, set Temperature to 292.8. In the Velocity Inlet window of phase-1 tag, select Species tag, set h2o as 0.008202.
20. Under the Boundary Conditions, click Inlet and select phase-2 tag, set Phase option as phase-2.
21. In the Velocity Inlet window of phase-2 tag, select Momentum tag, set Velocity Specification Method as 'Magnitude, Normal to Boundary', set Reference Frame as Absolute, set Velocity Magnitude as 0, set Granular Temperature as 0.0001.
22. In the Velocity Inlet window of phase-2 tag, select Thermal tab, set Temperature to 273. In the Velocity Inlet window of phase-2 tag, select Multiphase tab, set Volume Fraction to 0.
23. Under the Boundary Conditions, click Outlet and select phase-1 tag, set Phase option to phase-1.
24. In the Pressure Outlet window of phase-1 tab , select Thermal tab, set Backflow Total Temperature to 300.
25. In the Pressure Outlet window of the phase-1 tab, select the Species tab, set h2o to 0. Under the Boundary Conditions, click Outlet and select the phase-2 tab, set the Phase option to phase-2.
26. In the Pressure Outlet window of the phase-2 tab, select the Thermal tab, and set Backflow Total Temperature to 300.
27. In the Pressure Outlet window of the phase-2 tab, select the Multiphase tab, set Backflow Granular Temperature to 0.0001, set Volume Fraction Specification Method to Backflow Volume Fraction, and set Backflow Volume Fraction to 0.
28. 2.28 Under the Boundary Conditions, click Wall and select the cold-wall tab, set the Phase option to mixture.
29. In the wall window of the cold-wall tab, select the Momentum tab, set the Wall Motion option to Stationary Wall, set the Wall Roughness Models option to Standard, set the Roughness Height to 0, and set the Roughness Constant to 0.5.
30. In the wall window of the cold-wall tab, select Thermal tab, select Temperature as Thermal Conditions, set Temperature to 252.65, and set Material to steel.
31. Under the Solution, click Methods and open the Solution Methods window. In the Solution Methods window, select Phase Coupled SIMPLE as Pressure-Velocity Coupling Scheme, select Least Squares Cell Based as Gradient Spatial Discretization, select Second Order as Pressure Spatial Discretization, select First Order Upwind as Density Spatial Discretization, select First Order Upwind as Momentum Spatial Discretization, select First Order Upwind as Volume Fraction Spatial Discretization, select First Order Upwind as Turbulent Kinetic Energy Spatial Discretization, select First Order Upwind as Turbulent Dissipation Rate Spatial Discretization, select First Order Upwind as Energy Spatial Discretization, select First Order Upwind as phase-1 h2o Spatial Discretization, select First Order Implicit as Transient Formulation.
32. Under the Solution, click Controls and open the Solution Controls window. Set Pressure Under-Relaxation Factor to 0.4, set Density Under-Relaxation Factor to 1, set Body Forces Under-Relaxation Factor to 1, set Momentum Under-Relaxation Factor to 0.4, set Volume Fraction Under-Relaxation Factor to 0.4, set Granular Temperature Under-Relaxation Factor to 0.3, set Turbulent Kinetic Energy Under-Relaxation Factor to 0.3, set Turbulent Dissipation Rate Under-Relaxation Factor to 0.3, set Turbulent Viscosity Under-Relaxation Factor to 0.3, set Energy Under-Relaxation Factor to 0.4, set phase-1 h2o Under-Relaxation Factor to 0.4.
33. Under the Solution, click Initialization and open Solution Initialization window. Select Standard Initialization as the Initialization Method, select Relative to Cell Zone as Reference Frame.
34. In the Solution Initialization window, set Gauge Pressure to 0, set Turbulent Kinetic Energy to 0.00135, set Turbulent Dissipation Rate to 0.001143987, set phase-1 X Velocity to 0, set phase-1 Y Velocity to 0, set phase-1 h2o to 0.008202, set phase-2 X Velocity to 0, set phase-2 Y Velocity to 0, set phase-2 Volume Fraction to 0, set phase-2 Granular Temperature to 0.0001, set phase-2 Temperature to 273. Then, click Initialize.

3. Post-processing and data export configuration

1. Under the Results, click Contours to open the Contours window. In the Contours window, enable Filled, Node Values, Boundary Values, Global Range, and Auto Range options. Select Phases as the Contours of Type and select Volume fraction. Then, choose Phase 2 as the Phase, and click Save/Display.
2. Under the Calculation Activities, click Solution Animations to open the Animations Definition window.
3. In the Animations Definition window, set Record after every to 1, and choose time step. Select HSF File as the Storage Type. Select contour-1 under the Animation Object options, and then click OK.
4. Under the Results, choose the Surfaces option, then click the New Line/Rake Surface option to open the Line/Rake Surface window.
5. In the Line/Rake Surface window, enable Line, set x0 [m] to 0.21, set x1 [m] to 0.21, set y0 [m] to 0, set y1 [m] to 0.015, and then click Create.
6. Under the File tab, choose the During Calculation option under the Export. Then, click Solution Data to open the Automatic Export window.
7. In the Automatic Export window, choose ASCII as the File Type. Select Cell Center as Location, select Space as Delimiter. Set Export Data Every to 1 and then choose Time Step. Under the Surfaces option, choose line-1. Under Quantities option, select Density (phase-2) and Volume Fraction (phase-2). Then, click Browse, open the Select File window. Under the Select File window, click OK. Under the Automatic Export window, click OK.
8. Under the Calculation Activities, click Autosave (Every Flow Time) to open the Autosave window. Then set Save Data File Every [s] to 100, select Flow Time, select Only if Modified as the Save Associated Case Files Type, and click OK.
9. Under the Solution, click Run Calculation to open the Run Calculation window. Then select Fixed as the Time Advancement Type and User-Specified as the Time Advancement Method.
10. In the Run Calculation window, set Number of Time Steps to 7200, set Time Step Size to 1, set Max Iterations/Time Step to 20, set Reporting Interval to 1, and set Profile Update Interval to 2. Then, click Calculate.

Results

The proposed modified numerical model effectively captures the key features of frost formation. It is based on the fundamental mechanisms underlying frost growth and employs an Eulerian multiphase flow approach, coupled with the Lee phase change model. This approach allows the model to update the maximum frost volume fraction, thereby accounting for density variations throughout the frosting process. The simulations show that the model is rigorously validated by comparing its predictions of frost thickness, density, and ...

Discussion

This study develops a numerical model capable of simulating frost formation on a low-temperature, horizontal cold surface by dynamically adjusting the upper limit of the volume fraction in accordance with both time and operating conditions, thereby reproducing variations in frost density. Although the validation presented here is confined to conventional cold-surface temperatures, the model guarantees that the maximum deviation between simulated and experimental frost thickness does not exceed -20% (MARD = 8.97%), and th...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Fluent	ANSYS
ICEM	Dassault system
SpaceClaim	ANSYS