Wood Subjected to Hygro-Thermal and/or Mechanical Loads

The FV method was originally developed for fluid flow, heat and mass transfer calculations (Patankar, 1980), and later generalized for stress analysis in isotropic linear and non-linear bodies (Demirdžic & Muzaferija, 1994; Demirdžic et al., 1997; Demirdžic & Martinovic, 1993). For the purpose of the stress analysis in the wood, the method is modified to take into account the anisotropic nature of the wood and influence of the moisture content and the temperature on the deformation and stresses (Horman, 1999). Also, performance of the wood is found to be very sensitive to the moisture content and the temperature. Thus, it is of a great importance to be able to predict behavior of such materials under different hygrothermo-mechanical loads. In order to demonstrate the methods capabilities, a transient analysis of fields of temperature, moisture, and stresses and displacement in the wood subjected to hygro-thermal or mechanical loads is performed.


Introduction
The FV method was originally developed for fluid flow, heat and mass transfer calculations (Patankar, 1980), and later generalized for stress analysis in isotropic linear and non-linear bodies (Demirdžić & Muzaferija, 1994;Demirdžić et al., 1997;Demirdžić & Martinović, 1993). For the purpose of the stress analysis in the wood, the method is modified to take into account the anisotropic nature of the wood and influence of the moisture content and the temperature on the deformation and stresses (Horman, 1999). Also, performance of the wood is found to be very sensitive to the moisture content and the temperature. Thus, it is of a great importance to be able to predict behavior of such materials under different hygrothermo-mechanical loads. In order to demonstrate the methods capabilities, a transient analysis of fields of temperature, moisture, and stresses and displacement in the wood subjected to hygro-thermal or mechanical loads is performed.

Governing equations
The behaviour of an arbitrary part of a solid, porous body at any instant of time can be described by the following energy, mass and momentum balance equations which, when written in a Cartesian tensor notation, read: In these equation, t is time, i x is the Cartesian coordinate,  is the mass density, q c and m c are the specific heat and the specific moisture, T is the temperature, M is the www.intechopen.com moisture potential and i u is the displacement, q s and m s are the heat and mass source, i b is the body force, and j q , j m  and ij  are the heat and mass flux vector, and stress tensor components, respectively.

Constitutive relations
In order to close the system of Eqs.
(1)-(3) the constitutive relations for heat and mass flux based on the theory of Luikov (1966) which takes into account both the Soret and Duffort effect, together with the constitutive relation for a solid body are used:  for Eqs. (1) and (2) 11  11  12  31  11  11  11   22  12  22  23  22  22  22   33  31  23  33  33  33  33   12  44  12   23  55  23   31  66  31   000  000   000  000  00   0000  0  00000   AAA  T  M  AAA  T  M   AAA  T Note that the pair of constitutive Equations (4) and (5) can be extended to take into account the effect of the pressure gradient on the heat and mass transfer.
 for a thermo-elasto-plastic isotropic material for Eqs (3) and (1) In the case of elastic conditions, the expression within the brackets vanishes, and the constitutive realtion (11) reduces to the Duhamel-Neumann form of Hooke's law.

Initial and boundary conditions
In order to complete the mathematical model, initial and boundary conditions have to be specified. As initial conditions, the temperature, the moisture potential, and the displacement and velocity components have to be specified at all points of the solution domain.
For a wood heat treatment process, boundary conditions can be either of Dirichlet or Von Neuman type, i.e. temperature and/or heat flux and displacements and/or forces (surface tractions) have to be specified at all boundaries.
For a convective wood drying process, the following boundary conditions are normally appropriate: where q h and m h are the (convective) heat and mass transfer coefficients, respectively, si f is the surface traction, and all quantities are calculated at the solution domain boundary, except for those with subscript a which correspond to the ambient air.

Generic transport equation
Before the construction of a numerical algorithm is started, it is important to notice that the governing Eqs. (1)-(3) or Eqs. (1) and (3) when combined with constitutive Eqs. (7)-(9), or Eqs. (4) ( 0 j m   ) and (11) can be written in the form of the following generic transport equation: which can be integrated over an arbitrary solution domain V bounded by the surface A , with unit outer normal vector j n to yield: The generic variable  stands for T , M or i u .

Finite volume discretisation
As all numerical methods, the present one consists of time, space and equations discretisation. The time interval of interest is subdivided into a number of subintervals t  , not necessarily of the the same length. The space is discretised by a number of contiguous, non-overlapping hexahedral control volumes ( CV ), with the computational points at their centres (Fig.1). Then the integrals in generic Equation (17) are calculated by employing the midpint rule, the gradients are evaluated by assuming a linear variation of the dependent variable  between the computational points, and a fully implicit temporal scheme is employed. As a result a non-linear algebraic equation of the following form for each CV is obtained: where the coefficients K a and b are defined as: where the subscripts P and e denote values at the centre of the CV and at the centre of the east cell-face, respectively, e A is the area of the east cell face, V is the volume of the CV , e x  is the distance between points P and E , and all quantities refer to the current time level, except for those with the superscript o which refer to the previous, "old" time level.

Solution algorithm
After assembling Eqs. (18) for all CVs and for all transport equations, five (four in 2D case) sets of N mutually coupled non-linear algebraic equations are obtained, where N is the number of CVs. Those equations are solved by employing the following segregated iterative procedure.
First, all dependent variables are given their initial values. Then the boundary conditions which correspond to the first time step are applied, and the sets of equations for each individual dependent variable ( T , M , i u ) are linearised and temporarily decoupled by assuming that coefficient K a and source terms b are known (calculated by using depedent variable values from the previous iteration or the previous time step), resulting in a system of linear algbraic equations of the form: where p and q are typically of the order 3 10  , R  is a suitable normalisation factor and superscripts m and 1 m  denote values at two successive iterations.
In the next time step the whole procedure is repeated, except that the initial values are replaced by the values from the previous time step.
The present discretisation procedure ensures that the matrix A  has the folowing desirable properties: it is seven (five in 2D case) -diagonal, symmetric, positive definite and diagonally dominant, which makes Eq. (20) easily solvable by a number of iterative methods which retain the sparsity of the matrix A  . Note that it does not make sense to solve Eq. (20) to a tight tolerance since its coefficients and sources are only approximate (based on the values from the previous iteration/time step). Normally, reduction of the absolute residuals for one order of magnitude suffices.
The segregated solution strategy employed enables re-use of the same storage for the matrix A and vestor b for all depedent variables , thus requiring only 8N storage locations ( 6N in a 2D case). It is also important to mention that the fully implicit time differencing used, avoids stability-related time step restrictions. In principle, it allows any magnitude of the time step to be used, and in practice it is limited only by the required temporal accuracy.
When constitutive Eqs. (11) for a thermo-elasto-plastic isotropic material are applied, an elastic deformation is assumed at the beginning of iterations of each (load increment) time step (the expression within the brackets in Eqs. (11) is omitted). In the next iteration step in CVs in which the effective stress has reached the yield stress an elasto-plastic deformation is assumed and the expression within the brackets in Eqs. (11) is activated.
After each time step (load increment) displacements and stresses are updated adding displacement and stress increments in the current time step to the total displacements and total stresses from the previous time step. This procedure is repeated until the prescribed number of time steps (or load increments) is completed.

Application of the method
The method described in the previous sections has been applied to a number of both linear and non linear solid body deformation problems, few of which will be presented.

Numerical predictions of the wood drying process
The wood drying process is an important step in the manufacturing of wood products. During that process a non-uniform distribution of moisture content and temperature causes deformation and stresses in the wood and may result in a deformed and/or cracked end-product.
A wood drying process can be described as an unsteady process of heat, mass and momentum transfer in an orthotropic continuum with variable physical properties. The method solves a coupled set consisting of energy, moisture potential and momentum equations (1-3) with the constitutive relations (4-6).
Beech-wood beams (600x50x50 mm 3 ) are exposed to the uniform, unsteady flow of hot air in a laboratory dryer with an automatic control of the ambient air parameters (Horman, 1999).
The temperature and/or moisture dependent physical properties of the wood, obtained by fitting available experimental data, are given in Table 2. The others are considered constant and are given in Table 3. The timber is known to be cylindrically orthotropic. However, the wood samples used in this study are taken from the outer region of a cylindrical timber log and the rectilinear isotropy of samples is a reasonable assumption. At the beginning of the drying process the wood samples had a uniform distribution of temperature, moisture, displacement and velocity:  For the purpose of the numerical calculations the problem is considered to be a 2D plane strain problem. Due to the double symmetry, only one quarter of the cross-section is taken as the solution domain. For all calculations presented in this study, a uniform numerical mesh consisting of 20 x 20 CV was employed, while the time step was varied from 10 to 100 min (first seven time steps of 10 min, 31 time steps of 30 min, and finally 140 time steps of 100 min). These results are found to be grid and time independent by performing a systematic grid and time-step refinement (difference between the results on the 20 x 20 CV mesh differ from ones obtained on a 40 x 40 CV mesh for less than 1%, while the results obtained with 3 t   h practically coincide with results obtained with


During the initial phase of drying ( 0 2 t   h) the moisture content is above the fiber saturation point and the deformation is a consequence of the thermal stress only. Figure 2 shows the calculated fields at 70 t  min. One can see that an increase in temperature (Fig. 2a) causes the expansion of wood sample (Fig. 2b) and that the outer region is subjected to compressive and the inner region to extensive stresses (Fig. 2c, d). h) the deformation and stresses due to hygroscopic loads dominate. Figure 3 shows that at 108 t  h the moisture content has fallen below the fiber saturation point (Fig. 3a) and that this causes the shrinking of the wood sample (Fig. 3b). Around 100 t  h the stresses reach their maximum values and are extensive in the outer region and compressive in the interior of the sample (Fig. 3c, d). By comparing the values of stresses at 70 t  min and 108 t  h, it can be seen that the thermal stresses are around 200 times smaller than the stresses caused by the drop in the moisture content below the fiber saturation point.
If one plots the contours of the effective stress at 108 t  h, when it is at its maximum (Fig.  4), one can see that the effective stress is greater than the yield stress ( 10 y   MPa at 10% moisture; 6 y   MPa at 30%) only in a very narrow surface region (1mm deep), which indicates that the plastic defomation did not take place in the interior of the sample, and that the drying schedule is well designed.
At the end of the drying process ( 246 t  h), the moisture content in the sample varies from 11,1 to 14, 4 % (Fig. 5a), while Fig. 5b and 5c illustrate the anisotropy of the wood sample, the contraction is 1,3 mm in the x and 0,6 mm in the y direction, or 6,5% and 3,3% (axis x and y ).
In order to confirm the validity of the FV predictions, the calculated temperature, moisture and displacements are compared with experimental data (Horman, 1995., Institut für Holzphysik und mechanische Technologie des Holzes, Hamburg) at reference points (Fig. 6). Figures 7 and 8 show temperature and moisture content histories at two reference points. It can be seen a good agreement between calculations and experiment: maximum difference for both temperature and moisture was 8%, and the average difference was less than 2% (Martinović et al., 2001).  Figure 9 shows how the displacements at two points on the surface of the sample vary during the drying process. One can see very little deformation during the initial phase ( 1000  t min) and a considerable shrinking of the sample afterwards, and that predictions closely follow experimental data (maximum difference 15%, average difference 5%).

Numerical predictions of the wood heat treatment process
The prediction of temperature, stresses and displacements in logs during their thermal preparation in the veneers production (wood steaming) is an important step for designing satisfying heating regime of logs preparation, without damaging in wood. The equations governing heat and momentum balance (Eqs. (1) and (3)) with corresponding constitutive relations (Eq. 11) in thermo-elasto-plastic material are solved.
For a mathematical description of a thermo-elasto-plastic deformation of the body the incremental plasticity theory is applied. The problem is considered to be a 2D plane strain problem (Horman et al., 2003).
A beech log with a diameter of 0,42 m and length of 5,1 m was exposed to steam, which temperature history during the phases of heating up, through-heating and cooling down is in Fig. 10 depicted. For numerical calculations the heat transfer coefficient 7840 q h  W/m 2 K, and thermal and mechanical properties of the wood given in Table 4

Numerical analysis of stress and strain conditions of a three-dimensional furniture skeleton construction and its joints
At the design stage of some pieces of furniture, their complex skeleton construction is subjected to stress and strain analysis. That allows them to satisfy all the functional demands (comfort), aesthetic demands, but also the strength and stiffness both by their shape and their dimensions. To achieve that, it is necessary to carry out a numerical simulation of the stress of a complex construction.
The finite volume method is used in the calculation. Orthotropy of the wood material is accounted for by approximating it with an isotropic material whose elastic modulus E and Poisson's ratio ν are calculated by employing the least-square method.
The coefficients of the stiffness matrix ij A are given in Eqs. (10).
The physical model is angle 3D joint and skeleton construction chair (Fig. 13.) a) b) Fig. 13. a) The angle joint, b) the model of an examined chair www.intechopen.com Mechanical properties of wood, spruce, for temperature 20°C and moisture content 9,8 %, are given in Table 5. Mass density is 0,44 g/cm 3 .   (24) Fig. 15. a) Tangential stress at the plane xy , at the distance of 4,6 mm from the symmetry plane b) tangential stress at the plane xy and resulting stress at the planes xz (the place of osculation of the planes of the tenon) The highest value of the compressive stress xx  is in the symmetry plane, at undermost point of the tenon (~15,9 MPa), and the highest value of the tensile stress is at upper point of the tenon (~11,9 MPa). The place of the highest value of the effective stress is at the place of the highest compressive stress ( max 13, 5 eff   MPa). Tangential stress is presented on the plane xy at the distance of 4,6 mm from the symmetry plane. Figure 15a shows that the places of maximal stress ( max ~5 MPa) are at 35 x  mm. At the same plane xy and the plane xz , in the place of osculation of the planes of the tenon, tangential stress and resulting stress (normal yy  and tangential stress yx  ) are calculated, respectively and it are presented in Figure 15b.
The maximal stress is ~ 10 MPa and it can be seen at the under part of the tenon.
In the end, the stress -strain analysis is done for the symmetrical half of the loaded chair (Horman et al., 2010). Mass load of the horizontal underframe of the whole chair is 100 kg and of the vertical frame is 22 kg. Effective stress contours at the elements of the chair frame and at the joints of the highest stresses (

Conclusion
The presented finite volume method for solution of the problems of energy, mass and momentum balance in conjugation with heat and mass transfer in an anisotropic, elastoplastic, porous body is successful applied. Predictions of temperature, moisture content, strain and stress field in the wood drying as well as wood heat treatment process show high accurate results for course numerical grids due the second order accurate fully conservative spatial differencing scheme. The fully implicit unconditionally stable temporal differencing scheme enable large time steps during heat treatment processes. The applied finite volume discretisation procedure results in the diagonal dominant system of algebraic equations which are suitable for an iterative solution algorithm. The segregated iterative solution algorithm comprising the linearization and temporary decoupling of the system of equations for each dependent variable shows efficiency as well robustness solving highly nonlinear system of equations.