Three-dimensional simulation of isothermal wood drying of radiata pine using effective diffusion coefficients

C.H. Salinas^‡, Y.A. Gatica^†, H.A. Vargas^* and R.A. Ananias^§

^‡ Departamento de Ingeniería Mecánica, Universidad del Bío-Bío, Collao 1202, Concepción, Chile.
^† Doctoral Student, Ciencias e Industrias de la Madera, Universidad del Bío-Bío, Collao 1202, Concepción, Chile
^* Magíster Student, Matemática Aplicada, Universidad del Bío-Bío, Collao 1202, Concepción, Chile
^§ Departamento de Ingeniería en Maderas, Universidad del Bío-Bío, Collao 1202, Concepción, Chile. casali@ubiobio.cl, ygatica@ubiobio.cl, hvargasp@ubiobio.cl, ananias@ubiobio.cl

Abstract— The objective of the present work was to simulate three-dimensional drying of wood, based on the concept of the effective diffusion coefficient. For this, we used a conventional drying process of radiata pine in a controlled environment (dry and wet bulb temperatures of 44 and 36 °C, respectively). This experiment allowed us to gather data on: a) the spatial distribution of moisture content at different drying times, b) the drying curve, and c) the surface emission coefficient. A differential, partial, non-linear, second-order mathematical model was used, and the exponential functions of the effective diffusion coefficient characterized the three-dimensional orthotropic (radial, tangential, and longitudinal) transport of moisture content in the wood. This mathematical model was integrated numerically through the control volume finite element method, which contemplates: a) tetrahedral elements for discretization, b) implicit Euler method for the time differential, and c) Gauss-Seidel with successive over-relaxation method to resolve the linear equation system. We compared these results for the three-dimensional spatial distributions of moisture content after 22, 44, 66, and 88 (h) of drying, and the resulting drying curves were in good agreement with the experiments.

Keywords— Mathematical Model; Wood Drying; Radiata Pine; Effective Diffusion; CVFEM.

I INTRODUCTION

A detailed understanding of the wood drying process and the optimization of its technological aspects require to know the transitory distribution of moisture content (M) inside the wood (Keey et al., 2000). In this context, the present work shown an experimental and numerical methodology that allows postulating a tridimensional model based on effective diffusion coefficient for wide range of M like suggested by Hukka (1999) and analyzed by Chen (2007).

Transport models of M can be grouped according to the phenomenology (Ananías et al., 2009; Salinas et al., 2008) or physical aspects of the transport phenomena (Keey et al., 2000). The latter approach includes classic diffusive models such as those presented by Pang (1997), models based on the thermodynamic of irreversible processes as established by Luikov (1966), and models developed using Whittaker's multiphase approach (Whittaker, 1977).

In particular, the diffusive models stand out for their simplicity and ease of implementation. These models are widely accepted for M below the fiber saturation point (FSP), given the diffusive nature of the transport when M<FSP. Research in this line has been done by Smith and Langrish (2008), Hukka (1999), Pang (1997), and others (Zhan et al., 2007; and Pereira et al., 2011). Above the FSP, diffusive models lose validity due to the predominance of capillarity and permeability in detriment to the diffusive phenomenon (Keey et al., 2000). Thus, researchers have proposed models differentiated by drying stage, such as that proposed by Davis et al., (2002). Nonetheless, diffusive models can be used for M>FSP by incorporating what is known as the effective diffusion coefficient (EDC), as analyzed by Chen (2007) and applied to the simulation of the drying kinetic for M<50% by Defo et al., (2004) and for the entire range of M by Rozas et al., (2009).

The problem of determining the values of the diffusive coefficient (or functions of the type given by Comstock, 1963) can be approached inversely: that is, given a known distribution of M, the diffusion coefficients can be determined (Simpson and Liu 1997; Liu et al., 2001; Olek and Weres, 2007; Kang et al., 2009). This has resulted in several proposals of inverse strategies for obtaining the EDC. The present study is framed within this context; specifically, we propose a three-dimensional simulation of conventional drying of radiata pine for the entire range of M based on the EDC differentiated by orthotropic (radial, tangential, longitudinal) direction. In Gatica et al., (2011), the EDC of the radial and tangential directions are reported. Herein, these will be complemented by the EDC of the longitudinal direction. For purposes of validation, experimental data are generated in the form of three-dimensional distributions of M and drying curves. Based on these experimental data, we determine surface emission coefficients (S_E), which are required by the mathematical model as a boundary condition.

The resulting mathematical model is characterized by a partial, differential, second-order diffusive transport equation with variable, non-linear coefficients (exponential variation with M) exposed to Neumann boundary conditions. This mathematical model is integrated numerically through the control volume finite element method (CVFEM; Baliga and Patankar, 1980), characterized by a mesh based on elements in the form of tetrahedrons and an implicit formulation based on the Euler method. The systems of equations are linearized in the update manner and solved iteratively through the Gauss-Seidel method with SOR (system over-relation). Along with the numerical validations, the results of simulations applied to the conventional drying (44/36 °C) of radiata pine are also reported and compared with similar ones experimentally obtained, that is, three-dimensional spatial distributions of M at drying times of 22, 44, 66 and 88 h, and drying curves.

Specifically, the objectives of this work are: 1) to determine the EDC in the longitudinal direction, complementing the work of Gatica et al. (2011) in the radial and tangential directions; 2) to collect three-dimensional experimental data on drying that allows the validation of the proposed model; 3) to determine S_E coefficients for the three-dimensional case; and 4) to perform an effective three-dimensional simulation of the conventional isothermal drying process of radiata pine wood.

II. MATERIALS AND METHODS

A. Description of the experiment

The experiment consisted of drying radiata pine wood inside a climate chamber (see Fig. 1) using dry and wet bulb temperatures of 44 and 36 °C, respectively, and an average flow velocity of 1.6 m/s in the drying environment. Two runs (E1 and E2) were carried out, using seven samples of planed timber (P1 to P7) sized to 80x80x80 mm and oriented as shown in Fig. 2.

Figure 1. Diagram of experimental equipment: 1) Fan, 2) digital balance (A&D model GF-4000), 3) climate chamber (Heraeus Votsch model Virk 150), 4) data acquisition system (Fluke model Hydra II), 5) computer, 6) dried-oven (Memmert model U-15), 7) precision balance (Boeckel model BPB32), and 8) sawmill (Makita model LS1040).

Figure 2. Placement of the samples for assay E2.

Run E1 consisted of drying two wood samples (P1 and P2). Sample P1 was placed on an A&D scale (model DF4000) that was monitored by computer (experimental drying curve) for the continuous determination of variations in the mass during the drying process. With sample P2, 10 T-type thermocouples were installed in the wood and connected to a data acquisition system (Fluke Hydra II) and a computer in order to monitor the temperature inside the wood. The last one was made to ensure that isothermal condition will be effectively imposes. Run E1 concluded when the variation of mass of P1 tended to zero. The anhydrous mass of P1 was determined using the gravimetric method. For this, the sample was dried in a Memmert U-15 oven at 103 °C for 24 h, according to Chilean standards INN 1984 and 1986. Once the anhydrous mass (m₀) was known, we determined the drying curve, using for a transitory evaluation. This experimental drying curve was used to determine the total drying time at 168 h, and four characteristic drying times (22, 44, 66, 88 h) were selected to evaluate the spatial distributions of M (Run E2).

(1)

where is the moisture content (%), m_w (kg) is the moist mass, and m₀ is the anhydrous mass (kg).

For run E2, samples P3 to P7 were placed inside the climate chamber and subjected to the same drying regime as run E1 (Fig. 2). Sample P3 was connected to the A&D scale, and its loss of mass was monitored continuously; this gave rise to the experimental drying curve. Samples P4 to P7 were used to determine the three-dimensional distributions of M at the four times selected previously. For this, a sample was removed from the climate chamber at each characteristic time and sectioned into eight parts in each principal direction, generating 512 cubic subsamples (approximately 8x8x8 mm) (see After this, each subsample was weighed on a high-precision balance (Boeckel BPB32) and, once the anhydrous mass was obtained as described in E1 for P1, we determined the M of each subsample according to Eq. (1).

Figure 3. Sectioning of samples P4 to P7.

B. Effective diffusion coefficient (EDC) and surface emission (S_E)

The EDC was determined for conventional drying (44/36) of radiata pine in the longitudinal direction, as was done in the radial and tangential directions by Gatica et al, (2011). This was necessary for the proposed three-dimensional orthotropic transport model of M.

The surface emission coefficient (S_E), which models the convection of mass at the wood/environment drying interface, was determined by solving the following inverse problem: given the spatial distribution of M at four characteristic drying times, we determined the coefficients S_E such that the difference between the simulated (M_Sim) and experimental (M_Exp) values, defined by Eq. (2), were minimal.,

(2)

In general, the distribution of S_E was not uniform and had to be modeled. In this case, after an analysis of the experimental data, we opted for differentiated modeling by sides, as shown in the diagram of Fig. 4, parameterized by an average value .

Figure 4. Distribution of S_E.

III. MATHEMATICAL MODEL

The phenomenon of transitory three-dimensional transport of M in wood was modeled through the following differential, partial, non-linear, second-order equation for a non-stationary state, with coefficients of variable diffusion (D) as function of the M and the Neumann boundary condition:

(3)

with where M is the moisture content (%), D is the effective diffusion coefficient (m²/s), and t is time (s).

Considering a volume V of boundary S and external normal direction (n), the initial and boundary conditions were as follows:

	(4)

where M₀, M_S and M_∞ are initial, superficial and equilibrium moisture content (%), respectively.

The diffusion coefficient ( D ) in function of M of the exponential type, similar to that proposed by Hukka (1999) and differentiated by the principal direction (radial, tangential, longitudinal), is given by the following expression:

(5)

where , FSP=0.603-0.001T_d (Bramhall, 1979), T_d dry bulb temperature (K), a_i,b_i with i1,2,3; model parameters.

IV. NUMERICAL MODEL

The mathematical model, given by Eq. (3), was integrated according to the CVFEM. For this, Green's Theorem was applied to the aforementioned equation, resulting in:

(6)

where V is the domain of boundary S, whose external unitary normal vector is .

If the domain V consisted of nl finite volumes (FV) (V^l with l=1,nl) of boundary S^l or nk finite elements (FE) (V_k with k=1,nk) of boundary S_k, we could write:

(7)

In turn, each was composed of three surface polyhedral planes with four vertices, defined based on medians (edges), a centroid of areas (sides), and a centroid of the FE (volume), as shown in Fig. 5 and Table 1. For the case of FV centered on local vertex 1, with external normal unitary vectors equal to and , we have that .

Figure 5. Subdivision of the FE (tetrahedron).

Table 1. Values of S_i

Integration of the transitory term

When considering the local variation ∂M / ∂t in the centroid of the FV as predominant, we have for the transitory term:

(8)

And evaluating the temporal differential in an implicit first-order form (implicit Euler), we obtained the integrated and differentiated values of the transitory term in function of the actual time (t=m Δt) and the previous time (t=(m-1) Δt). That is:

(9)

with being the volume of FV l and (t = m Δt)

Integration of the diffusive term

The diffusive term was spatially integrated over the boundary of the FE k equal to with k=1,nk for the time t=m Δt (omitted for simplicity). Thus, the partial contributions of nk FE, in agreement with that shown in Fig. 5, to the FV centered in the local node l equal to 1, were equal to:

(10)

Similar expressions can be written for the contributions of the FE to the FV centered on local nodes 2, 3, and 4 (l₂, l₃, and l₄).

For the above integration, we defined a linear variation of M within the FE, that is:

(11)

where A, B, C and D are constants defined in function of nodal values of the dependent and independent variables: M_i and (x_i, y_i, z_i) with (i = 1,2,3,4), respectively. That is:

(12)

where:

From which we obtained:

(13)

According to the notation of Eq. , assuming the gradient, the normal and the coefficient D (evaluated in t=(m-1)Δt) as constants on the surface , the integral that represents the diffusive contributions, the surface to the centered FV l₁ is equal to:

(14)

with being the area of the boundary .

Likewise, the contributions of the complement of segments that make up the boundary of the FV V^l were determined.

In summary, the diffusive contributions of the nk FE to the FV centered on l were equal to:

(15)

Thus, when considering the transitory contributions given by (9), added to the diffusive terms given by (16), we get the equation of transitory three-dimensional discreet diffusion of M in t=mΔt for a generic volume l:

with l = 1,nl

(16)

where and .

Therefore, for each value l, we will get an algebraic equation formulated for the average value of M in each FV, configuring a system of nl x nl equations of the form .

V. RESULTS

Taking the experimental data expressed in drying curves and three-dimensional distributions of M as our base, we proceeded to simulate the drying using a three-dimensional diffusion equation with variable coefficients of M as the mathematical model, according to the mathematical model given by (3) or in its numerically integrated form expressed by (15). For this, we used the relations of EDC dependent on M determined by Gatica et al. (2011), complemented with those of the longitudinal direction contributed during the development of the present work. Furthermore, we determined the surface emission coefficients (S_E) according to the distribution shown in Fig. 4.

Table 2 summarizes the values of the parameters that determine these coefficients differentiated by the principal direction and drying stage (CMC: above and under critical moisture content). The values reported for S_E are similar to those published in the specialized literature (Siau, 1984; Pang, 1996). Likewise, the values of EDC for M around the FSP, shown in Fig. 6, are comparable to those published by Davis et al. (2002) for radiata pine and Rozas et al. (2009) for slash pine (Pinus elliottii).

Table 2. Parameters for EDC and S_E.

Figure 6. EDC for M<FSP.

In general it has been found that the EDC tangential was higher than radial, these results can be interpreted by consideration of the variations in the structure of radiate pine, is presumed that there is bigger M migration through the pit cavities on the radial plane of the cell-wall of the radiata pine tracheid's, compared with the rays M migration, this agrees with some results showed by Kang et al. (2008a,b). Of course, this aspect become very important data for industrial processes of drying like explain in Keey (2000).

The experimental results expressed by the drying curve (see Fig. 7) show the classic characteristics of conventional drying, which include the drying process of wood differentiated by three drying stages. These drying stages are delimited by the so-called CMC and FSP (M of around 60 and 30 %, respectively). Figure 7 shows the experimental data, the experimental movable average (obtained in a range of 0.5 h), and the numerical simulation through CVFEM. The experiments data of temperature shows small range of transitory condition (about 30 min) after which the isothermal condition was established with standard deviation less than 2 °C.

Figure 7. Experimental and simulated drying curves (uniform mesh of 80 FV).

Figures 8 and 9 show the spatial distributions (experimental and simulated) of the M within the wood sample at drying times of 22, 44, 66, and 88 h. In particular, Figure 8 shows the average one-dimensional spatial distributions evaluated in the direction of the flow of the drying environment (tangential). After 22 h of drying, the simulation moved away from the experimental M distribution, moving particularly towards the surface; however, as the drying progressed, the approximation improved notably. This could be due mainly to the model, which does not consider all the mechanisms at work in this drying stage. Nevertheless, the approximation is very good in terms of average M, as can be seen in Fig. 7. At this point, we can infer that the method was used correctly to determine the coefficients EDC and S_E. Nevertheless, a better approach to the M distributions in the initial drying phases has yet to be found. We believe that more experimental data are needed to improve accuracy of the results and use different kinds of EDC function in the first stage of drying process.

Figure 8. Average M at 22, 44, 66, and 88 h in the tangential direction.

Figure 9. Simulated (colored sheets) and experimental (symbols) transitory distribution of average M at 22, 44, 66, and 88 h in the radial and tangential planes.

Likewise, Fig. 9 shows the average two-dimensional spatial distributions evaluated every 22 h of drying in the radial and tangential planes. This figure reveals the same phenomenon as Fig. 8, that is, a poor simulation at the onset of drying with high and irregular experimental values of M within the wood.

Figure 10 corresponds to that presented in Figure 9: isolines of M, ie, spatial distributions of two-dimensional averages evaluated every 22 h in the radial and tangential planes. The experimental results are shown as isolines and the numerical simulations as colored sheets. This figure is useful for observing the transitory distributions of M according to the progress of the drying, highlighting the asymmetries of this due to natural and forced convection. In general, dry areas occur above the leading edge, favoring the evolution of the maximum M near the center towards the area below the trailing edge. Later, this tends to return to the center of the sample.

Figure 10. Average experimental (isolines) and simulated (levels) distributions of M at 22, 44, 66, and 88 h in the tangential-radial planes.

V. Conclusions

The three-dimensional numerical modeling done herein offers a satisfactory simulation (average deviation < 2%) of the average behavior of the M during conventional isothermal drying of radiata pine wood. This is based on a mathematical model of transport with EDC, exponentially dependent on M, and differentiated by the orthotropic direction (radial, tangential, and longitudinal) and drying stage (over and under the CMC).

The determination of the EDC, based on one-dimensional transport experiments of M, and the following determination of S_E using an inverse problem scheme are effective for the three-dimensional simulation of the drying process. This is particularly valid for the average values expressed by the drying curve. Nonetheless, large local deviations of the simulation, particularly in the initial stages of drying (M>CMC), which tend to disappear along with the drying, approach the FSP. Among others, phenomena such as cementation, which the present model does not incorporate, could be at work.

The values of EDC obtained for the longitudinal direction e^{-20.3+15.1 M} and e^{-19.77+10.05 M} for M above and below the CHC, respectively) are coherent with those published previously for the radial and tangential directions (Gatica et al., 2011), and with those published by Davis et al. (2002) and Rozas et al. (2009). Moreover, the uniform distribution by sides of the S_E based on its average value is very convenient for the search method and effective for reproducing the effects of asymmetry due to the conditions and flow of the drying environment. Unlike the case of the EDC, it is not necessary to differentiate values of S_E by drying stages, indicating the predominance of the flow conditions of the AS (constants in the present modeling) over the transport conditions within the wood.

The experimental results show that the analyzed drying process is differentiated in three stages: from green wood to the CMC, from the CMC to the FSP, and from the FSP to dry wood. The present modeling partially addresses this phenomenon by differentiating two stages of the EDC: before and after the M. In this sense, better approaches can be generated by modeling the differentiated process in three stages. This difference can be given at the level of the parameters of the function that approximates the EDC and/or the type of function used.

REFERENCES
1. Ananías, R.A., W. Broche, M. Alvear, C. Salinas and R.B. Keey, "Using an overall mass-transfer coefficient for prediction of drying of Chilean coigüe," Wood and Fiber Science, 41, 426-432 (2009).         [ Links ]
2. Baliga, B.R. and S.V. Patankar, "A new finite element formulation for convection- diffusion problems," Numerical Heat transfer, 3, 393-410 (1980).         [ Links ]
3. Bramhall, G., "Mathematical model for lumber drying. II. The model," Wood Science, 12, 22-31 (1979).         [ Links ]
4. Comstock, G.L., "Moisture diffusion coefficients in wood as calculated from adsorption, desorption and steady state data," Forest Products Journal, 13, 97-103 (1963).         [ Links ]
5. Chen, X.D., "Moisture diffusivity in food and biological materials," Drying Technology, 25, 1203-1213 (2007).         [ Links ]
6. Davis, C.P., C.G. Carrington and Z.F. Sum, "The influence of compression wood on the drying curves of pinus radiata dried in dehumidifier conditions," Drying Technology, 20, 2005-2026 (2002).         [ Links ]
7. Defo, M., Y. Fortin and A. Cloutier, "Modeling superheated steam vacuum drying of wood," Drying Technology, 22, 2231-2253 (2004).         [ Links ]
8. Gatica, Y., C. Salinas and R. Ananias, "Modelling Convectional One-Dimensional Drying of Radiata Pine Based on the Effective Diffusion Coefficient," Latin American Applied Research, 41, 183-189 (2011).         [ Links ]
9. Hukka, A., "The effective diffusion coefficient and mass transfer coefficient of nordic softwoods as calculated from direct drying experiments," Holzforschung, 53, 534-540 (1999).         [ Links ]
10. Kang, W., W.Y. Chung, C.D. Eon and H. Yeo, "The effect of opening on combined bound water and water vapor diffusion in wood," J. Wood Sci., 54, 343-348 (2008a).         [ Links ]
11. Kang, W., W.Y. Chung, C.D. Eom and H. Yeo, "Some considerations in heterogeneous nonisothermal transport models for wood: a numerical study," J. Wood Sci., 54, 267-277 (2008b).         [ Links ]
12. Kang, W., Y.H. Lee, W.Y. Chung and L.X. Xu, "Parameter estimation of moisture diffusivity in wood by an inverse method," J. Wood Sci., 55, 83-90 (2009).         [ Links ]
13. Keey, R., T. Langrish and J. Walker, Kiln drying lumber, Springer, N. York (2000).         [ Links ]
14. Luikov, S.A.V., Heat and mass transfer in capillary porous bodies, Pergamon Press, Oxford (1966).         [ Links ]
15. Liu, J., W. Simpson and S. Verrill, "An inverse moisture diffusion algorithm for the determination of diffusion coefficient," Drying Technology 19, 1555-1568 (2001).         [ Links ]
16. Olek, W. and J. Weres, "Effects of the method of identification of the diffusion coefficient on accuracy of modeling bound water transfer in wood," Transp Porous Med, 66, 135-144 (2007).         [ Links ]
17. Pang, S., "External heat and mass transfer coefficients for kiln drying of timber," Drying Technology, 14, 859-871 (1996).         [ Links ]
18. Pang, S., "Relationship between a diffusion model and a transport model for softwood drying," Wood and Fiber Science, 29, 58-67 (1997).         [ Links ]
19. Pereira, W., L. Duarte, C. Silva and P. Nascimento, "Optimization and simulation of drying processes using diffusion models: application to Wood drying using forced air at low temperature," Wood Sci Technol, 45, 787-800 (2011).         [ Links ]
20. Rozas, C., I. Tomaselli and E. Zanoelo, "Internal mass transfer coefficient during drying of softwood (Pinus elliotti) boards," Wood Science and Technology, 43, 361-373 (2009).         [ Links ]
21. Salinas, C., R.A. Ananías and P. Ruminot, "Modelación Fenomenológica del Secado de Madera de Pino Radiata a Alta Temperatura," Maderas. Ciencia y tecnología, 10, 207-217 (2008).         [ Links ]
22. Siau, J.F., Transport Processes in Wood. Springer, New York (1984).         [ Links ]
23. Simpson, W.T. and J.Y. Liu, "An optimization technique to determine red oak surface and internal moisture transfer coefficients during drying," Wood and Fiber Science, 29, 312-318 (1997).         [ Links ]
24. Smith, S. and T. Langrish, "Multicomponent solid modeling of continuos and intermittent drying of Pinus radiata sapwood below the fiber saturation point," Drying Technology, 26, 844-854 (2008).         [ Links ]
25. Whittaker, S., "Simultaneous heat, mass, and momentum transfer in porous media: A theory of drying," Advances in Heat Transfer, 13, 119- 203 (1977).         [ Links ]
26. Zhan, J., J. Gu and Y. Cai, "Analysis of moisture diffusivity of larch timber during convective drying condition by using Crank's method and Dincer's method," Journal of Forestry Research, 18, 199-203 (2007).         [ Links ]

Received: October 25, 2011.
Accepted: May 9, 2012.
Recommended by Subject Editor Cecil Coutinho.