Sage Journals: Discover world-class research

Abstract

Drying is a required operation in timber manufacturing industries. Most drying processes utilise convective kilns, which involve coupling between transient heat and mass (moisture content) transfer in solids and convection in a flowing medium. In the present work, the flow field in the kiln is determined using a system approach based on a general head loss equation, whereas the coupling between transient heat and mass transfer in the timber is accomplished by a solution of a local differential model across the timber stack. The scheme is applied to study the effects of distinct geometric kiln configurations and air flow temperatures upon drying rates and moisture content.

1. Introduction

The technological advancements in the last years have made available to industries several conventional and nonconventional wood drying methods such as radiofrequency heating, press, solar, and dehumidification drying strategies [1]. Notwithstanding, convective drying is the mostly used method in agricultural, food, paper, and timber processing industries [2]. In general applications, convective drying of moist objects (e.g., lumber manufacturing, food processing, etc.) conjugates heat and mass transfer both inside the solid and in the flowing medium. The understanding of the related physics of the problem is crucial for the process control and ensuing product quality.

The industrial importance of drying processes has been translated by the increasing number of works reporting numerical prediction of temperature distribution and moisture content in convective drying problems. In convective timber drying, temperature, moisture content, air velocity, and drying rates are the most important parameters to affect product quality; that is, an incorrect setting of global parameters may cause defects such as colour heterogeneities, warp, and splits amongst others. Therefore, development of simple and effective predictive models to evaluate the process parameters is of paramount importance in lumber manufacturing industries. Riepen and Paarhuis [3] present a numerical study on the optimization of flow velocity in convection kilns. The authors conclude that there is a strong relationship between drying rate and air velocity distribution inside the kiln. Kaya et al. [4 –6] devise 2D numerical models to predict heat and mass transfer during convective drying. The authors solve the flow problem to obtain the convection heat transfer coefficient in the boundary of the moist object considering a laminar flow regime. Analogy between thermal and concentration boundary layers is adopted to estimate the convection mass transfer coefficient on the surface, followed by solution of the transient diffusion problem (heat and moisture inside the object). In a similar direction, Mohan and Talukdar [2] present a 3D computational model to study the convective drying of a moist object. More recently, two works assessing convective drying of wood have been reported in the literature [7, 8]. Younsi et al. [7] utilise a 3D model to compute the turbulent flow inside convective kilns and evolution of the temperature and moisture content. Thibeault et al. [8] combine the Ansys CFX (heat and moisture transfer) and FESh++ (mechanical approximation) commercial codes to predict the thermo-hygro-mechanical behaviour of wood during convective drying.

Most literature works involving convective drying focus only on the local physics of some particular moist object (small scale). To the authors' best knowledge, there are no works in the literature based upon a system approach where the convective kiln is modelled in the real (full) scale and air flow is computed using a lumped approximation. It is expected that the kiln geometric configuration, as well as lumber stacking, and sticker number, geometry, and placement affect both the head loss and flow aerodynamics and, as a consequence, the flow operating condition (and, ultimately, the drying rate). Based upon the authors' experience in modelling and simulation of forced convection heat transfer in industrial problems (e.g., [10 –12]), a full scale, globally iterative coupled simulation (in conjunction with the fan performance curve to obtain the operation point) of air flow and heat and mass transfer is computationally intensive and, within the industrial environment, would require CFD advanced users to select proper model settings. Therefore, the present system approach aims at developing a simpler model and is yet able to obtain good approximations to the physical phenomena. Following this route, it is proposed a mixed numerical model (a global approach to determine flow velocities and a local differential model to compute evolution of the timber temperature and moisture content) to predict the effects of distinct kiln configurations and air flow temperature on the evolution of drying rates, moisture content, and temperatures in a full geometric scale.

2. System Approach: The Flow Operating Condition

The strategy used in this work to estimate air velocity across the timber stack is based upon determination of the flow operating point of the global system—the system approach. This technique requires (i) knowledge of the fan performance curves—provided by the manufacturer and (ii) computation of the system characteristic curve. The typical system characteristic curve correlates the total head loss (or pressure drop) and flow rate in the inlet/exit plenums and timber stack, accounting for friction losses and effects of flow separation. This section describes the procedure to obtain the characteristic curve of the system and the operating condition for the compact dry kiln depicted in Figures 1 and 2.

Figure 1:

Compact dry kiln geometry: lumber stack and sticker placement.

Figure 2:

Compact dry kiln geometry: lateral view.

Air flow between two consecutive timber layers is dictated by the pressure difference between the inlet and exit plenums. Ideally, the kiln design should provide uniform pressure in both plenums, making it possible to obtain the same velocity in all flow channels. However, the plenum configuration or inadequate sticker dimensions/placement may cause pressure differences and, consequently, variations in air velocity across the stack height, which may lead to undesirable, nonuniform drying rates. The present formulation follows the concept of Nijdam and Keey [13], who developed a differential head loss equation to compute channel velocities across timber stacks assuming that the kiln configuration is such to prevent boundary layer separation at the top corner. Aiming at a discrete solution for a known number of timber layers and stickers, the proposed approach determines the flow rate based upon an equivalent flow circuit and accounts for channel inlet and exit effects and flow through parallel channels between stickers.

The geometrical information is indicated in Figures 1 and 2. The main parameters are the stack height, H, width, W, and length, L; number of stickers and corresponding dimensions, N_P and L_P × H_S × W; number of lumber layers across the stack height, N_L; and maximum and minimum inlet and exit plenum widths, W_H and W₀. Such configuration causes air to flow between timber layers through parallel, N_S = N_P – 1 individual channels of cross-section A_S = H_S × L_S. The stack level is indicated by k = 1, 2, …, N_C (upward direction), in which N_C = N_L + 1. The geometry of the inlet and exit plenums, depicted in Figure 2, defines a continuous cross-section variation, A_Pⁱ(z) = A_P^o(z) = t(z) × L, so that t(z = 0) = W₀ and t(z = L) = W_H.

The pressure drop, Δp_c, between entrance and exit sections of a generic channel of length l_c can be readily derived from the mass conservation principle and the classical head loss equation [14], so that

Δ p = p_{in} - p_{out} = [\frac{ρ_{a} α_{t, out}}{2 A_{out}} - \frac{ρ_{a} α_{t, in}}{2 A_{in}} + \frac{ρ_{a}}{2 \bar{A}} (λ \frac{l_{c}}{D_{h}} + K_{in} + K_{out})] \times {\dot{Q}}_{}^{2} = R_{} {\dot{Q}}_{}^{2},

(1)

where subscripts “in” and “out” indicate entrance and exit sections, K is a loss coefficient (when separation takes place in contraction/expansion flows), α_t is the kinetic energy correction factor, ρ_a is the specific mass, A is the cross-sectional area, $\dot{Q} = {\bar{V}}_{1} A_{1} = {\bar{V}}_{2} A_{2} = \bar{V} \bar{A}$ is the volumetric flow rate, R is the total hydraulic resistance, and λ is the friction factor, approximated using Swamee and Jain's equation [15]:

\begin{matrix} λ = 0.25 {[\log (\frac{e / D_{h}}{3.7} + \frac{5.74}{R e_{D_{h}}^{0.9}})]}^{- 2}, \\ where R e_{D_{h}} = \frac{ρ_{a} \bar{V} D_{h}}{μ_{a}}, \end{matrix}

(2)

in which e is the channel roughness, μ_a is the air viscosity, and D_h is the hydraulic diameter.

It is relevant to notice that (1) is nonlinear, which prevents use of a direct analogy with Ohm's liner relation between voltage (pressure drop) and electrical current (flow rate). However, a similar approach can be used based upon the mass conservation and by assuming equality of pressure at the entrance/exit of the channels at a given stack level.

2.1. Equivalent Global Flow Circuit

Figures 1 and 2 show that air flows through the inlet plenum, splits into N_S parallel channels at every discrete level k of the timber stack, and exits through the exit plenum. The configuration of the equivalent flow circuit is shown in Figure 3, which highlights the flow through parallel channels for stack levels k and k + 1 (Figure 3(a)) and flow through the timber stack and inlet/exit plenum stretches (Figure 3(b)). In view of the present kiln geometry and flow circuit presented in Figure 3, application of (1) requires computation of the hydraulic resistances for the inlet/exit plenum stretches, R_kⁱ and R_k^e (between stack levels k and k + 1), and parallel channels, R_k^j (between stickers for stack level k), as

\begin{matrix} R_{k}^{i} = ρ_{a} \frac{α_{k}^{i}}{2 A_{P}^{k}} - ρ_{a} \frac{α_{t, k + 1}^{i}}{2 A_{P}^{k + 1}} + ρ_{a} λ_{P, k}^{i} \frac{Δ H_{P}}{{\bar{D}}_{h, k}^{P}} \frac{1}{2 {\bar{A}}_{P, k}}, \\ R_{k}^{e} = ρ_{a} \frac{α_{k + 1}^{e}}{2 A_{P}^{k + 1}} - ρ_{a} \frac{α_{t, k}^{e}}{2 A_{P}^{k}} + ρ_{a} λ_{P, k}^{e} \frac{Δ H_{P}}{{\bar{D}}_{h, k}^{P}} \frac{1}{2 {\bar{A}}_{P, k}}, \\ R_{k}^{j} = ρ_{a} \frac{1}{2 {\bar{A}}_{S}} (K_{S}^{1, j} + K_{S}^{2, j}) + ρ_{a} λ_{S, k}^{j} \frac{W}{D_{h}^{S}} \frac{1}{2 {\bar{A}}_{S}}, \end{matrix}

(3)

in which the subscripts P and S represent plenum and individual channel measures and k is the stack level across the timber stack height; the superscripts i and e indicate inlet and exit plenum stretches and j denotes an individual flow channel at a given stack level k. The hydraulic diameter for a flow channel, D_h^S, and the corresponding mean values for the inlet/exit plenum stretches, ${\bar{D}}_{h}^{P}$ , are

D_{h}^{S} = \frac{2 A_{S}}{H_{S} + L_{S}}, {\bar{D}}_{h}^{P} (z) = \frac{2 {\bar{A}}_{P} (z)}{t (z) + L},

(4)

where H_S, L_S, L, and t(z) are indicated in Figures 1 and 2.

Figure 3:

Equivalent flow circuit: (a) channels (between two timber layers) and (b) inlet/exit plenum stretches and channels (across the stack height).

2.2. Individual Channel Approximation

Figure 1 shows that insertion of stickers between two timber layers creates parallel channels with equal inlet and exit pressures, p_kⁱ and p_k^e, respectively, at every stack level k. Furthermore, the hydraulic resistance for each individual channel at a given level k is constant, R_k^j = R_k^S, (i.e., equal flow rate, ${\dot{Q}}_{k}^{j} = {\dot{Q}}_{k}^{S}$ , at each parallel channel). Therefore, solving (1) for the total flow rate, ${\dot{Q}}_{k}^{C}$ , in the stack level k,

\begin{matrix} {\dot{Q}}_{k}^{C} = {(\frac{Δ p_{k}}{R_{k}^{C}})}^{1 / 2}, \\ where {\dot{Q}}_{k}^{C} = {\dot{Q}}_{k}^{1} + \cdot + {\dot{Q}}_{k}^{N_{S}} = \sum_{j = 1}^{N_{S}} {\dot{Q}}_{k}^{j}, \\ {\dot{Q}}_{k}^{j} = {(\frac{Δ p_{k}}{R_{k}^{j}})}^{1 / 2}, \end{matrix}

(5)

one obtains ${\dot{Q}}_{k}^{C}$ and the hydraulic resistance, R_k^C, for level k as

\begin{matrix} {\dot{Q}}_{k}^{C} = {(\frac{Δ p_{k}}{R_{k}^{C}})}^{1 / 2} = \sum_{j = 1}^{N_{S}} {\dot{Q}}_{k}^{j} = \sum_{j = 1}^{N_{S}} {(\frac{Δ p_{k}}{R_{k}^{j}})}^{1 / 2} = N_{S} {(\frac{Δ p_{k}}{R_{k}^{S}})}^{1 / 2}, \\ R_{k}^{C} = \frac{R_{k}^{S}}{{N_{S}}^{2}} . \end{matrix}

(6)

2.3. Inlet and Exit Plenum Approximation

The inlet and exit plenums are split into stretches of equal height ΔH_P, as indicated in Figure 2, for which hydraulic resistances R_kⁱ and R_k^e are computed using (3). Therefore, in addition to R_kⁱ and R_k^e, the global circuit configuration includes also the hydraulic resistance for stack level k, R_k^C, as illustrated in Figure 3(b). The pressure drop in the channels and plenum stretches for two consecutive stack levels, k and k + 1, are, respectively,

\begin{matrix} Δ p_{k} = p_{k}^{i} - p_{k}^{e} = R_{k}^{C} {({\dot{Q}}_{k}^{C})}^{2}, \\ Δ p_{k + 1} = p_{k + 1}^{i} - p_{k + 1}^{e} = R_{k + 1}^{C} {({\dot{Q}}_{k}^{C})}^{2}, with \\ p_{k + 1}^{i} - p_{k}^{i} = R_{k}^{i} {({\dot{Q}}_{k}^{i})}^{2}, \\ p_{k}^{e} - p_{k + 1}^{e} = R_{k}^{e} ({\dot{Q}}_{k}^{e}), \end{matrix}

(7)

in which ${\dot{Q}}_{k}^{i} = {\dot{Q}}_{k}^{e} = {\dot{Q}}_{k}^{P}$ is the flow rate in the plenum stretch. By combining (7), it is possible to determine a recurrence equation for ${\dot{Q}}_{k + 1}^{C}$ and ${\dot{Q}}_{k + 1}^{P}$ as a function of the flow rate for the parallel flow channels and plenum stretches at stack level k as

\begin{matrix} {\dot{Q}}_{k + 1}^{C} = \sqrt{\frac{(R_{k}^{i} + R_{k}^{o})}{R_{k + 1}^{C}} {({\dot{Q}}_{k}^{P})}^{2} + \frac{R_{k}^{C}}{R_{k + 1}^{C}} {({\dot{Q}}_{k}^{C})}^{2}}, \\ {\dot{Q}}_{k + 1}^{P} = {\dot{Q}}_{k}^{P} + {\dot{Q}}_{k + 1}^{C} . \end{matrix}

(8)

The total flow rate, ${\dot{Q}}_{N_{C}}^{P}$ , and pressure drop, Δp_{N
_C} = p_{N
_C}ⁱ – p_{N
_C}^e, in the system are determined by acknowledging that the flow rate in level k = 1 (under the first timber layer) is the same as in the first inlet/exit plenum stretches, ${\dot{Q}}_{1}^{P} = {\dot{Q}}_{1}^{C}$ . The subsequent values across the stack height $(Δ p_{2}, {\dot{Q}}_{2}^{P})$ , $(Δ p_{3}, {\dot{Q}}_{3}^{P}), \dots$ , $(Δ p_{N_{C}}, {\dot{Q}}_{N_{C}}^{P})$ are determined recursively using (7) and (8). Therefore, the system characteristic curve can be readily determined by assuming different values of ${\dot{Q}}_{1}^{C}$ within a range $[0, {\dot{Q}}_{1, \max}^{C}]$ , from which the corresponding total flow rate, ${\dot{Q}}_{N_{C}}^{P}$ , and pressure drop, Δp_{N
_C}, are determined at the end of the recursive computation.

2.4. Fan Performance Curve and Operating Point

Conceptually, the operating point consists in the balancing flow condition between the mechanical energy provided by fans and the losses imposed by the flow topology in the dry kiln, represented by the point at which the fan performance curve at a given rotational speed intersects the system characteristic curve. Mathematically, the flow operating condition $(Δ p_{D}, {\dot{Q}}_{D})$ can be determined by the solution of the nonlinear expression $G (\dot{Q}, Δ p) - F (\dot{Q}, Δ p) = 0$ , where $G (\dot{Q}, Δ p)$ and $F (\dot{Q}, Δ p)$ are the system and fan characteristic equations obtained by application of nonlinear (polynomial) regression to the system and fan characteristic curves as

\begin{matrix} G (\dot{Q}, Δ p) = Δ p - \sum_{r = 0}^{n_{r}} γ_{r} {\dot{Q}}_{}^{(r)}, \\ F (\dot{Q}, Δ p) = Δ p - \sum_{p = 0}^{n_{p}} β_{p}^{} {\dot{Q}}_{}^{(p)}, \end{matrix}

(9)

in which (n_r,γ_r) and (n_p, $β_{p}$ ) are the maximum order of the regression polynomials and corresponding coefficients for the system and fan, respectively. A second-order polynomial, n_r = 2, for the system would suffice for a coefficient of determination R² ≅ 1.0.

The fan characteristic curves are provided by the manufacturer and should indicate the respective rotational speeds; otherwise the so-called fan laws [16] can be used to estimate the performance curves within a limited range of speeds. It is worthy to mention that the fan performance curve may present a nonmonotonic behaviour with respect to the flow rate, which might require a piecewise approximation to improve accuracy. Furthermore, additional care should also be exercised to avoid operation near the fan instability region, that is, where the system characteristic curve intersects the fan curve in more than one point or both curves present similar slopes.

At the operating point, the flow rate, ${\dot{Q}}_{D} = {\dot{Q}}_{N_{C}}^{P}$ , is computed at the entrance of the inlet plenum using the recurrence equation (8). Therefore, in addition to ${\dot{Q}}_{D}$ , the corresponding flow rates, ${\dot{Q}}_{D, k}^{S}$ , in each individual channel at stack level k are also available, making it possible to estimate the individual channel velocities ${\bar{V}}_{D, k}^{S} = {\dot{Q}}_{D, k}^{S} / A_{S}$ .

3. Heat and Mass Transfer

The transient heat and mass diffusion processes in a timber layer are modelled by a one-dimensional formulation. The global and local approaches are illustrated in Figures 1, 2 and 4. Figures 1 and 2, show the timber stack and lateral view of the kiln in a full scale, whereas the local configuration is presented in Figure 4, which portraits a detailed view of the local geometry with its main dimensions and variables.

Figure 4:

Local geometry: flow channels, where H_S, W, and L_S are the channel height, length, and depth (the latter is not represented in the picture: see Figure 1), respectively, and timber layers with thickness H_T.

Assuming that W ≫ H_T (timber layers are thin), the mathematical model comprising the energy and chemical element (moisture) conservation principles reduces, respectively, to

\frac{\partial T}{\partial t} = α_{T} \frac{\partial^{2} T}{\partial z^{2}}, \frac{\partial M}{\partial t} = \frac{\partial}{\partial z} (D_{M} \frac{\partial M}{\partial z}),

(10)

where z is the direction normal to the gas flow, T is the timber temperature, α_T = μ_Tc_T/κ_T is the thermal diffusivity of the wood, M represents the moisture content in the solid (kg of moisture/kg of dry wood), and D_M is the diffusion coefficient of moisture in the timber. The diffusion coefficient, D_M, is temperature dependent, being given by the Arrhenius law

D_{M} = D_{Γ} \exp (- \frac{D_{∊}}{T}),

(11)

in which D_∊ and D_Γ are experimental constants. The present work assumes D_Γ = 4.66 × 10⁻⁵ m²/s and D_∊ = 3.771 × 10³ K, corresponding to the drying wood process in a wide range of temperatures [17].

Computation of the global drying curve requires solution of (10)-(11) for each individual timber layer owing to possible differences in the flow velocities across the timber stack. In addition, the boundary conditions at the solid/gas interface and the initial condition for transient heat and mass transfer problems must also be specified. The energy balance at both interfaces of a timber layer of thickness H_T leads to

\begin{matrix} {{- κ}_{T} \frac{\partial T}{\partial z} |}_{z = 0} = h_{c} [T (z = 0) - T_{a}], \\ {{- κ}_{T} \frac{\partial T}{\partial z} |}_{z = H_{T}} = h_{c} [T (z = H_{T}) - T_{a}], \end{matrix}

(12)

where κ_T is the thermal conductivity of the solid (timber layer), T_a is the air flow temperature, T(z = 0) and T(z = H_T) are temperatures at the solid/gas interfaces, and h_c is the convective heat transfer coefficient. The latter is determined by assuming fully developed flow of the gas in the channels so that [18]

\frac{h_{c} D_{h}}{κ_{a}} = 8.23, \frac{h_{c} D_{h}}{κ_{a}} = 0.023 R {e_{D_{h}}}^{0.8} P r^{0.4},

(13)

for laminar and turbulent flow regimes, respectively. In (13), D_h is the hydraulic diameter, k_a is the thermal conductivity of the gas, and $R e_{D_{h}} = ρ_{a} {\bar{V}}_{D, k}^{S} D_{h} / μ_{a}$ and Pr = μ_ac_a/κ_a are the Reynolds and Prandtl numbers, respectively. Therefore, the individual mean channel velocity at the operating point, ${\bar{V}}_{D, k}^{S}$ , (detailed in the preceding section) and definition of the air temperature make it possible to obtain the heat transfer coefficient, h_T (computed individually for each flow channel, k).

Application of the mass balance for the chemical element (moisture) at the air/timber interfaces leads to

\begin{matrix} - D_{M} {\frac{\partial M}{\partial z} |}_{z = 0} = h_{M} [M (z = 0) - M_{a}], \\ - D_{M} {\frac{\partial M}{\partial z} |}_{z = H_{T}} = h_{M} [M (z = H_{T}) - M_{a}], \end{matrix}

(14)

where M_a is the moisture content of the drying air, M(z = 0) and M(z = H_T) are the moisture content in the timber/air interfaces, and h_M is the convective mass transfer coefficient. The coefficient h_M is determined by using the analogy between thermal and concentration boundary layers, δ_t/δ_c ≈ Le^1/3 [18], that is,

h_{M} = h_{c} (\frac{D_{AB} L e^{1 / 3}}{k_{a}}),

(15)

where Le = α_a/D_AB is the Lewis number and D_AB is the diffusion coefficient of the chemical component A (moisture) in B (air). It is important to state that the procedure adopted in the present work for modelling conjugated heat and mass diffusion has also been adopted by other authors [2, 4 –6] within similar framework.

The heat and mass conservation equations (10) are fully coupled through the mass diffusion coefficient, D_M, and mass and heat convective coefficients, h_M and h_c, thereby requiring a numerical solution. The finite volume method was used in this work owing to its local conservative character and discretisation/implementation simplicity. The space derivatives are discretized by central second-order formulae whereas the time derivative is approached by the second-order accurate Cranck-Nicholson scheme. The resulting systems of algebraic equations are solved by the well-known TDMA algorithm. The procedure is classical in the literature and for more details the reader is referred to [19]. The numerical approximations for the local problems were validated against analytical solutions for a one-dimensional plate subjected to convective transfer on both sides [9]. The largest local errors for the nondimensional parameters were smaller than 1.8% (very early in the process), substantially decreasing as the process advances. Further discussions and validation of the local approximation are presented in [12].

4. Numerical Examples

4.1. Validation of the Flow Rate Model

The present system approach combines a global solution to obtain flow rates and velocities inside the channels and a local approximation to heat and mass transfer for timber layers. The critical aspect of the present model is the accurate determination of the channel velocities across the timber stack. Validation is, therefore, performed for the flow solution against a CFD simulation of a reduced size dry kiln using the Ansys CFX commercial package.

The simulations are performed for two plenum sizes, W_H = W₀ = 0.090 and 0.180 m, and 12 timber layers. The entrance velocity at the inlet plenum is assumed V_{N
_C}^P = 4.2 m/s for both cases. Table 1 indicates the geometrical data and flow information. It is important to note that, in this example, only the hydrodynamic solution is determined, especially the channel velocities across the stack height.

Table 1:

Dimensions and operational parameters for a reduced size dry kiln.

Parameter/measure	Symbol	Value

Stack
Height	H	0.185 m
Length	L	1.000 m
Width	W	0.300 m
Sticker
Number	N _P	2
Width	L _P	0.1 × 10⁻³ m
Height	H _S	5 × 10⁻³ m
Timber
Thickness	H _T	10 × 10⁻³ m
Head loss
Loss coefficient—inlet	$K_{S}^{1}$	0.5
Loss coefficient—exit	$K_{S}^{2}$	1.0
Inlet/exit kinetic factors	$α_{t}^{i} = α_{t}^{e}$	1.0
Plenum roughness	e _P	2.0 × 10⁻⁶ m
Timber roughness	e _S	3.0 × 10⁻³ m
Air
Temperature	T _a	373 K
Specific mass +	$ρ_{a} = \frac{p_{a}}{286.9 T_{a}}$ [kg/m³]
Viscosity*	$μ_{a} = \frac{1.45 \times 1 0^{- 6} {T_{a}}^{3 / 2}}{110 + T_{a}}$ [Pa s]

The assumption of ideal gas is used in the hydrodynamic solution [9], in which p_a and T_a are the kiln average pressure and temperature, respectively.

Viscosity is computed using Sutherland's law [9].

The velocity distribution for plenum sizes (a) W_H = W₀ = 0.090 m and (b) W_H = W₀ = 0.180 m is presented in Figure 5. The same entrance velocity at the inlet plenum (located at z = H in Figure 2) used in both cases gives rise to substantially different flow rates across the timber stack. In case (a), the total flow rate is ${\dot{Q}}_{(a)}^{P} = 0.378$ m³/s, whereas the larger plenum of case (b) leads to a higher flow rate, ${\dot{Q}}_{(b)}^{P} = 0.756$ m³/s. Since the stack height, lumber thickness, and channel height are the same for both settings, case (b) provides channel velocities twice as high as case (a); that is, ${\bar{V}}_{(b)} = 11.99$ m/s and ${\bar{V}}_{(a)} = 5.94$ m/s, respectively. Figure 5 shows that differences between the Ansys CFX and present approach are acceptably small, with average values 3.62% and 1.64% for cases (a) and (b), respectively.

Figure 5:

Channel velocities across the stack height for plenum velocity V_{N
_C}^P = 4.2 m/s and plenum sizes (a) W_H = W₀ = 0.090 m and (b) W_H = W₀ = 0.180 m.

4.2. Full Scale Simulation of a Dry Kiln

It has been well established in the literature that the air temperature and velocity dictate the drying rate. Therefore, the ideal compact dry kiln would have to combine not only uniform drying rates, but also high lumber volume and compact dimensions (e.g., by designing narrow inlet and exit plenums). Within this framework, this example focuses on two issues: (a) the aerodynamic behaviour of the dry kiln and (b) evolution of the lumber moisture content and drying rates.

The timber stack comprises 30 lumber layers (L × W = 2.40 × 1.60 m) of thickness H_T = 25.4 × 10⁻³ m, as depicted in Figures 1 and 2. The other dry kiln dimensions and process parameters are presented in Table 2. The simulations were performed for plenum sizes W_H = W₀ = 0.10, 0.15, 0.20, 0.30, 0.50, and 0.80 m. The fan performance curve is described using a fourth-order polynomial in $F (\dot{Q}, Δ p)$ of equation (9) with a coefficient of determination R² = 0.99987, as indicated in Table 3. In the present case, the operating region is located after the instability zone, thereby ensuring smooth operating conditions.

Table 2:

Dry kiln dimensions and operational parameters.

Parameter/measure	Symbol	Value

Stack
Height	H	2.19 m
Length	L	2.40 m
Width	W	1.60 m
Sticker
Number	N _P	3
Width	L _P	46 × 10⁻³ m
Height	H _S	46 × 10⁻³ m
Timber
Initial temperature	T _T ⁰	25°C (298 K)
Initial moisture content	M ₀	12 kg/kg
Thickness	H _T	25.4 × 10⁻³ m
Thermal conductivity	κ_T	0.18 Wm/s
Specific heat	c _T	1883 J kg/K
Specific mass	ρ_T	825 kg/m³
Head loss
Loss coefficient—inlet	K _S ¹	0.5
Loss coefficient—exit	K _S ²	1.0
Inlet/exit kinetic factors	$α_{t}^{i} = α_{t}^{e}$	1.0
Plenum roughness	e _P	2.0 × 10⁻⁶ m
Timber roughness	e _S	3.0 × 10⁻³ m
Air
Air moisture	M _a	Dry air
Moisture/air mass diffusivity	D _AB	2.6 × 10⁻⁵ m²/s
Thermal conductivity	κ_a	2.889 × 10⁻² W/mK
Specific heat	c _a	1013.2 J kg/K
Specific mass	$ρ_{a} = \frac{p_{a}}{286.9 T_{a}}$ [kg/m³]
Viscosity	$μ_{a} = \frac{1.45 \times 1 0^{- 6} {T_{a}}^{3 / 2}}{110 + T_{a}}$ [Pa s]

Table 3:

Fan performance curve for a rotational speed ω = 1740 rpm: fourth order polynomial regression, n_p = 4, in $F (\dot{Q}, Δ p)$ of equation (9).

Flow rate, $\dot{Q}$ [m³/s]
Coefficient β_p	Value

β₀	74.46793666109
β₁	−9.741128005452
β₂	−2.599301152791
β₃	0.686533446875
β₄	−0.041653532981

4.2.1. Hydrodynamic Solution: Velocity Distribution across the Timber Stack

This section highlights the effect of the plenum size on the hydrodynamic performance of the compact dry kiln using the formulation discussed in the previous sections. The system approach provides the operating point corresponding to the intersection between the fan performance curve and the system characteristic curve. Figure 6 presents the operating conditions for different plenum sizes, in which $\dot{Q}$ is the global flow rate and Δp is the pressure difference between the entrance of the inlet plenum and exit of the exit plenum. Narrow plenums (e.g., in this example, W_H = W₀ = 0.10 m) impose significant head losses leading to smaller global flow rates and higher pressure drops. As the plenum size increases, the wall effect decreases, permitting balance flow conditions at smaller pressure drops and higher flow rates, converging to ${\dot{Q} |}_{\infty} \approx 9.02$ m³/s and ${Δ p |}_{\infty} \approx 1.04$ mm H₂O for the present compact kiln. Theoretically, when W_H → ∞ and W₀ → ∞, the pressure in the inlet and exit plenums is uniform and Δp is due only to the pressure drop in the flow channels (flow between two consecutive lumber layers) across the timber stack.

Figure 6:

System approach: fan and system characteristic curves and operating points for plenum sizes W_H = W₀ = 0.10, 0.15, 0.20, 0.30, 0.50, and 0.80 m.

The operating point represents a global measure of the hydrodynamic balance in the dry kiln. For each such point depicted in Figure 6, the global flow rate, $\dot{Q}$ , encompasses individual contributions by each flow channel, as illustrated in Figures 3(a) and 3(b). As discussed in the previous paragraph, the dimensions of the inlet and exit plenums affect head losses, which in turn cause different flow rates and, consequently, air velocities, V(z), across the timber stack. Therefore, narrow plenums entail higher head losses, causing air to flow primarily across the top channels and producing highly inhomogeneous flow rates and velocity distributions across the timber stack. Conversely, the virtually uniform pressure distributions in large plenums not only lead to higher global flow rates, but also make it possible uniform velocity distributions across the timber stack (highly desirable to ensure uniform drying rates). Figure 7 shows the velocity distribution in all 31 channels across the lumber stack, in which the differences caused by narrow plenums can be clearly observed. The simulations show that the relative differences between the maximum and minimum velocities increase from ∼1% to ∼77% for plenum sizes W_H = W₀ = 0.8 m and W_H = W₀ = 0.1 m, respectively.

Figure 7:

Velocity distribution across the timber stack for plenum sizes W_H = W₀ = 0.10, 0.15, 0.20, 0.30, 0.50, and 0.80 m.

4.2.2. Lumber Moisture Content and Drying Rates

The physical model dictated by the governing equations and corresponding boundary conditions establishes that, in general, higher velocities yield higher convective and mass transfer coefficients at the lumber-air interface. In addition, the strong temperature dependency of the moisture diffusivity leads to an equally strong correlation between the air flow temperature and moisture content evolution along the drying process. Therefore, a physical understanding of the dry kiln behaviour regarding both effects is essential to defining drying strategies within practical industrial applications. This section is focused on such physical aspects, in which the effects of the air flow temperature and plenum size are investigated. The timber and other kiln/process parameters are defined in Table 2. The time step used in the simulations to solve the heat and mass conservation equations is Δt = 0.1 s.

The previous section features the effect of the plenum size upon the velocity distribution across the timber stack. The ensuing effect upon the drying process for a narrow plenum is shown in Figures 8 and 9 which presents the average moisture content, $\bar{M}$ , and the moisture content at z = H_T/2 (centre of the timber layer) for timber layers 1 (bottom of the stack) and 30 (top of the stack) corresponding to air flow temperatures 60°C, 90°C, and 120°C (333 K, 363 K, and 393 K). One could highlight two aspects regarding higher air flow temperatures: (i) drying rates increase and (ii) the effects of an inhomogeneous velocity distribution across the timber stack are magnified. Higher air flow temperatures cause the timber temperature and moisture diffusivity to increase, as (15) indicates, leading to higher moisture flow towards the timber-air interface with consequent increase in drying rates. Such higher moisture flow amplifies the sensitivity of the convective heat and mass transfer coefficients with regard to local velocities, leading in turn to greater differences of the drying rate across the timber stack. For instance, Figure 8 shows that the differences between the average moisture content for the 1st and 30th timber layers after 8 hours of drying time increases from 2% to 15% for air flow temperatures ranging from 60°C to 120°C (333 K to 393 K).

Figure 8:

Average moisture content for a plenum W_H = 0.1 m.

Figure 9:

Moisture content at z = H_T/2 for a plenum W_H = 0.1 m.

The heat and mass transfer coupled effects are featured in Figure 10 which shows the differences in heating and drying rates for the 1st and 30th timber layers for a narrow plenum. The aforementioned higher influence of the air flow temperature is expressed in higher heating and drying rates, that is, variations of the average temperature and moisture content within a single timber layer. Furthermore, the differences in air velocities across the timber stack (see Figure 7) cause substantial changes in heating and drying rates, especially early in the process.

Figure 10:

Heating rate, $d \bar{T} / d t$ , versus dehumidification rate, $| d \bar{M} / d t |$ , for air flow temperatures 60°C (333 K) and 120°C (393 K).

It is worthy to notice that, as the plenum size increases, differences in air velocity decrease, leading to homogeneous heating and drying rates across the timber stack. Figure 11 presents $d \bar{T} / d t \times | d \bar{M} / d t |$ for the 1st and 30th timber layers and plenum sizes W_H = 0.1, 0.2, and 0.8 m, in which iso-time curves are also shown. The larger thermal diffusivity causes timber layers to heat at higher rates than the corresponding dehumidification rates, that is, after 2 hours of drying time all timber layers reach approximately the air flow temperature, whereas drying rates, $| d \bar{M} / d t |$ , are still at higher levels. Figure 11 also shows that, for larger plenums, for example, W_H = 0.8 m, heating and drying rates across the timber stack are virtually uniform (the curves for 1st and 30th layers are almost coincident in the figure).

Figure 11:

Heating rate, $d \bar{T} / d t$ , versus dehumidification rate, $| d \bar{M} / d t |$ , for air flow temperature 120°C (393 K) and plenum sizes W_H = 0.1, 0.2, and 0.8 m.

The dynamics of the drying process is summarised in Figure 12 which shows $\bar{T} \times \bar{M}$ for larger plenums (uniform heating and drying rates across the timber stack) and air flow temperatures 60°C, 90°C, and 120°C (333 K, 363 K, and 393 K). For the sake of greater accuracy, the iso-time curves are computed using intermediate air flow temperatures in the range between 60°C and 120°C. Firstly, as discussed early in this section, it is possible to observe that higher air flow temperatures lead to higher drying rates; for example, after 24 hours of drying time the average moisture contents are 0.16, 1.50, and 4.62 kg/kg for air flow temperatures 120°C, 90°C, and 60°C, respectively. It is also interesting to notice that, in all cases, the timber layers quickly reach temperature equilibrium with air flow (around 4 hours of drying time).

Figure 12:

Evolution of the average timber temperature, $\bar{T}$ , and moisture content, $\bar{M}$ , for plenum W_H = 0.8 m.

5. Conclusions

Convective drying is one of the most widely used drying strategies in industry. The process is performed in dry kilns, in which heated air flows through a lumber stack comprising timber layers separated by stickers. The physics of the problem encompasses the hydrodynamic description of the air flow, conjugated heat, and moisture transfer within timber layers and convective heat and moisture transfer at lumber-air interface. In this work, the actual air flow distribution is determined using a system approach, which provides the kiln operating flow condition (the point at which the fan performance curve at a given rotational speed intersects the system characteristic curve). The air flow velocities across the timber stack at the operating point are used to compute the convective heat and moisture coefficients at the timber surface, which in turn are used to solve the fully coupled heat and moisture transfer problems.

This class of problems requires definition of geometric and operational dry kiln parameters able to provide both high drying rates and uniform moisture contents across the timber stack. The former is justified by economic reasons, whereas the latter is required with the objective of ensuring good quality of the final product (avoid warping, decolouration, etc.). The most critical aspect of this problem is the correct computation of the channel velocities across the stack height. Therefore, in order to ensure that the present approach is acceptable, channel velocity solutions for a reduced size dry kiln are compared with the Ansys CFX commercial code. With the objective of evaluating the uniformity of the drying process across the stack height, the effects of the kiln geometric configuration (plenum sizes) and air flow temperatures are investigated. The fully coupled simulations of a compact dry kiln have shown that higher air flow temperatures increase drying rates (a positive aspect); however, the effects of a nonuniform velocity distribution are magnified, that is, inhomogeneous drying rates across the timber stack (a negative aspect). Therefore, an industrial solution must seek a compromise between both effects, which may eventually lead to changing the sticker number and dimensions (to improve velocity uniformity), establish progressive air heating (to reduce the effect of inhomogeneous velocity distributions) and even replace the air circulating system (providing a different fan characteristic curve).

Footnotes

Nomenclature

Conflict of Interests

The authors affirm that they have no financial support or conflict of interests with manufacturers of any commercial software mentioned in the paper.

References

Shmulsky

and Jones

P. D.

, Forest Products & Wood Science, Wiley, Chichester, UK, 6th edition, 2011.

Mohan

V. P. C.

and Talukdar

, “Three dimensional numerical modeling of simultaneous heat and moisture transfer in a moist object subjected to convective drying,” International Journal of Heat and Mass Transfer, vol. 53, no. 21–22, pp. 4638–4650, 2010.

Riepen

and Paarhuis

, “Analysis and optimisation of the airflow distribution in convection kilns,” in Proceedings of the 1st Workshop on State of the Art for Kiln Drying, Edinburgh, UK, October 1999.

Kaya

Aydin

, and Dincer

, “Numerical modeling of heat and mass transfer during forced convection drying of rectangular moist objects,” International Journal of Heat and Mass Transfer, vol. 49, no. 17–18, pp. 3094–3103, 2006.

Kaya

Aydin

, and Dincer

, “Numerical modeling of forced-convection drying of cylindrical moist objects,” Numerical Heat Transfer Part A, vol. 51, no. 9, pp. 843–854, 2007.

Kaya

Aydin

, and Dineer

, “Heat and mass transfer modeling of recirculating flows during air drying of moist objects for various dryer configurations,” Numerical Heat Transfer Part A, vol. 53, no. 1, pp. 18–34, 2008.

Younsi

Kocaefe

Poncsak

, and Kocaefe

, “Computational and experimental analysis of high temperature thermal treatment of wood based on ThermoWood technology,” International Communications in Heat and Mass Transfer, vol. 37, no. 1, pp. 21–28, 2010.

Thibeault

Marceau

Younsi

, and Kocaefe

, “Numerical and experimental validation of thermo-hygro-mechanical behaviour of wood during drying process,” International Communications in Heat and Mass Transfer, vol. 37, no. 7, pp. 756–760, 2010.

Hirsch

, Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics, vol. 1, Butterworth-Heinemann, Oxford, UK, 2nd edition, 2007.

10.

Zdanski

P. S. B.

Ortega

M. A.

, and Fico

N. G. C. R.

Jr. , “Heat transfer studies in the flow over shallow cavities,” Journal of Heat Transfer, vol. 127, no. 7, pp. 699–712, 2005.

11.

Zdanski

P. S. B.

Vaz

, and Dias

A. P. C.

, “Forced convection heat transfer of polymer melt flow inside channels with con-traction/expansion sections,” International Communications in Heat and Mass Transfer, vol. 38, no. 10, pp. 1335–1339, 2011.

12.

Zdanski

P. S. B.

Vaz

Jr. Cerqueira

R. F. L.

, “Numerical simulation of the coupling between transient heat and mass diffusion process,” in Proceedings of the 7th National Congress of Mechanical Engineering, São Luis, Brazil, July 2012, (Portuguese).

13.

Nijdam

J. J.

and Keey

R. B.

, “The influence of kiln geometry on flow maldistribution across timber stacks in kilns,” Drying Technology, vol. 18, no. 8, pp. 1865–1877, 2000.

14.

Post

, Applied Computational Fluid Mechanics, Jones and Barlett, Sudbury, Mass, USA, 2011.

15.

Swamee

P. K.

and Jain

A. K.

, “Explicit equations for pipe-flow problems,” Journal of the Hydraulics Division, vol. 102, no. 5, pp. 657–664, 1976.

16.

Thumann

and Mehta

, Handbook of Energy Engineering, Fairmont Press, Lilburn, Ga, USA, 6th edition, 2008.

17.

Musch

H. E.

Barton

G. W.

Langrish

T. A. G.

, and Brooke

A. S.

, “Nonlinear model predictive control of timber drying,” Computers and Chemical Engineering, vol. 22, no. 3, pp. 415–425, 1998.

18.

Bergman

T. L.

Lavine

A. S.

Incropera

F. P.

, and DeWitt

D. P.

, Fundamentals of Heat and Mass Transfer, Wiley, Chichester, UK, 7th edition, 2011.

19.

Patankar

S. V.

, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC, USA, 1980.

Conjugated Heat and Mass Transfer in Convective Drying in Compact Wood Kilns: A System Approach

Abstract

1. Introduction

2. System Approach: The Flow Operating Condition

2.1. Equivalent Global Flow Circuit

2.2. Individual Channel Approximation

2.3. Inlet and Exit Plenum Approximation

2.4. Fan Performance Curve and Operating Point

3. Heat and Mass Transfer

4. Numerical Examples

4.1. Validation of the Flow Rate Model

4.2. Full Scale Simulation of a Dry Kiln

4.2.1. Hydrodynamic Solution: Velocity Distribution across the Timber Stack

4.2.2. Lumber Moisture Content and Drying Rates

5. Conclusions

Footnotes

Nomenclature

Conflict of Interests

References