Influence of the geometry on the numerical simulation of the cooling kinetics of cucumbers

In this paper, the effect of the geometric representation of cucumbers on the numerical simulation of its cooling kinetics is studied. It is supposed that the diffusion model with boundary condition of the third kind satisfactorily describes the cooling, and that the thermo-physical parameters are constant during the process. The geometries used to represent the cucumber are: infinite cylinder, f inite cylinder, and ellipsoid. The diffusion equation was solved through the finite volume method, with a fully implicit formulation, using cylindrical and generalized coordinates. The convective heat transfer coefficient and the thermal diffusivity were determined through optimization, using the inverse method. The best model in the representation of the cucumber’s shape was the ellipsoid, but the time demanded in its optimization was about 66 times greater than the time for the infinite cylinder. Additional key words: convective heat transfer coefficient; Cucumis sativus; cylindrical coordinates; diffusion; finite volume method; generalized coordinates; thermal diffusivity.


Introduction
In order to describe the drying and cooling kinetics of foodstuffs, several authors use diffusion models (Dincer, 1996;Gastón et al., 2002;Doymaz, 2005;Amendola and Queiroz, 2007;Resende et al., 2007;Hacihafizoglu et al., 2008;Silva et al., 2008;Cuesta and Lamúa, 2009).If the thermo-physical parameters and the dimensions of the solid are considered constant, the diffusion equation has an analytical solution for several shapes such as a sphere, a cylinder and a slab, among others.Such solutions can be obtained, for instance, in Luikov (1968) and Crank (1992).In order to solve the diffusion equation, the simplification of the shape is made by many authors (e.g., Dincer, 1996;Doymaz, 2005;Amendola and Queiroz, 2007;Resende et al., 2007;Silva et al., 2008), and in these works, good results are reported in the determination of the cooling and drying parameters.Naturally, the approaching of the geometry of the solid to a simpler shape introduces an error in the calculated values for the thermo-physical parameters and also in the description of the drying or cooling kinetics.
There are analytical solutions in the literature considering the shrinkage and the variable thermo-physical parameters (Ruiz-López and García-Alvarado, 2007) but, in general, for a solid with any shape the diffusion equation should be numerically solved.One of the methods used to solve the diffusion equation numerically is the finite volume method, with fully implicit formulation (Patankar, 1980).Wu et al. (2004) used this method to describe the heat and mass transfer inside a single rice kernel during its drying process.In this work, the authors used generalized coordinates due to the irregular shape of the grain.
The discrepancy on the drying kinetics description due to the shape of the solid was studied by Hacihafizoglu et al. (2008) using experimental data for rice.The geometries used to approach the grain's shape were: sphere, finite cylinder, and prolate spheroid.In this study the water diffusivity and the volume of the rice grain were considered constant and this consideration made it possible to use analytical expressions as solutions of the diffusion equation, regarding the boundary condition of the first kind.The authors found that a prolate spheroid represents better the drying kinetics of rice.
The geometry effect on the water diffusivity was studied in PROINTA-Isla Verde and broom wheat cultivars (Gastón et al., 2003).The authors compared the spherical and axisymmetric ellipsoidal geometries applying the finite element method to solve numerically the diffusion equation.In this work, the researchers considered the water diffusivity as having a constant value and the boundary condition of the first type.Although the geometric effect has been studied in several drying processes, no works were found in the literature with the objective of studying this effect on the cooling of fruits and vegetables, particularly using the finite volume method in arbitrary domain.
This paper uses numerical solutions for the diffusion equation to determine thermo-physical parameters and to compare cooling kinetics of cucumbers, considering three shapes for the vegetable: infinite cylinder, finite cylinder and ellipsoid.In order to accomplish that, the diffusion equation will be numerically solved through the finite volume method, using a fully implicit formulation.For all the shapes, the thermal diffusivity and the convective heat transfer coefficient will be determined through optimization, using the inverse method.

Material and methods
The mathematical models for the numerical solution of the diffusion equation presuppose the following hypotheses: -The solid must be homogeneous and isotropic.
-The spatial distribution of the temperature must be initially uniform.
-The only mechanism of transport of heat within the solid is diffusion.
-The process parameters are constant during the diffusion.
-The boundary condition is of the third kind.
-The cooling medium temperature must be constant.
-There must be no phase change in the product during the process.
-The mass loss and the source term due to the respiration heat are negligible during the cooling of the cucumber.

Solution of the diffusion equation: cylindrical coordinates
In cylindrical coordinates, the heat conduction equation in an infinite cylinder (length L c much greater than the radius R c ) is given in the following way (Luikov, 1968;Crank, 1992): Influence of the geometry on the description of the cooling of cucumbers [1] where Φ is the temperature, r defines a radial position inside the infinite cylinder, Γ Φ is the thermal conductivity, λ = ρc p , in which ρ is the density and c p is the specific heat, while S is the source term.In order to discretize Eq. [1], a uniform mesh will be created, with N control volumes, in the circular section of the infinite cylinder, as is observed in Figure 1.
The diffusion equation for the model will be solved through the finite volume method (Patankar, 1980), with a fully implicit formulation.A similar problem was numerically solved for the boundary condition of the first kind by Silva et al. (2008).Thus, for an internal control volume (from control volume number 2 up to N), and for the control volume number 1 (in the centre), the discretized equations can be obtained in Silva et al. (2008).
The convective boundary condition is expressed by the following expression: [2] where h is the convective heat transfer coefficient, Φ b is the temperature at the boundary, and Φ ϱ is the temperature of the cooling medium.
[1] with respect to space (2π r p ∆r L c ) and time (∆t), and taking into account Eq. [2], we obtain the following result for the control volume P at the boundary: [3] in which [4a-c] where the superscript 0 means «former time» t and its absence means «current time» t + ∆t.On the other hand, the subscripts «w» and «e» mean west and east (see Fig. 1).This solution will be used for the description of the cooling kinetics of the cucumber represented by an infinite cylinder.The source term S will be considered zero and Gauss-Seidel Method is applied for all geometries considered.

Solution of the diffusion equation: generalized coordinates
For a solid with arbitrary shape, the diffusion equation can be solved using structured meshes (Thompson et al., 1985;Maliska, 2004) and generalized coordinates (Tannehill et al., 1997;Maliska, 2004).The diffusion equation in this system of coordinates, assuming variable thermo-physical parameters, can be written in the following way (Maliska, 2004;Wu et al., 2004;Silva et al., 2009a): [5] where , η, and γ are the axes in the generalized coordinates system (transformed domain), is the time and J is Jacobian of the Cartesian (x,y,z)-generalized (, η, γ) coordinates transformation.In Eq. [5], Φ is the temperature, λ and Γ Φ are process parameters referring to the thermo-physical properties: λ ϵ ρc p , in which ρ is the density and cp is the specific heat, while Γ Φ is the thermal conductivity.If the solid is obtained by the revolution of a vertical plane area (η or xy) around y axis, Jacobian J is given by the determinant [6] In addition, due to the symmetry, terms involving ‫‪γ‬ץ/ץ‬ and the coefficients α 31 (or α 13 ), α 32 (or α 23 ) and α 33 are zero.The coefficients α 11 , α 22 , α 12 and α 21 are given by the expressions: In Eqs.[7], the terms of the type «v u » represent the partial derivative of «v» with respect to «u».Admitting the simplifications above, Eq. [5] becomes: [8] + The discretized equation for an internal control volume is given by Silva et al. (2009a), where a similar numerical solution is presented for the boundary condition of the first kind.For a control volume at the boundary (for instance, east boundary), using the finite volume method, with a fully implicit formulation, the following result is obtained for a control volume P, subjected to the third kind boundary condition, by integrating Eq. [8] with respect to space and time (Silva et al., 2010a): [9] where: [10a-g] Again, the superscript means «former time» t and its absence means «current time» t + ∆t.On the other hand, the subscripts «e», «w», «n», «s» mean east, west, north and south, while «se» and «ne» mean south-east and north-east.In a general way, the factor «f» in Eqs.[10] is given by the following expression: [11] where ∆n represents the distance from the nodal point P up to an interface of the control volume.
The aforementioned system of equations can also be solved for the dimensionless temperature, which is defined in the following way: [12] where Φ 0 is the initial temperature (K), and is supposed uniform.In this case, for t = 0, Φ* = 1 and for t → ϱ, Φ* = 0.

Optimization
An optimization algorithm using the inverse method was developed considering the following requirements: -Minimization of the chi-square regarding the fit process of a simulated curve to the experimental data; -Use of Levenberg-Marquardt's algorithm (Press et al., 1996;Silva et al., 2006), with corrections in sequence of the parameters.
The expression for the chi-square regarding the fit of a simulated curve to the experimental data is given by (Taylor, 1997): [13] where Φ exp i is the experimental value of Φ at time t i , Φ sim i is the correspondent simulated value, N p is the number of experimental data, and 1/σ 2 i is the statistical weight of the experimental point «i».In the absence of information, those weights in Eq. [13] can be made equal to the unit.Obviously, the chi-square depends on Φ sim , which depends on the process parameters Γ Φ and h, and these can be determined by the minimization of χ 2 .If the cooling involves an interval of temperature which is not much large, the parameters Γ Φ and h can be considered constant.Then, the minimization of the objective function is accomplished in cycles involving the following steps (Silva et al., 2010a,b): 1) Inform the initial values for the parameters «Γ Φ » and «h».Solve the diffusion equation and determine the chi-square.
2) Inform the value for the correction of «h».
3) Correct the parameter «h», maintaining the parameter «Γ Φ » with constant value.Solve the diffusion equation and calculate the chi-square.
4) Compare the latest calculated value of the chisquare with the previous one.If the latest value is smaller, return to the step 2; otherwise, decrease the last correction of the value of «h» and proceed to step 5. 5) Inform the value for the correction of «Γ Φ ». 6) Correct the parameter «Γ Φ », maintaining the parameter «h» with constant value.Solve the diffusion equation and calculate the chi-square.
7) Compare the latest calculated value of the chisquare with the previous one.If the latest value is smaller, return to the step 5; otherwise, decrease the last correction of value of «Γ Φ » and proceed to step 8. 8) Begin a new cycle coming back to the step 2 until the stipulated convergence for the parameters «Γ Φ » and «h» is reached.
The algorithm presented above is based on Levenberg-Marquardt algorithm (Press et al., 1996) in the sense that the corrections of the parameters are variable.In each cycle, the value of the correction of each parameter can be initially modest, compatible with the tolerance of convergence imposed to the problem.Then, for a given cycle, in each return to the step 2 or 5, the value of the new correction can be multiplied by the factor 2. If the modest correction initially informed in step 2 or 5 does not minimize the objective function, in the next cycle its value can be multiplied by the factor -1. On the other hand, initial values for the parameters can be estimated through obtained values for similar products available in the literature.

Cooling of cucumbers: experimental data
Experimental data obtained by Dincer (1996) for the cooling of a cucumber are explored in the present paper.The temperatures were measured with a thermocouple placed in the centre of the cucumber.As the uncertainties of the temperatures are not available, the statistical weights were made equal to 1 in Eq. [13].The dimensions of the cucumber used in the experiment were: radios R c = 0.019 m and height L c = 0.160 m.The decimal moisture content of the cucumber was X = 0.96 (wb), and the mass loss during the cooling process was found to be negligible.The initial temperature of the cucumber was 295 K, while the temperature of the cooling air was 277 K.The cooling air velocity was kept at 2 m s -1 and its relative humidity was 80%.
In the present paper, the specific heat was estimated from the expression given by c p = 1.381 + 2.930X (Sweat, 1986;ASHRAE, 1993), with c p given in kJ kg -1 K -1 when the moisture content is decimal in wet basis.Using this latest information, the value for the specific heat of cucumbers was estimated: c p = 4,190 J kg -1 K -1 .The density of the cucumber was determined by Fasina and Fleming (2001) for a vegetable with the same moisture content of the cucumber currently analyzed: ρ = 959 kg m -3 .For comparison with the results to be obtained in this paper; an estimate of the thermal conductivity can be obtained from expression Γ Φ = 0.148 + 0.493X (ASHRAE, 1993), which for the cucumber of this experiment gives Γ Φ = 0.621 W m -1 K -1 .

Developed software and statistical indicators
The tests about the influence of the geometry on the numerical simulations were performed in a computer Intel Pentium IV with 2 GB (RAM).The compilation of the developed source codes was made in the studio Compaq Visual Fortran (CVF) 6.6.0Professional Edition, using the programming option QuickWin Application, and the platform was Windows Vista.On the other hand, in all the optimization processes, the relative tolerance on the parameters was 1.0 × 10 -4 .Besides the chi-square which will be used as a statistical indicator, the coeff icient of determination R 2 (Taylor, 1997) was also utilized to determine the quality of the fit.

Results
In this paper, it was assumed that the location of the thermocouple within the cucumber is exactly in the centre of the vegetable as in Dincer (1996), from where the experimental data was obtained.

Initial values: considerations
The best way to save time, when simulating the cooling kinetics for the three geometries, is to initially perform the optimization for the inf inite cylinder, which is a faster process (one-dimensional) than the others, determining the parameters Γ Φ and h.Then the obtained values should be used as initial values for the finite cylinder, which involves an orthogonal structured mesh, and that makes the solution of the diffusion equation faster than for the case of a non-orthogonal structured mesh.Thus, the last obtained values for the parameters should be used as initial values in the solution for the ellipsoidal geometry, which is a slower process, since it involves a non-orthogonal structured mesh.However, in order to compare the times of execution of the simulations, it was established 1,000 time steps for all the geometries, and the following initial values for the parameters: Γ Φ = 0.4 W m -1 K -1 and h = 20 W m -2 K -1 .For the finite cylinder, the following values for the parameters were obtained: Γ Φ = 0.5765 W m -1 K -1 , and h = 25.9899W m -2 s -1 .
The distribution of the dimensionless temperature at t = 402 s is shown through the contour plot on the quarter-rectangle which represents the cucumber, and can be observed in Figure 6.
In Figure 6, note that for the f inite cylinder, the difference between the dimensionless temperatures in the center and at the surface, along the radial direction, depends on axial position.In that figure, it is possible to observe that, if the cucumber is considered a finite cylinder, cooling occurs more quickly at the extremities, in the corners in contact with cooling air.

Ellipsoid
For the dimensions of the cucumber, the quarterellipse generating the ellipsoid can be drawn in a paper sheet.If adequate points are chosen at the boundaries of the quarter ellipse, such figure can be digitized, and the non-orthogonal structured mesh can be generated through the software 2D Grid Generation as shown in Figure 7.
After the refinement of the generated initial mesh, the non-orthogonal structured mesh 20 × 40 is obtained and it is shown through Figure 8, which is presented with a rotation of 90°counterclockwise with respect to the Figure 7.As in the case of the finite cylinder, it should be observed that Figure 8 indicates the boundary conditions, taking into consideration the symmetry condition which was imposed to the problem.On the other hand, after the optimization process, the results shown in Figure 9 are obtained.
For the ellipsoid, the following values for the parameters were found after the optimization process: the thermal conductivity is given by Γ Φ = 0.6431 W m -1 K -1 , and the convective heat transfer coefficient is given by h = 23.1744W m -2 s -1 .
The distribution of the dimensionless temperature at t = 402 s is shown through the contour plot on the quarter-ellipse generating the ellipsoid, and it can be observed through Figure 10.
A look in Figure 10 indicates that, for the ellipsoid, the difference between the dimensionless temperatures in the center and at the surface, along the radial direction, depends on axial position.In Figure 10, it is possible to observe that, for the ellipsoidal geometry, cooling occurs more quickly at the axial extremities of the cucumber.
The results obtained in this article for the three studied geometries can be summarized through Table 1.
In order to estimate the uncertainties of the determined parameters for the ellipsoid, 50 virtual experiments were performed as presented by Le Niliot and Lefèvre (2004) and also by Mariani et al. (2009).To this end, the experimental measurements were disrupted with 50 different Gaussian error distributions with zero mean value and standard deviation given by σ = 2σ Φ* (95%) in which σ Φ* is the standard deviation of the simulation.
The average values of the parameters for the ellipsoid and also the uncertainties of them are given in Table 2.
Student's t-test is a statistical indicator that gives the probability P(t) of a parameter be zero, despite its value obtained.The values of «t» and P(t) are given in Table 3 for the parameters obtained through the virtual experiments.
As can be seen in Table 3, values of Γ Φ and «h» are significant due to the fact that for all them P(t) is very close to zero.
Initially, through the results for h and Γ Φ , it is possible to conclude that the infinite and finite cylinder geometries are almost equivalent, and that could be foreseen due to the values of the height and radius of the product.

Discussion
It was used for the three examined cases: 1) the same optimizer, 2) the same method to solve the system of equations, 3) the same initial values in the optimizations processes, 4) the same tolerances for convergence.
Influence of the geometry on the description of the cooling of cucumbers 249

‫ץ‬Φ* ‫ץ‬x
Experimental Numerical Thus, it is possible to conclude that the results in Table 1 are different only due to the different geometries.
In Table 1 it can be noticed that the thermal conductivity value increases when the considered geometry approaches the real shape of the cucumber.The discrepancy between the greater and the smaller values is about 13%.On the other hand, the convective heat transfer coefficient decreases when the considered geometry approaches the real shape of the cucumber.The discrepancy between the greater and the smaller values is also about 13%.Analyzing the statistical indicators presented in Table 1, it can be noticed that the simulation of the cooling kinetics of the cucumber using an ellipsoidal mesh produces the best results.The worst simulation is that regarding the cucumber as an infinite cylinder, but such simulation gives results close to that obtained for finite cylinder as observed by Erdogdu and Turhan (2006), for a cylinder with such dimensions.On the other hand, it should be observed that the optimization process regarding the ellipsoidal mesh is about 66 times slower than in the case of the infinite cylinder.
The focused discussion in this article is referring to the influence of the geometry used in the numerical simulation of the cooling of cucumbers.In this sense, it can be observed through the statistical indicators that, the closer the geometry is to the real shape of the cucumber, the better the obtained results are in the description of the cooling process.Similar results were obtained by Silva et al. (2009b) who studied the influence of the geometry on the numerical simulation of the isothermal drying kinetics of bananas, considering boundary condition of the first kind.Gastón et al. (2003) studied the geometry effect on water diffusivity estimation in broom wheat cultivars, with similar results to the obtained in this article.However, in none of these articles was made an analysis in the sense to answer if the obtained parameters were significant.
The compatibility between the values of the parameters determined through the experimental data and the corresponding values of the parameters obtained from the virtual experiments is worthy of notice.For Γ Φ , the module of the percent discrepancy is 0.98%, while for «h» is 1.00%.
Analyzing the statistical indicators of the fits, it was observed that the more the shape approaches the real shape of the cucumber, the better are the results obtained for the fit of the simulated curve to the experimental data of the cooling kinetics.However, the time demanded in the optimization process of an ellipsoid is substantially larger than the time for the infinite cylinder (about 66 times).Thus, the decision about the geometric model which represents the cucumber should take into consideration the desired cost-benefit relation to be obtained in the simulation.
The value of the thermal conductivity obtained in the present paper for the ellipsoid geometry agrees with the value available in the literature for an agricultural product with moisture content of 96% (wb), and the discrepancy between them is only 3.5%.This information makes it possible to conclude that the obtained convective heat transfer coefficient should be correct; considering a similar discrepancy.Γ Φ (W m -1 K -1 ) 151.0 0.0 h (W m -2 K -1 ) 227.2 0.0 nação de Aperfeiçoamento de Pessoal de Nível Superior) and CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico).
Figure 1.a) Uniform mesh (R = 0.019 m): N control volumes with thickness ∆r, and its nodal points; b) Control volume P and its neighbours to west (W) and to east (E).

Figure 2 .
Figure 2. Cooling kinetics in the centre of the cucumber, considered as an infinite cylinder.

Figure 3 .
Figure 3. Contour plot representing the distribution of the dimensionless temperature on the circular section of the cucumber (R c = 0.019 m), considered as an infinite cylinder, in the instant t = 402 s.
Generation available on the internet at http://zeus.df.ufcg.edu.br/labfit/gg.htm.The points at the boundary were defined in a file to be read by the software which generated the two-dimensional mesh 20 × 40 on the quarter-rectangle, as shown in Figure4.It should be observed that Figure4indicates the boundary conditions, taking into consideration the symmetry condition that was imposed to the problem.On the other hand, after the optimization process, the results shown in Figure5are obtained.

Figure 4 .Figure 5 .
Figure 4. Orthogonal structured mesh 20 × 40 generating the half finite cylinder which represents the cucumber: revolution around the axis x.

Figure 6 .
Figure 6.Contour plot representing the distribution of the dimensionless temperature on the quarter-rectangle generating the half finite cylinder which represents the cucumber, in t = 402 s.

Figure 9 .
Figure 9. Cooling kinetics in the centre of the cucumber, considered as an ellipsoid.

Figure 10 .Figure 8 .
Figure 10.Contour plot representing the distribution of the dimensionless temperature on the quarter-ellipse generating the ellipsoid which represents the cucumber, in t = 402 s.

Table 1 .
Summary of the optimization processes for the three geometries

Table 2 .
Statistical treatment for the 50 virtual results using ellipsoidal geometry