Modeling and Simulation of Forced-air Cooling of Strawberries Using Variable Convective Coefficient

In the forced-air cooling process of fruits occurs, besides the convective heat transfer, the mass transfer by evaporation. The energy need in the evaporation is taken from fruit that has its temperature lowered. In this study it has been proposed the use of empirical correlations for calculating the convective heat transfer coefficient as a function of surface temperature of the strawberry during the cooling process. The aim of this variation of the convective coefficient is to compensate the effect of evaporation in the heat transfer process. Linear and exponential correlations are tested, both with two adjustable parameters. The simulations are performed using experimental conditions reported in the literature for the cooling of strawberries. The results confirm the suitability of the proposed methodology.


INTRODUCTION
Brazil is the third largest producer of fruits, total of, in 2008, approximately 43 million tons (IBRAF, 2010).The total volume of strawberries produced in the country is around 100 thousand tons, of which just over 50% is still consumed fresh.The remainder is used as ingredient by food industries or to produce frozen pulp.In general, strawberries sold in the form of juices or jellies get their values raised by up to 1000% (IBRAF, 2010).
There is a great waste of the fruit produced, due to its high perishability, so it is very important to develop technologies that will extend the life of these foods, so that they can be offered over long distances and periods between harvests.
Among the various post-harvest techniques, the fast cooling of fresh fruit for their optimum temperature storage is essential to maintain quality and increase service life (FERRUA & SINGH, 2009a).
The forced-air cooling consists on passing cold air through the stack of packages containing fruit.During the cooling, the convective coefficient of heat transfer is established depending on the size of packaging, its opening area of distribution, the characteristics of the fruits (temperature, water content, specific heat and geometric shapes), and the characteristics of the cooling air (temperature and velocity) (THOMPSON et al., 1998).
Temperature control also prevents the development of pathogenic and spoilage-causing microorganisms.To better use the technology of fast cooling, experimental analysis or aid of mathematical tool are necessary.The cooling is given as completed once the center of the strawberry reaches 7 / 8 of cooling, obtained through the difference between the initial temperature and the one wanted to store the product (PIROZZI, 2002).
Particularly, the problem of predicting the cooling time of the products and characterization of the process parameters of heat transfer has been studied over many years.The common goal of these studies is to develop a simple method of prediction, which requires a minimum amount of data and a suitable computer program.On the other hand, the heat transfer convective coefficient has a very important role in processes involving convection, one of the most common causes of miscalculating the temperature inside the solid (TERUEL et al., 2001).The coefficient convective is the object of study of AMENDOLA & PIROZZI (2005) and AMENDOLA ( 2009), for cooling of strawberries and figs subjected to fast cooling.Correlation for heat transfer convective coefficients in forced-air cooling of tomatoes and oranges are presented in KUMAR et al. (2008).
The mathematical models used to study the cooling of fruits, in general, express transient conduction problems, with many geometry and boundary conditions simplified, which in some cases can be solved by analytical solutions and in other cases by numerical methods.Mathematical modeling and numerical simulation of strawberries forced-air cooling through the use of numerical methods are tools that can create conditions that will serve as support in decision making, both in the use of this technology, and the planning, implementation, manufacture and operation of equipment used in forced-air cooling (PIROZZI, 2002).
In the process of fast cooling of strawberries in cardboard box, there is water loss of the fruit by evaporation.Under the experimental conditions applied in the PIROZZI (2002) and AMENDOLA & PIROZZI (2005) simulations and in the present study, this loss is about 1 to 2% of the total mass.Early in the cooling process, the cold air in contact with the strawberry's surface is heated and has its relative humidity decreased as the saturation humidity increases with temperature.As the surface of the strawberry has free water to the same initial temperature, the process of water evaporation is favored.In the process, the strawberry's surface temperature is reduced, which leads to a decrease in evaporation.The process of water evaporation consumes energy (latent heat of vaporization) removed from the strawberries.The joint simulation of the processes of heat and mass transfer is not usually applied, since the goal is to determine the cooling time and not the amount of water evaporated.However, if it is used, the computational cost involved and the complexity of the model are far superior.In FERRUA & SINGH (2009b), FERRUA & SINGH (2009c) and FERRUA & SINGH (2009d) complete modeling is performed, considering the heat and mass transfer, and momentum in the process of a forced-air cooling in a strawberry box.These authors use a commercial package CFD (Computational Fluid Dynamics) for the numerical solution of the model.It should be noted that the CFD analysis and the complexity of modeling involves a very high computational cost.This paper proposes the modeling of heat transfer during cooling using a convective coefficient which varies with temperature, in order to compensate for this effect of heat loss by evaporation.It is noteworthy that the existence of these two types of mechanisms involved in the fast cooling by forced air (heat transfer by convection and energy loss by evaporation) is reported in FERRUA & SINGH (2009a).The idea of the proposed approach is the development of a less complex and low computational cost, able to be simulated on any personal computer, unlike the models based on CFD analysis.

Mathematical Model Applied
The equation that rules the process of one-dimensional heat transfer (temperature as a function of radius and time) in a strawberry, considered spherical, with no heat generation, is described below: In which, T is the temperature; r is the radial position; t is the time;  is the thermal diffusivity.
The initial condition of eq. ( 1) is the temperature before cooling the strawberry: In which, T m0 is the initial temperature of the strawberry.
Since the equation is of second order, it has two boundary conditions.The first is radial symmetry: The second boundary condition is the heat flux on the strawberry's surface: In which, h is the convective coefficient of heat transfer; T amb the ambient temperature (cooling chamber) during cooling; k is the thermal conductivity.
Considering that the convective coefficient is constant during the cooling process, the proposed problem has analytical solution (INCROPERA et al. 2008).However, the strategy of this coefficient of variation with temperature, adopted in this paper, promotes a non-linearity in the boundary condition (eq.4) that leaves the model with no analytical solution.This way, the use of numerical method is necessary and a program in Matlab language is built.
As numerical solutions have truncation errors, the program built was validated by simulating the problem with constant h and the results are compared with the analytical solution.

Analytic Solution of the Problem with constant h
In order for this solution to be generic and not specific to each set of values of , k and h, the dimensioning of the variables of the problem becomes needed.The following equations provide dimensionless used in the analytical solution.
In which, In which, y is the dimensionless radius; R the strawberry's radius In which, Bi is the Biot number (dimensionless widely used in the characterization of heat transfer processes).
In which, Fo is the Fourier number (dimensionless for the time).
From these dimensionless values, the analytical solution of the problem with constant h is an infinite series (INCROPERA et al. 2008): In which, C n are series coefficients calculated by the equation: It is observed that the center of the strawberry (y = 0), eq. ( 9) has a singularity and its value should be calculated by the limit of this expression, the result is displayed: To find the times of 1 / 2 and 7 / 8 cooling, it must determine the values of Fo (from which can be determined the values of time) that make eq.( 11) equal to the desired value.That is,  (0,Fo)=1/2 to t 1/2 e  (0,Fo)=7/8 to t 7/8 .

Numerical Solution of the Problem with variable h
As already discussed, the nonlinearities in the boundary condition (eq.4) leave the model of heat transfer in the strawberry with no analytical solution.Several methods could be used in numerical solution of this model, whereas this study was chosen by the method of finite differences.
The spatial direction (r ranging from zero to R), was discretized.It was used an uniform mesh with n internal points, and the temperature at the center of strawberry given by T 0 and the surface by T n+1 .For each of the internal points, central finite difference approximations of second order were used in the determination of the derivatives: In which, r is the distance between two points of discretization (r = R/n); i is the discretization index (i = 1 the n for the internal points).
The temperature at the center of the strawberry, T 0 , is calculated by discretizing the symmetry condition given by eq. ( 3), by finite differences ahead of the second order: The temperature at the center of the strawberry, T n+1 , is calculated by boundary condition, given by eq. ( 4), discretized by finite differences behind of the second order: It is noteworthy that the first order approximations for the derivatives in the contours could have been applied, as used in PIROZZI & AMENDOLA ( 2005), but the convergence of the mesh would be slower (more discretization points are needed).
After the discretization of all spatial derivatives, the system of partial differential equations becomes a system of ordinary differential equations (ODEs), the independent variable time.The discretization of time derivatives by finite differences, as adopted, is possible, but in general these approximations are of first order which also increases the numerical error.In this study, it was opted for the integration of these ODEs, using a MATLAB function implementing the Runge-Kutta Method of the fourth order.

Validation Procedures and Empirical Functions Tested
In order to demonstrate the convergence of the mesh, it was simulated the experimental conditions used in PIROZZI (2002) and presented in Table 1, considering convective coefficient constant reported by these authors and comparing the result with the analytical solution.The numerical solution reported by PIROZZI ( 2002) is also shown for comparison.TABLE 1. Experimental conditions and used properties (PIROZZI, 2002).After the demonstration of the convergence of the mesh and validation of the method used, it was adopted two strategies to change the convective coefficient as a function of temperature in order to incorporate the effect of evaporation in the process of heat transfer.The empirical equations tested are based on linear variation and the exponential variation of convective coefficient with temperature on the surface of the strawberry as the eq.( 15) and eq.( 16).
In which, h 0 is a parameter that gives the value of h when T n+1 tends 0 o C;   is a parameter whose value controls the effect T n+1 in h.
It is observed that both equations have two parameters (h 0 and ) whose values are estimated to minimize the sum of squared difference between the numerical and experimental values of times 1/2 and 7/8 of cooling reported in PIROZZI ( 2002) and PIROZZI & AMENDOLA (2005).In the minimization process, it was used the implementation of Quasi-Newton's method of MATLAB.

Validation of the Numerical Method
As discussed, to validate the numerical method proposed, it was held comparison of the results obtained for times of 1/2 and 7/8 of cooling with the values obtained by analytical solution.In the analytical solution, it was applied 10 terms with this condition and the truncation error is less than 10-20, for all simulations.The simulated experimental conditions are described in Table 1.The results obtained are presented in Table 2, together with the simulated results presented in PIROZZI (2002).
It can be seen in this Table 2 that the simulated data in this study differ with regard to the analytical solution up to 0.04 minutes or 0.15% (t 1/2 in condition 3), indicating that the numerical solution converged.There is, however, that the numerical results presented in PIROZZI ( 2002) and PIROZZI & AMENDOLA (2005) differ by up to 5.58 minutes (t 7/8 in condition 6), the maximum relative error of 14.44% (t 7/8 in condition 4).It is noteworthy that in these studies, the analytical solution is not used in the validation being used only the convergence of the mesh space (increasing the number of discretization points in space until the solution does not vary significantly).Another important point is that the authors discretize time with a first order approximation with step size (t used in the integration) constant, which results in higher integration errors.These errors were minimized in this study since the integration is performed using a MATLAB function which uses the Runge-Kutta of fourth order with variable step size in order to limit the integration error.The h value used for each experimental condition was estimated at PIROZZI (2002) and notes that there is a big difference between the simulated and experimental of cooling times for all conditions.It is not possible to set a constant value of h that promotes a good estimate of both cooling times.To show this impossibility, it presents the value of relative error (less experimental simulated) as a function of h for conditions 3 and 4 (see FIGURE 1).FIGURE 1. Variation with the h value of the relative error of t 1/2 and t 7/8 using (a) experimental condition 3 and (b) experimental condition 4.
It can be seen in FIGURE 1, that for the third experimental condition, the value of h should be between 12.1 and 19.5 5 W.m -². 0C -1 , which are, respectively, the optimal values of this parameter for determining the times t 7/8 and t 1/2 (values of h that zero out the relative errors for each of these simulated temperatures).As for the four experimental conditions, these optimal values are, respectively, 21.9 and 124.3 W.m -². 0C -1 .In both cases, the value of h estimated by PIROZZI (2002) is between these two values, which is expected since it has an optimization problem with two goals and it is necessary to establish a compromise between them (a goal is to reset the error to the value of t 7/8 and the other corresponds to the value of t 1/2 ).It is also observed that in both conditions, the value of h should be higher in order to reproduce the value of t 1/2 , indicating that early in the process there is a greater heat loss to the strawberry, corroborating the assumption that the energy required for evaporation, while the surface of the strawberry is not yet cooled, has influence in the process.
Table 3 presents the results of the strategy to offset the effect of evaporation in the heat transfer process using a linear function of h with the temperature on the strawberry surface (eq.16).It is noteworthy that in the process of estimating the parameters, the restriction of positive h 0 was applied in order not to be obtained negative values of h for surface temperatures close to 0° C, since this would not have physical meaning.It can be seen in this Table 3 that for the experimental conditions 1, 3 and 6, both times t 1/2 t 7/8 was reproduced with an error less than 0.01 minutes (0.05%).For other conditions, the errors are much lower than those reported for h constant (Table 2), but it was not possible to cancel them completely, especially for condition 2. This way, the linear variation of h with the surface temperature of the strawberry was not effective for all conditions tested.
The results obtained using the h exponential variation with temperature (eq.17) are presented in Table 4.It can be seen in Table 4, that all values of the cooling times were obtained with a maximum error of 0.01 min (0.05%).It is also observed that the highest values of h 0 and  were obtained for conditions that cannot be efficiently simulated using a linear variation of h with temperature.These higher values promote higher convective coefficients at the beginning of the cooling, and also being higher their decrease rate with the temperature of the strawberry surface.
While the parameters of the curves have been estimated taking into account only the time 1/2 of cooling and 7/8 of cooling, it was observed that the simulated results for the temperature variation in the center of the strawberry were closer to the experimental than those presented in PIROZZI ( 2002 It can be seen in FIGURE 2, that although the coefficients of the function for calculation of h as a function of temperature has considered only two experimental points (t 1/2 e t 7/8 ), the simulated results for the temperature at the center of the strawberry is closer to the experimental values than the results presented in PIROZZI (2002) and PIROZZI & AMENDOLA (2005), whose estimation of h considers all points of the curve.This result confirms the suitability of using h varies exponentially with the temperature of the strawberry surface.

CONCLUSION
In this study it was propose the use of empirical functions for the variation of the convective coefficient of heat transfer (h) as a function of surface temperature in the simulation of the cooling of strawberries considered spherical.The functions proposed for h make non-linear the boundary conditions of the problem of heat conduction in one-dimensional spherical coordinates.In this way a numerical solution of the problem is needed.All spatial derivatives are discretized using finite difference approximations of second order.The system of ordinary differential equations, resulting from discretization of the spatial direction, is integrated in time using Matlab implementation of the Runge-Kutta method of fourth order.The validation of the numerical method is performed by solving the problem with the h constant, because for this situation an analytical solution exists.
The validation results obtained indicate that the numerical procedure introduces errors less than 0.15%.Furthermore, it was found that the numerical results of PIROZZI ( 2002) and PIROZZI & AMENDOLA (2005) are not converged.It was also observed that using the h constant, it is not possible to reproduce both t 1/2 and t 7/8 , and the results indicate that h should be greater at the beginning of the cooling process.
The best correlation for the h variation with the strawberry surface temperature was the exponential function, being the experimental cooling values t 1/2 and t 7/8 reproduced with errors less than 0.05%.Applying this strategy, there is a significant improvement compared with literature data, the reproduction of the experimental values of the temperature at the center of the strawberry as a function of cooling time.