CETvol87 DOI: 10.3303/CET2187036 Paper Received: 7 October 2020; Revised: 21 February 2021; Accepted: 8 April 2021 Please cite this article as: Bassani A., Duserm Garrido G., Giuberti G., Dordoni R., Spigno G., 2021, Comprehensive Mathematical Model for Freezing Time Prediction of Finite Object, Chemical Engineering Transactions, 87, 211-216 DOI:10.3303/CET2187036 CHEMICAL ENGINEERING TRANSACTIONS VOL. 87, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Laura Piazza, Mauro Moresi, Francesco Donsì Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-85-3; ISSN 2283-9216 Comprehensive Mathematical Model for Freezing Time Prediction of Finite Object Andrea Bassani*, Guillermo Duserm Garrido, Gianluca Giuberti, Roberta Dordoni, Giorgia Spigno Università Cattolica del Sacro Cuore, Department for Sustainbale Food Process (DiSTAS), Via Emilia Parmense 84, 29122 Piacenza, Italy andrea.bassani@unicatt.it Frozen food sees a continuous increase in consumption thanks to the capability to preserve the organoleptic properties and, at the same time, to increase the shelf-life of the food itself. Therefore, a proper design of the freezing process is crucial and is strictly related to an accurate evaluation of freezing time because this establishes the minimum residence time of the product in a continuous freezer. For this reason, several mathematical models have been proposed and investigated for predicting the freezing time, starting from empirical models (lower accuracy and computational time) up to computational fluid dynamics simulation (higher accuracy and computational time). An excellent compromise between accuracy and computational effort seems to be a model that combines empirical laws for property evaluation and heat diffusion equation solved in one dimension. This model can numerically be solved using the method of lines in which spatial derivatives are discretized by the finite difference method and the resulting system of ordinary differential equations is integrated using an appropriate solver. This work aimed to fill some gaps to develop a comprehensive and more accurate model for freezing time prediction. Indeed, the key idea is to validate a model that could be used to optimize the refrigeration process for energy-saving and be the base for a design of experiment in case of lack of experimental data. The newly developed model includes the evaluation of freezing time for finite food shapes because, in some cases, it is not possible to assume a characteristic direction for heat flux. The property of unfrozen and frozen food (e.g. density, thermal conductivity, apparent specific heat, etc.) are evaluated basing on the principal constituent of food (e.g. water, fiber, protein, etc.), while heat transfer coefficient is evaluated using empirical equations, depending on adimensional numbers. In this way, it is possible to be flexible and not strictly related to the evaluation of the properties of food for each different case. The proposed model was validated with different sets of experimental data related to beef and strawberry pulp freezing. For these cases, the R2 index is equal to 0.982 and 0.996 respectively, showing a good agreement between experimental and predicted data. Finally, a case study of spinach cubes freezing was provided to show the potentiality of the novel model. 1. Introduction Nowadays frozen products are experiencing a continuous increase in consumption, both in Europe and in Italy, up to an annual rate of 2.7% for a market value of over 100 billion euros. IIAS (Italian institute of frozen food) reports a national consumption of 849,900 tons in 2019 showing a plus of 1.3% compared to 2018. This corresponds to an annual individual consumption of more than 14 kg, for a total market value estimated around 4,400/4,700 million euros. Vegetables are the leading foods of the sector along with other products such as seafood, recipes, desserts, fruit. In 2019, vegetables represented more than 43% of the retail sector and both simple vegetables such as peas, followed by spinach, green beans, and potatoes (IIAS, 2019). Indeed, the frozen foods are products capable of maintaining the nutritional value and organoleptic properties over time, thus reducing, as a secondary effect, the food waste resulting from the deterioration of the product itself (Erickson and Hung, 2012). For these reasons, the cold chain and its management assume a critical role for the maintenance of process efficiency and product quality during the processing, freezing, and preservation steps, up to the final product consumption. In this regard, different aims can be pointed out for 211 proper management of the cold chain like the increasing of the sustainability of production, the improvement of the quality of the frozen product, and the optimization of process energy consumption (Taoukis et al., 2016). In order to meet some of these aims, recently different research projects, both at Europen and national level, were financed. For instance, the project CoACh (Cold management in Agro-food Chains: solutions for process digitalization), aims of develop a prototype of an integrated system for the optimized management of the cold chain in a process of transformation of plant products and its validation in an industrial environment. In this framework, the aim of this work is to develop and validate a suitable mathematical model for the freezing process. This mathematical model will be used inside the CoACh project with optimization algorithms to manage the production of air used as cooling fluid most efficiently. In particular, the model will be able to predict the effects of variations on the air produced, for example in terms of relative humidity or temperature, on the product itself, always ensuring that the freezing characteristics are maintained. Moreover, this model can be used for the optimization of the freezing time reducing energy consumption. Therefore, since the model, based on mass and energy balances, has to be included in an optimization process, must require a low computational effort, while maintaining good accuracy in the evaluation of key parameters. In the following paragraphs, several validated mathematical models for the simulation of the refrigeration process will be briefly presented and useful laws for the estimation of the properties involved will be illustrated. This review process was useful to select a suitable model for the future energy optimization of the process foreseen by the CoACh project. However, this model has some shortcomings, illustrated below, that will be addressed in this work. Therefore, several simulations for the new model validation are needed and will be presented. Finally, the potential advantages of the new model will be briefly discussed. 2. Materials and Methods As previously mentioned, in this section are briefly presented the different freezing mathematical models investigated with regards to the specifics required by the CoACh project (i.e. high accuracy and low computational effort). The selected one is deeply described showing the involved equations. Then, the food and freezing properties evaluation were presented. 2.1 Freezing Model In literature, several mathematical models have been investigated and validated to simulate the refrigeration process (Zhao and Takhar, 2017) and can be summarized in three categories. The first one included the models that are based on empirical formulas for the evaluation of both the freezing time and the different associated properties such as ice fraction, bound water fraction, thermal property of food (e.g. specific heat, density, etc.), and the convective heat transfer coefficient (Pham, 1986). These models have a computational effort near to zero, but it is possible to evaluate only one key parameter: the freezing time. On the contrary, the second category presents models that solve the equation of heat diffusion in 3-D using computational fluid dynamics (CFD). This allows to calculate, through the resolution of partial differential equations, the related properties such as the convective heat transfer coefficient or the initial freezing temperature. However, these models require high computational times, making them unsuitable for process optimization (Cevoli et al. 2018). The last kind of models exploits some of the peculiarities of each of the two previous categories. In particular, empirical formulas are used for the evaluation of properties (e.g, food thermal property, ice fraction, heat convective coefficient), while the equation of heat diffusion is solved in one dimension, allowing the evaluation of the temperature trend inside the food (Ferreira and Rojas, 2019). These models are reasonably accurate and, at the same time, have a low computational effort, making them suitable for optimization purposes. In particular, the partial differential equation was solved using the method of lines (MOL), which consists in the discretization of the spatial derivative and solve the reaming ordinary differential equation system with a proper solver. However, these models present some shortcomings, such as the impossibility to evaluate the temperature variation for food treated with finite shape, even if simple, such as a parallelepiped or a cylinder (e.g grilled eggplants or French fries). For this reason, in this work a novel model that combines the one proposed by Ferreira and Rojas (2019) and the EHTD (Equivalent Heat Transfer Dimensionality) shape factor approach (Pham, 2014) is proposed. This led to the following equation to be discretized and solved: = (1) Where and are time and spatial coordinates, while , , are density, thermal conductivity and apparent specific heat of food evaluated at temperature . The latter is a specific heat that includes the latent heat effect (Pham, 2014). is the shape factor, which allows to take into account different food shapes like 212 cylinder, rectangular brick and also ellipsoid for irregular shape and can be evaluated as reported by Pham (2014). Equation (1) can be discretized as: ( ) = ( ) ( − 1) − 2 ( ) + ( + 1)( ) ( )(∆ ) (2) Where indicates the number of nodes which represent the spatial discretization. For sure, a higher number of nodes means higher accuracy but also higher computational time. For these reasons it is necessary to select an optimal value, that, for this study, is equal to 20. The related boundary and initial conditions are: = − ≤ ≤ + = 0 (3) − ( ) = 0 = 0 ≥ 0 (4) − ( ) = ℎ( − ) = ≥ 0 (5) Where is the distance between the center of the food and the surface, is the air or coolant temperature, and ℎ is the convective heat transfer coefficient. For the reasons presented above, this model is selected to be investigated and validated for future process energy optimization. It is important to underline that equation (1) was developed starting from the case of infinite plate shape but can be easily translated also for cylinder and sphere starting from the equation presented by Ferreira et al. (2016) and Ferreira (2017). The resulting ordinary differential equation system (ODE) is solved using Matlab 2019b, which is a multi-paradigm numerical computing environment and fourth- generation programming language developed by MathWorks. 2.2 Food and Freezing Properties Once the model was selected, another key aspect was the evaluation of the food properties (e.g., specific heat, conductivity, density, etc.). These properties depend strongly on the chemical composition of the food and its temperature. Therefore, it is difficult to produce an experimental database for all possible food conditions and compositions. For this reason, in this work different types of food were defined basing on their main constituents that are water, protein, fat, carbohydrates, fiber and ash. From these constituents, it is possible to predict, with certain reliability, the thermal properties of the food as a function of temperature (Becker and Fricke, 2004). It is certainly limiting to define food as a mixture of its main constituents, but this allows to be flexible to energetically optimize the process and to potentially reduce freezing times. The definition of the food according to its constituents can also allow a rapid characterization of the food itself through the use of widely validated analytical methods (Deng et al., 2017). In addition, it could be possible to adapt the freezing times according to the different batches of the same food to be frozen and not only according to the different types of food. In other words, it could allow the possibility to study the effects of the different compositions of the same food on the final freezing time and on the energy consumed. Regarding the other property, the following equations were selected for evaluating ice fraction ( ), bounded water ( ), initial freezing temperature ( ) and apparent specific heat ( ) (Pham, 2014): = ( − ) 1 − (6) = −4.66 − 46.4 (7) = 0.45 + 0.3 (8) ≤ = + + − ( − ) − (9) Where is the mass fraction, the subscripts are related to the main constituent (e.g. is water) and subscript 0 means at initial conditions. The temperatures are in Celsius degree. 213 Finally, the convective heat transfer coefficient can be evaluated form the Nusselt number using the following relation for laminar and turbulent flows (Zigunov et al. 2018) = 0.664( ) . ( ) / < 10 (10) = 0.664( ) . ( ) / + 0.037( ) . ( ) 1 − . ≥ 10 = 10 (11) 3. Results and discussions The novel model, as reported in section 2.1, was developed starting from the model presented by Ferreira (Ferreira, 2017) including the shape factor to simulate the freezing process even for food that must be represented with finite geometric shapes. For this reason, it is necessary to validate the model basing on different sets of experimental data used also for validating the previous EHTD shape factor model (Pham, 2014). It is useful to point out that the data used for validating the model proposed by Ferreira was not used both because is based on tylose gel, that has a fixed properties, that cannot be evaluated starting from their main constituents, and because the shape of tylose gel presents a main direction for heat diffusion (Ferreira, 2019). De Michelis and Calvelo (1983) measured different freezing times in blocks of beef frozen in molds having different shapes and sizes. Briefly, the molds were made of aluminum or acrylic of different thicknesses and they are immersed in cold air with controlled temperature. The block to be frozen was supported by a structure with air blowers, which forced a highly turbulent stream on each face. The results of model simulations are reported in Figure 1 in which have been further compared both with the predicted freezing time using the Pham method (Pham, 1986) and using the model presented directly by De Michelis and Cavelo (1983), that is based on a shape factor approach. It is useful to underline that, as reported, the new model estimates food properties starting from its initial composition, defined in terms of its main constituents (i.e. water, proteins, fats, carbohydrates, fibers, ashes). However, since the composition of beef is not available for the simulations, it has been estimated starting from a typical composition of beef (Pham et al., 1994) as reported in Table 1. Table 1: Composition of foods investigated in this work Beef – De Michelis and Cavelo (1983) Strawberry pulp - Salvadori et al. (1996) Spinach Water 0.6340 0.8975 0.9010 Protein 0.3064 0.0265 0.0350 Fat 0.0528 0.0063 0.0070 Carbohydrate 0.0015 0.0402 0.0290 Fiber 0.0000 0.0198 0.0190 Ash 0.0053 0.0097 0.0090 Figure 1 shows a good agreement between experimental and simulated data, confirmed by the evaluation of R2 coefficient that is equal to 0.9816, 0.9821, and 0.9953 for the new model, Pham model and De Michelis model respectively. The slight differences could be explained by the assumption of the food composition. The new model was further validated with another set of experimental data, provided by Salvadori et al. (1996), who studied the freezing of strawberry pulp. Likewise, it is necessary to estimate the initial composition of the strawberry pulp because only the water content was provided. The composition is reported in Table 1. Briefly, the tests were performed using a container (0.145x0.500x0.520 m) with three different configurations. The first (P1) with all its later insulated with polystyrene, the second (P2) with only the lateral of 0.520 m insulated and the third (P3) with no insulation. Table 2 reports the operating conditions for each configuration and the results of the simulation. The latter show a good agreement considering the assumption of the initial composition of the strawberry pulp and considering the intrinsic error of the evaluation of its properties starting from its constituents (Becker and Fricke, 2004). After validation, a first simulation was done in order to show the potentiality of this new approach taking as reference food the “Spinacio Cubello”, which is frozen spinach with a cube shape. The latter is kindly provided by Orogel, which is an Italian industry specialized in the production of frozen foods and a partner of the CoACh project. The spinach composition, reported in Table 2, is assumed based on the property reported by Cevoli et al. 2018 and on the CREA database (2020). The principal dimension of the spinach cube is equal to 0.04 m. The operating conditions of the refrigeration process have been initially hypothesized and in particular, the air temperature and the initial temperature of the spinach have been assumed equal to -40°C and 20°C respectively. 214 Figure 1: Scatter diagram – Experimental and predicted freezing time from De Michelis and Calvelo (1983) Table 2: Experimental conditions and freezing time (Salvadori et al., 1996) and predicted freezing time Tinlet [°C] Tout [°C] Tair [°C] Freezing time – experimental [hr] Freezing time – Predicted [hr] P1 26.30 -18.00 -35.00 37.72 41.09 P2 26.30 -18.00 -35.00 33.10 35.64 P3 26.30 -18.00 -35.00 28.30 31.00 Refrigerant air has been assumed in a turbulent flow. The results obtained from the simulation are reported in Figure 2. In particular, are shown the variation of temperature as a function of time in the geometric center of spinach cub and on the surface (Figure 2-A) and as a function of the distance from the geometric center at a fixed instant of time (Figure 2-B). Figure 2: Temperature variation inside "Spinacio Cubello" as function of time (A) and distance from the center (B) The key idea is that, with this model, it could be possible to evaluate the properties of food at different times and positions with reasonable computing efforts. For instance, the enthalpy of food could be predicted in order to evaluate the amount of energy need for the freezing process and, in this way, optimize the process itself from an energy point of view. For sure, this model cannot give the same in-deep information if compared to CFD models. Indeed, the latter can predict, for example, the formation of hot and cold spots that could lead to the formation of ice crystals of different sizes directly influencing the final quality of the frozen food (Ullah et al., 2014). However, as already mentioned, the model is suitable to perform an energy optimization basing on the operating conditions (e.g., refrigerant temperature or final food temperature) and then perform a CFD simulation at optimal conditions in order to have a complete and detailed picture of the freezing process. 215 4. Conclusions In this work, a comprehensive mathematical model for the freezing process was developed and validated starting from the one proposed by Ferreria and Rojas (2019). Indeed, the latter presented some gaps like not considering the finite shapes of food. The developed model is efficient, flexible and, at the same time, accurate if the food considered can be assimilated to a defined geometric figure (e.g., cylinder sphere or parallelepiped). This model can be used in an effort to reduce the freezing time of the food and, in the framework of the CoACh project, to energetically optimize the process or the production of the cooling fluid. In this regard, it is interesting to highlight how the freezing process can be included in energy optimization management with the next steps, such as transport and storage. From this point of view, perhaps preferring not to modify the process conditions during the freezing phase, the frozen product can also be used as an energy accumulator to be exploited in the following phases, always guaranteeing the quality of the food. Although there are several studies in this regard, this type of approach can still be improved both in terms of empirical properties and in terms of the geometric shapes allowed to define a food. For instance, it could be possible to include some empirical laws to evaluate the formation and the dimensions of ice crystal (Kiani and Sun, (2011)) which is a crucial parameter to determine the quality of frozen food. Acknowledgments The research is part of the CoACh project funded by the European Funds of the Emilia-Romagna Region - POR FESR 2014-2020. References Becker B. R., Fricke B. A., 2004, Heat transfer coefficients for forced-air cooling and freezing of selected foods, International Journal of Refrigeration, 27(5), 540-551. Cevoli C., Fabbri A., Tylewicz U., Rocculi P., 2018, Finite element model to study the thawing of packed frozen vegetables as influenced by working environment temperature, Biosystems Engineering, 170, 1-11. CREA, 2020, Consiglio per la ricerca in agricoltura e l’analisi dell’economia agraria, , accessed 13.12.2020. De Michelis A., Calvelo, A., 1983, Freezing time predictions for brick and cylindrical‐shaped foods, Journal of Food Science, 48(3), 909-913. Deng Q., Penner M. H., Zhao Y., 2011, Chemical composition of dietary fiber and polyphenols of five different varieties of wine grape pomace skins, Food Research International, 44(9), 2712-2720. Erickson M. C., Hung Y. C, 2012, Quality in frozen food, Springer Science & Business Media, London, UK. Ferreira S. R., Rojas L. O. A., Souza D. F. S., Oliveira J. A, 2016, Freezing time of an infinite cylinder and sphere using the method of lines, International Journal of Refrigeration, 68, 37-49. Ferreira S. R., 2017, Freezing time of a slab using the method of lines, International Journal of Refrigeration, 75, 77-94. Ferreira S. R., Rojas L. O. A, 2019, Freezing times using time derivative of temperature on surface of foods, International Journal of Refrigeration, 98, 436-443. IIAS, 2019, Istituto Italiano Alimenti Surgelati, , accessed 13.12.2020 Kiani H., Sun D. W., 2011, Water crystallization and its importance to freezing of foods: A review, Trends in Food Science & Technology, 22(8), 407-426. Pham Q. T., 1986, Simplified equation for predicting the freezing time of foodstuffs, International Journal of Food Science & Technology, 21(2), 209-219. Pham Q. T., Wee H. K., Kemp R. M., Lindsay D. T., 1994, Determination of the enthalpy of foods by an adiabatic calorimeter, Journal of Food Engineering, 21(2), 137-156. Pham Q. T., 2014, Food freezing and thawing calculations, Springer, London UK. Taoukis P. S., Gogou E., Tsironi T., Giannoglou M., Dermesonlouoglou E., Katsaros G., 2016, Food cold chain management and optimization, In Emerging and Traditional Technologies for Safe, Healthy and Quality Food, (pp. 285-309), Springer, Cham. Ullah J., Takhar P. S., Sablani S. S., 2014, Effect of temperature fluctuations on ice-crystal growth in frozen potatoes during storage, LWT-Food Science and Technology, 59(2), 1186-1190. Zhao Y., Takhar P. S., 2017, Freezing of foods: mathematical and experimental aspects, Food engineering reviews, 9(1), 1-12. Zigunov F., Zinani F., da Silva J., 2018, Numerical simulations of food freezing times using thermophysical food property models, Estudos Tecnológicos em Engenharia, 12(1), 1-13. 216