Drying of Clay Slabs during the Falling Rate Period: Optimization and Simulation of the Process Using Diffusion Models

Mass transfer in clay slabs during drying is described using diffusion models with equilibrium (model 1) and convective (model 2) boundary conditions. Drying experiments using hot air were performed for clay slabs with initial moisture content of 0.11 (dry basis, db) at 50, 60, 70, 80 and 90 °C. For this initial moisture content, drying occurred in falling rate. For the slabs, the one-dimensional solutions referring to models 1 and 2 were coupled to optimizers to determine the process parameters in each temperature. Thus, equations can be obtained to express these parameters as a function of the drying air temperature. The analyses of drying processes have indicated that there is a good agreement between the experimental results and the corresponding simulations using the model 2, for all temperatures.


Introduction
In order to describe a drying process, a mathematical model is normally used.For ceramic products, for instance, several drying models are available in the literature.Among these models the following can be cited: empirical models (Sander et al., 2003); semi-empirical models (Sander et al., 2003); models based on Darcy law (Chemkhi & Zagrouba, 2008;Mihoub & Bellagi, 2012) and diffusion models (Sander et al., 2003;Chemkhi & Zagrouba, 2005;Farias et al., 2012;Silva et al., 2012aSilva et al., , 2012b)).Some of these models will be highlighted in the following.Sander et al. (2003) studied three drying models for clay slabs: two empirical models and one diffusion model with constant diffusivity.The initial moisture content ranged from 0.22 to 0.28 (db).Drying experiments were carried out at several temperatures (from 30 to 70 °C); and it was observed two periods for the process.In this work, the liquid diffusion model was applied to describe the whole process.According to the authors, a good agreement between experimental results and the results predicted by the diffusion model has been found.On the other hand, Silva et al. (2012a) described the drying process of roof tiles, shaped from red clay, using diffusion models.In their study, the researchers used the first 25 terms of the series that represents the analytical solution of the diffusion equation with boundary condition of the first and third kind.To determine the process parameters, the authors coupled the analytical solutions to optimizers based on the inverse method.They concluded that the boundary condition of the third kind is better than the first kind in the description of roof tiles drying.On the other hand, the researchers pointed out the limitations of the models, which use constant process parameters, and suggested a numerical solution for a more rigorous study.Silva et al. (2012b) described water migration in clay slabs during drying, emphasizing two distinct periods: constant and falling rate.For the first period, the process was described by an empirical equation and, for the second one, drying was described through the solution of the diffusion equation with boundary condition of the third kind.To eliminate some restrictions found in the literature, a three-dimensional numerical approach was used to describe the process.This diffusion model enables to eliminate some common restrictions such as constant process parameters and representation one-dimensional of the geometry, among others.However, the time spent in optimizations processes can be considered very high.In this context, the proposal of a quick model to describe drying of clay slabs deserves consideration.
The main objective of this article is to describe drying of clay slabs during the falling rate period, using one-dimensional diffusion models, including optimization and simulation of the process.To this end, optimizers were coupled to the one-dimensional analytical solutions with boundary condition of the first and third kind.Then, using experimental datasets, the process parameters were obtained and the results were analyzed.

Diffusion Equation
As is known, the diffusion equation can be used to study drying of clay materials during the falling period (Silva et al., 2012b).For infinite slabs, the one-dimensional diffusion equation is written as (Luikov, 1968;Crank, 1992) where M is the moisture content in dry basis (db); D is the effective mass diffusivity (m 2 •s -1 ); t is the time (s) and x is the Cartesian coordinate of position (m).Figure 1 presents an infinite slab, highlighting its thickness L x in the direction of the main mass flux, where the distributions of moisture will be analyzed.
Figure 1.Infinite slab highlighting the thickness L x where the moisture distribution will be analyzed In this article, the assumptions to solve Equation (1) are: (1) diffusion is the only transport mechanism inside the slab; (2) the thickness L x does not vary during drying; (3) the initial moisture distribution is uniform; (4) the infinite slab is considered homogeneous and isotropic; (5) the process parameters are constant during drying, (6) the origin of the axis x is located at the central point of the slab.

Boundary Condition of the Third Kind
For an infinite slab, the convective boundary condition is defined in the following way: in which h is the convective mass transfer coefficient (m•s -1 ); M(x,t) is the moisture content (dry basis, db) in a position x at time t; M eq is the equilibrium moisture content (db); L x is the thickness of the infinite slab.
For an infinite slab with the uniform initial moisture content M 0 and the boundary condition defined by Equation (2), the analytical solution M(x,t) of Equation ( 1) is given by (Luikov, 1968;Crank, 1992): where the parameter A n is given by sin sin( ) and μ n are the roots of the characteristic equation for the infinite slab: The parameter Bi is the mass transfer Biot number, defined as follows: The expression for the average moisture content ( ) M t at time t is given by: where the parameter B n is given by ( ) Equation ( 5) is a transcendental equation which can be solved for a specified mass transfer Biot number.An auxiliary program was written in Fortran, using the bisection method, and the first 16 roots of Equation ( 5) were calculated for 469 values of mass transfer Biot numbers from Bi = 0 (which corresponds to an infinite resistance of the water flux at the surface) to Bi = 200 (which practically corresponds to an equilibrium boundary condition).

Boundary Condition of the First Kind
For the equilibrium boundary condition, the solutions of Equation ( 1) are also given by Equations ( 3) and ( 7).However, for an infinite mass transfer Biot number (that is the equilibrium boundary condition) Equation ( 5) is given by which implies in: The coefficients A n defined by Equation ( 4) are now known and the local moisture content at any instant t, ( , ) M x t , can be calculated by Equation (3).On the other hand, when Bi   Equation ( 8) becomes: With the coefficients B n calculated by Equation ( 11), the average moisture content at instant t, ( ) M t , can be calculated by Equation (7).

Optimizer for the Boundary Condition of the First Kind
In order to determine the effective mass diffusivity, the optimization algorithm developed by Silva et al. (2009) and (2012a) will be briefly reviewed here.The objective function is defined by the chi-square obtained through the fit of the analytical solution to the experimental points (Bevington & Robinson, 1992;Taylor, 1997) exp where is the i th experimental point of the average moisture content; ana i M (D) is the average moisture content given as a function of D at the same point i and is calculated from Equation ( 7) supposing boundary condition of the first kind; 2 i 1 / σ is the statistical weight of the experimental average moisture content at the point i and N p is the number of experimental points.The chi-square depends only on the effective mass diffusivity, since the mass transfer Biot number is known for the prescribed boundary condition ( Bi   ).If 2 i 1 / σ was not obtained from the experiment, and is therefore unknown, the same value for the statistical weight, for instance 2 i 1 / σ 1  , should be attributed to all experimental points.The algorithm to determine D is given in the following: Initially, ana i M (D ) is calculated according to Equation (7), using a value stipulated for the number of terms N (for instance, 200) and an initial value of D (for instance, 10 -20 ).For the boundary condition of the first kind, the values for n  and n B are given by Equation ( 10) and ( 11), respectively.The values of ana i M (D ) are used to calculate χ 2 according to Equation (12).Then, the value of D is doubled giving a new value for ana i M (D ) , which is now used to calculate a new χ 2 .The value of χ 2 obtained in the first time is compared with the second value obtained for χ 2 .If the statement the second χ 2 is lower than the first χ 2 is satisfied, the program will continue the previously described procedure, otherwise the program will finish this procedure because a minimum value for chi-square is found.The last three values for D and χ 2 are then recorded.The recorded D interval, which contains the minimal value of χ 2 , is then divided into n equal parts, and each part of this interval is defined as a step.Thus, the χ 2 of each step is determined, and a more refined interval that contain the minimum is determined.This procedure is repeated until a convergence criterion is satisfied.

Optimizer for the Boundary Condition of the Third Kind
The optimizer presented in Section 2.3 can be applied for the convective boundary condition if Equation ( 12) is adequately adapted.In this case, the objective function for the minimization becomes (Silva et al., 2010): in which is the average moisture content calculated by Equation ( 7) with n  and n B given, respectively, by Equation ( 5) and ( 8).In this case, Equation ( 7) is a function of D and Bi.Consequently, the objective function given by Equation ( 13) depends on D and also on Bi or, equivalently, on h.However, for a specified Biot number, Equation ( 13) is only a function of D. Thereby, the optimizer presented in Section 2.3 can be applied for the minimization of Equation ( 13), determining an optimum value of D for a specified Bi.The minimization procedure can be repeated for all 469 Biot numbers established in the Item 2.2.1 to cover the domain from 0 to 200.After a scan of all Biot numbers, the lowest chi-square among the 469 minimums is identified.

Model 1
This model involves the analytical solution of the diffusion equation with boundary condition of the first kind, and this solution will be coupled to the optimizer presented in Section 2.3.The number of terms of Equation ( 7) was established as 200.

Model 2
This model involves the analytical solution with boundary condition of the third kind coupled to the optimizer presented in Section 2.4.The number of terms of the analytical solution given by Equation ( 7) was established as 16.

General Considerations
Since D and h are determined by the optimization, Equation ( 3) is used to determine the distribution of the moisture content inside the infinite slab as a function of the position x for a given instant t.Moreover, Equation ( 7) is used to simulate the drying kinetics of the slabs.As mentioned earlier, for an optimization process involving boundary condition of the first kind, the number of terms stipulated for the summation of Equation ( 7) was 200, while for boundary condition of the third kind this number was equal to 16.A detailed study about the number of terms of the summation for the two boundary conditions and the cut-off errors is presented by Silva et al. (2012a).Such study shows that the larger the Biot number, the greater must be the number of terms of the series so that the truncation errors can be considered negligible, particularly for the initial instants of drying.
The two optimizers coupled to Equation (7) were developed in a computer Intel Pentium IV with 2 GB (RAM).The source codes were compiled by Compaq Visual Fortran (CVF) 6.6.0Professional Edition, using the programming option QuickWin Application under the Windows Vista platform.For the two boundary conditions analyzed in this article, the convergence criterion stipulated for the determination of the diffusivity was 1×10 -15 .On the other hand, the determination coefficient R 2 and the chi-square χ 2 were used as statistical indicators for the analysis of the results.
During drying, as the thermal diffusivity is much greater than the mass diffusivity, it was considered that the drying process occurs at isothermal conditions.Thus, the effective mass diffusivity as a function of the drying air temperature can be expressed through an Arrhenius-type equation (Nastaj & Witkiewicz, 2009) exp ( . ) In Equation ( 14), D 0 is a pre-exponential factor (m 2 •s -1 ), E a is the activation energy (J•mol -1 ) and R is the universal gas constant (8.314J•mol -1 K -1 ).The convective mass transfer coefficient can be also given as a function of the temperature of the drying air by an Arrhenius-type equation: where a and b are fitting parameters.

Experimental Datasets
Initially, the clay (41.80 % SiO 2 ; 22.09 % Al 2 O 3 ; 10.43 % Fe 2 O 3 ; 3.79 % MgO; 2.69 % K 2 O; 1.25 % CaO and others) was dried at 110 ºC, de-agglomerated in a ball mill and sieved through 80-mesh sieve (180m).Then, water was added to the powder to form a homogeneous mixture with a moisture content of 0.11 (dry basis) approximately.After, the obtained mass was kept standing for 24 hours to improve the plasticity.The clay slabs were obtained by extrusion, and the mass shaped as slab was cut by hand.The temperature and the relative humidity of the room in which the kiln was placed were on average, respectively, 25 o C and 75 %.
The behavior of the clay samples during drying was investigated by monitoring and recording the changes in weight of the test samples during drying using a kiln at the temperatures of 50, 60, 70, 80 and 90 o C.Among the five experimental datasets, three of them were used to determine the process parameters D and h by optimization, and the general information about these three datasets are given in Table 1.
Table 1.Drying air temperature (T), initial moisture content (M 0 ), equilibrium moisture content (M eq ) and thickness (L x ) of the clay slabs used in the optimizations 50.0 0.1116 0.0162 6.05x10 -3 60.0 0.1086 0.0105 6.02x10 -3 90.0 0.1046 0.0024 5.92x10 -3 The moisture contents were measured by the gravimetric method, and the process took place until the mass reached its equilibrium value.At the end of each drying, the kiln temperature was set at 105 ºC and the slabs remained there for 24 h, enabling the measurement of dry matter.The thickness L x of the slabs was measured with a digital caliper.
With the process parameters calculated by optimization at 50, 60 and 90 o C, Equation ( 14) and ( 15) can be used to estimate D and h at other temperatures.In order to evaluate these results, two other experiments at the temperatures of 70 and 80 o C were performed, and their general information is available in Table 2.
Table 2. Drying air temperature (T), initial moisture content (M 0 ), equilibrium moisture content (M eq ) and thickness (L x ) of the clay slabs used to test the expressions obtained for the process parameters 70.0 0.1090 0.0077 5.96x10 -3 80.0 0.1129 0.0048 6.13x10 -3 The experimental datasets for the five drying kinetics in falling rate are presented in Figure 2.

Model 1: Results
Using the optimizer presented in Section 2.3 and the experimental datasets characterized in Table 1 and given in Figure 2, the results obtained for the model 1 are presented in Table 3 for the temperatures of 50, 60 and 90 o C.
With the effective mass diffusivities obtained through the model 1, the drying kinetics for the temperatures of 50, 60 and 90 o C are shown in Figure 3.

Optimizations
For the boundary condition of the third kind, the process parameters are determined by optimization at 50, 60 and 90 o C, and the results are presented in Table 4.
The results obtained in Table 4 for the process parameters enable to present the drying kinetics at 50, 60 and 90 o C, as can be seen in Figure 4.

Arrhenius-Type Equations
The statistical indicators supposing convective boundary condition (Table 4) are better than those obtained supposing equilibrium boundary condition (Table 3).Thus, using nonlinear regression to fit Equation ( 14) to the data (T,D) given in Table 4, the following result is obtained: Table 3. Process parameters and statistical indicators obtained by optimization using model 1 Comparing Equation ( 16) with Equation ( 14), the activation energy is obtained: E a = 14.9 kJ•mol -1 .On the other hand, Figure 5 shows the behavior between the effective mass diffusivity and the drying air temperature.h (m min -1 ) 4.92×10 -4 90 °C Bi 14.0 Using nonlinear regression to fit Equation ( 15) to the data (T, h) given in Table 4, the following relationship is found: 2 1914 h 9.57 10 exp (T 273.15 Figure 6 shows the relationship between h and T in the interval 50-90 o C.  16) and ( 17) Since the Equations ( 16) and ( 17) are determined, the optimization process is no longer needed to calculate D and h.This calculation can be performed directly from Equations ( 16) and ( 17).In order to confirm this affirmation, the pairs D and h were calculated through Equations ( 16) and ( 17) for the drying air temperatures of 70 and 80 °C.The following results were obtained: D=7.84×10 -8 m 2 min -1 and h=3.62×10 -4 m min -1 for T = 70 o C; and D=9.09×10 -8 m 2 min -1 and h=4.24×10 -4 m min -1 for T = 80 o C. Thus, using Equation ( 7), the simulations of the drying kinetics can be presented together with the experimental datasets characterized in Table 2 and given in Figure 2 (70 and 80 °C).The results can be observed in Figure 7.
Figure 7. Simulations of the drying kinetics using parameters determined by Equations ( 16) and ( 17 With the process parameters D and h known for the five drying air temperatures, the moisture content distribution within the slab can be determined at a given instant t, using Equation (3).As an example, for t = 20 min, Figure 8 presents the moisture distribution for all drying air temperatures studied herein.

Discussion
While the use of a three-dimensional numerical solution coupled to an optimizer can take several hours to determine D and h (Silva et al., 2012b), with a three-dimensional analytical solution and the optimization algorithms presented here this time is reduced to few minutes (Silva et al., 2012c).On the other hand, the use of these algorithms coupled to one-dimensional analytical solutions reduces the optimization time to few seconds (Silva et al., 2009).The simplicity and the time reduction justify the large quantity of works in the literature that describe water migration in porous solids using approximated geometries, generally one-dimensional, and analytical solutions together with some type of optimization algorithm (Cunningham et al., 2007;Ruiz-López & García-Alvarado, 2007;Hacihafizoglu et al., 2008;Mariani et al., 2008).
According to Silva et al. (2012a), the use of only the first term of series that represents the analytical solution can lead to significant cut-off errors, particularly in the description of the initial instants of the process.Therefore, in this article, it was used 16 terms of the series that represents the solution for the boundary condition of the third kind, and the obtained results can be considered good.
With respect to geometry, the samples are parallelepipeds, but the thickness is much smaller than the other two dimensions.In this sense, the one-dimensional models are reasonable, despite these models overestimate the effective mass diffusivity due to the consideration of the flux only in the two largest surfaces of the parallelepiped, ignoring the fluxes in the other four smaller areas.Nevertheless, the results obtained are reasonable, and the optimization times are small (only a few seconds).
An inspection in Tables 3 and 4 indicates that the results obtained with the boundary condition of the third kind presents better statistical indicators than the results obtained with the boundary condition of the first kind.However, this last boundary condition is also reasonable.A new inspection in Table 4 indicates the reason for this fact.The mass transfer Biot numbers for the three temperatures (50, 60 and 90 o C) are approximately 14, and this number can be considered very high, indicating only a little resistance at the surfaces of the slabs.
According to Silva et al. (2012b), once the optimizations have been used for some temperatures of the drying air, enabling the obtaining of Equations ( 16) and ( 17), these optimizations can be substituted by the idea presented here.Equations ( 16) and ( 17) can be used, instead the optimizations, to describe a drying kinetics at other temperature between 50 and 90 o C.This means that the parameters D and h can be calculated with no optimization process, and then the simulations of the drying kinetics can be performed.
Once the parameters of the drying kinetics have been calculated, Equation (3) makes it possible to determine the moisture content at any point within the slab, at an instant previously stipulated.For instance, the contour plots with the distribution of the moisture content at instant t = 20.0 min are shown in Figure 8, for all temperatures of the drying air.As is known, information on how the moisture content (and temperature) is different in different positions within the slabs is important because such differences generate stresses (Collard et al., 1992;Mihoubi & Bellagi, 2012) that can damage the product.

Conclusion
The statistical indicators enable to conclude that the model 2 is better than model 1 in the description of drying of clay slabs.The presented optimization algorithm has permitted to determine the process parameters D and h for several drying temperatures in the range 50-90 o C.
It was possible to determine Arrhenius-type expressions to calculate the parameters D and h as a function of the drying air temperature.The values calculated from the Arrhenius equations make it possible to simulate the drying process in a chosen temperature within the interval 50-90 o C, with no need of new optimizations.The comparison between these simulations and experimental datasets enables to affirm that the methodology presented in this article is a useful resource in substitution to the optimization processes.
Model 2 permits to predicts the moisture distribution in a given instant, enabling to study stresses that can damage the product.

Figure 2 .
Figure 2. Experimental datasets for the falling rate period with drying air temperature T at: (a) 50 o C; (b) 60 o C; (c) 70 o C; (d) 80 o C; (e) 90 o C

Figure 3 .
Figure 3. Drying kinetics of clay slabs using the model 1 at temperatures (a) T = 50 o C; (b) T = 60 o C; (c) T = 90 o C

Figure 4 .
Figure 4. Drying kinetics of clay slabs using the model 2 at temperatures (a) T = 50 o C; (b) T = 60 o C; (c) T = 90 o C

Figure 5 .
Figure 5. Effective mass diffusivity for model 2 as a function of the drying air temperature

Figure 6 .
Figure 6.Convective mass transfer coefficient for model 2 as a function of the drying air temperature Figure 7. Simulations of the drying kinetics using parameters determined by Equations (16) and (17) with the drying air temperatures at: (a) T = 70 o C; (b) T = 80 o C