Numerical validation of ice formation on a lattice structure evaporator

This paper presents a comparative analysis between experimental and numerical results of the ice formation process on the surface of the evaporator in a Sparkling Water Dispenser. The experimental setup comprises an R290 refrigeration cycle that is fully equipped with thermocouples, pressure transducers, a wattmeter, and a Coriolis-meter. The evaporator, which is specifically designed to enhance the ice formation speed, is situated within a water-filled tank, where the ice formation process takes place. Due to the phase change phenomenon, which involves the interface between two phases moving, solving transient heat transfer problems involving solidifications can be inherently challenging. Analytical solutions are only possible under certain simplified circumstances. In cases where exact solutions are not available, semi-analytic, approximate, and numerical methods can be utilized to address phase-change problems. A numerical model of the solidification process based on the energy equation and conjugate heat transfer was developed using COMSOL Multiphysics. The software was found to be effective in simulating the physical processes associated with heat transfer through conduction and convection, as well as the behaviour of phase change. The results showed a very good agreement between experimental and numerical results.


Introduction
Refrigeration systems are widely used in many industries and applications.Among the various applications of refrigeration systems, we also find water dispensers.These devices extract potable water from a source (usually from the water supply or a reservoir), cool it to the setpoint temperature (typically below 10 °C), and dispense it.Given the physical properties of water and the specific requirements of this type of device (which must be small and consume low electrical power), an auxiliary system is needed to handle intermittent peaks of thermal power.Therefore, these devices are generally associated with a thermal storage unit, which mitigates load peaks and prevents the refrigeration circuit from being oversized and activated whenever there is a demand for chilled water.Thermal storage plays a central role in the operation of a Sparkling Water Dispenser (SWD).Due to its availability, low cost, and high energy density (333 kJ/kg), ice water storage is the most widely used technology in SWDs.
Another critical component of a refrigeration system is the evaporator.In these devices, the evaporator is immersed inside the storage unit and charges the storage by solidifying the water.Ice formation can significantly reduce the heat transfer rate and increase the system's energy consumption.The ice formation process is complex and depends on various factors, such as the evaporator geometry, ambient and refrigerant temperatures, and the flow rate of the refrigerant.To better understand this process, researchers have developed various experimental and numerical models to investigate the ice formation mechanism and optimize the design of the heat exchanger [1][2][3][4].Given objective to optimize the heat exchanger geometry, the utilization of a numerical model offers a significant advantage in a cost-effective and time-efficient manner.The phase change process has been studied for a variety of geometries, the most frequent of which being shell-and-tube systems, encapsulated systems, and ice-on-coil systems.Ice-on-coil energy storage tanks are widely used in phase change commercial cooling applications [5][6][7][8][9].Numerous studies on ice formation on the coil heat exchangers have been conducted during the past decay using commercial CFD programs.Afsharpanah et al. [9] numerically examined a cuboid-shaped ice container with serpentine tubes and plates by evaluating various dimensionless parameters regarding the flow and geometric aspects.In numerical research, Hamzeh & Miansari [10] looked explored the quantity and configuration of refrigerated tubes as well as the impact of the fin dimensions on ice development.A novel ice-on-coil cold storage system with coil tube was investigated in experimental and numerical research by Mousavi Ajarostaghi et al.,[11] by proposing a new evaporator concerning better cooling process and uniform ice production.In all mentioned research, Ansys Fluent was used for the numerical simulation of the ice formation process.COMSOL Multiphysics is one such software that also can be used for this purpose.However, it is important to validate the simulation results with experimental data to ensure their accuracy.
In this paper, we will compare the experimental data obtained from lattice structure evaporators with the simulation results obtained using COMSOL Multiphysics.We will also discuss the agreement between the experimental and simulated results.

Refrigerant Loop
Figure 1 shows the sketch of the experimental apparatus.The facility consists of a propane (R290) refrigeration cycle composed of a reciprocating compressor (COM), a water-heated evaporator (EV), an air-cooled condenser (CO), and an expansion device (CT).An oil filter (OF) is located before the capillary tube.The evaporator, consisting of a copper finned-tube structure, is placed inside a tank of water kept in motion by an agitator.The temperature of the hot source is kept constant in the laboratory environment, while the cold source's temperature evolves during ice formation on the evaporator's surface.
Temperatures at each point of interest are measured with calibrated T-type thermocouples with a measurement uncertainty of ±0.2 °C.Thermocouples are located at the inlet of the capillary tube (T1), at the outlet of the capillary tube (T2), at the outlet of the evaporator (T3), at the outlet of the compressor (T4), and in front of the condenser fan (T5).Thermocouples T1, T2, T3, and T4 are placed inside the copper tube in direct contact with the refrigerant fluid.The pressure at the sides of each component is measured with absolute pressure transducers with an uncertainty of ±0.025 bar.The pressure transducers (KELLER-PAA-23 SY Ei) are placed at the inlet of the capillary tube (P1), at the outlet of the capillary tube (P2), at the outlet of the evaporator (P3) and the outlet of the compressor (P4).A Coriolimeter (Bronkhorst -MINI Cori-Flow, CF) placed at the condenser outlet measures the R290 mass flow rate with an uncertainty of ±0.2%.Finally, a digital energy meter (Socomec COUNTIS E03/E04, WM) measures the electrical power absorbed by the compressor.

Evaporator and Thermal Storage
The evaporator is made of a spiral-shaped copper tube with an 8 mm outer diameter installed on a square base with 19 cm on each side and a height of 21 cm.The heat exchanger tube has seven turns and a longitudinal copper fin that runs its entire length except for the corners.The fin has a thickness of 1 mm and a length of 24 mm.In addition, 80 copper bars with a 3 mm thickness, spanning the entire height of the heat exchanger, are positioned perpendicular to the fins.The thermal storage, a water-filled tank containing 5.4 litres of water, holds the evaporator.A 10 cm-long radial agitator is in the center.Six T-type thermocouples, arranged between the fourth and fifth turns starting from the top, are used to instrument the heat exchanger.These are positioned between two fins and two rods and installed into a 3D-printed support, as shown in figure 2. Thermocouples are calibrated and of the same type as those used in the refrigerant loop, with an uncertainty of 0.2 °C.

Figure 2.
Six thermocouples location on two different surfaces.

Experimental Procedures
The experimental procedure unfolds as follows.The storage is at 27 °C, then, the refrigeration cycle and the agitator inside the storage are turned on, starting the sensible cooling phase.When the average temperatures measured by the thermocouples are close to 0°C, the agitator is switched off, and the ice formation begins.The test ends when 4 kg of ice has been formed.The water level and the energy balance on the evaporator assess the amount of ice.Three experiments were conducted one day apart making sure that the same environmental temperature of 22 °C was maintained

Numerical simulation
The most suitable method for simulating our case is the conjugated heat transfer approach, as it accurately represents the physics involved.This approach considers the heat transfer occurring at the interface between a copper tube carrying the heat transfer fluid (refrigerant) in motion and the stationary storage fluid (water), which undergoes a phase change in both processes.As shown in figure 3a, the study was carried out using a 3D geometry, and transient analysis under conduction and convection heat transfer was used to account for the time dependence of the problem as well as for the heat transfer and fluid flow in the entire system.The propane temperature was determined based on the pressure that was observed throughout the test.

Grid Mesh
In accordance with the physics and geometry concerns, the computational grid was created using COMSOL Meshing (figure 3b).A sensitivity analysis was done to choose the best grid while taking accuracy and computing cost into account.To do this, the simulation was run with varying cell sizes coarser (1.9 million cells), coarse (3.4 million cells), normal (6.8 million cells), and fine (9.6 million cells).By comparing the temperature variations at five different times, the findings in the normal and fine grid meshes can be shown to have a small difference.As a result, all simulations were done using a normal grid mesh.

Governing Equations
The buoyancy effects are taken into account using the Boussinesq approximation: The continuity: ∇.  ⃗ = 0 (2) The momentum equation: The energy equation: Where ℎ  and ℎ  are sensible and latent enthalpy.In these equations, ,   and  are calculated based on the properties of two phases: Where,  1→2 is the latent heat energy, and  is the thermal diffusivity.

Boundary Conditions, Initial Conditions
For boundary conditions, the following presumptions are made: -The outer, top, and bottom boundaries are considered as perfectly insulated surfaces.
-Due to geometry symmetry, one-eighth of geometry was simulated, and two surfaces were considered as symmetric.-No-slip conditions was considered on all walls.
-The temperature is based on the measured pressure and mass flow rate inside the evaporation tube used for the convective heat flux boundary condition for the inner surface of the tube.The temperature difference between tube rows was neglected.
The initial condition for simulations is a temperature field of 27 °C for all domains.All the properties of the material were used in simulations considered as temperature variants.Some parameters related to the simulation are presented in Table 2.

Solver Setting
Conjugated heat transfer, turbulent, and laminar flow models are employed as the physics for the simulation.With a total duration of 3500 seconds and a time step of 1 second, a time-dependent approach is used.Three variables, pressure (P), velocity (V), and temperature (T), are used to solve the physics.The solver employed by COMSOL is coupled direct solver PARDISO.

Results and Discussion
The simulation is implemented by dividing it into two stages (0-1250 and 1251-3500).The first stage involves modeling the turbulent flow physics interface, considering the experimental conditions with a working impeller.In the second stage, the simulation focuses on the laminar flow to replicate the conditions without the impeller.The results of the numerical modeling of the ice-formation process in the geometry under the study and its comparison with the experimental data are presented in this section.Figure 4a shows the comparison between numerical and experimental data for the six temperatures along the process and the absolute error between experimental and numerical results represented in figure 4b.The maximum absolute error for six thermocouples is below two degrees.It is clear that the simulation results are in good agreement with the experimental data.The velocity profile and streamlines at surface 1 are presented in figure 5a.The temperature contours for the surface one for different times have been illustrated in figure 5b.The observation reveals that as time passes, a growing region surrounding the heat exchanger attains the temperature required for solidification.Subsequently, this area solidifies, and the temperatures within these regions even fall below the phase change temperature.To examine how ice is distributed within the storage, one can utilize the contours of the liquid fraction.Figure 5c shows phase fraction contours at different times.These findings are in excellent accord with the outlines of the temperature shown in figure 5b.

Conclusion
This study used a 3D numerical simulation to examine the ability of COMSOL Multiphysics in simulating the ice formation process in a find-and-tube heat exchanger.The numerical results were validated with experimental results of an evaporator part of a refrigeration cycle.The findings of this study demonstrate that the various physical phenomena involved in water flow, heat transfer through conduction and convection, and phase change can be effectively simulated using COMSOL Multiphysics software.Future endeavors will involve investigating the complete refrigeration cycle, encompassing ice formation and ice melting using COMSOL Multiphysics and performing an experimental validation of the numerical results.

Figure 1 .
Figure 1.Schematic of the experimental setup

Figure 4 .
Figure 4. Comparison between experimental and numerical results for six thermocouples, a) temperature, b) absolute errors.

Figure 5 .
Figure 5. Numerical results for, a) velocity profile, b) temperature profiles at different times, and c) liquid fraction at different times (blue as liquid and white as solid).

[ 1 ]
Jannesari H and Abdollahi N 2017 Experimental and numerical study of thin ring and annular fin effects on improving the ice formation in ice-on-coil thermal storage systems Appl.Energy.189 369-384 [2] Mađerić D Čarija Z Pavković B and Delač B 2021 Experimental and numerical study on water ice forming on pipe columns in a limited-volume storage Appl.Therm.Eng.194 117080 [3] Li Y Yan Z Yang C Guo B Yuan H Zhao J and Mei N 2017 Study of a Coil Heat Exchanger with an Ice Storage System Energies 10 1982 [4] Yang T Sun Q and Wennersten R 2018 The impact of refrigerant inlet temperature on the ice storage process in an ice-on-coil storage plate Energy.Procedia.145 82-87 [5] Saraceno L Boccardi G Celata G P Lazzarini R Trinchieri R 2011 Development of two heat transfer correlations for a scraped surface heat exchanger in an ice-cream machine Appl.Therm.Eng. 31 4106-4112 [6] Wang Baolong, Li Xianting, Zhang Maoyong, Yang Xudong 2011 Experimental Investigation of Discharge Performance and Temperature Distribution of an External Melt Ice-on-Coil Ice Storage Tank HVAC&R Res. 9 291-308 [7] Ezan M A Erek A 2012 Solidification and Melting Periods of an Ice-on-Coil Latent Heat Thermal Energy Storage System J. Heat.Transf.134 062301 [8] López-Navarro A Biosca-Taronger J Torregrosa-Jaime B Martínez-Galván I Corberán J M Esteban-Matías J C Payá J 2013 Experimental investigation of the temperatures and performance of a commercial ice-storage tank Int.J. Refrig.36 1310-1318 [9] Mousavi Ajarostaghi S S Poncet S Sedighi K Amiri L 2023 Solidification analysis in an ice-on-coil ice storage system: Experimental and numerical approaches J. Energy Storage 65 107291 [10] Afsharpanah F Pakzad K Mousavi Ajarostaghi S S and Arıcı M 2022 Assessment of the charging performance in a cold thermal energy storage container with two rows of serpentine tubes and extended surfaces J. Energy.Storage.51 104464 [11] Hamzeh H A and Miansari M 2020 Numerical study of tube arrangement and fin effects on improving the ice formation in ice-on-coil thermal storage systems Int.Commun.Heat.Mass.Transf.113 104520

Table 1 .
Table 1 summarizes the main characteristics and accuracy of the measuring devices.Characteristics and accuracy of the measurement elements.