Use of a Combined Constant Rate and Diffusion Model to Simulate Kiln-Drying of Pinus radiata Timber

This paper presents the use of a combined constant drying-rate and diffusion model to simulate the drying of Pinus radiata timber under kiln-drying conditions. The constant drying-rate and diffusion coefficients of the model, which control the drying rate of individual pieces of timber, were determined from calibrating the model against the experimental drying curves obtained under the kiln-drying conditions. The experimental drying curves were obtained from the gravimetric measurements of the moisture content of timber during kiln drying. Statistical relationships were developed for the constant drying-rate and the diffusion coefficients of the model as functions of kiln temperature and the dry basis density of timber. To determine the effects of variability of timber, a simulation scheme was developed based on the model, the probability distribution of the density of timber, the equations for the constant drying-rate coefficient and the diffusion coefficient. The model and the associated simulation method provides a simple way to estimate the drying time of a stack of timber using parameters determined from experimental results for the specific timber kiln.


List of Symbols t
= the elapsed drying time, s. M e = the equilibrium moisture content, kg kg -1 .M i = the initial moisture content, dry basis, kg kg -1 .D x and D y = the diffusion coefficients in the x and y directions, m 2 s -1 .a and b = the dimensions of the timber, m.

ϕ
= the relative humidity, decimal.T db = the dry bulb temperature, °C.THICK = the timber thickness, m.MC MAX = the moisture content at total saturation, decimal.ρ = the basic density of the board, kg m -3 .CR = the constant drying rate coefficient, s -1 .MC = the dry weight basis moisture content, decimal.

Introduction
Timber drying models have been extensively researched in order to develop a model that can accurately predict the drying of timber from high moisture contents (in excess of 100 %, dry basis) to levels of approximately 10 %.The movement of moisture through wood is a complex process.Drying requires a lower concentration of water vapor in the air surrounding the wood than in the wood at the wood/air interface.When this occurs moisture can be evaporated from the wood surface.The loss of moisture from the surface of wood causes moisture to move from the wetter interior to the drying surfaces.In the initial stages of drying, the free water replaces moisture evaporated from the surface, as less energy is required to enable moisture movement.Once all of the free water has been removed, the fiber saturation point (FSP) had been reached (usually at moisture contents of 25-30 %, dry basis).
The FSP represents the point at which the cell walls are saturated with water, but no water is present in the cell lumen or intercellular spaces.
To enable bound water movement the hydrogen bonds holding the moisture in the cell walls need to be broken (Pratt and Turner 1986).Numerous timber-drying models have been developed during the last 30 years and these have generally been classified into three categories: empirically based, diffusion based, and combined heat and mass transfer based.Empirical models are based on the assumption that the rate of drying is proportional to the instantaneous difference between the moisture content of the material and the equivalent equilibrium moisture content (M e ) of the drying air.The diffusion models have been widely used to describe the movement of moisture in wood at moisture contents below the FSP.Difficulties arise in the measurement and selection of the diffusion coefficient.The heat and mass transfer models are based on the well-known Luikov equations and also contain experimentally determined coefficients.The models are further being extended as the knowledge of wood anatomy, fluid movement through porous media, and mass and heat transfer in cellular material are integrated.Kamke and Vanek (1995) investigated 12 diffusion and heat and mass transfer drying models and found that the high variability recorded between the models and experimental measurements was largely due to coefficient selection.The more sophisticated models did not perform any better than the simpler models if the physical property data was inadequate.Therefore, a simpler model can adequately perform well in most practical situations encountered in the kiln drying of timber.The aim of this research is to develop such a simple model that can be used to predict drying times for various timber dimensions and kiln conditions.The model can be used to predict the drying time for a stack of timber (with different densities and moisture contents) provided that expressions for the constant drying rate and diffusion coefficients have been determined for the particular kiln.

Constant and Falling Rate Drying Periods of Timber
In this research the drying of timber was considered to consist of two periods.During the early stages of drying the material consists of so much water that liquid surfaces exist and drying proceeds at a constant rate (Henderson and Perry 1976, Moyne and Basilico 1985, Waananen et al. 1993).Constant drying rates are achieved when surface free water is maintained and the only resistance to mass transfer is external to the wood.Since the moisture source for the surface is internal moisture, constant drying rates can only be maintained if there is sufficient moisture transport to keep the surface moisture content above the FSP.If this level is not maintained then some of the resistance to mass transfer becomes internal and neither the drying rate nor the surface temperature remains constant (Kayihan 1987) and drying proceeds to the falling rate period.
Drying in the falling drying-rate period involves the movement of moisture to the surface of the material and the removal of the moisture from the surface.Due to the complex physical nature of wood and the wood-water bond, the application of classical diffusion theory to the moisture movement problem involves a certain degree of approximation (Moschler and Martin 1968).For this study the falling drying-rate period was modeled using an analytical solution of Fick's second law (Equation 1) which determines the moisture content at a point (x,y) in 2-dimensional space at some time t (Crank 1979).
The analytical solution assumes isothermal air conditions (temperature, velocity and relative humidity) during the kiln cycle; 2-dimensional mass transfer; moisture transfer is mainly by transient diffusion through the wood; the diffusion coefficient is constant; and, the evaporation rate is proportional to the difference between the actual moisture content on the surface and the moisture content required to maintain equilibrium with the surrounding atmosphere (Mounji et al. 1991).
The diffusion coefficient is an important coefficient in the use of diffusion based models.It has been reported (Comstock 1963, Luikov 1966, Simpson 1993, Chen et al. 1994, Pang 1997) that the diffusion coefficient is influenced by the drying temperature, density and moisture content of timber.The diffusion coefficient of water in cellophane and wood substance was shown by Stamm (1959) to increase with temperature in proportion to the increase in vapor pressure of water.Stamm and Nelson (1961) also observed that the diffusion coefficient decreased with increasing wood density.The exact nature of the relationship between moisture content and the diffusion coefficient is still not known.Simpson (1993) recorded an exponential relationship for the diffusion coefficient in the moisture content range of 5 to 30 %.Other factors affecting the diffusion coefficient that are yet to be quantified are the species (specific gravity) and the growth ring orientation.
Several assumptions were made in the development of our model: 1) The constant drying-rate period was assumed to exist until such a time that the slope of the diffusion model (using an initial moisture content of 40 %) was equal to the slope of the drying curve during the constant drying rate period (Fig. 1).
The assumption that the initial moisture content for the diffusion model was equal to 40 % was based on this being approximately the average moisture content of timber when the surface moisture content reached the FSP, approximately 30 % (Kininmonth and Whitehouse 1991, Walker 1993, Pang 1995).
2) The diffusion coefficient in the analytical solution was assumed to be constant in the x and y directions.Literature has suggested that the ratios of radial and tangential diffusion coefficients vary for different tree species (Luikov 1966).The radial diffusion coefficient of New Zealand Pinus radiata has been estimated to be approximately 1.4 times the tangential diffusion coefficient (Kininmonth and Whitehouse 1991).
3) The summation portion of Equation 1 was simplified to include the terms representing the summation from i = 0 to 3 and j = 0 to 3, to reduce computing time.This decision was based on determining the effect of additional terms on the average moisture content values obtained from the model.After 2 hours a difference of 0.3 % was observed between the average moisture content determined from the model using i = 0 to 3 and j = 0 to 3 and that using i = 0 to 20 and j = 0 to 20.The only time when the number of terms appears to be significant is in the determination of the initial average moisture content.4) An empirical relationship was required for M e .
M e for a material depends on the physical nature of the material and the temperature and humidity of the surrounding atmosphere.M e values are usually published in table form, however to enable a versatile model an empirical relationship was required.The relationship that Pang (1995) used to determine the M e of Pinus radiata timber in his research was used in this work. where: 5) During kiln drying the humidity inside the drying chamber is determined by controlling the dry and wet bulb temperatures.To enable the M e to be calculated from Equation 2 the relative humidity was required to be determined from the dry bulb and wet bulb temperatures using a series of psychometric equations (ASHRAE ... 1981).

Development of a Single Board Drying Model
Initially the model was developed to allow constant drying rate and diffusion coefficients to be determined for a series of boards dried under various kiln conditions.These coefficients were determined by trial and error to adjust the resulting drying curves to represent those measured using gravimetric techniques.The intention of this was to determine values of both coefficients for different dimension timber dried at various kiln settings and to apply statistical techniques to obtain relationships for the coefficients for the use in the model.Constant drying rate and diffusion coefficients were determined for a total of 80 boards at dry bulb temperatures from 90 to 140 °C, with timber thicknesses of 25, 40 and 50 mm as show in Table 1.The table of sample boards shows that gravimetric results could not be collected for some thickness timber at some of the kiln settings.The absence of measurements particularly at a kiln setting of 140/90 may result in a lack of accuracy when predicting coefficients for these temperatures.The 50 mm timber also dried at 140/90 was also all low density timber which have higher diffusion and constant rate coefficients.
The average moisture content of a particular board at any time can be determined using Equations 3 and 4. The decision on which equation to use requires determining whether drying is proceeding in the constant rate or falling rate period.This decision is based assumption 1 (that drying proceeds at a constant rate until the slope of the drying curve in the falling rate is equal to the constant rate coefficient).If drying is proceeding in the constant rate period then the moisture content is found from: If drying is proceeding in the falling rate period then the average moisture content is determined by integrating Equation 1 over the dimensions of the timber and dividing this by the cross sectional area of the timber.This can be represented by:

Results
The trial and error fitting of constant drying-rate and diffusion coefficients gives an accurate indication of the drying curve for individual boards at temperatures from 70 to 140 degrees Celsius.Drying curves, such as those shown in Figs 2  and 3, were obtained at several kiln temperatures, timber dimensions and initial moisture contents.These results show that this type of model can be used to accurately represent the drying curves for differing dimension timber drying at different kiln settings.

Examination of the Diffusion and Constant Drying-rate Coefficients
An empirical relationship for the diffusion and constant drying rate coefficients using the board density, thickness and drying temperature was developed using regression analysis.Using the single board model, the diffusion and constant drying rate coefficients were determined from gravimetric data for 80 pieces of Pinus radiata timber.The coefficients were determined by calibrating the model to obtain a drying curve that matched the gravimetric measurements that had previously been obtained.The coefficients were determined for three timber thicknesses (25, 40 and 50 mm), four discrete drying conditions (90/ 60, 110/70, 130/80 and 140/90) and each board had a different density (ranging from 346 to 548 kg m -3 ).Initially a correlation matrix was calculated to examine whether any of the variables were highly correlated with each other.This analysis showed that both the dry bulb and wet bulb temperatures were highly correlated, as a result the wet bulb temperature was omitted from the analysis.

Diffusion Coefficient
A General Linear Model (GLM) in statistical analysis was used to assess the significance of  the main effects and interactions of thickness, dry bulb temperature and density in the diffusion coefficient.Initially all of the terms of interest and the interaction terms were included in the analysis.The analysis of variance showed that the only significant interaction (p < 5 %) was the T DB ×THICK (dry bulb temperature × thickness) term.This indicated that the other interaction terms could be omitted from the analysis as they had little effect on the parameter determination.
The significance of the T DB ×THICK interaction term was of particular interest as no other literature suggests a relationship between timber thickness and the diffusion coefficient.It must be noted that a full selection of the temperature and thickness combinations was not available.
The significance of the thickness (and interaction term) in the model appears to be an important term although the physical situation does not justify this inclusion.It was decided that the interaction term should not be included in the model until a greater number of samples had been obtained to further investigate the effect.Completing the analysis without the interaction term gave a model where the thickness term became non-significant (p > 5 %) hence reducing the model to including T DB and the board density.The regression analysis for the diffusion coefficient using the dry bulb temperature and the density had an R 2 of 49.9 %.

Constant Rate Coefficient
A similar process was used to determine the model for the constant rate coefficient as mentioned above for the diffusion coefficient.The inclusion of all of the interaction terms in the initial GLM analysis showed that there were no significant interactions between T DB , thickness and density.A regression procedure of the main terms was used to obtain a model to predict the constant rate coefficient.The regression analysis showed that all of the terms were highly significant (p < 1 %) and the model had an R 2 of 76.0 %.It must be emphasized that the deterministic relationships obtained in this section were for a specific timber kiln.This kiln is a small experimental kiln that will differ markedly from larger industrial kilns.This work does show that, as reported previously, the timber density and drying temperature can be used to predict the diffusion coefficient.What has not been examined previously is the possible effect of the timber thickness on the diffusion coefficient.The analysis also showed that the constant drying-rate coefficient could be predicted from the dry bulb temperature, timber density and thickness with a good degree of accuracy.

Development of a Monte-Carlo Simulation Model
The Monte-Carlo model was an extension of the single board drying model in that the drying of a charge of timber (in this case 200 boards) was simulated.Parameters such as the board density and initial moisture content were generated using probability distributions determined from the data collected for this purpose.The relationships for the constant drying rate and diffusion coefficients determined from the single board model were used to model the drying of each board in the charge.The Monte-Carlo simulation determines the moisture content of each board at each time interval (usually a 15 minute increment) and evaluates whether 90 % of the boards are deemed dry (having average moisture contents within 2 % of the target moisture content, 12 % dry weight basis).On the completion of drying the moisture contents of each board are written to a file and these can be further analyzed to obtain the average drying rate and the distribution of moisture content within the charge at the completion of drying.
The densities (basic densities) of the individual boards were assumed to be normally distributed.The mean was assumed to be 450 kg m -3 , based on the measurements in the single boards.
A standard deviation was 30 kg m -3 which represented approximately 99 % of the timber densities being in the range of 350-550 kg m -3 (±3 standard deviations).
The initial moisture content of each board was determined from the density of the generated piece of timber.The initial moisture content was determined by calculating the maximum moisture content (moisture content at total saturation) using the Equation 7. The maximum moisture content was then adjusted to represent the moisture loss before the timber is put into the kiln using a uniformly generated value between 10 and 50 %.These values that simulated moisture loss (during felling, sawing, transportation and stacking) prior to the boards being placed in the kiln.

Results
The Monte-Carlo simulation was conducted at kiln settings of 90/60, 110/70, 120/70, 130/80 and 140/90 with 25, 40 and 50 mm timber.Dur-ing each simulation the moisture content distribution on the 200 boards was examined (Fig. 4).As expected the variation of the moisture contents reduced as drying proceeded and at the completion of drying 75 % of the boards had moisture contents of between 8.5 and 12.0 %.
At the completion of drying, the moisture contents of a charge of 50 mm, Pinus radiata dried at 110/70, were approximately normally distributed (with a mean of 11.0 % and a standard deviation of 1.7).While this result suggests that there is evidence of over dried boards, schedules for kiln drying of Pinus radiata include an equalisation phase to minimize the final moisture content variation (Producing ... 1996).During the equalisation period the wet bulb temperature is raised (to achieve an M e of 2 percent below the target moisture content) to restrict over drying of the fastest drying boards while maintaining the drying of the slower drying pieces.For example, the fastest option for drying in at 90/60 in an accelerated conventional kiln with a target moisture content of 10 % is to maintain the dry bulb at 90 °C and raise the wet bulb to 80 °C during equalisation.
A summary of the time taken for 90 % of the timber to have a moisture content of < 14 % is provided in Table 2.The table shows that the drying times differ considerably with timber thickness.The drying times obtained from the model were consistent with the drying times obtained  from another set of boards which were not used to determine the diffusion and constant drying-rate coefficients, but dried under the same drying conditions.This analysis shows that the use of a combined constant rate and diffusion model can be used to predict the drying times for different thickness timber at various kiln settings.

Conclusions
Due to the complex nature of the physical processes occurring as timber dries and the high variability in wood structure, modeling timber drying is a difficult task.The model presented in this paper used gravimetric measurements to determine relationships between density and drying temperature to predict diffusion and constant drying-rate coefficients.A Monte-Carlo simulation was conducted at several kiln settings to examine the coefficients and how differences in timber density effect the moisture content variation at the completion of drying.The simulations showed that this method of modelling the drying of timber was accurate.The differences in the predicted drying times were small and resulted from: -temperature changes during the heat up period for the kiln not being represented; -loss of heat from the kiln while gravimetric samples are removed; -estimation of the diffusion coefficient had a moderate R 2 and further factors may need to be examined; and -gravimetric measurements were not available for each board thickness and kiln setting.

Discussion
The gravimetric measurements used to determine the relationships for the constant rate and diffusion coefficients were obtained for a 8 m 3 conventional kiln.During these drying experiments no measure was made of the drying quality.Drying at these temperatures will result in some loss in timber quality which is why industrial drying includes equalisation and reconditioning phases in kiln schedules to minimise this loss.These changes in the kiln conditions could not be incorporated into this model due to the assumptions required with the analytic solution used to model the falling rate period.
The simulations showed that this method of modeling the drying time of timber was accurate, however improvements could be made by: -analyzing more boards using the single board model to improve the coefficient prediction; -examining the inclusion of more variables in the coefficient determination, such as fan direction and velocity and timber dimensions; -examining the use of a diffusion equation that allows changes in the kiln conditions during the kiln cycle to be included to account for changes during the heat up, drying, equalisation and reconditioning phases of the kiln cycle; -examining the use of the simulation method in different types of kilns.

Fig. 1 .
Fig. 1. Outline of the method used to model the drying of a piece of timber, showing the constant dryingrate and falling drying-rate periods.

Fig. 3 .
Fig. 3. Measured and predicted drying curves for 25 and 40 mm Pinus radiata timber drying in a kiln at 130/80.

Fig. 4 .
Fig. 4. Drying curve showing the predicted median and quartiles of the 200 50-mm Pinus radiata boards at a kiln setting of 110/70.

Table 1 .
Number and density range of samples used to determine relationships for the diffusion and constant rate coefficients.
Note: NA measurements were unavailable at these settings.

Table 2 .
Predicted average drying times in hours for Pinus radiata timber using different kiln settings.Experimental values are in parentheses.
Note: NA Measurements were unavailable at these kiln settings