MODELING AND TESTING OF AN ICE BANK FOR MILK COOLING AFTER MILKING

Based on mass and energy balance equations and data collected in specialized bibliographies, it was proposed a methodology for modeling of an ice bank for thermal storage of energy at low temperature to be used in milk cooling. The data obtained by the equations arrangement were used in the design and installation of a real system, which was tested to collect data and to compare it with the estimated. The data obtained in terms of prediction of ice formation, temperature of the ice bank, water and capacity of the selected refrigeration system prove the adequacy of the equations arrangement, showing that it can be used as a tool for sizing solid ice bank in smooth pipe coils, with good approximation.


INTRODUCTION
Refrigeration equipment are getting more popular in the agriculture area because reduction in temperature is an efficient method to preserve perishable products (Senthilkumar et al., 2015), keeping its characteristics of appearance and taste for a longer period, which allows these products to be delivered to the big consumer centers, far from the production places (Rezaeinejad & Nazarian, 2013).
The cold reduces the microbial and enzymatic activities of a product. Therefore, the faster a product is cooled after harvesting, the better is its quality and thus its shelf life (França et al., 2015;Siqueira et al., 2014). This is the case of milk, which should be cooled immediately after milking in a period of less than two hours (Taffarel et al., 2013;Porcionato et al., 2008).
Fast cooling means reducing the temperature of a product in a small period. It means that a certain amount of thermal load must be taken out of the product/place, which means that a powerful equipment is required when compared with refrigeration equipment used only to maintain the temperature (Neves Filho, 2001).
In this case, ice accumulation is an interesting option because it allows storing energy in ice banks while thermal demand is low. This thermal energy is to be used during periods of peak demand (Njoku et al., 2014;Oró et al., 2013). With ice bank, it is possible to distribute the electric demand of energy during the day, using the low thermal load demand to produce ice to be used when the thermal load demand is high. In this case, a lower capacity equipment is required (Yi & Dong, 2015). The conventional refrigeration system (direct expansion of refrigerant fluid) requires a powerful equipment in this stage to remove the thermal load instantaneously. If the milk storage tank is undersized, there will be a burst of starts and stops of the compressor, which is not recommended because it increases the energy consumption and decreases the working life of the equipment (Rutberg et al., 2013). In this case, there is higher expenditures for maintenance and repairs.
Comparing thermal storage with cold water, ice bank is advantageous because it requires smaller storage volumes (Ismail et al., 2016a) because of its capacity of storing energy in the latent form (Shinde & Suresh, 2014), which is eighty times higher than the liquid cold water heat capacity.
Despite of largely use in fast cooling processes in the agro industrial area, it was not found Brazilian bibliography about a specific method to size an ice bank. The few bibliographies found are not clear and are based on complex equations with multivariable and no ready solution. It makes difficult to apply these equations in the engineering area.
Despite of what it seems, it is not easy to size an ice bank. The heat transfer process between the coil and the cold water happens in unsteady-state, where there is variation in the ice thickness during the thermal storage process (Ismail et al., 2016b). Therefore, it is not only a matter of getting the product of cooling capacity by the latent heat of fusion of ice.
This study was written because of the difficulty in finding data for the design of an ice bank, in order to provide input to future designers and researchers that face the same dilemma. It shows a simplified model for sizing ice banks, comparing model data with data from a real system built based on it.

MATERIAL AND METHODS
The research was carried out in the Thermodynamic and Energy Laboratory of the Faculty of Agricultural Engineering of the State University of Campinas -Unicamp. It was divided in two parts: the first consisted in dimensioning the ice bank and the second consisted in building the cooling system (ice bank) and testing itto obtain data and its comparison with the estimated by the arrangement of the equations (model).
The purpose of the ice bank is the production of thermal capacitance for fast cooling of milk after milking on dairy farms. For example, the thermal cooling load was based on 600 liters of milk per cycle at 30 ° C, resulting in about 74,400 kJ, which were accumulated in the gaps between each milking, in periods of 10 hours.
For the ice bank thermal reservoir, an existing tank in the laboratory was employed. The tank was made of timber and insulated with expanded polystyrene 50 mm thickness, with internal dimensions of 1.00 x 0.90 x 0.80 m. Considering the calculated thermal load and the stored energy in the liquid water surrounding ice, a mass of 186 kg of ice was determined to cool the 600 liters of milk.
To size the evaporator (coil where the ice is accumulated) it was taken in account the ice forming process. Modeling based on energy balance to the ice accumulator tank was applied to size the evaporator (Ismail et al., 2014). Considering the energy flow ( Figure  1): the infiltration heat rate due to heat transfer from environment to the tank ( I Q • ), the heat rate from water to the ice surface ( a Q • ) and the heat rate removed from the ice surface by the refrigerant flowing through the evaporator coil (  For scaling purpose, it was considered that ice formation occurs uniformly all around the evaporator coil tubes. In the real system, it can not occur by heat exchange phenomena. Thus, by adjusting the heat transfer equations to the case (Borgnakke & Sonntag, 2009), whereas the ice surface temperature is equal to 0 ° C, the rate of heat removed of the ice surface by the refrigerant fluid was calculated by [eq.(1)].
) .( . The global heat transfer coefficient between the ice surface and the refrigerant fluid is proportional to the inverse of the sum of the thermal resistances to heat transfer (ice resistance, pipe wall resistance and the internal resistance of the tube). It was calculated by [eq.(2)]. The coefficients adopted in [eq.(2)], suggested by Neves Filho (2001), were: -heat transfer coefficient between the refrigerant and the inner wall of the coil i h =1740W m -2 ºC -1 (forced convection); -Thermal conductivity of the pipe (copper) t The heat transfer rate between the water and the ice surface was calculated by the [eq.(3)], considering the surface temperature of the ice equal to 0 °C.
Where, a Q • -heat transfer rate between the water and ice (kW); e h -external heat transfer coefficient between the water and the ice surface (kW m -2 ºC -1 ), a T -water temperature from the ice thermal storage tank (ºC).
According to Neves Filho (2001), considering that there was no agitation in the tank during the ice formation, the value of 81.2 W m -2 ºC -1 for natural convection was adopted for the external heat transfer coefficient between the water and the ice surface.
Considering that during freezing the heat flow ( causes the ice formation, the amount of ice accumulated in the time interval dt can be calculated using [eq.(4)]. The term ( 2 / dD ) represents the ice layer formed over the time interval dt .
In case the of the water's temperature during the ice formation process, considering the balance between the heat removed from the water by the ice surface and the heat from the external environment transmitted through the tank's walls (infiltration heat), the variation of the water temperature in a time interval dt was calculated by [eq. (5)].
dT is positive and the water temperature will increase. Otherwise, a dT is negative and the water temperature will decrease, which happens when the refrigeration system is operating and accumulating ice.
The infiltration heat rate was obtained by the [eq. (6) The global heat transfer coefficient between the inner wall of the tank and the environment was calculated using [eq. (7)]. Considering that the ice bank should operate in continuous regime, with successive cycles of ice formation and use of thermal load, it was established that the initial temperature of the water in the thermal storage process would be 3.5 °C in order that the end temperature of the cooled milk is 4 °C. For the temperature of the external environment, it adopted the value of 24 °C, which is the average temperature of dry bulb for Campinas / SP region.
Adding the eqs (1) and (3) in the corresponding terms in [eq.(4)] and applying a numerical solution with small finite differences (Ismail et al., 2015;Ismail et al., 2016a;Ismail et al., 2016b, Kandula, 2011, the differential time to occur a variation of ice diameter considering constant D during this process is given by: Jaboticabal, v.38, n.4, p.510-517, jul./ago. 2018 In the same time differential, replacing the corresponding terms in [eq. (5)] by eqs (6) and (3), the change in water temperature applying numerical solution is determined by: The ice mass formed in the t  interval is given by the following equation: It was adopted to work with D  variation of 0.001 m. With constant D , it was determined a t  time interval and the water temperature variation to form the first ice layer corresponding to D  . In the next interval, ice diameter, water temperature and ice and water mass were obtained by the following expressions: And so on for each subsequent interval, the new values were calculated.
It was considered that in the beginning there was no ice in the freezing process. Therefore, To simplify calculation, the equations used were inserted in an Excel spreadsheet. The design was based on the determination of the total length of the coil to the accumulation of ice mass ice M within the stipulated period of 10 hours. The length and diameter of the tube were inserted as input variables. Once the diameter was selected, assigning values for tube length, output data were obtained: the mass of ice formed and the total time for its formation. The pipe length for which the pre-established conditions (mass of accumulated ice and time for accumulation) were met was adopted.
As suggested by Neves Filho (2001), as the evaporator is a heat exchanger that uses refrigerant and water, there must be a temperature difference of at least 5 °C between the refrigerant and the water. In this case, since the ice is at 0 °C, the evaporation temperature adopted was -5 °C. To simplify the calculation, this assumption was considered throughout the whole process.
Based on the data obtained by the calculation, the sizing and assembling of the ice coil (evaporator) were done as well as the selection of refrigeration components (compressor, condenser and expansion valve) to build the ice bank.
According to the results obtained by the model, considering the diameter of the selected copper tube (15.9 mm), the length of the evaporator was 45 m. To minimize pressure loss and excessive refrigerant superheat at the evaporator outlet, resulting in increased formation of ice next to the expansion valve, the length was dividedin6 horizontal coils of 7.5 meters each, arranged one on another and supported by a metal base (Figure 2). To avoid coupling between ice cylinders, which would affect the heat transfer, the tubes have 100 mm of space between each other, considering that the model accused a final diameter of 84.9 mm of ice.
For the length of the evaporator according to the modeling, if running the equipment 9.98 hours the produced ice mass is 219 kg. For the 186 kg of ice estimated to cool the milk, the time predicted by the model to run the process is 8.24 hours, which increases the downtime to the equipment. Except the evaporator, which was designed and assembled, the other components were selected by the technical catalogs, considering the data provided by the model and the previously established operating conditions.
The compressor selected was a Block III model from Bitzer with R22 refrigerant fluid and refrigeration capacity varying between 4.63 and 2.47 kW for evaporating temperature between -5 and -10 °C and condensation temperature varying from 30 to 50°C.
Considering the cooling capacity and the compressor power, to be able to reject the sum of the heat removed from the evaporator and the compression work, it was determined that the capacity of the condenser should be 5.30 kW. Therefore, a brazed plate heat exchanger model WP 4 20 from Apema was selected.
To meet the cooling capacity of the compressor and the evaporator, a thermostatic expansion valve with external equalization model TEX2-1.0 from Danfoss was selected. The valve is able to provide 4.8 kW of cooling capacity in the initial conditions.
Once the ice bank was designed for dairy farms, where there is often interruption in electricity supply, an alternative to activate the refrigeration equipment was considered. Thus, the system is constructed with an electric motor and a combustion engine. The combustion engine received a conversion kit for gas, which allows to use biogas produced from the manure of the stabled animals in the facilities. Besides, the system was designed to use the rejected heat from the refrigeration system to heat water used to clean the tanks and milking equipment. No coincidence a water condenser was chosen.
After the acquisition of the remaining components of the refrigeration system (compressor, condenser, expansion valve and control valves), the refrigeration equipment was put together and instrumented for tests aiming confirmation of the model parameters, as operation time, ice diameter, ice mass accumulated and the final temperature of the ice bank's water.
For temperature measurement, thermocouples "T" type were installed in the inlet and outlet of each refrigeration component as well as in the water from the ice bank. To measure the suction and discharge pressures, Wika pressure transmitters with a range of 0 to 40 bar were employed. For the mass flow measurement, a Coriolis flow meter from Danfoss was used, with a working range from 0 to 1000 kg h -1 . Everything was connected to an acquiring and recording data system consisting of a board A/D CAD 12/32-32 and a MCS-1000 signal conditioner, both from Lynx brand.
During the tests, the ice diameter was measured using a digital caliper, in 10-minute intervals. Measurements were taken in the center of the evaporator, in the first layer.
The mass of accumulated ice during the tests was obtained multiplying the ice density by the volume of ice calculated according to the length of the ice serpentine tube, the tube diameter and the ice diameter measured by the caliper (Equation 17). Using the temperature and pressure measured data and thermodynamic tables for R22, the enthalpy of the system was applied in Equations 18 and 19 to determine respectively the heat removal rate from the evaporator (ice bank) and the rejected heat in the evaporator by the R22 refrigerant. Tests with the combustion engine were carried out using compressed natural gas.
The diameter of the ice formed in the ice bank was the main variable calculated by the "model", where its relative mean error (P) was verified by [eq.(20)].
Where, no-number of experimental observations; Y -experimental value observed; and Ŷvalue estimated by the model.

RESULTS AND DISCUSSION
In Figure 3, it is shown the ice cylinders formed around the evaporator coil after a few minutes of the completion of a test and depletion of the water tank. As also observed by Habeebullah (2007), there is a variation in the ice diameter formed over the coil, which follows the direction of the coolant flow, with the largest diameter at the coil inlet, where the biphasic mixture contains more refrigerant in liquid state and, therefore, the capacity to remove heat is larger. In the coil outlet, with the saturation in the capacity of removing heat from the refrigerant, which is already in the superheated vapor state and with a smaller heat transfer coefficient, the heat transfer is lower and consequently there is a smaller formation of ice. As the ice diameter measurements were performed in the middle of the first coil of the evaporator, which was very close to the tank surface, the data obtained reflects an average for the process. Maybe even a little less, once the coils that were below had greater formation of ice. It happened because these coils were more submerged and became more isolated from the heat infiltration from the external environment.
Analyzing the refrigerant temperature variation in the evaporator inlet (evaporation temperature) and in the evaporator outlet (Figure 4) it is possible to understand the thickness variation. The refrigerant superheat in the evaporator outlet leads to a lower global heat transfer coefficient and therefore a decrease in the heat transfer between the refrigerant and the ice surface, as observed by Ismail et al. (2016a) and Habeebullah (2007). Different from what was proposed in the model, the evaporation temperature (evaporator inlet) varied between -6 and -9 °C, remaining stable from the middle of the operating time. Figure 5 shows the curves of the observed and calculated ice diameters as functions of operating time. Following the established conditions of sizing, all tests started with the water temperature of the ice bank close to 4 °C, with small variation. Engenharia Agrícola, Jaboticabal, v.38, n.4, p.510-517, jul./ago. 2018 As in Ismail et al. (2014), Ismail et al. (2015), Ismail et al. (2016a) and Ismail et al. (2016b),the behavior of the ice diameter curve as a function of the operating time for the test approached the predicted by the model with a P value of 5.18. According to Mohapatra & Rao (2005), models with P value above 10% are unappropriated to represent its process, where P value indicate the deviation of the observed data from the curve obtained by the estimated process Mohapatra & Rao (2005), Therefore, the obtained P value indicates the equations arrangement used is appropriate to predict the ice diameter. Figure 6 shows ice forming through operation time, where it was observed a linear behavior in accordance with the behavior obtained by Habeebullah (2007). Figure 7 shows the temperature of the water in the ice bank, which had a similar behavior as observed by Habeebullah (2007). Despite variation in the beginning and through the process, the final temperature was always around 0°C as predicted by calculation.  Analyzing the refrigerant temperature variation in the evaporator inlet (evaporation temperature) and in the evaporator outlet (Figure 4) it is possible to understand the thickness variation. The refrigerant superheat in the evaporator outlet leads to a lower global heat transfer coefficient and therefore a decrease in the heat transfer between the refrigerant and the ice surface, as observed by Ismail et al. (2016a) and Habeebullah (2007).
Despite the lower correspondence for ice mass and for water temperature, the calculated data is still close to the mean observed values and within the standard deviation.
Even that it was observed in certain periods of the process a difference between observed and calculated values, in the end of the process (the end period of operation) the minimum required values (186 kg and 77 mm of diameter) were achieved, with mean values of 79.8 mm for ice diameter and 199.8 kg for ice mass.
The trials with the combustion driven motor showed the lower means of ice diameter and ice mass Engenharia Agrícola, Jaboticabal, v.38, n.4, p.510-517, jul./ago. 2018 (76.4 mm and 184.6 kg respectively). When using the electric driven motor, the means were 82.8 mm and 214.8 kg for the ice diameter and ice mass respectively, where these values are closer to the estimated. The difference was caused by the difference in the rotation due to the pulleys used in the combustion driven motor, where the rotation of the compressor was slightly lower, which affected its performance, as observed by Jordan et al. (2016).

CONCLUSIONS
The modeling used based on the arrangement of the energy balance equations and on the heat transfer showed good approximationto predict parameters related to ice forming. Therefore, it is possible to use the modeling in ice bank dimensioning.
The operational conditions affected the obtained results, and therefore the adequacy of the used modeling, as the case of the compressor rotation.
Using more realistic parameters would enhance the modeling adjust used to predict the ice formation, as the case of the evaporation temperature.