John Eriksson, Sigurdur Ormarsson and Hans Petersson Department of Structural Mechanics
Chalmers University of Technology, Göteborg, Sweden e−mail: johne@sm.chalmers.se
ABSTRACT Introduction
Wood is a complex and moisture-sensitive material. This often causes problems in industry. In the sawmill industry, for example, moisture-related warping of sawn timber during drying is a serious problem.
Accordingly, in analysing distortion of timber boards or of wood products such as doors, window frames, and the like it is of substantial interest to have through knowledge of the distribution of the moisture content (MC) and of its variation over time. Simulations here need to take account of the fibre orientation, the orthotropy of the wood and variations in the material parameters as a function of moisture and temperature. This suggests use of full 3D non-linear diffusion model for analysing the moisture flow below the saturation point. Figure 1 shows the material directions l, r, t in wood. Most of simulation tools currently available are based on calculations of average MC or of one-dimensional moisture flow [1]. The influence of temperature, although it can be important, is not studied here.
Figure 1. The local and the global coordinate systems.
Theory
The non-linear diffusion equation for the moisture content w can be written as [2]
Q
where the source term Q, or the general moisture rate per unit volume and second, is added to the equation.
∇
∇
∇
∇w is the moisture gradient vector and Dw the diffusivity matrix. For anisotropic material the matrix Dw
generally has nine independent coefficients. In this study, the wood involved is assumed to be an orthotropic material and the Dw matrix to be given by
where the subscripts l, r and t stand for the longitudinal, radial and tangential directions, respectively.
The constitutive equation for moisture flow below the fibre saturation point is w
where q is the moisture flux vector. The amount of moisture passing through a unit surface area per unit of time is denoted as qn = qTn, where n is the outward normal to the surface, see [3] for details. The transformation between the orthotropic directions and the global coordinate system is given by
where l, r, t and i, j, k are the unit vectors along the axes of the local and the global coordinate systems, respectively. The transformation matrix A consists of nine direction cosines that take into account the effect of annual rings, spiral grain and the conicity of the log [4]. When the flow and the moisture gradient vectors are transformed to the local coordinate system, the following relationship is obtained:
q
Combining Eqs. (1.3)–(1.6) yields to the relationship between the global and the local diffusivity matrix.
w T w(w) AD (w)A
D = (1.7)
The moisture content can be approximated by w(x,y,z,t)=N(x,y,z)⋅a(t) where a is a vector containing all the moisture values at the nodal points and a&is the time derivative of the vector a.Multiplying Eq. (1.1) by a weight-function vector chosen according to the Galerkin method, and then integrating and applying the Green-Gauss theorem yields
a V and S refer to the volume and the boundary surface, respectively, of the body in question. Equation (1.8) can in brief form be written as
f a C
Ka+ &= (1.9)
If convection at the boundary is considered in Eq. (1.8), K and f need to be modified. The flux there can be expressed by the convection boundary condition qn = α(w-w ), where ∞ α is a convection coefficient dependent on a variety of conditions. The terms to be added due to the convection condition are
∫
αThe method chosen to solve the equation system (1.9) or (1.10) is a standard finite difference scheme. The time domain is divided up by different time points t1, t2, t3, …,tn, on the time axis. Various types of finite difference schemes with differing characteristics can be obtained, considering either Eq. (1.8) or Eq. (1.10) at time t = tn + β∆tn and choosing the parameter β so that 1 ≥ β ≥ 0. The moisture at time step n+1 is
In the examples given below β = 0.5 is selected which means that the scheme is unconditionally stable.
Implementation
Thus far, only 2D cases have been considered. Implementation is performed analogous with CALFEM, in a Matlab environment, using bilinear elements.
D(w)
The annual rings of the sawn timber boards studied are represented as ellipses, see Fig.3. The shape of the ellipses can be described by γ, the ratio of ar and br.
Figure 3. Material direction model in 2D, pith location and the lengths ar and br defining one annual ring.
Verification example
Verification can be achieved by comparing a case of transient moisture flow with a known analytical solution. The case studied is a rectangle with a constant initial MC of 1 and a boundary MC of 0. The solution to this problem can be found by use of eigenfunctions. A check of the analytical solution by use of Parseval’s identity was performed, see [7] for details. The material was assumed to be isotropic with a diffusivity of 7.78·10-7, and the dimensions of the rectangle were taken to be 0.1x0.05. A MC-profile at y = 0.025, t = 150 (∆t = 0.5) shows convergence with the FE solution, see Fig.4.
Figure 4. Analytical and FE solution of a transient moisture flow problem, (a) 140 elements (b) 480 elements.
Comparison for a 2D example of solutions using linear versus nonlinear diffusion coefficients.
The case studied was a 100x50 mm2 board with the diffusion coefficients [1] shown in Fig.5. The boundary condition was set to an MC of 8% and the initial condition to an MC of 27%.
Figure 5. Influence of moisture content on the diffusion coefficients.
For the MC range for which no measured data were available, a least square extrapolation fit was performed. The influence of pith location was also studied, the MC of two different boards dried in the same manner and equal in their material parameters but with differing pith locations being compared.
r
-- FE solution, 140 elements
-* Analytical solution (a) -* Analytical solution
x 10-1 x 10-1
(a)
(c)
Figure 6. Moisture content of boards subjected to a sudden drop in boundary MC and dried for 39 hours (a) Linear coefficients and pith in the middle of the cross-section
(b) Linear coefficients and pith at the centre of the lowest edge (c) Nonlinear coefficients and pith in the middle of the cross-section (d) Nonlinear coefficients and pith at the centre of the lowest edge.
In the simulations presented, circular annual rings were assumed (γ =1). The results shown in Fig.6(a)-(d) indicate pith position to be of importance for the MC profile and nonlinearity to affect the MC profile considerably. Table 1 summarises the results of the simulations performed.
When the MC of a node close to the boundary is plotted versus time for the four cases just considered, the node could be seen to dry faster at the beginning in cases (c) and (d) but to dry faster in cases (a) and (b) as time progresses. This result is expected on the basis of the variation in diffusivity parameters that is shown in Fig.5. If γ = 0.5, that is if the annual ring representation is elliptical and if all the other prerequisites are the same as in Fig.6(c), the average MC at time 39 hours increases, see Tab.1.
Table 1. Average MC for two pith locations as solved using a linear and a nonlinear technique and an example of elliptical annual rings as solved by a linear technique (the boundary nodes not being considered in the average).
Case Fig.6(a) Fig.6(b) Fig.6(c) Fig.6(d) Fig.6(a) but γ = 0.5
Average MC [%] 12.8 13.2 15.2 15.2 13.8
REFERENCES
[1] Anders Rosenkilde. Mätningar av fuktkvotsgradienter och ytfenomen vid virkestorkning. Licentiate thesis, Report TRITA-BYMA 1996:1, Stockholm, Byggnadsmaterial KTH (1996).
[2] J. Arfvidsson, J. Claesson. Isothermal moisture flow in building materials: modelling, measurements and calculations based on Kirchhoff’s potential. Building and environment 35 (2000) 519-536.
[3] N.S. Ottosen and H. Petersson. Introduction to the Finite Element Method. Prentice-Hall, New York, (1992).
[4] Sigurdur Ormarsson. Numerical analysis of moisture-related distortions in sawn timber. Doctoral thesis, Publication 99:7, Göteborg, Structural Mechanics, CTH (1999).
[5] Alf Samuelsson, Nils-Erik Wiberg. Finite element method. Studentliteratur, Lund, (1998).
[6] O.C.Zieniewicz, R.L.Taylor. The finite element method. volume 2, fifth edition, solid mechanics, pages 413-423. Butterworth Heinemann, (2000).
[7] Gunnar Sparr, Annika Sparr. Kontinuerliga system. Studentliteratur, Lund, (1999, 2000).
(b)
(d)