Simulações numéricas de tempo de congelamento de alimentos utilizando modelos de propriedades termofísicas de alimentos

Resumo This work aims to evaluate the performance of a numerical modeling for prediction of food freezing times by comparing numerical with experimental results. A one-dimensional heat diffusion model including temperature dependent thermophysical properties, sudden variations of thermophysical properties during phase change and surface mass transfer was approximated using the finite differences method. The thermophysical properties of foods were modeled as functions of food composition, temperature and phase, using Choi and Okos correlations (Choi & Okos, 1986). Simple geometries such as thin slab, long cylinder and sphere were modeled, and the numerical results were compared with experimental data obtained in a pilot scale freezing tunnel. Numerical simulations were performed for some selected foods, namely sausages, potatoes, hamburger and cheese, in different geometries and sizes. The boundary conditions of the freezing surfaces were of heat and mass convection. The heat transfer coefficients were taken after usual correlations for these geometries. The mass transfer modeling was done using mass convection correlations, considering a known surface wetness. It was found that mass transfer due to moisture evaporation, thermodependency of properties and an accurate estimate for the heat transfer coefficient were crucial elements for the correct prediction of freezing curves. The numerical results were only able to predict the experimental freezing curve by adjusting the theoretical value of the heat transfer coefficient by a factor, varying from 0.7 to 1.3 in most cases, with some Este trabalho tem o objetivo de avaliar a performance do uso de métodos numéricos para previsão do tempo de congelamento de produtos alimentícios através da comparação de resultados numéricos e experimentais. Um modelo de difusão de calor unidimensional que inclui propriedades termofísicas variáveis com a temperatura e grandes variações das mesmas durante a mudança de fase, bem como a transferência de massa por evaporação de água da superfície, foi construído através de um método de diferenças finitas. As propriedades termofísicas dos alimentos foram modeladas como funções da composiçao do alimento, sua temperatura e estado físico, através do emprego das correlações de Choi e Okos (Choi & Okos, 1986). Geometrias simplificadas, tais como placa plana, cilindro longo e esfera foram modeladas, e os resultados numéricos foram comparados com dados experimentais extraídos de um túnel de congelamento piloto. Simulações numéricas foram realizadas para alguns alimentos selecionados, como salsichas, batatas, hambúrguer e queijo prato, em diferentes tamanhos e geometrias. As condições de contorno empregadas nas superfícies do alimento foram de convecção de calor e massa. Os coeficientes de transferência de calor foram inicialmente estimados utilizando as correlações típicas para estas geometrias. A transferência de massa foi modelada através de correlações apropriadas, dada uma umidade superficial calculada iterativamente. Foi observado que a transferência de massa através da evaporação da água, a dependência térmica das propriedades físicas do alimento, bem como uma correta estimativa para o coeficiente de transferência de calor foram elementos cruciais para gerar previsões adequadas das curvas de congelamento dos Numerical simulations of food freezing times using thermophysical food property models Fernando Zigunov, Flávia Zinani, Janice da Silva 2 Estudos Tecnológicos em Engenharia, vol. 12, n. 1, p. 1-13, jul/dez 2018 outliers up to 2.4. This means that although the heat conduction inside the food itself seemed to generate reasonable food freezing rates, the convection coefficients produced experimentally seemed to vary wildly from the ones predicted in the theoretical relationships. Therefore, one should be very keen of the magnitude of the convection coefficient while performing predictions for a given food freezing application problem. produtos. Os resultados numéricos somente geraram predições das curvas de congelamento dos produtos após o ajuste do coeficiente de convecção através de um fator de correção, que variou de 0.7 a 1.3 na maioria dos casos, podendo chegar até 2.4. Embora a condução de calor dentro do alimento pareça ter sido modelada com taxas de resfriamento adequadas, os coeficientes de convecção produzidos experimentalmente foram muito discrepantes em alguns casos, em relação aos coeficientes obtidos através das correlações utilizadas. Portanto, é necessário um cuidado especial com os valores do coeficiente de convecção empregados durante as simulações numéricas quando previsões deste tipo forem necessárias para a aplicação em problemas de congelamento de alimentos.

This work aims to evaluate the performance of a numerical modeling for prediction of food freezing times by comparing numerical with experimental results. A one-dimensional heat diffusion model including temperature dependent thermophysical properties, sudden variations of thermophysical properties during phase change and surface mass transfer was approximated using the finite differences method. The thermophysical properties of foods were modeled as functions of food composition, temperature and phase, using Choi and Okos correlations (Choi & Okos, 1986). Simple geometries such as thin slab, long cylinder and sphere were modeled, and the numerical results were compared with experimental data obtained in a pilot scale freezing tunnel. Numerical simulations were performed for some selected foods, namely sausages, potatoes, hamburger and cheese, in different geometries and sizes. The boundary conditions of the freezing surfaces were of heat and mass convection. The heat transfer coefficients were taken after usual correlations for these geometries. The mass transfer modeling was done using mass convection correlations, considering a known surface wetness. It was found that mass transfer due to moisture evaporation, thermodependency of properties and an accurate estimate for the heat transfer coefficient were crucial elements for the correct prediction of freezing curves. The numerical results were only able to predict the experimental freezing curve by adjusting the theoretical value of the heat transfer coefficient by a factor, varying from 0.7 to 1.3 in most cases, with some Este trabalho tem o objetivo de avaliar a performance do uso de métodos numéricos para previsão do tempo de congelamento de produtos alimentícios através da comparação de resultados numéricos e experimentais. Um modelo de difusão de calor unidimensional que inclui propriedades termofísicas variáveis com a temperatura e grandes variações das mesmas durante a mudança de fase, bem como a transferência de massa por evaporação de água da superfície, foi construído através de um método de diferenças finitas. As propriedades termofísicas dos alimentos foram modeladas como funções da composiçao do alimento, sua temperatura e estado físico, através do emprego das correlações de Choi e Okos (Choi & Okos, 1986). Geometrias simplificadas, tais como placa plana, cilindro longo e esfera foram modeladas, e os resultados numéricos foram comparados com dados experimentais extraídos de um túnel de congelamento piloto. Simulações numéricas foram realizadas para alguns alimentos selecionados, como salsichas, batatas, hambúrguer e queijo prato, em diferentes tamanhos e geometrias. As condições de contorno empregadas nas superfícies do alimento foram de convecção de calor e massa. Os coeficientes de transferência de calor foram inicialmente estimados utilizando as correlações típicas para estas geometrias. A transferência de massa foi modelada através de correlações apropriadas, dada uma umidade superficial calculada iterativamente. Foi observado que a transferência de massa através da evaporação da água, a dependência térmica das propriedades físicas do alimento, bem como uma correta estimativa para o coeficiente de transferência de calor foram elementos cruciais para gerar previsões adequadas das curvas de congelamento dos outliers up to 2.4. This means that although the heat conduction inside the food itself seemed to generate reasonable food freezing rates, the convection coefficients produced experimentally seemed to vary wildly from the ones predicted in the theoretical relationships. Therefore, one should be very keen of the magnitude of the convection coefficient while performing predictions for a given food freezing application problem.

Introduction
Food freezing as a conservation process plays an important role in the global food industry. Even though it is so widely used, the complexity involved in the process makes it very difficult for one to predict the food freezing times accurately. Nevertheless, the time of freezing is perhaps the most important variable to know before building a food freezing unit. In continuous freezing processes, usually performed in spiral or carton freezers, the size of the freezing equipment necessary to perform a certain freezing task is proportional to the freezing time. Therefore, underestimating freezing time can lead to insufficient production capability of the equipment, whereas overestimating it may result in unnecessary costs in building materials, physical space and energy.
A major issue when dealing with modeling of food thermal processes is dealing with the wildly varying food thermophysical properties, which are dependent on food temperature and composition.
These properties experience large and sudden changes during freezing, mainly because the formation of ice crystals. As foods do not present a freezing point but a threshold temperature at which freezing begins, the problems get more complex to solve. At this point, thermophysical properties experiment steep changes. Choi and Okos (1986) have modeled these changes, developing models for food properties as functions of temperature and composition. Their models are indicated by ASHRAE (2015) to be used in the design of freezing systems. As the use of realistic models for thermophysical properties is a key point Numerical simulations of food freezing times using thermophysical food property models Fernando Zigunov, Flávia Zinani, Janice da Silva 3 Estudos Tecnológicos em Engenharia, vol. 12, n. 1, p. 1-13, jul/dez 2018 to achieve reliable solutions, the use of analytical methods to determine freezing times is highly impaired.
Empirical correlations, such as those developed by Plank (apud Stoeker, 2004) and Pham (1986) produce results that might not be much accurate as they do not consider all the mechanisms involved in the process of freezing. The results produced by these two correlations will be compared with the numerical method presented here.
In the present work, a numerical model to simulate heat conduction during food freezing was developed. It is based on the finite difference method, which employs the equations of Choi and Okos (1986) to model thermophysical properties. An experimental test bench, consisting of a pilot scale freezing tunnel, was also built in order to validate the numerical results and elucidate some issues about the physics of the process. The numerical model was designed to solve one-dimensional transient conduction problems; therefore, it could only be applied to geometries that have one-dimensional temperature profiles: thin slab, long cylinder and sphere. Even though these geometries might be uncommon in carton freezers, they are frequently found in linear or spiral conveyor belt individual quick freezing (IQF), such as hamburgers, meatballs, sausages, fish fillets, etc.
While a food product undergoes a freezing process, several phenomena happen at various scales.
The refrigerating medium (air, in this case) gains heat from the surface of the product by convection, and humidity present in the surface of the product diffuses to the air. The water present in the surface of the product evaporates or sublimates, causing surface dehydration and weight loss in unpackaged products (Campañone et. al., 2005). This phase change that happens in the surface water causes an increase in the heat flux, since the latent heat needed to cause phase transition to vapor is given by the product.
Inside the product, heat is transferred mainly by thermal conduction. While the product freezes, its apparent specific heat increases greatly, as demonstrated experimentally by Pham et. al. (1994). Since the transient freezing front is local, the product experiences a local sudden decrease of thermal diffusivity.
The result is that the freezing front advances slowly, and the temperature profile at any given time is relatively complex.
Concerning the microscale phenomena, the aqueous solutions present in the food products are subject to a freezing temperature shift depending on their concentration. Also, during freezing the ice phase gradually separates from the aqueous phase, and in this process the aqueous phase concentrates, increasing the freezing temperature shift (Corry, 1987apud Erickson, 1997. Due to this phenomenon, the aqueous solution freezes gradually, releasing its latent heat over a wide range of temperatures. As Pham (2006) reviews, there are several approaches to solve a mathematical grid that computes the complexities involved with variable thermophysical properties, such as apparent specific heat methods and enthalpy methods. In this work, an apparent specific heat method was used to compute the grid, due to its ease of implementation. The drawback is the larger amount of computational effort to Numerical simulations of food freezing times using thermophysical food property models Fernando Zigunov, Flávia Zinani, Janice da Silva 4 Estudos Tecnológicos em Engenharia, vol. 12, n. 1, p. 1-13, jul/dez 2018 solve a given problem accurately. Since the problems examined here would be in the simpler 1D grid scenario, the computational time was still acceptable.

Numerical model
In order to predict freezing times, a mathematical model for the transient temperature field inside a food sample subjected to convective flow was built. The problem domain consisted of the food sample, so the heat transfer is governed by the general heat equation with null source term (Bergman et al., 2011) and the boundary conditions of heat convection and heat flux, in cases in which water evaporation was considered. Only one-dimensional problems were considered. Food samples consisting of flat plate, long cylinder and sphere geometries were employed, and the heat equation in rectangular, cylindrical or spherical coordinate system was employed in each case, respectively: where x is the position (radial distance in the case of cylindrical and spherical coordinate systems), k is the thermal conductivity, T is the temperature,  is the mass density, cp is the specific heat and t is the time. The coefficient m is equal to 0 for rectangular, 1 for cylindrical and 2 for spherical coordinate systems, respectively (Bergman et al., 2011).
The finite difference method, with a Crank-Nicholson scheme, was employed to solve Eq. 1 for the temperature field in space and time (Chapra and Canale, 2010). An algorithm was implemented using Microsoft Excel® VBA to solve the resulting system iteratively, with a time step of 1/300 of the total experiment time. The termophysical properties , and were calculated at each time step and position, for they were assumed to be functions of temperature and food composition during the freezing process, as explained in item 2.2.
The boundary conditions implemented were Neumann in the center (zero heat flux) and in the surface (known heat flux in the surface, calculated iteratively from convection parameters or evaporation parameters, using 10 iterations per time step). In order to obtain the heat flux in the surface, two individual heat flux components were assumed: one due to temperature difference and another due to water concentration difference (mass convection). The latter was included after the execution of some experiments, which showed that the wet surface effects were of great importance in freezing phenomena.
To obtain the convection heat transfer coefficient, the following usual correlations were employed:  In order to obtain the heat flux due to combined heat and mass transfer, the enthalpy potential model derived by Stoecker and Jones (1983) was employed. In this model, the total heat flux through a wet surface, q", is not only a function of the heat convection coefficient, h, and the temperature difference between the surface and the fluid. The heat flux is given by where , is the enthalpy of the air, given the surface temperature and humidity, ,∞ is the enthalpy of the surrounding air and , is the moist air specific heat, which can be calculated as in which , is the dry air specific heat, , the water specific heat and w the local water content.

Food properties
The correlations by Choi and Okos (1986) were employed to model thermophysical properties of foods in the present work. These correlations are based on the assumption that the food pieces are homogeneous mixtures of six basic components: water, ice, carbohydrates, fat, proteins and ashes. The correlations predict thermal conductivity, specific heat, models for properties of each of these components as a function of temperature, a series of equations are employed to get the apparent thermal properties of a given food. ASHRAE (2015) as well as Fricke & Becker (2001) review the models for predicting food thermophysical properties. These references agree in assuming foods as homogeneous mixtures of six basic components: water, ice, carbohydrates, fat, proteins and ashes. By using the Choi & Okos (1986) thermophysical models for properties of each of these components as a function of temperature, a series of equations are employed to get the apparent thermal properties of a given food.
The graphs shown in Figure 1 show two sample property curves for the apparent specific heat and enthalpy of a hamburger that is composed of 66.3% of water, 19.7% of proteins, 11.4% of fat and 2.6% of ashes, with an initial freezing temperature of -1.9 ºC, for the mentioned model. This data was used as an input to the numerical model.

Experimental setup and procedure
A pilot-scale freezing tunnel, capable of achieving internal temperatures of -30 ºC and air speeds up to 4.5 m/s was used to gather experimental data to compare with the simulation results.
The experimental setup is shown in Figure 2.  Table 1. The test conditions tried to reproduce the usual application parameters in industrial IQF freezers.

Data processing
The data gathered from each experiment was a set of temperature measurements, collected from the beginning of the experiment until the product reached approximately -20 ºC. The temperature measurements were acquired with a frequency of 0.2 Hz, and the measurement uncertainty of the PT-100 probes used was determined to be ±1 ºC. This set of measurements was plotted into a curve of temperature versus time. Afterwards, the temperature measurement set was loaded into the mathematical solver, along with the experimental parameters described on Tab. 1.
The mathematical solver is a VBA program developed inside a Microsoft Excel ® Workbook.
This platform was mainly chosen because of its large availability in the industry. It runs the algorithm described in Section 2.1, based on product and process features inputted by the user, and plots a curve Numerical simulations of food freezing times using thermophysical food property models Fernando Zigunov, Flávia Zinani, Janice da Silva 8 Estudos Tecnológicos em Engenharia, vol. 12, n. 1, p. 1-13, jul/dez 2018 that represents the predicted temperature drop curve for that product, based on the thermophysical properties model given by ASHRAE. Varying from center to surface

Comparison between experimental and numerical results
Representative temperature evolutions through time of test subjects are shown in Figure 3.
After each experiment, the product was cut and the sensor position was measured in relation to the surface, and the measured values are shown below each graph. The predictions from the simulations were far from the results obtained from the experiments. In Figure 4 a chart showing the predicted freezing times (until -18 ºC as measured by the sensor) versus the times obtained experimentally is shown to summarize the data obtained. The model proposed here predicted within ±15% only 55% of the samples, with prediction errors varying between -25% and +96%.   In order to test this hypothesis and quantify the deviation of the experimental convection coefficient from the convection coefficient predicted by the correlations, a multiplier parameter ℎ , shown in Eq. (8), was inserted in the mathematical solver to change the effective coefficient ℎ in the simulation. With this free parameter, a curve fit was possible and the results obtained were within R 2 >0.9 for all experiments. Nevertheless, when this parameter was changed, the curve fit had a very large fidelity to the experimental data, as displayed in Figure 6.

Comparison of the numerical solver predictions against other models in literature
Two correlations present in the literature were used as references to compare the quality of the predictions of the numerical solver developed in this work. The predictions produced by these correlations are compared with the experimental freezing times in Figure 7. The first correlation, shown as blue triangles is Planck's correlation. It dates from 1944, and was retrieved from Stoecker (2004).
The second correlation, shown as red circles, is designed by Pham (1986), and builds upon Planck's correlation by separating the pre-freezing, freezing and post-freezing times.  The method presented here, without any convection coefficient adjustment, performed similarly as the correlation by Pham, both predicting 25 out of 40 experiments within ±30%. It seems that Pham's correlation tends to overestimate more the freezing time than the method presented here.
On the other hand, only 11 experiments out of 40 were predicted by Planck's correlation within ±30%. Also, it underpredicted the freezing time of 31 experiments (77.5%), which is highly undesirable in an equipment sizing scenario.

Final remarks
The results of this experimental work show that the prediction of a food product's freezing time is not a straightforward task, where simple correlations will likely yield poor results. Based on the data obtained, it is possible to affirm that the transient conduction problem inside a food product can be satisfactorily modeled by using the food property models presented by Choi and Okos (1986) and the basic heat diffusion equations. It seems that other authors also struggle to get good predictions with their numerical methods for this application due to lack of the convection coefficient accuracy. Evans et. al. (1991) for example, reports prediction inaccuracies up to 24.6% on a boiling fluid boundary with a FDM scheme. Authors that have reported good fit between the predicted and experimental freezing times often have in some way measured the convection coefficient at the boundary (Perussello et. al., 2011;Wilson and Singh, 1987).
However, from a technology application point of view the fact that the freezing time is such a strong function of the convection coefficient at the boundary makes it fundamental to perform good predictions on the latter. Therefore, simply solving the transient heat equations for a food freezing might be of limited application with regards to accuracy when the convection coefficient at the boundary is highly variable or hard to determine, as is the case for many industrial applications.