Model of moisture transport in wood, taking into account diffusion and the accompanying adsorption of water vapour through the skeleton

The paper presents a model of moisture transport in wood taking into account diffusion and the accompanying adsorption of water vapour through the skeleton. A two-parameter form of the source term was proposed, depending on the distance of the current mass concentration of bound water from the equilibrium state. The tests on cubic samples with a side of 2 cm were carried out which allowed to determine the coefficients of the proposed model on the basis of the reverse method. The tests were performed for pine, larch, oak and ash in all directions of orthotropy. Also the tests on thin samples were performed to verify the source term.


Introduction
The Wood is widely used in construction as a structural and finishing material. The properties of wood such as strength, dimensional stability and durability are strongly dependent on its humidity. This means that the spatial distribution of moisture and its changes versus time play a key role in predicting its behaviour. The most important from the point of view of material engineering is the transport of moisture in the range of the moisture content of the wood corresponding to dry conditions and the moisture content occurring at full saturation of the fibres, i.e. in the mass moisture range from 0% to 30% [1,2]. In this range of moisture in wood there is water in the form of water vapour and water bound by surface forces, and changes in its content lead to shrinkage or swelling of the wood [3].
One of the first models for predicting moisture distributions in the hygroscopic range under isothermal conditions are based on the balance of total moisture diffusion. These models, although they do not fully explain the physical problem and the results obtained by means of the diffusion equation cannot always be matched to experimental data, which was analysed by Shi [4], are attractive due to the small number of parameters needed to describe the phenomenon. On the other hand, if the diffusion coefficient is assumed to depend on the moisture content, a very good fit of the model to experimental data is obtained [5,6].
Recently, more advanced models have also been developed, which analyse the combined flow of vapour and bound water [7][8][9][10][11][12]. These models assume that vapour and bound water are transported by diffusion. However, the form of the diffusion coefficient is not entirely clear. Siau [13] assumes that the resistance to vapour flow is influenced by both cell lumens and cell walls. Therefore, the coefficient can only be considered as an apparent diffusion coefficient. In fact, it is assumed that the vapour passes through the cell walls. On the one hand, it is adsorbed on the cell wall and then desorbed on the other. However, this is only one theory of the mechanism of vapour transmission. According to Dinwoodie [14], the main vapour transport route leads through pits in cell walls. A diffusion coefficient independent of moisture content is then to be expected. This case was assumed and confirmed by the following analyses.
The issue of the transport of bound water itself is also not fully explained. According to Kärger and Valiullin [15] diffusion in mesopores is the result of the interplay of several mechanisms. Thus, at partial pore-space saturation one may distinguish between surface diffusion and the process of mass transfer through the pore gas phase. On the other hand Hozjan and Svensson [11] claim that the diffusion of bound water in the porous structure is very slow compared to the diffusion of vapour. Therefore, for moisture transport in larger wood samples, the transport of bound water is minor importance and may be omitted [11,16].
Another important issue is the possibility of surface resistance on boundary surfaces of wooden samples. This implies the adoption of more complex boundary conditions, e.g. boundary conditions of the third type instead of the first. This issue was considered experimentally by Rosen [17] analysing the effect of air velocity on water vapour adsorption for both longitudinal and transverse directions. He noted that for air velocities above 3 m/s, this effect is negligible and can be ignored.
The last issue, and probably the most important one, is the form of the source term. Its form has a very large impact on the results obtained from the model. A wrong one may give the apparent impression that, for example, the vapour diffusion coefficient is variable or that the transport of bound water is of great or minor importance.
The literature review shows that the topics related to moisture transport in the hygroscopic range are not fully explained and should be further analysed. This article proposes a model of moisture transport in wood that takes into account the diffusion and accompanying adsorption of water vapour through the skeleton. A new form of the source term has been proposed, which is an original contribution to modelling moisture flows in wood compared to those found in the literature. Its form was verified on thin pine samples for the relative humidity range from 30% to 50% and from 50% to 70%. The constancy of water vapour diffusion coefficients in the model for each of the orthotropic directions, boundary conditions of the first type for water vapour and boundary conditions with a source of moisture sorption for bound water were adopted. On the basis of measurements of changes in the mass of samples by minimising the target function, material coefficients were determined for the proposed model. The tests were performed on samples of the following dimensions 2 cm × 2 cm × 2 cm of pine, larch, oak and ash in the relative humidity range from 25% to 45%. The proposed model allowed to describe the process of moisture diffusion occurring in wood with a good fit for experimental research taking into account only 3 material parameters, i.e. the diffusion coefficient and two parameters describing the source term. Such a small number of model parameters allows to optimise the reverse methods of determining coefficients in wood and wood-based materials.

Transport equations
The moisture flow in wood in the hygroscopic range is mainly carried out by water vapour diffusion and sorption at the phase boundary between the air in the pores of the material and its skeleton. There is also a slight diffusion of bound water. However, it is much slower compared to the diffusion of water vapour in a porous structure [11] and thus negligible. For a one-dimensional case, without the influence of diffusion of bound water, the description of the phenomenon comes down to solving the following system of transport equations =̇. For equation (1) a uniform initial condition is adopted and boundary conditions of the first type where 0 is the water vapour density corresponding to the initial relative air humidity in which the wood was conditioned [kg/m3], is the water vapour density corresponding to the relative air humidity of the surrounding air [kg/m3], is the sample thickness [m]. For equation (2) a uniform initial condition is adopted and it is assumed that there is a boundary condition on the boundary of the material with the same source of moisture sorption as inside it. In equation (6) 0 is the mass concentration of bound water corresponding to the initial relative air humidity in pores [kg/kg]. Such boundary conditions are correct for ambient air velocities greater than 3 m/s [17], which in the studies carried out in this paper was met.
The change in average concentration of bound water can be determined directly by measuring the mass of water that enters the material. It is also possible to determine the spatial distribution of the local concentration in the sample. Isotopic or NMR methods are used [10,18]. The average concentration is linked to the local concentration by a relationship

Moisture equilibrium
Moisture adsorption is the process responsible for the bounding of water vapour molecules on the inner surface of the material. The amount of absorbed moisture is greater the higher the relative air humidity and less the higher the temperature [19]. By studying the changes in the absorbed moisture with respect to relative air humidity, a so-called sorption isotherm is obtained, which describes how much the material is capable of absorbing water at a constant relative humidity and a constant temperature to move to a state of equilibrium. This equilibrium is determined by comparing the chemical potential of water vapour and bound water. For vapour, this potential is well defined on the basis of thermodynamic considerations and has the form (e.g. in [19]) in which is the gaseous constant of water [J/(kg K], is temperature [K] and is the saturation pressure of water vapour [Pa]. Several forms of this potential have already been proposed for water in porous materials [20][21][22], but it does not seem that any of them is particularly preferred. This article proposes a new form in which is the maximum hygroscopic mass concentration of bound water [kg/kg] while and are material constants expressed in [J/kg] and [-] respectively. By comparing equations (8) and (9) and assuming that, for a state of equilibrium, the concentration becomes an equilibrium concentration an expression is obtained for the mass concentration at which the equilibrium between moisture balance between the material and the environment is achieved

Source term
As sorption stops at equilibrium, its rate should be a function of the difference between the equilibrium and the current mass concentration of bound water where is the coefficient describing the rate of adsorption [kg/(m3s]. However, the study [8] showed that the adsorption rate cannot be constant. This article adopts the original form of a source term similar to the function describing the chemical potential of associated water, with the difference that the maximum concentration was replaced by an equilibrium concentration where 0 is the coefficient of absorption rate [kg/(m3s] and is the material coefficient [-]. This form of function may be interpreted as depending on the difference between the equilibrium and the current concentration with a variable adsorption coefficient or as depending on the distance of the current bound water concentration to the equilibrium state.

Water vapour diffusion coefficient
The coefficient is the effective coefficient of water vapour diffusion, which may be related to the coefficient of water vapour diffusion in air according to the relationship = (13) in which is the dimensionless coefficient taking into account the resistance of the wood due to its porous structure and is the coefficient of water vapour diffusion in air. The coefficient adopted was as in the article [23] given as [24] ( ) = 2,16 · 10 −5 ( 273, 15 )

Verification of the adopted forms of sorption isotherm and the source term
The verification of the functions describing the sorption isotherm and source term was carried out on wood chips obtained from pine wood (Fig. 1). They were about 0.1 mm thick. This allowed, in the case of the determination of the sorption curve, to significantly shorten the testing time in relation to the testing of samples of several centimetres and, in the case of the source term, to minimise the influence of diffusion on the adsorption process. This type of testing is also used for other fibrous materials [25].

Figure 1. Thin wood samples used in the tests
Initially, the samples were dried at 60°. After drying, they were moved to a climatic chamber with constant climatic conditions = 25° and = 30%. After the mass of the samples was stabilised, the relative humidity in the chamber was changed to = 50%, and then the samples were weighed several times for 5 hours. After the mass of the samples was stabilised again (about a week), the relative humidity in the chamber was changed again to = 70% and the procedure was as in the previous stage. Once the mass was stabilised and the samples were weighed, they were placed in a exsiccator with relative air humidity ≈ 98% and temperature = 25°. This made it possible to determine the close to maximum concentration that can be achieved by the wood still in the hygroscopic range.
On the basis of the measurements of the change of masses, the concentrations of bound water at sequent moments were calculated according to the formula where and are the wet and dry masses of the sample. They were used to determine the parameters of the equation describing the sorption isotherm and the parameters of the model of moisture adsorption by wood.
The coefficients of equation (10) were determined by the Levenberg-Marquardt algorithm implemented in the Matlab package. The sorption isotherms for the determined parameters are shown in Figure 2.

Figure 2. Sorption isotherm of the tested wood
A very good fit of the model curve to the data obtained from the measurements confirmed that equation (9) describes the potential of bound water very well. In addition, it can be seen that in the range = 20 ÷ 60% the dependence of the mass concentration on relative air humidity is almost linear, which has been used in the calculation (shown in section 6) to determine the equilibrium concentration.
At this stage of research, a verification of the adopted source term was also carried out. For such thin samples, the water vapour density is almost immediately equalised in the entire pore volume. With a constant relative humidity of the surrounding air, we obtain, therefore, a homogeneous in the whole volume and constant in time distribution of the equilibrium mass concentration of bound water = ( ) = . This means that equation (2), taking into account (12), is reduced to the form of There is not analytical solution for the above equation. The solution was obtained using the finite difference method (FDM) with an explicit scheme. The coefficients of the equation were determined so that the values obtained according to equation (16) were as close as possible to those obtained during the experiment. This was achieved by finding the minimum of the function by means of a domain search. This function was in the form of where is the number of measurements taken, , and , are the mass concentrations measured and calculated according to the model, and is the vector of sought parameters ℎ 0 and . The graph of the change of the mass concentration of bound water versus time, obtained on the basis of the determined parameters, is shown in Figure 3.  Fig. 3 we can see a very good match between the values measured and obtained from the model in the humidity range from 30% to 70%. The material tested differed in terms of structure from the samples used in the main testing, but the results obtained confirm that the adopted form of the source term describes well the process of moisture adsorption by wood.

Description of main experiment
The main tests were performed on cubes with the following dimensions 2 cm × 2 cm × 2 cm, sourced from 4 types of wood (pine, larch, oak and ash), 36 pieces for each type. The samples were cut so that the growth rings are orthogonally aligned to their sides (Fig. 4). Initially, the samples were dried at 60° to determine the mass of dry material. They were then transferred to a climatic chamber with constant climatic conditions = 25° and = 25%. Once the masses were stabilised, the samples were insulated on 4 sides in such a way as to force a one-dimensional moisture flow. In this way 12 samples were obtained each of the three anatomical directions: radial, tangential and longitudinal ( Figure 5) and weighed again to calculate the mass of insulation used.

Figure 5. Insulated wood samples
The wooden cubes prepared in this way were placed in the climatic chamber where, after the masses were stabilised again, the air humidity was increased to = 45%. The change in the mass of the samples on the first day was measured every several hours. The time between sequent successive measurements was gradually extended to 4 days in the final stage of measurements. The measurements of mass change took 39 days (total approx. 3 months). The measured masses were used to determine the average mass concentrations of bound water in the samples according to the formula where and are the current mass and mass of dry sample respectively. Then, for each direction, 9 cubes were selected for which the results were the least divergent and the average concentration of the samples in the series was calculated for each time moment according to the formula where is the number of samples in the series. For selected cubes the average apparent density of dry wood was calculated . The density of the cellular substance for all types of wood is almost identical and is about 1500 kg/m 3 (assumed as in [26] citing as in [27]). This allows to calculate the approximate porosity of the tested wood from the formula The average densities and approximate porosities of the tested wood depending on its species are shown in Table 1. The data obtained in this way made it possible to determine the coefficients of the model proposed in the article.

Numerical solution
By discretising equations (1) and (2), using the implicit and centred finite difference scheme, we obtained where ∆ is the time increase, ∆ is the spatial increase, ( ) is the index referring to the current moment in time for which the saturation values are known, ( + 1) is the index referring to the future moment in time for which the saturation values are searching and ( − 1, ) and ( + 1) are the indices referring to spatial nodes. In order to obtain a solution to equations (21) and (22), it is more advantageous to present in a matrix form where The condition for the end of the calculation is in the form ‖ ( +1 )‖ < where is a small value that determines the permissible error or assuming in advance the number of required iterations.

Optimization procedure
The proposed model of moisture transport in wood has three parameters: two source term parameters , and the water vapour diffusion coefficient . However, the tests lasted 39 days, which is too little time for the samples to reach equilibrium. The final concentration (equilibrium concentration) that the samples would achieve after a very long time was also required. So eventually, the number of unknowns rose to four parameters.
The model coefficients were determined in such a way that the differences between the bound water mass concentrations determined during the experiment and calculated using the model were as small as possible. The target function was minimised by the domain method. In the case of a single series, the target function was in the form of where is the number of measuring moments, ̅ and ̅ are the mass concentrations measured and calculated according to the model in time and is a vector for the sought parameters: , , = , ∆ ∞ = − ( -last measured value).
Matching errors were calculated to evaluate the fitting of model results to the measured ones. Local error according to the formula and the global matching error according to the formula used in the article [28] = √ ∑ ( ̅ , − ̅ , ) The results for the pine samples are shown in Table 3 and in the form of graphs in Figure 6.  The diagrams in Figure 6 show a very good fit of the model to the results obtained from measurements for all three directions of orthotropy. However, it can be seen that the coefficients of the source term shown in Table 2 differ slightly (although they should be the same because the source term is independent of direction). The reasons include measurement errors and heterogeneity of wooden samples. In order to determine the parameters and which would describe with a good approximation the moisture flow in all directions, the next step was to look for the best fit of the results for the three directions simultaneously. Unfortunately in this case there are 8 coefficients to be determined. For this reason, this stage was divided into two parts. In the first one, for pre-determined diffusion coefficients and final mass concentrations, the parameters of the source term were determined by the domain method. In the second part, for the parameters set and , the remaining coefficients for the best fit were sought. These calculations were repeated several times until the search values stopped changing. As shown in Table 3, for the results obtained in this way, the matching errors hardly changed compared to step one. This is due to the fact that the first stage already gives a very good approximation of the final results. The algorithm for determining coefficients of model is also shown in Figure 7.

Results and Discussion
As a result of calculations carried out in accordance with the scheme presented in section 4, the parameter values of the proposed model of moisture transport in wood were obtained. They are presented in Tables 5÷8. For greater transparency and easier comparison of results, the diffusion coefficients are replaced by the vapour resistance factor for the successive orthotropic directions calculated according to the transformed formula (13)  (30) where is the number of samples in the series. All the diagrams show a very good match between the curves obtained from the model and the measured values. This is confirmed by very small values of matching errors: local (max( ) = 0,77% ÷ 1,71% ) and global ( = 0,37% ÷ 0,97% ). In addition, the model's matching to the measurements is similar for the entire course of the test, which indicates that the model describes the process well for moisture levels far from and close to equilibrium. The results of several researchers presented by Time [30] show that for spruce at = 25% the values of diffusion resistance in the longitudinal and transverse directions should be in the ranges ∈ 〈1,7; 5,0〉 and ∈ 〈11; 150〉. Similarly Krabbenhoft and Damkilde [8] citing for Siau [13] show that ∈ 〈1,7; 2,8〉 and ∈ 〈56; 156〉 , although for higher humidity. Diffusion resistances obtained from the model presented in Tables 5÷8 do not differ much from the values determined by these researchers.
The results of the calculations also confirm the correct selection of the form of the source term. The coefficients and calculated for the three directions separately (Table 2) are very similar. A small difference is due to the heterogeneity of individual samples. In turn, the values calculated for the three directions simultaneously also describe the process very well for each direction. An incorrect selection of the form of the mass source would take effect of different values of parameters and for each direction or in the impossibility of obtaining, with constant diffusion coefficients, a satisfactory match between the calculation and the measurement. All this confirms the suitability of the proposed model for modelling moisture transport in wood.

Conclusions
The model of moisture transport in wood proposed in the paper describes very well the diffusion process in the range of moisture present in the studies. This is confirmed by very small matching errors not exceeding for local 2% and for global 1%. The source term adopted in the model, although describing a very non-linear adsorption process and in a large humidity range, has only two parameters. The small number of unknowns makes it easy to determine the water vapour diffusion coefficients in wood from the reverse method. The main difference between this model and most of the literature is that the diffusion coefficients are constant and independent of the humidity of the wood and the surrounding air. However, the main tests were performed for a fairly narrow humidity range, and although tests performed on thin samples give a high probability of model accuracy for higher humidity ranges, it is necessary to extend the tests to a larger humidity range to allow more far-reaching conclusions to be drawn.