ABSTRACT

This article compares diffusion models used to describe seedless grape drying at low temperature. The models were analyzed, assuming the following characteristics of the drying process: boundary conditions of the first and the third kind; constant and variable volume, V; constant and variable effective mass diffusivity, D; constant convective mass transfer coefficient, h. Solutions of the diffusion equation (analytical and numerical) were used to determine D and h for experimental data of seedless grape drying. Comparison of simulations of drying kinetics indicates that the best model should consider: 1) shrinkage; 2) convective boundary condition; 3) variable effective mass diffusivity. For the analyzed experimental dataset, the best function to represent the effective mass diffusivity is a hyperbolic cosine. In this case, the statistical indicators of the simulation can be considered excellent (the determination coefficient is R² = 0.9999 and the chi-square is χ² = 3.241 × 10^–4).

1. INTRODUCTION

After harvest, the time of fruit conservation under natural conditions is limited to only a few days. Grapes, for instance, belong to the most perishable fruits. They are very susceptible to microbial decay and moisture loss. In general, two mechanisms are used to prolong the consuming time of grapes: cooling and drying. In the second case, the life time of the product is much higher than in the first case, beyond resulting in a food very appreciated in several parts of the world: raisins. The drying process of grapes under natural conditions [1] is generally slow due to the resistance of its skin. In order to increase the drying rate, commonly a pretreatment is accomplished before drying. Several pretreatment techniques involving the use of chemical additives are described in the literature [2-6]. However, the demand for natural foods, without adding chemicals, is increasing in many parts of the world. Therefore, alternatives to chemical additives are used as pretreatment for the production of raisins, for instance: abrasion method [7] and dipping in hot water [6, 8].

For describing the drying process, a mathematical model must be used. Several drying models are available in the literature: empirical models [1-5], and diffusion models, which are very common to describe drying kinetics of grapes [1,3,4-9]. In order to describe the drying kinetics through a diffusion model, the initial and boundary conditions must be known. In general, the appropriate boundary condition is of the third kind [6,7,9], but the first kind one is also found in several research works [1, 3-6,8], especially if some pretreatment is used before the drying process. For some products, like rice, drying does not substantially reduce the volume of the product and the shrinkage can be discarded [10,11]. However, in products such as banana [12] and grape [6,8], the shrinkage is so significant that its effect should not be discarded. But if shrinkage is taken into account, its effect on the diffusivity should also be considered: the effective mass diffusivity should have a variable value throughout the process [6,8].

If a diffusion model is used to describe the drying of a product, the diffusion equation must be solved. Some research works present analytical solutions for the diffusion equation, particularly if the boundary condition is of the first kind [1,3-7]. If an analytical solution is proposed for the boundary condition of the third kind, the series which represents the solution is normally expressed by only the first term, and the process parameters are determined by regression [13,14]. However, this procedure only works well with experimental data referring to Fourier numbers grater than 0.2. For the boundary condition of the third kind, numerical solutions are frequently found in the literature, if the complete drying kinetic must be studied [9,15-20].

In order to use the aforementioned solutions for the description of seedless grape drying, the process parameters must be known. These parameters can be determined through optimization processes, normally by using the inverse method [21,22]. Once the process parameters are known, the simulation of the drying kinetics can be performed, and the moisture content can be determined at an instant t in any position r within the grape.

Since several models are available in the literature to describe the drying of grapes, this article investigates what is the most appropriate one. To this end, five diffusion models, involving the aspects discussed above, were used to describe seedless grape drying. Then, the results for each model were compared with each other and the best model was determined.

2. MATERIAL AND METHODS

The solutions of the diffusion equation to describe seedless grape drying presupposes the following hypotheses: 1) grapes are considered as spheres; 2) grapes are considered as homogeneous and isotropic; 3) the moisture distribution inside the grapes must present radial symmetry and must be initially uniform; 4) the convective mass transfer coefficient is constant; 5) diffusion is the only mass transport mechanism inside the grapes.

2.1. Diffusion Equation and Boundary Condition

Given the hypotheses above, the one-dimensional moisture diffusion equation can be written in spherical coordinates as

(1)

where M is the moisture content (db), r defines the position inside the sphere (m), D is the effective mass diffusivity (m·s^–1), and t is the time (s).

The boundary condition of the third kind may be expressed by the equation:

(2)

where h is the convective mass transfer coefficient (m·s^–1), is the moisture content at the boundary (db), and is the equilibrium moisture content (db).

2.2. Analytical Solution

The analytical solution presented in this article refers to a spherical geometry governed by Eqs.1 and 2, where a constant value for the effective mass diffusivity D is assumed. For a sphere with radius R and an initial moisture content M₀, the solution M(r,t) of Eq.1 with boundary condition given by Eq.2 is [23,24]:

(3)

with

(4)

where are roots of a characteristic equation to be defined later. In Eq.3, M(r,t) is the moisture content (dry basis) at the position r from the centre of the sphere at time t.

The average value of the moisture content at time t is defined by

(5)

The solution of the diffusion equation for the average value of a spherical body is obtained by substituting Eq.3 into Eq.5, and the result is:

(6)

where

(7)

Here, Bi is the Biot number given by

(8)

and are the roots of the characteristic equation for the sphere

(9)

If only the experimental data corresponding to a Fourier number Fo = Dt/R² > 0.2 are analyzed, the drying kinetic is well described by the first term of the sum in Eq.6 [14]. Using the dimensionless moisture content,

(10)

Eq.6 reduces for this situation to

(11)

where is given by Eq.7, with n = 1, and

(12)

Then, the parameters B₁ and A₁ can be determined by fit Eq.11 to experimental data. By substitution of Eq.9 into Eq.7, with n = 1, a transcendental equation is obtained, and the root, corresponding to n = 1, can be determined. Thus, the Biot number is also determined through Eq.9, and the effective mass diffusivity by Eq.12. Once Bi and D are known, Eq.8 can be used to calculate the convective mass transfer coefficient, h.

2.3. Numerical Solution

Eq.1 was also numerically solved by using the finite volume method with fully implicit formulation [20,25]. Figure 1 presents the uniform one-dimensional mesh of the corresponding sphere. The control volumes have a thickness (m) and the control volume number “i” has a nodal point “P”.

Integration of Eq.1 over space () and time () gives the following result for the control volume P:

(13)

where the superscript “0” means “former time t” and its absence means “current time t + Dt”. The indexes “P”, “e” and “w” refer, respectively, to the nodal point, and east and west interfaces of the control volume.

Discretizing Eq.2 gives

(14)

and therefore

(15)

where the subscript b refers to the boundary.

For an internal control volume (control volume from 2 up to N – 1), Eq.13 results in:

(16)

where

(17)

For the control volume 1, due to the symmetry condition (flux zero at the centre), Eq.13 becomes:

(18)

with

(19)

For the control volume N (boundary), the combination of Eqs.13 and 14 results in:

(20)

with

(21)

Considering the numerical solution, the system of Eqs.16, 18 and 20 can be solved for each time step, for instance, by the TDMA method (Press et al., 1996). Note that the coefficients A are calculated only once if the effective mass diffusivity D is constant, whereas B is calculated in each time step because its value depends on, the moisture content of the control volume P at the beginning of each time step. Otherwise, if the effective mass diffusivity is variable, the coefficients A also must be calculated in each time step, due to the nonlinearities caused by the variation of D. In this case, errors due to these nonlinearities can be eliminated by adequate time refinement. As the volume of the fruit is variable during the process, the position r of each nodal point must be recalculated in each time step. Once is numerically determined, the moisture content M_b (at the boundary) can be calculated from Eq.15 at a given time t.

The average value of M at a given time t may be calculated by Eq.5, written in the discretized form [20]:

(22)

where V_j are the volumes of the control volume “j” (m³) and

(23)

is the volume of the sphere (m³).

2.4. Effective Mass Diffusivity

In case of the numerical solution, the process parameter D may be calculated at nodal points from an appropriate relation between D and the moisture content M (or dimensionless moisture content X):

(24)

where “a” and “b” are parameters which fit the numerical solution to the experimental data, and they are determined by optimization.

At the interfaces between control volumes, for example “e” (Figure 1), a harmonic mean given by the following expression must be used to determine [20, 25]:

(25)

which is valid for uniform grids. Note that Eq.25 is also valid for a constant effective mass diffusivity, with a value D. In this case, ,; and Eq.25 becomes.

2.5. Optimization

The expression for the chi-square involving the fit of a simulated curve to the experimental data was used as objective function, and is given by [26,27]

(26)

where is the moisture content measured at the experimental point “i” (db), is the correspondent simulated moisture content (db), N_p is the number of experimental points, is the statistical weight referring to the point “i”. If the statistical weights are not known, they are set equal to 1. In Eq.26, the chi-square depends on, which depends on D and h. If h can be considered as constant and the effective mass diffusivity is given by Eq.24, the parameters a, b and h can be determined through the minimization of the objective function, which can be accomplished in cycles involving the following steps:

Step 1. Inform the initial values for the parameters “a”, “b” and “h”. Solve the diffusion equation and determine the chi-square;

Step 2. Inform the value for the correction of “h”;

Step 3. Correct the parameter “h”, keeping the parameters “a” and “b” constant. Solve the diffusion equation and calculate the chi-square;

Step 4. Compare the last calculated value of the chisquare with the previous one. If the last value is smaller, return to the step 2; otherwise, decrease the last correction of the value of “h” and proceed with step 5;

Step 5. Inform the value for the correction of “a”;

Step 6. Correct the parameter “a”, keeping the parameters “b” and “h” constant. Solve the diffusion equation and calculate the chi-square;

Step 7. Compare the last calculated value of the chisquare with the previous one. If the last value is smaller, return to step 5; otherwise, decrease the last correction of the value of “a” and proceed with step 8;

Step 8. Inform the value for the correction of “b”;

Step 9. Correct the parameter “b”, keeping the parameters “a” and “h” constant. Solve the diffusion equation and calculate the chi-square;

Step 10. Compare the last calculated value of the chisquare with the previous one. If the last value is smaller, return to the step 8; otherwise, decrease the last correction of the value of “b” and proceed with step 11;

Step 11. Begin a new cycle going back to step 2 until the stipulated convergence for the parameters “a”, “b” and “h” is reached.

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 step 2, 5 or 8, the value of the new correction can be multiplied by the factor 2. If the initially informed modest correction does not minimize the objective function, its value can be multiplied by the factor –1 in the next cycle. Note that the steps 8, 9 and 10 are not necessary when the effective mass diffusivity is supposed to be constant. Initial values for the parameters can be estimated through values obtained from similar products available in the literature, or through empirical correlations.

2.6. Software

The software used to determine h and D through numerical solution was developed in an Intel Pentium IV computer with 1 GB RAM. The program was compiled in Compaq Visual Fortran (CVF) 6.6.0 Professional Edition, using the programming option QuickWin Application, in a Windows XP platform. The developed software was also used to draw contour plots at specified times, showing the moisture content distribution within the sphere which represents the grape. Besides these characteristics, the software also draws the numerically simulated curve, which was fitted to the experimental data. The statistical indicators used in the analyses of the obtained results were chi-square, and determination coefficient R² [26,27]. LAB Fit Curve Fitting Software V. 7.2.48 (http://zeus.df.ufcg.edu.br/labfit) was used for the statistical treatment of the data.

2.7. Experimental Data

Experimental data obtained by Esmaiili et al. [6] for the drying of seedless sultana grapes (Vitis vinifera L.) with hot air are explored in the present paper. In order to increase the skin’s permeability to moisture, a pretreatment to be described in the following was used. The grapes were dipped for 15 s in hot water at 95˚C. The temperature of the drying air was 50˚C, its relative humidity was 10%, and its velocity was kept at 1.5 m·s^–1. The mean radius of the grapes was 6.65 × 10⁻³ m at the beginning of the process, the initial moisture content was 3.25 (db), and the equilibrium moisture content was 0.17 (db). During the drying process, the radius decreased according to the experimentally obtained expression

(27)

in which r is obtained in meter.

The dimensionless moisture content data were digitized by using xyExtract Digitizer from http://zeus.df.ufcg.edu.br/labfit/index_xyExtract.htm.

3. RESULTS AND DISCUSSION

Five analytical and numerical models were analyzed to describe drying of seedless grapes, and the results are reported. In each numerical simulation, the domain was divided into 100 control volumes. In order to get an idea about the dispersion of the experimental points in relation to the corresponding simulated values, the following expression was used to define the error at a point “i”:

(28)

3.1. Model 1: Constant Volume and Diffusivity with Equilibrium Boundary Condition

In model 1, the effective mass diffusivity was determined by the optimization algorithm proposed in section 2.5, coupled to the numerical solution presented in section 2. The mass transfer coefficient was kept constant at h = 1 × 10^{+10 m}·s^–1, which corresponds to a boundary condition of the first kind. In this case, the objective function has a single minimum, as noted by Silva et al. [21], and the optimization algorithm is not sensitive to the initial value of D. For the experimental data of seedless grape drying, for example, using several significantly different initial values of D, and imposing the relative tolerance of 1 × 10^–4 with 1000 time steps, the optimization delivers the same value for the effective mass diffusivity, as can be seen in Table 1.

Despite of the fact that the last initial value of D is 5 × 10⁵ times greater than the first, the same effective mass diffusivity was obtained in all optimizations. This is consistent with the results of Silva et al. [21], which observed that there is a single value for the diffusivity that minimizes the objective function for the boundary condition of the first kind. The drying simulation considering constant volume and diffusivity for the boundary condition of the first kind is shown in Figure 2.

The statistical indicators of this simulation are: and. Using Eq.28, each error can be calculated. The graph of the error distribution is shown in Figure 3, which also shows the possibility to fit a polynomial of degree 5 to the errors.

Figure 3 shows that this model presents a biased fit, with a high value for the average error, once the expected value is zero, and therefore the model should be discarded as an option to describe the drying kinetics of seedless grapes for the experimental conditions investigated in this paper.

3.2. Model 2: Constant Diffusivity, Variable Volume with Equilibrium Boundary Condition

Using the numerical solution presented in this paper and imposing a relative tolerance of 1 × 10^–4 for the convergence of the effective mass diffusivity, with the time of drying divided into 1000 steps, m²·s^–1 is obtained from model 2. A simulation using this value for D results in the graph shown in Figure 4.

The statistical indicators for model 2 are:

and. For this simulation, again using Eq.28, each error can be calculated. The

graph of the error distribution is shown in Figure 5, which also shows the possibility to fit a polynomial of degree 5 to the errors.

Model 2 also delivers a biased fit for the error distribution, with a high value for the average error. Comparison with model 1 indicates that the inclusion of the shrinkage has improved the results, but as can be observed, the equilibrium boundary condition is not really adequate to describe the drying kinetics of seedless grapes for the experimental data analyzed in this article.

3.3. Model 3: Constant Volume and Diffusivity with Convective Boundary Condition

For this model, the analytical solution presented in this paper was explored. In order to guarantee Fo > 0.2, Eq.11 was fitted to the experimental data without the first points, giving B₁ = 0.8792, A₁ = 9.822 × 10^–6 s^–1, and consequently, Bi = 4.115, D = 7.11 × 10^{–11 m2}·s^–1, h = 4.40 × 10^{–8 m}·s^–1. A graph of this fit is shown in Figure 6(a), while Figure 6(b) shows the drying simulation including all experimental data points.

The error distribution is shown in Figure 7(a), together with the average value of the errors, while Figure 7(b) shows the fit of a polynomial of degree 5 to the errors.

Figure 7 shows that model 3 also yields a biased fit of the error distribution, with a high value for the average error. The advantage of this model is its simplicity, but it does not consider the strong shrinkage of the grapes. In addition, the first experimental points cannot be properly described by the model. On the other hand, the obtained results can be used as initial values for other optimization processes.