# Modelling of airflow and wood drying inside a kiln: a comprehensive approach.

Abstract

The quality and uniformity of kiln-drying is strongly influenced by the airflow distribution inside the kiln and by the moisture and heat transport inside the wood. A new comprehensive mathematical model that combines both processes is presented here. A method for obtaining numerical solutions for this wood drying model is described and results on different grids are compared. The wood drying model predictions are compared with experimental data. A solution procedure for coupling airflow and wood drying sub-models is described. Results for a test case provide an opportunity to evaluate model performance and to suggest many practical applications. This paper is part of a larger effort to develop computer models for wood kilns that can be applied to industrial processes. The model has a large potential to become a useful tool for kiln design and a valuable part of a training program for kiln operators.

**********

Wood kilns are designed to remove moisture from the lumber pieces by supplying heat to the lumber to evaporate liquid, which then is carried away by the air circulation. Kiln-drying is influenced by heat and mass transfer between the airflow and wood as well as by the complex moisture transport processes happening inside each piece of lumber. Kiln designers and operators face a difficult task of balancing product quality with a demand for the reduction of drying time. When limited experimental data are available, it is quite natural that mathematical modelling would become an integral part of the research and design processes. For the past two decades an extensive amount of work has been done to improve our understanding of the processes inside a wood kiln and to develop models simulating kiln-drying. Most of the previous research focused on modelling either airflow inside the kiln (Hua et al. 2001) or wood drying processes inside the lumber (Kamke and Vanek 1994). The previous research achieved success in th e development of separate models and prepared the way for the development of a model that combines both airflow and wood drying sub-models into one comprehensive model. The purpose of this paper is to present that combined model, to show its application for a test case and to suggest some of its potential for future simulations.

Physical and mathematical models

Airflow model

A schematic diagram of a typical wood kiln is shown in Figure 1. Airflow in the kiln is created by axial fans located above the wood load. Flow passes around a top load baffle and enters into the entrance plenum, then passes through wood packages, providing the lumber with heat and taking moisture Out, until reaching the middle plenum. From the middle plenum, air goes into the second stack of lumber and out to the outlet plenum; from there it flows back to the fans. Heating coils can be located at the top and in the middle; vents, which allow moist air to leave the kiln and fresh air to enter it, are usually located at the top on the suction and pressure sides. Fans are usually reversible and motors can be two-speed or be equipped with variable speed drives. The average velocity in the channels between wood pieces usually varies from 2 to 6 m/sec., and the distance between wood rows from 16 to 38 mm. For that velocity range, density can be assumed constant and body forces such as gravity can be neglected. As most industrial kilns are very long and fans are spaced at regular intervals along the length, the air profile can stabilize before entering the plenum, and a two-dimensional approximation is justified. For turbulence modelling, the standard k-[epsilon] model with wall function treatment is chosen for its reliability and known robust performance. This airflow model has been developed and extensively used by the University of British Columbia to analyze the influence of the kiln geometry on the airflow pattern. Some of the findings regarding different plenum arrangement, various sticker thickness, and flow between an array of boards were presented previously (Hua at al. 2001).

Wood drying model

During the last 10 years, much effort was directed toward the development of wood drying models. As a result, there has been significant progress in the understanding of water movement during the drying process. An extensive overview of drying models has been made by Bian (2001). For implementation into the comprehensive kiln model, the "three variable" model based on an air and water mass balance, together with the energy balance developed by Perre and Turner (1996), was chosen. The main model variables are moisture content (MC), temperature, and gaseous pressure. The MC is defined as a ratio of the weight of all water in the wood to the ovendry weight of the wood. Three phases (gas, liquid, and solid) are considered, and the gas phase is separated into vapor and air parts. The solid and gas phases are always continuous, while the liquid phase is continuous only until it reaches the fiber saturation point, that is the point where the only liquid left in the wood is that held hygroscopically by the cell walls , and called bound water. Thermodynamic equilibrium is assumed, and the gas phase behaves like an ideal gas. The mass balance equations are written for all three phases.

[partial]([[rho].sub.l] + [[rho].sub.b])/[partial]t + [nabla]([[rho].sub.l]x[V.sub.l])+

[nabla]([[rho].sub.b]x[V.sub.b])=-<M>-

<[M.sub.b]> liquid phase [1]

[partial][[rho].sub.a]/[partial]t + [nabla]([[rho].sub.a]x[V.sub.a])=0

air part of gas phase [2]

[partial][[rho].sub.v]/[partial]t + [nabla]([[rho].sub.v]x[V.sub.v])=<M> +

<[M.sub.b]> vapor part of gas phase [3]

where:

[rho] = density (kg/[m.sup.3])

V = velocity (m/sec.)

<M> = mass flux (kg/[[m.sup.3]*sec.])

l = liquid water

b = bound water

a = air

v = vapor

An overbar implies averaging over the corresponding phase.

The enthalpy conservation is written:

[partial]([rho]H)/[partial]t + [nabla]([rho]xVxH)=

[nabla]x([K.sub.eff][nabla]T) [4]

where:

H = enthalphy (J)

T = temperature (K)

[K.sub.eff] = effective conductivity tensor (W/[m*K])

Coupling between airflow and wood drying models

As was mentioned earlier, past modelling work done for wood kilns focused mainly on either one of the two main aspects of kiln performance: airflow distribution or wood drying. This separate approach neglects the enormous influence that one of the processes has on the other. Non-uniformity of the airflow makes drying conditions different for different places inside the kiln. Due to different wall stresses on the wood air interface, heat and mass transfer varies from one wood piece to another and along the piece itself. Wood drying can influence the airflow by changing its temperature and composition. Moreover, moisture removal from different wood pieces is affected by the location of the piece in the stack and by neighboring pieces: the drying rate of a board upstream would influence the drying of a downstream board, because air conditions will change along the passage. For a more complete understanding of these processes inside the wood kiln during drying, these interactions must be taken into consideration. From a practical point of view, a comprehensive model could be used directly by kiln design engineers, because such a model is able to consider the whole kiln-drying process, providing answers on various "what if" questions. The comprehensive approach makes it easier to analyze real industrial problems, bridging the gap between theory and practice.

There are several options for coupling the two models. One is to have a complete transient airflow calculation including temperature and vapor concentration running simultaneously with wood drying calculations. Although this is not impossible, given today's computer power, it is highly uneconomical due to the computer time and resources required. Airflow parameters such as velocities, pressure, and turbulence change insignificantly during drying at a constant fan speed, but their solution still demands a lot of computational time. The logical alternative would be a second approach in which the airflow is calculated before drying starts. Once done, velocities, pressure, and turbulence parameters are fixed at their steady state values. The only parameters that change will be temperature and vapor density. Obviously, this brings a significant simplification to the problem without seriously compromising the main principles. But simplifications need not stop there. It is known that for the accurate solution of the continuity and momentum equations for the airflow, a fine mesh is required; while, at the same time, the convection-diffusion problem can tolerate a much coarser grid. This discussion leads to a third option in which the temperature and vapor density will be solved on a coarser grid, ensuring that most of the physics are captured for the least amount of computational cost.

For physically correct coupling, it is important to treat accurately the interface between the wood and the airflow. A typical interface is shown in Figure 2. Heat from the air is transferred to the wood; part of it is used for water evaporation from the wood surface and another part is conducted inside. The mass flux of water vapor into the air consists of liquid evaporated from the surface and vapor flux coming from inside the wood. The boundary conditions, which result from continuity and conservation requirements, are shown below.

[[rho].sub.v] [\.sub.w] = [[rho].sub.v] [\.sub.a], T[\.sub.w] = T[\.sub.a], P[\.sub.w] = P[\.sub.a] [5]

[m.sub.v] [\.sub.w] + [m.sub.ev .j] = [m.sub.v] [\.sub.a], q[\.sub.w] + [m.sub.ev.l] x

[DELTA][h.sub.vap] = q[\.sub.a] [6]

where:

[[rho].sub.v] = vapor density (kg/[m.sup.3])

P = pressure (Pa)

[m.sub.v] = vapor mass flux (kg/[[m.sup.2]*sec.])

[m.sub.ev.l] = mass flux of evaporated liquid (kg/[[m.sup.2].sec])

[DELTA][h.sub.vap] = latent heat of evaporation (J/kg)

q = heat flux (W/[m.sup.2])

w = wood side of the interface

a = air side of the interface

To determine the heat and mass transfer on the air side, simple equations are used:

[m.sub.v] [\.sub.a] = [h.sub.m] x ([[rho].sub.v] [\.sub.w]- [[rho].sub.v]), q[\.sub.a] =

h x (T[\.sub.w] - T) [7]

where:

[h.sub.m] = local mass transfer coefficient (m/sec.)

h = local heat transfer coefficient (W/[[m.sup.2]*K])

Values without a subscript are outside flow values.

Because the flow between wood pieces is essentially a turbulent flow over a flat plate, the Chilton-Colburn analogy between mass and momentum transfer can be employed for obtaining the local mass transfer coefficient from the airflow calculations.

[Sh.sub.x]/[Re.sub.x] x Sc - 1/3 = [C.sub.[f.sub.x]]/2 [8]

where:

[Sh.sub.x] = [h.sub.m] x x/D =

local Sherwood number

[Re.sub.x] = V x x/v =

local Reynolds number

Sc = v/D = Schmidt number

[C.sub.[f.sub.x]] = [[tau].sub.w]/1/2[rho]x[V.sup.2] =

The skin friction coefficient

Analogously, for local heat transfer coefficient, the Reynolds analogy can be used:

[St.sub.x] x [Pr.sup.2/3] = [C.sub.f]/2 [9]

where:

[St.sub.x] = h/[rho]x[c.sub.p]xV =

local Stanton number

Pr = [c.sub.p] x v x [rho]/k = Prandtl number

x = the distance from the board edge (m)

D = water vapor diffusivity through air ([m.sup.2]/sec.)

[nu] = kinematic viscosity of air ([m.sup.2]/sec.)

k = air thermal conductivity (J/[m.K])

[c.sub.p] = specific air heat capacity (J[kg.K])

[[tau].sub.w] = wall shear stress (Pa), which is calculated from the turbulence parameters based on wall functions

Numerical solution

Method for solving initial airflow equations

For solving the airflow Reynolds averaged Navier-Stokes equations, a coordinate transformation with domain segmentation is employed. Decomposition of the kiln into several computational blocks allows different types of grids with various densities to be used in different kiln parts. A sectional pressure correction procedure is used to accelerate convergence. A detailed description of the method can be found in previous study (Nowak 1998).

Method for solving wood drying equations

After some rearrangement and application of phenomenological laws, it is possible to use Equations [1] to [4] to obtain a system of partial differential equations for three primary variables: MC, temperature, and gaseous pressure. As a wood board has a rectangular shape, it is logical to use an orthogonal grid for the numerical solution. In the drying process, high gradients of drying variables occur near the exchange surfaces; in order to resolve these properly, a small mesh size near the wood--air boundary is required. An additional reason to have a sufficiently small mesh near boundaries lies in the implementation of boundary conditions. Properties at the interface are assumed to be the same as at the bordering cell centers. To justify this approach, small boundary cells are necessary. At the same time, in the middle of the wood piece, a much bigger grid size can be tolerated. Under the uniform grid approach, an excessively refined mesh has to be used everywhere, which make computation expensive and ineffi cient.

A control volume cell centered technique together with implicit Crank-Nicolson time discretization is used to obtain a system of linear algebraic equations. Central difference approximations are used for values at the cell interfaces. Coefficients at the cell interface are computed as weighted averages from neighboring cell values. After discretization, a block-structured non-symmetric matrix was obtained. This matrix is extremely sparse; it has zero elements everywhere except along five diagonals and on these diagonals it contains small 3 x 3 matrices. For a coupled solution of this system, the Generalized Conjugate Residuals method (GCR) (Wesseling 2001) is used. Although this method was originally developed for single matrices, its generalization for coupled systems of equations with block-structured matrices is straightforward. GCR, as with all of the methods from the conjugate gradient family, requires preconditioning for efficient performance. Here the ILU(0) (Incomplete lower upper factorization) preco nditioning (Ferziger and Peric 1999) was used and proved to be very efficient for this application. Implementation of this method allows more flexibility in the choice of model input parameters, which is necessary, especially when the wood drying model is combined with the airflow.

Method for solving transport equations in the air side

To obtain the solution for moisture and temperature distributions in the airflow, the conservation equations were discretized on the coarse grid with a control volume approach and fully implicit time advancement. First order up-wind approximations for cell interface values were used because they provide a high degree of robustness (Ferziger and Peric 1999). As some of the computational domains on the air side can be two-dimensional, care must be taken when interpolating values of the flow parameters from the fine grid, on which the airflow equations were solved, to the coarse grid, on which the vapor density and temperature equations were solved. Mass conservation has to be preserved under such interpolation, otherwise unphysical solutions may emerge. In this model, the velocity along the x-axis and all turbulence parameters are interpolated by means of a mass flux weighted scheme, and, to obtain velocity along the y-axis, the continuity equation was solved on the coarse grid. Systems resulting from discretiz ation of one-dimensional equations can be simply solved by marching methods, while for linear systems obtained from two-dimensional equations, the Bi-CGSTAB method was used (Van der Vorst 1992).

Method for providing coupling between airflow and wood drying

The coupling procedure is designed to achieve energy and vapor mass conservation inside the kiln, keeping computational cost at an acceptable level. Calculations start from the airflow model and flow parameters are calculated on the fine grid without heat and mass transfer. After convergence is reached, heat and mass transfer coefficients are obtained from Chilton-Colburn (Eq. [8]) and Reynolds (Eq. [9]) analogies and a new coarse grid is generated for the airflow domain. After solution initialization, the heat and mass transfer equations in the air side are solved simultaneously with the wood drying model equations and boundary conditions (Eq. [7]) applied at the wood air interface. Inlet and outlet plenums are solved with a one-dimensional approximation. Two different approaches are analyzed for the middle plenum solution: full mixing approximation or solution of the 2D-conservation equations. There are many control options that can be implemented for the modelling of the top inlet area: some of the moistur e can be removed, heat can be added to the airflow, and different time schedules can be studied. It is easy to see that many computations are independent of each other during a single time step. This leads to the idea of using parallel computing to further reduce computational time. The present computer program is parallelized by using the Message Passing Interface (MPI) that allows a cluster of workstations to be used as a single parallel computing system with distributed memory (Gropp at al. 1994). With this approach, wood rows are distributed evenly between available computers that communicate between each other at the end of every iteration, sending heat and mass fluxes through their boundaries and receiving back updated outside flow conditions.

Results and discussion

Wood drying model verification

For a computational investigation of the wood drying, a rectangular wood board 0.2 by 0.05 m is chosen. All of the wood properties are taken from the work by Perre and Degiovanni (1990). At first, the model behavior with different types of grids is studied. The dry-bulb temperature is kept at 80[degrees]C, wet-bulb at 60[degrees]C, initial temperature is 20[degrees]C and MC is 0.99. Uniform grids 30 x 14, 160 x 40 and 240 x 60 are compared with a non-uniform grid 30 x 14 (Fig. 3). For comparison, the MC values in the middle of the wood piece are plotted against time. The results shown in Figure 4 confirm the effectiveness of the nonuniform grid. It is easy to see that the non-uniform grid allows the use of substantially less grid cells with the same accuracy. To validate the wood drying model, its predictions are compared with the results of the experiments conducted by Forintek Canada Corporation (Oliveira 2000). Two cases are investigated: in the first, the initial MC in the piece of lumber is 0.66, and in the second, 0.78. Dry- and wet-bulb temperatures are set to the same values measured during the experiment. A plot of the average MC change with time is presented in Figures 5 and 6. Relatively good agreement between calculated and measured values is achieved. Although there are some quantitative deviations, the qualitative behavior of the model is in the good agreement with the experiment.

Sample coupling case results and analysis

The kiln configuration for this numerical experiment contains many of the features of an industrial kiln. Its height is 6.65 m, width is 4.7 m, inlet and outlet plenums width is 1.04 m, and middle plenum width is 0.76 m. It contains four wood packages divided into front and back stacks, one package sitting on top of the other. Each package contains 290 lumber pieces of 2 by 6's (38.1 by 139.7mm) assembled into 29 rows with a distance between rows (sticker thickness) of 19 mm. Recent research (Sun 2001) shows that the influence of small gaps between the boards on heat and mass transfer coefficients is insignificant. This supports the assumption that all boards in the same row have exactly the same height and no spaces exist between them. Based on that, each row of boards is modelled as the same piece that can be divided in separate boards during the analysis stage. The velocity at the inlet is assumed to be 2m/sec., which corresponds to an average velocity value of about 4mlsec. between wood pieces. Dry- and w et-bulb temperatures are kept constant and equal to 80[degrees]C and 60[degrees]C, respectively. Initial MC is 1.0. Pressure is assumed to be equal to its atmospheric value, an acceptable assumption because all pressure variations in the airflow side, obtained from the airflow solution, are substantially smaller than the pressure variations inside the wood during drying.

For a full kiln, modelling is simplified by dividing its geometry into computational sub-domains (Fig.7). The kiln is separated into top inlet area, top outlet area, inlet, middle and outlet plenums, and separate air channels between rows of wood pieces. In spite of the simplifications, all main physical features of the flow inside the kiln are retained. The grid for the airflow solution is shown on the left part of Figure 8. It should be adequate to provide all desired details of the flow. For the coupled solution, two choices of grid are investigated: in the first one, shown in the middle of Figure 8, the air side is one dimensional throughout all sub-domains, in the second one (right side of Fig. 8), a two-dimensional grid is generated in the middle plenum. Both cases are calculated to compare their performance. During drying, data for MC, temperature, and pressure inside the wood pieces and temperature and vapor density in the air are collected.

The results showing the average MC distribution around the kiln are compared in Figure 9. Several important observations can be made from this plot. First, it is easy to see a reduced drying rate at the top front package (rows 1 to 29). This can be explained by lower velocity and transfer coefficients in the upper channels, due to the kiln geometry. Second, less MC at the very top and the very bottom row of boards can be noticed. This occurs because the air passing through the channels above the top and below the bottom boards has lumber only on one side, so this portion of air can maintain a higher temperature and lower MC compared with the air coming through other channels. From this figure, it can be clearly seen that proper resolution of the middle plenum is critical for obtaining realistic results for the back packages, especially for their top section, although at the same time it has almost no influence on the front packages. This can be explained with the help of temperature distribution contours arou nd the kiln, for the case where the 2D solution was obtained in the middle plenum. (Fig. 10) Hot air from the big gap between the packages goes up and creates highly non-uniform air inlet conditions for the top back package. From this analysis, it can be concluded that the only case when the middle plenum solution can be neglected is when there is no gap or little gap between top and bottom packages.

The average MC distribution from wood piece to wood piece along the row is presented in Figure 11. Lumber pieces are numbered from left to right, that is in the same direction as the airflow In every row of lumber, lower MC occurs in the boards that are closer to the airflow entrance. This is a consequence of higher transfer coefficients and higher difference between wood and air properties at the beginning of the air channel. The non-uniformity of kiln-drying inside the row of lumber is shown in Figure 12. The moisture distribution inside some of the wood rows is shown as a surface in a space where x- and y-axes represent the width of stack and the board thickness of lumber and the z-axis corresponds to the MC. Not only the non-uniformity between the rows and between different boards can be analyzed with this model, as has already been discussed, but also non-uniformity within the board itself can be analyzed. In every wood board, a moisture gradient from the center to the surface can be observed. This gradi ent is responsible for stresses, occurring in the wood during drying, which can cause wood defects such as checking and warping. Modelling can provide critical kiln locations as well as critical times where and when high moisture gradients are taking place. After determination of the problem, parametric studies for different time schedules and wood arrangements can be conducted and practical solutions can be obtained.

Obviously, all of these comments apply to the specific kiln configuration that is chosen for the numerical investigation and are intended not as general design recommendations but as an example of modelling capability.

Conclusions

A comprehensive wood kiln model has been presented and numerical methods showing an efficient strategy for the coupled kiln solution have been described. Experimental data for wood drying were compared with the model. A test kiln-drying case was investigated on two different grids and results were analyzed. The capability of the model to simulate heat and mass transfer processes inside the real industrial kiln was confirmed. The kiln model can be an important tool for the design of new kilns, for the analysis of problems encountered in existing kilns, and optimization of their operation. Computational methods can considerably reduce costly experimental time and increase the understanding of drying processes by allowing researchers to look inside the kiln and "inside" the wood piece when a kiln is in operation. Kiln performance can be analyzed by changing initial conditions and drying schedules, reversing the airflow, implementing different control sequences, varying positions of control sensors, and by trying many other alternatives. Model results can be incorporated into kiln control systems, providing relationships between values, measured at some kiln locations, and important kiln performance parameters. This model constitutes a substantial advance towards the development of virtual kiln simulators, which can be used for kiln design and for training kiln operators.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

Literature cited

Bian, Z. 2001. Airflow and wood models for woodkilns. M.A.Sc. thesis. Univ. of British Columbia, Vancouver, BC, Canada.

Gropp, W., E. Lusk, and A. Skjellum. 1994. Using MPI: Portable Parallel Programming with the Message-Passing Interface. MIT Press, Cambridge, MA.

Hua, L., E. Bibeau, P. He, I. Gartshore, M. Salcudean, Z. Bian, and S. Chow. 2001. Modelling of airflow in wood kilos. Forest Prod. J. 51(6):74-81.

Ferziger, J. and M. Peric. 1999. Computational Methods for Fluid Dynamics. Springer-Verlag, NY. 389 pp.

Kamke, F.A. and M. Vanek. 1994. Comparison of wood drying models. In: 4th IUFRO Inter. Conf. on Wood Drying. IUFRO, Vienna, Austria. pp. 1-21

Nowak, P. 1998. GEO-MGFD: A system of computer programs for steady-state flow of variable density fluid in the complex three dimensional regions using multi-grid and multi-segment techniques. Technical rept. Dept. of Mechanical Engineering, Univ. of British Columbia, Vancouver, BC, Canada.

Oliveira, L. 2000. Wood Drying Scientist, Forintek Canada Corporation. Personal communication.

Perre, P. and A. Degiovanni, 1990. Simulation par volumes finis des Transferts Couples en Milieu Poreux Anisotropes: sechange du bois a basse et a haute temperture. (Control volume formulation of simultaneous transfers in anisotropic porous media: simulations of softwood drying at low and high temperature.) Int. J. Heat Mass Transfer 33(1l):2463-2478. (In French, English abstract).

_____ and I. Turner. 1996. The use of macroscopic equations to simulate heat and mass transfer in porous media: some possibilities illustrated by a wide range of configurations that emphasise the role of internal pressure. In: Mathematical Modelling and Numerical Techniques in Drying Technology. I. Turner and A.S. Mujumdar, eds. Marcel Dekker, NY.

Sun, Z.F. 2001. Numerical simulation of flow in an array of in-line blunt boards: mass transfer and flow patterns. Chemical Engineering Sci. 56(2001): 1883-1896.

Van der Vorst, HA. 1992, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 13(2):631-644.

Wesseling, P. 2001. Principles of Computational Fluid Dynamics. Springer-Verlag, Berlin. 654 pp.

The authors are, respectively, Research Engineer, Graduate Student, Professor, and Professor, Dept. of Mechanical Engineering, Univ. of British Columbia, 2324 Main Mall, Vancouver, BC, Canada V6T 1Z4; and Wood Drying Scientist, Forintek Canada Corp., 2665 East Mall, Vancouver, BC, Canada V6T lW5. The authors would like to acknowledge financial help from the Natural Sciences and Engineering Research Council of Canada and from Forintek Canada Corp. This paper was received for publication in January 2002. Article No. 9432. [C]Forest Products Society 2003. Forest Prod. J. 53(4):46-54.

The quality and uniformity of kiln-drying is strongly influenced by the airflow distribution inside the kiln and by the moisture and heat transport inside the wood. A new comprehensive mathematical model that combines both processes is presented here. A method for obtaining numerical solutions for this wood drying model is described and results on different grids are compared. The wood drying model predictions are compared with experimental data. A solution procedure for coupling airflow and wood drying sub-models is described. Results for a test case provide an opportunity to evaluate model performance and to suggest many practical applications. This paper is part of a larger effort to develop computer models for wood kilns that can be applied to industrial processes. The model has a large potential to become a useful tool for kiln design and a valuable part of a training program for kiln operators.

**********

Wood kilns are designed to remove moisture from the lumber pieces by supplying heat to the lumber to evaporate liquid, which then is carried away by the air circulation. Kiln-drying is influenced by heat and mass transfer between the airflow and wood as well as by the complex moisture transport processes happening inside each piece of lumber. Kiln designers and operators face a difficult task of balancing product quality with a demand for the reduction of drying time. When limited experimental data are available, it is quite natural that mathematical modelling would become an integral part of the research and design processes. For the past two decades an extensive amount of work has been done to improve our understanding of the processes inside a wood kiln and to develop models simulating kiln-drying. Most of the previous research focused on modelling either airflow inside the kiln (Hua et al. 2001) or wood drying processes inside the lumber (Kamke and Vanek 1994). The previous research achieved success in th e development of separate models and prepared the way for the development of a model that combines both airflow and wood drying sub-models into one comprehensive model. The purpose of this paper is to present that combined model, to show its application for a test case and to suggest some of its potential for future simulations.

Physical and mathematical models

Airflow model

A schematic diagram of a typical wood kiln is shown in Figure 1. Airflow in the kiln is created by axial fans located above the wood load. Flow passes around a top load baffle and enters into the entrance plenum, then passes through wood packages, providing the lumber with heat and taking moisture Out, until reaching the middle plenum. From the middle plenum, air goes into the second stack of lumber and out to the outlet plenum; from there it flows back to the fans. Heating coils can be located at the top and in the middle; vents, which allow moist air to leave the kiln and fresh air to enter it, are usually located at the top on the suction and pressure sides. Fans are usually reversible and motors can be two-speed or be equipped with variable speed drives. The average velocity in the channels between wood pieces usually varies from 2 to 6 m/sec., and the distance between wood rows from 16 to 38 mm. For that velocity range, density can be assumed constant and body forces such as gravity can be neglected. As most industrial kilns are very long and fans are spaced at regular intervals along the length, the air profile can stabilize before entering the plenum, and a two-dimensional approximation is justified. For turbulence modelling, the standard k-[epsilon] model with wall function treatment is chosen for its reliability and known robust performance. This airflow model has been developed and extensively used by the University of British Columbia to analyze the influence of the kiln geometry on the airflow pattern. Some of the findings regarding different plenum arrangement, various sticker thickness, and flow between an array of boards were presented previously (Hua at al. 2001).

Wood drying model

During the last 10 years, much effort was directed toward the development of wood drying models. As a result, there has been significant progress in the understanding of water movement during the drying process. An extensive overview of drying models has been made by Bian (2001). For implementation into the comprehensive kiln model, the "three variable" model based on an air and water mass balance, together with the energy balance developed by Perre and Turner (1996), was chosen. The main model variables are moisture content (MC), temperature, and gaseous pressure. The MC is defined as a ratio of the weight of all water in the wood to the ovendry weight of the wood. Three phases (gas, liquid, and solid) are considered, and the gas phase is separated into vapor and air parts. The solid and gas phases are always continuous, while the liquid phase is continuous only until it reaches the fiber saturation point, that is the point where the only liquid left in the wood is that held hygroscopically by the cell walls , and called bound water. Thermodynamic equilibrium is assumed, and the gas phase behaves like an ideal gas. The mass balance equations are written for all three phases.

[partial]([[rho].sub.l] + [[rho].sub.b])/[partial]t + [nabla]([[rho].sub.l]x[V.sub.l])+

[nabla]([[rho].sub.b]x[V.sub.b])=-<M>-

<[M.sub.b]> liquid phase [1]

[partial][[rho].sub.a]/[partial]t + [nabla]([[rho].sub.a]x[V.sub.a])=0

air part of gas phase [2]

[partial][[rho].sub.v]/[partial]t + [nabla]([[rho].sub.v]x[V.sub.v])=<M> +

<[M.sub.b]> vapor part of gas phase [3]

where:

[rho] = density (kg/[m.sup.3])

V = velocity (m/sec.)

<M> = mass flux (kg/[[m.sup.3]*sec.])

l = liquid water

b = bound water

a = air

v = vapor

An overbar implies averaging over the corresponding phase.

The enthalpy conservation is written:

[partial]([rho]H)/[partial]t + [nabla]([rho]xVxH)=

[nabla]x([K.sub.eff][nabla]T) [4]

where:

H = enthalphy (J)

T = temperature (K)

[K.sub.eff] = effective conductivity tensor (W/[m*K])

Coupling between airflow and wood drying models

As was mentioned earlier, past modelling work done for wood kilns focused mainly on either one of the two main aspects of kiln performance: airflow distribution or wood drying. This separate approach neglects the enormous influence that one of the processes has on the other. Non-uniformity of the airflow makes drying conditions different for different places inside the kiln. Due to different wall stresses on the wood air interface, heat and mass transfer varies from one wood piece to another and along the piece itself. Wood drying can influence the airflow by changing its temperature and composition. Moreover, moisture removal from different wood pieces is affected by the location of the piece in the stack and by neighboring pieces: the drying rate of a board upstream would influence the drying of a downstream board, because air conditions will change along the passage. For a more complete understanding of these processes inside the wood kiln during drying, these interactions must be taken into consideration. From a practical point of view, a comprehensive model could be used directly by kiln design engineers, because such a model is able to consider the whole kiln-drying process, providing answers on various "what if" questions. The comprehensive approach makes it easier to analyze real industrial problems, bridging the gap between theory and practice.

There are several options for coupling the two models. One is to have a complete transient airflow calculation including temperature and vapor concentration running simultaneously with wood drying calculations. Although this is not impossible, given today's computer power, it is highly uneconomical due to the computer time and resources required. Airflow parameters such as velocities, pressure, and turbulence change insignificantly during drying at a constant fan speed, but their solution still demands a lot of computational time. The logical alternative would be a second approach in which the airflow is calculated before drying starts. Once done, velocities, pressure, and turbulence parameters are fixed at their steady state values. The only parameters that change will be temperature and vapor density. Obviously, this brings a significant simplification to the problem without seriously compromising the main principles. But simplifications need not stop there. It is known that for the accurate solution of the continuity and momentum equations for the airflow, a fine mesh is required; while, at the same time, the convection-diffusion problem can tolerate a much coarser grid. This discussion leads to a third option in which the temperature and vapor density will be solved on a coarser grid, ensuring that most of the physics are captured for the least amount of computational cost.

For physically correct coupling, it is important to treat accurately the interface between the wood and the airflow. A typical interface is shown in Figure 2. Heat from the air is transferred to the wood; part of it is used for water evaporation from the wood surface and another part is conducted inside. The mass flux of water vapor into the air consists of liquid evaporated from the surface and vapor flux coming from inside the wood. The boundary conditions, which result from continuity and conservation requirements, are shown below.

[[rho].sub.v] [\.sub.w] = [[rho].sub.v] [\.sub.a], T[\.sub.w] = T[\.sub.a], P[\.sub.w] = P[\.sub.a] [5]

[m.sub.v] [\.sub.w] + [m.sub.ev .j] = [m.sub.v] [\.sub.a], q[\.sub.w] + [m.sub.ev.l] x

[DELTA][h.sub.vap] = q[\.sub.a] [6]

where:

[[rho].sub.v] = vapor density (kg/[m.sup.3])

P = pressure (Pa)

[m.sub.v] = vapor mass flux (kg/[[m.sup.2]*sec.])

[m.sub.ev.l] = mass flux of evaporated liquid (kg/[[m.sup.2].sec])

[DELTA][h.sub.vap] = latent heat of evaporation (J/kg)

q = heat flux (W/[m.sup.2])

w = wood side of the interface

a = air side of the interface

To determine the heat and mass transfer on the air side, simple equations are used:

[m.sub.v] [\.sub.a] = [h.sub.m] x ([[rho].sub.v] [\.sub.w]- [[rho].sub.v]), q[\.sub.a] =

h x (T[\.sub.w] - T) [7]

where:

[h.sub.m] = local mass transfer coefficient (m/sec.)

h = local heat transfer coefficient (W/[[m.sup.2]*K])

Values without a subscript are outside flow values.

Because the flow between wood pieces is essentially a turbulent flow over a flat plate, the Chilton-Colburn analogy between mass and momentum transfer can be employed for obtaining the local mass transfer coefficient from the airflow calculations.

[Sh.sub.x]/[Re.sub.x] x Sc - 1/3 = [C.sub.[f.sub.x]]/2 [8]

where:

[Sh.sub.x] = [h.sub.m] x x/D =

local Sherwood number

[Re.sub.x] = V x x/v =

local Reynolds number

Sc = v/D = Schmidt number

[C.sub.[f.sub.x]] = [[tau].sub.w]/1/2[rho]x[V.sup.2] =

The skin friction coefficient

Analogously, for local heat transfer coefficient, the Reynolds analogy can be used:

[St.sub.x] x [Pr.sup.2/3] = [C.sub.f]/2 [9]

where:

[St.sub.x] = h/[rho]x[c.sub.p]xV =

local Stanton number

Pr = [c.sub.p] x v x [rho]/k = Prandtl number

x = the distance from the board edge (m)

D = water vapor diffusivity through air ([m.sup.2]/sec.)

[nu] = kinematic viscosity of air ([m.sup.2]/sec.)

k = air thermal conductivity (J/[m.K])

[c.sub.p] = specific air heat capacity (J[kg.K])

[[tau].sub.w] = wall shear stress (Pa), which is calculated from the turbulence parameters based on wall functions

Numerical solution

Method for solving initial airflow equations

For solving the airflow Reynolds averaged Navier-Stokes equations, a coordinate transformation with domain segmentation is employed. Decomposition of the kiln into several computational blocks allows different types of grids with various densities to be used in different kiln parts. A sectional pressure correction procedure is used to accelerate convergence. A detailed description of the method can be found in previous study (Nowak 1998).

Method for solving wood drying equations

After some rearrangement and application of phenomenological laws, it is possible to use Equations [1] to [4] to obtain a system of partial differential equations for three primary variables: MC, temperature, and gaseous pressure. As a wood board has a rectangular shape, it is logical to use an orthogonal grid for the numerical solution. In the drying process, high gradients of drying variables occur near the exchange surfaces; in order to resolve these properly, a small mesh size near the wood--air boundary is required. An additional reason to have a sufficiently small mesh near boundaries lies in the implementation of boundary conditions. Properties at the interface are assumed to be the same as at the bordering cell centers. To justify this approach, small boundary cells are necessary. At the same time, in the middle of the wood piece, a much bigger grid size can be tolerated. Under the uniform grid approach, an excessively refined mesh has to be used everywhere, which make computation expensive and ineffi cient.

A control volume cell centered technique together with implicit Crank-Nicolson time discretization is used to obtain a system of linear algebraic equations. Central difference approximations are used for values at the cell interfaces. Coefficients at the cell interface are computed as weighted averages from neighboring cell values. After discretization, a block-structured non-symmetric matrix was obtained. This matrix is extremely sparse; it has zero elements everywhere except along five diagonals and on these diagonals it contains small 3 x 3 matrices. For a coupled solution of this system, the Generalized Conjugate Residuals method (GCR) (Wesseling 2001) is used. Although this method was originally developed for single matrices, its generalization for coupled systems of equations with block-structured matrices is straightforward. GCR, as with all of the methods from the conjugate gradient family, requires preconditioning for efficient performance. Here the ILU(0) (Incomplete lower upper factorization) preco nditioning (Ferziger and Peric 1999) was used and proved to be very efficient for this application. Implementation of this method allows more flexibility in the choice of model input parameters, which is necessary, especially when the wood drying model is combined with the airflow.

Method for solving transport equations in the air side

To obtain the solution for moisture and temperature distributions in the airflow, the conservation equations were discretized on the coarse grid with a control volume approach and fully implicit time advancement. First order up-wind approximations for cell interface values were used because they provide a high degree of robustness (Ferziger and Peric 1999). As some of the computational domains on the air side can be two-dimensional, care must be taken when interpolating values of the flow parameters from the fine grid, on which the airflow equations were solved, to the coarse grid, on which the vapor density and temperature equations were solved. Mass conservation has to be preserved under such interpolation, otherwise unphysical solutions may emerge. In this model, the velocity along the x-axis and all turbulence parameters are interpolated by means of a mass flux weighted scheme, and, to obtain velocity along the y-axis, the continuity equation was solved on the coarse grid. Systems resulting from discretiz ation of one-dimensional equations can be simply solved by marching methods, while for linear systems obtained from two-dimensional equations, the Bi-CGSTAB method was used (Van der Vorst 1992).

Method for providing coupling between airflow and wood drying

The coupling procedure is designed to achieve energy and vapor mass conservation inside the kiln, keeping computational cost at an acceptable level. Calculations start from the airflow model and flow parameters are calculated on the fine grid without heat and mass transfer. After convergence is reached, heat and mass transfer coefficients are obtained from Chilton-Colburn (Eq. [8]) and Reynolds (Eq. [9]) analogies and a new coarse grid is generated for the airflow domain. After solution initialization, the heat and mass transfer equations in the air side are solved simultaneously with the wood drying model equations and boundary conditions (Eq. [7]) applied at the wood air interface. Inlet and outlet plenums are solved with a one-dimensional approximation. Two different approaches are analyzed for the middle plenum solution: full mixing approximation or solution of the 2D-conservation equations. There are many control options that can be implemented for the modelling of the top inlet area: some of the moistur e can be removed, heat can be added to the airflow, and different time schedules can be studied. It is easy to see that many computations are independent of each other during a single time step. This leads to the idea of using parallel computing to further reduce computational time. The present computer program is parallelized by using the Message Passing Interface (MPI) that allows a cluster of workstations to be used as a single parallel computing system with distributed memory (Gropp at al. 1994). With this approach, wood rows are distributed evenly between available computers that communicate between each other at the end of every iteration, sending heat and mass fluxes through their boundaries and receiving back updated outside flow conditions.

Results and discussion

Wood drying model verification

For a computational investigation of the wood drying, a rectangular wood board 0.2 by 0.05 m is chosen. All of the wood properties are taken from the work by Perre and Degiovanni (1990). At first, the model behavior with different types of grids is studied. The dry-bulb temperature is kept at 80[degrees]C, wet-bulb at 60[degrees]C, initial temperature is 20[degrees]C and MC is 0.99. Uniform grids 30 x 14, 160 x 40 and 240 x 60 are compared with a non-uniform grid 30 x 14 (Fig. 3). For comparison, the MC values in the middle of the wood piece are plotted against time. The results shown in Figure 4 confirm the effectiveness of the nonuniform grid. It is easy to see that the non-uniform grid allows the use of substantially less grid cells with the same accuracy. To validate the wood drying model, its predictions are compared with the results of the experiments conducted by Forintek Canada Corporation (Oliveira 2000). Two cases are investigated: in the first, the initial MC in the piece of lumber is 0.66, and in the second, 0.78. Dry- and wet-bulb temperatures are set to the same values measured during the experiment. A plot of the average MC change with time is presented in Figures 5 and 6. Relatively good agreement between calculated and measured values is achieved. Although there are some quantitative deviations, the qualitative behavior of the model is in the good agreement with the experiment.

Sample coupling case results and analysis

The kiln configuration for this numerical experiment contains many of the features of an industrial kiln. Its height is 6.65 m, width is 4.7 m, inlet and outlet plenums width is 1.04 m, and middle plenum width is 0.76 m. It contains four wood packages divided into front and back stacks, one package sitting on top of the other. Each package contains 290 lumber pieces of 2 by 6's (38.1 by 139.7mm) assembled into 29 rows with a distance between rows (sticker thickness) of 19 mm. Recent research (Sun 2001) shows that the influence of small gaps between the boards on heat and mass transfer coefficients is insignificant. This supports the assumption that all boards in the same row have exactly the same height and no spaces exist between them. Based on that, each row of boards is modelled as the same piece that can be divided in separate boards during the analysis stage. The velocity at the inlet is assumed to be 2m/sec., which corresponds to an average velocity value of about 4mlsec. between wood pieces. Dry- and w et-bulb temperatures are kept constant and equal to 80[degrees]C and 60[degrees]C, respectively. Initial MC is 1.0. Pressure is assumed to be equal to its atmospheric value, an acceptable assumption because all pressure variations in the airflow side, obtained from the airflow solution, are substantially smaller than the pressure variations inside the wood during drying.

For a full kiln, modelling is simplified by dividing its geometry into computational sub-domains (Fig.7). The kiln is separated into top inlet area, top outlet area, inlet, middle and outlet plenums, and separate air channels between rows of wood pieces. In spite of the simplifications, all main physical features of the flow inside the kiln are retained. The grid for the airflow solution is shown on the left part of Figure 8. It should be adequate to provide all desired details of the flow. For the coupled solution, two choices of grid are investigated: in the first one, shown in the middle of Figure 8, the air side is one dimensional throughout all sub-domains, in the second one (right side of Fig. 8), a two-dimensional grid is generated in the middle plenum. Both cases are calculated to compare their performance. During drying, data for MC, temperature, and pressure inside the wood pieces and temperature and vapor density in the air are collected.

The results showing the average MC distribution around the kiln are compared in Figure 9. Several important observations can be made from this plot. First, it is easy to see a reduced drying rate at the top front package (rows 1 to 29). This can be explained by lower velocity and transfer coefficients in the upper channels, due to the kiln geometry. Second, less MC at the very top and the very bottom row of boards can be noticed. This occurs because the air passing through the channels above the top and below the bottom boards has lumber only on one side, so this portion of air can maintain a higher temperature and lower MC compared with the air coming through other channels. From this figure, it can be clearly seen that proper resolution of the middle plenum is critical for obtaining realistic results for the back packages, especially for their top section, although at the same time it has almost no influence on the front packages. This can be explained with the help of temperature distribution contours arou nd the kiln, for the case where the 2D solution was obtained in the middle plenum. (Fig. 10) Hot air from the big gap between the packages goes up and creates highly non-uniform air inlet conditions for the top back package. From this analysis, it can be concluded that the only case when the middle plenum solution can be neglected is when there is no gap or little gap between top and bottom packages.

The average MC distribution from wood piece to wood piece along the row is presented in Figure 11. Lumber pieces are numbered from left to right, that is in the same direction as the airflow In every row of lumber, lower MC occurs in the boards that are closer to the airflow entrance. This is a consequence of higher transfer coefficients and higher difference between wood and air properties at the beginning of the air channel. The non-uniformity of kiln-drying inside the row of lumber is shown in Figure 12. The moisture distribution inside some of the wood rows is shown as a surface in a space where x- and y-axes represent the width of stack and the board thickness of lumber and the z-axis corresponds to the MC. Not only the non-uniformity between the rows and between different boards can be analyzed with this model, as has already been discussed, but also non-uniformity within the board itself can be analyzed. In every wood board, a moisture gradient from the center to the surface can be observed. This gradi ent is responsible for stresses, occurring in the wood during drying, which can cause wood defects such as checking and warping. Modelling can provide critical kiln locations as well as critical times where and when high moisture gradients are taking place. After determination of the problem, parametric studies for different time schedules and wood arrangements can be conducted and practical solutions can be obtained.

Obviously, all of these comments apply to the specific kiln configuration that is chosen for the numerical investigation and are intended not as general design recommendations but as an example of modelling capability.

Conclusions

A comprehensive wood kiln model has been presented and numerical methods showing an efficient strategy for the coupled kiln solution have been described. Experimental data for wood drying were compared with the model. A test kiln-drying case was investigated on two different grids and results were analyzed. The capability of the model to simulate heat and mass transfer processes inside the real industrial kiln was confirmed. The kiln model can be an important tool for the design of new kilns, for the analysis of problems encountered in existing kilns, and optimization of their operation. Computational methods can considerably reduce costly experimental time and increase the understanding of drying processes by allowing researchers to look inside the kiln and "inside" the wood piece when a kiln is in operation. Kiln performance can be analyzed by changing initial conditions and drying schedules, reversing the airflow, implementing different control sequences, varying positions of control sensors, and by trying many other alternatives. Model results can be incorporated into kiln control systems, providing relationships between values, measured at some kiln locations, and important kiln performance parameters. This model constitutes a substantial advance towards the development of virtual kiln simulators, which can be used for kiln design and for training kiln operators.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

Literature cited

Bian, Z. 2001. Airflow and wood models for woodkilns. M.A.Sc. thesis. Univ. of British Columbia, Vancouver, BC, Canada.

Gropp, W., E. Lusk, and A. Skjellum. 1994. Using MPI: Portable Parallel Programming with the Message-Passing Interface. MIT Press, Cambridge, MA.

Hua, L., E. Bibeau, P. He, I. Gartshore, M. Salcudean, Z. Bian, and S. Chow. 2001. Modelling of airflow in wood kilos. Forest Prod. J. 51(6):74-81.

Ferziger, J. and M. Peric. 1999. Computational Methods for Fluid Dynamics. Springer-Verlag, NY. 389 pp.

Kamke, F.A. and M. Vanek. 1994. Comparison of wood drying models. In: 4th IUFRO Inter. Conf. on Wood Drying. IUFRO, Vienna, Austria. pp. 1-21

Nowak, P. 1998. GEO-MGFD: A system of computer programs for steady-state flow of variable density fluid in the complex three dimensional regions using multi-grid and multi-segment techniques. Technical rept. Dept. of Mechanical Engineering, Univ. of British Columbia, Vancouver, BC, Canada.

Oliveira, L. 2000. Wood Drying Scientist, Forintek Canada Corporation. Personal communication.

Perre, P. and A. Degiovanni, 1990. Simulation par volumes finis des Transferts Couples en Milieu Poreux Anisotropes: sechange du bois a basse et a haute temperture. (Control volume formulation of simultaneous transfers in anisotropic porous media: simulations of softwood drying at low and high temperature.) Int. J. Heat Mass Transfer 33(1l):2463-2478. (In French, English abstract).

_____ and I. Turner. 1996. The use of macroscopic equations to simulate heat and mass transfer in porous media: some possibilities illustrated by a wide range of configurations that emphasise the role of internal pressure. In: Mathematical Modelling and Numerical Techniques in Drying Technology. I. Turner and A.S. Mujumdar, eds. Marcel Dekker, NY.

Sun, Z.F. 2001. Numerical simulation of flow in an array of in-line blunt boards: mass transfer and flow patterns. Chemical Engineering Sci. 56(2001): 1883-1896.

Van der Vorst, HA. 1992, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 13(2):631-644.

Wesseling, P. 2001. Principles of Computational Fluid Dynamics. Springer-Verlag, Berlin. 654 pp.

The authors are, respectively, Research Engineer, Graduate Student, Professor, and Professor, Dept. of Mechanical Engineering, Univ. of British Columbia, 2324 Main Mall, Vancouver, BC, Canada V6T 1Z4; and Wood Drying Scientist, Forintek Canada Corp., 2665 East Mall, Vancouver, BC, Canada V6T lW5. The authors would like to acknowledge financial help from the Natural Sciences and Engineering Research Council of Canada and from Forintek Canada Corp. This paper was received for publication in January 2002. Article No. 9432. [C]Forest Products Society 2003. Forest Prod. J. 53(4):46-54.

Printer friendly Cite/link Email Feedback | |

Author: | Pougatch, Konstantin; Bian, Zhengbing; Gartshore, Ian; Salcudean, Martha; Oliveira, Luiz |
---|---|

Publication: | Forest Products Journal |

Geographic Code: | 1USA |

Date: | May 1, 2003 |

Words: | 4731 |

Previous Article: | A rapid method for fluoride analysis of treated wood. (Technical Note). |

Next Article: | Effect of drilled holes on the bending strength of large dimension Douglas-fir lumber. |

Topics: |