Developing sustainable solutions for thermal fluid systems requires reproducible methods that integrate computational modelling with advanced materials design. This study establishes a CFD-guided workflow for predicting and mitigating hot corrosion in Sulphur recovery unit (SRU) carbon steel piping under blended acid gas–sour water stripper (SWS) conditions. Bare steel simulations revealed severe wall loss (~29.6 mm/year at 1373 K), while an internally applied nanocomposite hybrid multilayer (NCHH) coating reduced corrosion flux by three orders of magnitude (~0.007 mm/year) without disrupting hydrodynamic stability. The multilayer architecture of cerium-modified epoxy, hafnia nanoparticles, self-healing capsules, and a hydrophobic topcoat demonstrated adaptive protection by suppressing anodic dissolution, sealing micro- cracks, and minimizing electrolyte ingress. By coupling CFD outputs with electrochemical modelling, the workflow provides scalable design rules for sustainable coating systems that extend pipeline service life, reduce maintenance costs, and lower environmental impact in energy-intensive industrial operations.

CFD Guided Methods for Sustainable Nanocomposite Multilayer Coatings: Adaptive Hot Corrosion Control in Sulphur Recovery Unit Carbon Steel Piping

This research is supported by (XLYC2203184), (U23A20657) and (LJ222410153082).

1. Physical model and mesh
<ol>
	<li>Open SpaceClaim, select the Sketch tab, and choose the Rectangle option under the Create function.</li>
	<li>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.</li>
	<li>Open ICEM, go to the File tab, select Geometry tab, then choose the Open Geometry and import the 2D geometry model.</li>
	<li>Select the Boundary of the 2D Model, open the Create Part function under the Parts tab, and assign names to the different boundaries.</li>
	<li>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.</li>
	<li>In the Initialize Blocks tab, select 2D Surface Blocking as the type and check Initialize with Settings.</li>
	<li>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.</li>
	<li>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.</li>
	<li>Under the Merge Blocks Across Curves tab, select All for the Method, and set Ignore Size to 0.0.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>Right-click on the Pre-Mesh, then select Convert to Unstructured Mesh from the context menu.</li>
	<li>Go to the Output Mesh tab and click Solver Setup. In the Solver Setup, select ANSYS Fluent as the solver and click Apply.</li>
	<li>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.</li>
</ol>
2. Frost formation simulation software operation
<ol>
	<li>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.</li>
	<li>In the Solver settings, choose Pressure-Based Type, Absolute Velocity Formulation, and Transient Time.</li>
	<li>Set Gravitational Acceleration in the Y-direction to -9.81. Under Models, click Energy and enable the Energy Equation.</li>
	<li>Under Models, enable Viscous and select k-epsilon (2 eqn). Under the k-epsilon Model section, select Standard.</li>
	<li>For Near-Wall Treatment, choose Standard Wall Functions. Under the Turbulence Multiphase Model, choose Mixture.</li>
	<li>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.</li>
	<li>Under User-Defined Functions, set Turbulent Viscosity mixture, phase-1, phase-2 to none. Under Models, click Species and enable Species Transport.</li>
	<li>In the Species Transport Model window, check Diffusion Energy Source and select phase-1.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>In the Phases setup window of the Phases tab, select Phase Property as Granular Temperature Model, set Diameter to 0.0001.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>In the Phases Interaction tab of the Multiphase Model window, select the Interfacial Area tab. Then choose the ia-symmetric option.</li>
	<li>Under the Boundary Conditions, click Inlet and select the phase-1 tab, set the Phase option to phase-1.</li>
	<li>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.</li>
	<li>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.</li>
	<li>Under the Boundary Conditions, click Inlet and select phase-2 tag, set Phase option as phase-2.</li>
	<li>In the Velocity Inlet window of phase-2 tag, select Momentum tag, set Velocity Specification Method as &#39;Magnitude, Normal to Boundary&#39;, set Reference Frame as Absolute, set Velocity Magnitude as 0, set Granular Temperature as 0.0001.</li>
	<li>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.</li>
	<li>Under the Boundary Conditions, click Outlet and select phase-1 tag, set Phase option to phase-1.</li>
	<li>In the Pressure Outlet window of phase-1 tab , select Thermal tab, set Backflow Total Temperature to 300.</li>
	<li>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.</li>
	<li>In the Pressure Outlet window of the phase-2 tab, select the Thermal tab, and set Backflow Total Temperature to 300.</li>
	<li>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.</li>
	<li>2.28 Under the Boundary Conditions, click Wall and select the cold-wall tab, set the Phase option to mixture.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
</ol>
3. Post-processing and data export configuration
<ol>
	<li>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.</li>
	<li>Under the Calculation Activities, click Solution Animations to open the Animations Definition window.</li>
	<li>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.</li>
	<li>Under the Results, choose the Surfaces option, then click the New Line/Rake Surface option to open the Line/Rake Surface window.</li>
	<li>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.</li>
	<li>Under the File tab, choose the During Calculation option under the Export. Then, click Solution Data to open the Automatic Export window.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
	<li>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.</li>
</ol>

<ol><li>Saygin, A., Basol, A. M., Arik, M. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((An+experimental+study+on+the+frost+formation+over+a+flat+plate%3A+Effect+of+frosting+on+heat+transfer%5BTitle%5D)+AND+%22Exp+Therm+Fluid+Sci%22%5BJournal%5D)">An experimental study on the frost formation over a flat plate: Effect of frosting on heat transfer.</a> Exp Therm Fluid Sci. 144, 110862(2023).</li>
<li>Fang, X., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((A+new+frictional+pressure+drop+correlation+based+on+flow+patterns+for+hydrocarbon+refrigerants+condensation+flow%5BTitle%5D)+AND+%22Int+J+Refrig%22%5BJournal%5D)">A new frictional pressure drop correlation based on flow patterns for hydrocarbon refrigerants condensation flow.</a> Int J Refrig. 170, 214-223 (2025).</li>
<li>Rong, X., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Experimental+study+on+a+multi-evaporator+mutual+defrosting+system+for+air+source+heat+pumps%5BTitle%5D)+AND+%22Appl+Energy%22%5BJournal%5D)">Experimental study on a multi-evaporator mutual defrosting system for air source heat pumps.</a> Appl Energy. 332, 120528(2023).</li>
<li>Jia, Y., Xu, X., Li, Y., Liang, X., Yao, M. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Experimental+studies+on+frost+and+defrost+of+fine+tube+bundles+under+coolant+temperature+between+%E2%88%9220+and+%E2%88%925%E2%80%AF%C2%B0C%5BTitle%5D)+AND+%22Int+J+Heat+Mass+Transf%22%5BJournal%5D)">Experimental studies on frost and defrost of fine tube bundles under coolant temperature between −20 and −5 °C.</a> Int J Heat Mass Transf. 116, 617-620 (2018).</li>
<li>Song, M., Dang, C. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Review+on+the+measurement+and+calculation+of+frost+characteristics%5BTitle%5D)+AND+%22Int.+J.+Heat+Mass+Transf%22%5BJournal%5D)">Review on the measurement and calculation of frost characteristics.</a> Int. J. Heat Mass Transf. 124, 586-614 (2018).</li>
<li>Lee, J., Lee, K. -S. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((The+behavior+of+frost+layer+growth+under+conditions+favorable+for+desublimation%5BTitle%5D)+AND+%22Int+J+Heat+Mass+Transf%22%5BJournal%5D)">The behavior of frost layer growth under conditions favorable for desublimation.</a> Int J Heat Mass Transf. 120, 259-266 (2018).</li>
<li>Lee, Y. B., Ro, S. T. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Frost+formation+on+a+vertical+plate+in+simultaneously+developing+flow%5BTitle%5D)+AND+%22Exp+Therm+Fluid+Sci%22%5BJournal%5D)">Frost formation on a vertical plate in simultaneously developing flow.</a> Exp Therm Fluid Sci. 26 (8), 939-945 (2002).</li>
<li>Mengjie, S., Shangwen, L., Hosseini, S. H., Xiaoyan, L., Zhihua, W. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((An+experimental+study+on+the+effect+of+horizontal+cold+plate+surface+temperature+on+frosting+characteristics+under+natural+convection%5BTitle%5D)+AND+%22Appl+Therm+Eng%22%5BJournal%5D)">An experimental study on the effect of horizontal cold plate surface temperature on frosting characteristics under natural convection.</a> Appl Therm Eng. 211, 118416(2022).</li>
<li>Jeong, H., Byun, S., Kim, D. R., Lee, K. S. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Frost+growth+mechanism+and+its+behavior+under+ultra-low+temperature+conditions%5BTitle%5D)+AND+%22Int+J+Heat+Mass+Transf%22%5BJournal%5D)">Frost growth mechanism and its behavior under ultra-low temperature conditions.</a> Int J Heat Mass Transf. 169, 120941(2021).</li>
<li>Barzanoni, Y., Noorshams, O., Basirat Tabrizi, H., Damangir, E. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Experimental+investigation+of+frost+formation+on+a+horizontal+cold+cylinder+under+cross+flow%5BTitle%5D)+AND+%22Int+J+Refrig%22%5BJournal%5D)">Experimental investigation of frost formation on a horizontal cold cylinder under cross flow.</a> Int J Refrig. 34 (4), 1174-1180 (2011).</li>
<li>Jones, B. W., Parker, J. D. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Frost+formation+with+varying+environmental+parameters%5BTitle%5D)+AND+%22J+Heat+Transf%22%5BJournal%5D)">Frost formation with varying environmental parameters.</a> J Heat Transf. 97 (2), 255-259 (1975).</li>
<li>Cui, J., Li, W. Z., Liu, Y., Jiang, Z. Y. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((A+new+time-+and+space-dependent+model+for+predicting+frost+formation%5BTitle%5D)+AND+%22Appl+Therm+Eng%22%5BJournal%5D)">A new time- and space-dependent model for predicting frost formation.</a> Appl Therm Eng. 31 (4), 447-457 (2011).</li>
<li>Lenic, K., Trp, A., Frankovic, B. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Transient+two-dimensional+model+of+frost+formation+on+a+fin-and-tube+heat+exchanger%5BTitle%5D)+AND+%22Int+J+Heat+Mass+Transf%22%5BJournal%5D)">Transient two-dimensional model of frost formation on a fin-and-tube heat exchanger.</a> Int J Heat Mass Transf. 52 (1-2), 22-32 (2009).</li>
<li>Yu, J., Huo, R., Shen, H., Li, X., Zhu, Z. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((A+simulation+study+on+the+condensation+flow+and+thermal+control+characteristics+of+mixed+refrigerant+in+a+dimpled+tube%5BTitle%5D)+AND+%22Appl+Therm+Eng%22%5BJournal%5D)">A simulation study on the condensation flow and thermal control characteristics of mixed refrigerant in a dimpled tube.</a> Appl Therm Eng. 231, 120889(2023).</li>
<li>Yu, J., Jiang, Y., Cai, W., Li, X., Zhu, Z. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Condensation+flow+patterns+and+heat+transfer+correction+for+zeotropic+hydrocarbon+mixtures+in+a+helically+coiled+tube%5BTitle%5D)+AND+%22Int+J+Heat+Mass+Transf%22%5BJournal%5D)">Condensation flow patterns and heat transfer correction for zeotropic hydrocarbon mixtures in a helically coiled tube.</a> Int J Heat Mass Transf. 143, 118500(2019).</li>
<li>You, Y., Wang, S., Lv, W., Chen, Y., Gross, U. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((A+CFD+model+of+frost+formation+based+on+dynamic+meshes+technique+via+secondary+development+of+ANSYS+Fluent%5BTitle%5D)+AND+%22Int+J+Heat+Fluid+Flow%22%5BJournal%5D)">A CFD model of frost formation based on dynamic meshes technique via secondary development of ANSYS Fluent.</a> Int J Heat Fluid Flow. 89, 108807(2021).</li>
<li>Cai, W., Fang, X., Li, S., Qiu, G. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((A+modified+CFD+model+for+frosting+on+a+horizontal+plate%5BTitle%5D)+AND+%22Int+J+Heat+Mass+Transf%22%5BJournal%5D)">A modified CFD model for frosting on a horizontal plate.</a> Int J Heat Mass Transf. 229, 125726(2024).</li>
<li>Boyina, K. S., et al. <a target="_blank" href="http://www.ncbi.nlm.nih.gov/pubmed?term=((Condensation+frosting+on+meter-scale+superhydrophobic+and+superhydrophilic+heat+exchangers%5BTitle%5D)+AND+%22Int+J+Heat+Mass+Transf%22%5BJournal%5D)">Condensation frosting on meter-scale superhydrophobic and superhydrophilic heat exchangers.</a> Int J Heat Mass Transf. 145, 118694(2019).</li>
</ol>

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

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 that the maximum deviation in frost density is within 29.72% (MARD = 16.06%). Moreover, the simulated frost distribution closely matches experimental observations. These results are of great significance for the construction of high-accuracy numerical frost models.
While numerous numerical frost models have been proposed in the literature11,12,13,14,15,16, the present work distinguishes itself by validating a broader range of parameters and achieving smaller discrepancies between simulation and experiment. Crucially, it incorporates the dynamic evolution of frost density with respect to both operating conditions and time, thereby more faithfully representing the actual frosting process. Consequently, the model presented here exhibits superior theoretical rigor and general applicability.
Given that the ultimate objective extends beyond predicting frosting on horizontal cold surfaces -- encompassing three-dimensional and inclined geometries in practical applications -- further investigation integrating additional physical models is required. Equally important is the acquisition of extensive experimental data to validate numerical models under varied physical assumptions. Such efforts will deepen our understanding of how different parameters influence the frosting process and will inform the design of heat exchangers and the development of anti- and defrosting strategies in related fields.
Although surface wettability was not explicitly addressed in the present study, it is recognized as a factor that can affect frost nucleation behavior18. Future investigations may consider incorporating wettability-related parameters to extend the model&#39;s applicability to a broader range of surface conditions.

Frost formation is a common phenomenon observed across many different fields. The accumulation of frost on heat exchanger surfaces significantly impairs heat transfer efficiency1, obstructs fluid flow, and disrupts the overall performance of heat exchangers2, ultimately hindering their normal operation3. Therefore, understanding the mechanisms and behavior of frost formation is critical for addressing this problem in refrigeration systems4. 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 surface5,6,7,8,9,10. Numerous experimental findings have indicated that lower incoming air temperatures tend to result in thicker frost layers6, while higher humidity levels contribute to the formation of denser frost layers7. 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 density8. 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 hill9. 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 surfaces10. In addition, several studies11,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 theory11. 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 theory12. Their predictions of frost thickness showed a deviation of less than 13% from the experimental data provided by Lenic et al.13. In parallel, CFD studies on condensation in structured tubes have shown that geometric features like dimples14 or helical pitches15 enhance local heat and mass transfer. Recently, You et al.16developed 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&#39;s potential in resolving complex phase change phenomena, offering valuable insights for frost formation modeling.
In conclusion, a significant number of studies5,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 studies6,7,8,9,10,11,12,13 primarily validate models using frost thickness, limiting their broader applicability16. 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&#39;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.

<table><tbody><tr><td>Fluent</td><td>ANSYS</td><td></td><td></td></tr><tr><td>ICEM</td><td>Dassault system</td><td></td><td></td></tr><tr><td>SpaceClaim</td><td>ANSYS</td><td></td><td></td></tr></tbody></table>

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.

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 distribution with experimental data from various studies. Furthermore, the frost morphology predicted by the model aligns closely with experimental observations reported in the literature.
Figure 1 provides a schematic of the simulated physical model. The computational domain includes a fully developed inlet region at the front, with the upper wall set as adiabatic. The velocity inlet is defined on the left boundary, while the right boundary serves as the pressure outlet. This setup ensures realistic boundary conditions for frost formation under convective flow.
Figure 2 illustrates the simulated volume fraction distribution of frost at 3600 s. The simulation conditions are summarized in Table 1, where the wall temperature was 252.65 K, the air temperature was 292.95 K, the inlet velocity was 0.6 m/s, and the relative humidity was 58%. The maximum volume fraction reaches 0.07. The frost layer is evenly distributed, with a slightly higher accumulation at the center compared to the sides, reflecting the symmetry of flow and thermal conditions. Notably, the simulation results near the inlet show excellent agreement with experimental observations, supporting the accuracy of the model in capturing local frost growth behavior.
Table 1 shows all tested cases and their respective simulation conditions. The cold wall temperature ranges from 252.65 K to 269.15 K, the inlet air temperature from 275.15 K to 292.95 K, the inlet velocity from 0.6 m/s to 0.92 m/s, and the relative humidity from 58% to 85%.
Comparison with previous work: Figure 3 illustrates a comparison between the simulation results for frost thickness and the experimental data17, under the conditions of a cold wall temperature of 269.15 K, inlet humid air temperature of 289.15 K, inlet velocity of 0.7 m/s, and relative humidity of 80%. The simulation results show good agreement with experimental data, particularly in capturing the temporal evolution of frost growth. A rapid increase in frost thickness is observed during the initial stage, followed by a gradual growth phase. This consistency validates the accuracy of the proposed simulation model in predicting key characteristics of frost formation, including frost density and surface temperature, and confirms its capability to capture the underlying physical processes with reasonable fidelity. As the frosting process progresses, the deviation between the simulated frost thickness and the measured data increases. The analysis suggests that although this study accounts for the time-dependent variation in density, the overall trend still deviates from the actual results, which in turn affects the accuracy of the simulated frost thickness compared to the experimental data.
Successful simulations: Figure 4 presents a comparison between the simulated and experimental frost layer thicknesses17, with the horizontal axis denoting the experimental values and the vertical axis indicating the simulation results. The overall deviation remains within 20%, with a maximum deviation of -20% and a mean absolute relative deviation (MARD) of 8.97%. These results demonstrate that under optimal conditions, the simulation method yields reliable predictions with relatively small errors. The good agreement between simulation and experiment underlines the robustness of the modeling approach in capturing frost formation dynamics when boundary conditions and parameters are well constrained.
Sub-optimal results: Figure 5 compares the simulated and experimental frost density values17, with the horizontal axis representing the experimental data and the vertical axis indicating the corresponding simulation results. Overall, the deviations remain within 30%, with a maximum deviation of 29.72% and a mean absolute relative deviation (MARD) of 16.06%. However, more pronounced discrepancies are evident in certain sub-optimal simulations, particularly in the prediction of frost density. These larger deviations highlight the sensitivity of the model to computational factors such as boundary conditions, numerical schemes, and uncertainties in input parameters, underscoring both the complexity of the frost formation process and the need for further model refinement.
In this study, both thickness and density were simultaneously simulated to demonstrate the reliability of the proposed method. Our results highlight that density also deserves careful consideration. As shown in Figure 3 for one representative case, the deviation between simulation and experimental thickness increases over time. For Figure 4 and Figure 5, an analysis to distinguish markers for different test cases was conducted in our previous work17. With the present manuscript the intention was to provide more detailed operational procedures and new findings.
Analysis and suggestions for improvement: While the results from the optimal simulations support the reliability of the CFD technique, the larger deviations observed in the frost density calculations point to potential areas for model improvement. Factors such as the refinement of mesh resolution, better representation of material properties, and more precise input conditions may reduce errors. Further investigation into the sensitivity of the model to different parameters and environmental factors could help improve the overall predictive capability of the CFD simulations.
<img alt="Fluid flow diagram; velocity inlet to pressure outlet, adiabatic and cold walls, thermal analysis." class="xfigimg" src="/files/ftp_upload/68133/68133fig1.jpg" data-altupdatedat="1762509930000" data-alt="Figure 1"> 
Figure 1: Schematic diagram of the simulated physical model. The front section is the fully developed region, the upper wall is an adiabatic wall, the left side is designated as the velocity inlet, while the right side is the pressure outlet. H is the height of the physical model; L1 is the length of the fully developed inlet region; L2 is the length of the cooling wall; and L3 is the total length of the physical model. <a href="https://www.jove.com/files/ftp_upload/68133/68133fig1large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
<img alt="Fluid dynamics simulation; concentration distribution; color chart; volume fraction visualization." class="xfigimg" src="/files/ftp_upload/68133/68133fig2.jpg" data-altupdatedat="1762509930000" data-alt="Figure 2"> 
Figure 2: Volume fraction distribution contour at 3600 s of frost formation. The simulation conditions are summarized in Table 1, where the wall temperature was 252.65 K, the air temperature was 292.95 K, the inlet velocity was 0.6 m/s, and the relative humidity was 58%. The maximum volume fraction reaches 0.07. <a href="https://www.jove.com/files/ftp_upload/68133/68133fig2large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
<img alt="Frost growth measurement graph; frost thickness vs. time; experimental data comparison; conditions." class="xfigimg" src="/files/ftp_upload/68133/68133fig3.jpg" data-altupdatedat="1762509930000" data-alt="Figure 3"> 
Figure 3: Comparison between the simulation results for frost thickness and the experimental data. The experimental conditions were a cold wall temperature of 269.15 K, an inlet humid air temperature of 289.15 K, an inlet velocity of 0.7 m/s, and a humidity of 80%. Here, Twall is the cold wall temperature. Tair is the inlet humid air temperature, and RH is the relative humidity. It can be seen that the simulation results closely match the experimental data. This figure has been modified from17. <a href="https://www.jove.com/files/ftp_upload/68133/68133fig3large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
<img alt="Frost thickness correlation graph; experimental vs. simulated data, with error lines (+15%, -20%)." class="xfigimg" src="/files/ftp_upload/68133/68133fig4.jpg" data-altupdatedat="1762509930000" data-alt="Figure 4"> 
Figure 4: Comparison between the simulated frost thickness results and the experimental data. The horizontal axis represents the experimental thickness, and the vertical axis represents the simulated thickness. The continuous line represents the ideal agreement line (y = x), where the simulated frost thickness is equal to the experimental value. The two dotted lines indicate the ± deviation bands, corresponding to +15% and -20% of the experimental frost thickness. The red dots denote the paired data points obtained from the comparison between the simulated results and the experimental measurements. <a href="https://www.jove.com/files/ftp_upload/68133/68133fig4large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
<img alt="Scatter plot of frost density analysis: simulated vs. experimental density with error margins." class="xfigimg" src="/files/ftp_upload/68133/68133fig5.jpg" data-altupdatedat="1762509930000" data-alt="Figure 5"> 
Figure 5: Comparison between the simulated frost density results and the experimental data. The horizontal axis represents the experimental frost density, and the vertical axis represents the simulated frost density. The continuous line represents the ideal agreement line (y = x), where the simulated frost density is equal to the experimental value. The two dotted lines indicate the ± deviation bands, corresponding to +30% and -25% of the experimental frost density. The red dots denote the paired data points obtained from the comparison between the simulated results and the experimental measurements. <a href="https://www.jove.com/files/ftp_upload/68133/68133fig5large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
<table border="1" fo:keep-together.within-page="1" fo:keep-with-next.within-page="always">
	<tbody>
		<tr>
			<td>Name</td>
			<td>Density [kg/m2]</td>
			<td>Cp (Specific Heat)[J/(kg K)]</td>
			<td>Thermal Conductivity [W/(m K)]</td>
			<td>Viscosity [kg/(m s)]</td>
			<td>Pure Solvent Melting Heat [J/kg]</td>
			<td>Solidus Temperature [K]</td>
			<td>Liquidus Temperature [K]</td>
		</tr>
		<tr>
			<td>steel</td>
			<td>8030</td>
			<td>502.48</td>
			<td>16.27</td>
			<td>-</td>
			<td>-</td>
			<td>-</td>
			<td>-</td>
		</tr>
		<tr>
			<td>air</td>
			<td>1.225</td>
			<td>1006.43</td>
			<td>0.0242</td>
			<td>1.79E-05</td>
			<td>-</td>
			<td>-</td>
			<td>-</td>
		</tr>
		<tr>
			<td>water</td>
			<td>0.999972</td>
			<td>4182</td>
			<td>0.6</td>
			<td>0.001003</td>
			<td>334000</td>
			<td>273.15</td>
			<td>273.15</td>
		</tr>
	</tbody>
</table>
Table 1: Tested cases and their respective simulation conditions. The cold wall temperature ranges from 252.65 K to 269.15 K, the inlet air temperature from 275.15 K to 292.95 K, the inlet velocity from 0.6 m/s to 0.92 m/s, and the relative humidity from 58% to 85%. Twall is the cold wall temperature, Tair is the inlet humid air temperature, and RH is the relative humidity. The simulation results closely match the experimental data.

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.

cfd-simulation-of-frost-on-horizontal-cold-surfaces

CFD Simulation of Frost on Horizontal Cold Surfaces

Frost formation is a common phenomenon observed across many different fields, including refrigeration, construction, and natural gas processing. However, ...

Mustafa Demirci

Thermal and fluid systems are integral energy sources in today's technology-driven world. Enhancing their performance is crucial for optimizing energy consumption and reducing environmental impact. Efficiency and sustainability studies, utilizing experimental analysis and computational fluid dynamics (CFD), are essential for developing innovative solutions to improve these systems. This Methods Collection aims to showcase advanced approaches for boosting the efficiency, sustainability, and overall performance of thermal and fluid systems.

JoVE Engineering encompasses a broad range of experimental and instrumental techniques utilized in physics research. Investigations in this area strive to address and answer a broad range of scientific questions, such as device mechanisms and efficiencies, using physical tools. This approach often requires a combination of specialties, and research in this area tends to be interdisciplinary with contributions from mechanical, electrical, and chemical engineers.

Engineering

Experimental and Computational Fluid Dynamics (CFD) Studies on Thermal and Fluid Systems 

Thermal and fluid systems are integral energy sources in today's technology-driven world. Curated by Mustafa Demirci. Explore Engineering techniques on JoVE.

Guest Editor

Collection Overview

Articles

CFD Simulation of Frost on Horizontal Cold Surfaces

Abstracts