1. Introduction
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 moisture content 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] assumed 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. In this case, the diffusion coefficient should depend on the moisture content, e.g., in [
14]. However, this is only one theory of the mechanism of vapour transmission. According to Dinwoodie [
15], 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 [
16] 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,
17]. However, it should take into account that this assumption is valid only for relative air humidity below 65% [
18].
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 [
19] 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 three 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.
2. Model of Moisture Transport in Wood
2.1. 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 as mentioned in the introduction. 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:
in which
is the porosity of the material [m
3/m
3],
is the water vapour density [kg/m
3],
is the water vapour diffusion coefficient in the material in one of the orthotropic directions [m
2/s],
is the source term related to the moisture sorption process [kg/m
3s],
is the moisture content (mass of water divided by mass of dry material) [kg/kg] and
is the density of dry wood [kg/m
3].
For Equations (1) and (2) uniform initial conditions are adopted:
and boundary conditions of the first type for Equation (1):
where
is the water vapour density corresponding to the initial relative air humidity in which the wood was conditioned [kg/m
3],
is the moisture content corresponding to the initial relative air humidity in pores [kg/kg],
is the water vapour density corresponding to the relative air humidity of the surrounding air [kg/m
3],
is the sample thickness [m].
Such boundary conditions are correct for ambient air velocities greater than 3 m/s [
18], which in the studies carried out in this paper was met (3.5–4 m/s).
For Equation (2), no boundary conditions are defined because the model does not assume bound water transport. Its changes are related only to the adsorption of water vapour, which takes into account the source term.
The change in average moisture content 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 moisture content in the sample. Isotopic or NMR methods are used [
10,
20]. The average moisture content is linked to the local moisture content by a relationship:
2.2. 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 [
21]. 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 [
21]):
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 [
22,
23,
24], but it does not seem that any of them is particularly preferred. This article proposes a new form:
in which
is the maximum moisture content corresponding to fibre saturation point [kg/kg] while
and
are material constants expressed in [kg/J] and [-] respectively. By comparing Equations (8) and (9) and assuming that, for a state of equilibrium, the moisture content
becomes an equilibrium moisture content
an expression is obtained for the moisture content at which the equilibrium between moisture balance between the material and the environment is achieved:
2.3. Source Term
As sorption stops at equilibrium, its rate should be a function of the difference between the equilibrium and the current moisture content:
where
is the coefficient describing the rate of adsorption [kg/(m
3s].
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 moisture content was replaced by an equilibrium moisture content:
where
is the coefficient of absorption rate [kg/(m
3s] and
is the material coefficient [-]. This form of function may be interpreted as depending on the difference between the equilibrium and the current moisture content with a variable adsorption coefficient or as depending on the distance of the current moisture content to the equilibrium state.
2.4. 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:
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 [
25] given as [
26]:
3. 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 with density 459.6 kg/m
3 (
Figure 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 [
27].
Initially, the samples were dried at . After drying, they were moved to a climatic chamber with constant climatic conditions and . After the mass of the samples was stabilized (the mass changes of the samples were less than 0.002 g per 24 h), the relative humidity in the chamber was changed to and then the samples were weighed several times for 5 h. After the mass of the samples was stabilised again (about a week), the relative humidity in the chamber was changed again to 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 and temperature . This made it possible to determine the close to maximum moisture content that can be achieved by the wood still in the hygroscopic range.
On the basis of the measurements of the change of masses, the moisture content 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.
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
the dependence of the moisture content on relative air humidity is almost linear, which has been used in the calculation (shown in
Section 6) to determine the equilibrium moisture content.
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 moisture content
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 moisture contents measured and calculated according to the model, and
is the vector of sought parameters
and
.
The graph of the change of the moisture content versus time, obtained on the basis of the determined parameters, is shown in
Figure 3.
In
Figure 3 we can see a very good match between the values measured and obtained from the model in the humidity range from
to
. 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.
4. 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 (
Figure 4).
Initially, the samples were dried at
to determine the mass of dry material. They were then transferred to a climatic chamber with constant climatic conditions
and
. Once the masses were stabilized (the mass changes of samples were less than 0.002 g per 7 days), the samples were insulated on four 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.
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
. 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 moisture content 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 (only these were taken into account in the calculations). For these nine cubes the average concentration of the samples 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 wood substance for all types of wood is almost identical and depending on the research method [
28] oscillate around 1500 kg/m
3 (assumed as in [
29] citing as in [
30]). 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.
5. 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,
is the index referring to the future moment in time for which the saturation values are searching and
and
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
,
,
, but
, , are matrixes of coefficients depending on state variables (in the matrix
boundary conditions are taken into account),
is the vector of state variables in all spatial nodes for the
-th known point in time at which the initial conditions are taken into account,
is the saturation vector in all spatial nodes for the
-th search moment in time.
The value of the sought solution
is the argument
, , . Therefore it is necessary to approximate the solution. This approximation was performed using the Newton-Raphson iterative method. The values of the variable in the
-th iteration in the
-th time step has a form [
31]:
where:
The condition for the end of the calculation is in the form where is a small value that determines the permissible error or assuming in advance the number of required iterations.
Based on the above algorithm, in the Matlab environment an original program solving the system of Equations (1)–(6) was written.
6. 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 moisture content (equilibrium moisture content) 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 moisture content 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 moisture contents 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 [
32]:
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 moisture content, 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.
In this way, calculations were made for all four sample types and the final results are presented in the next section.
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
Table 4,
Table 5,
Table 6 and
Table 7. 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):
The results in graphical form are presented in
Figure 8,
Figure 9,
Figure 10 and
Figure 11. A model curve of the changes of moisture content
versus time and its measured values together with the standard deviation was plotted. The obtained results were also compared with the results for the multi-Fickian model (model m-F) presented in [
10] for parameters:
[1/s],
[-],
[1/s],
[-],
[-],
[-],
[-]. Standard deviation was calculated according to the formula:
where
is the number of samples in the series. Maximum and minimum values of standard deviation are presented in
Table 8.
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 () and global (). 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 [
33] show that for spruce at
the values of diffusion resistance in the longitudinal and transverse directions should be in the ranges
and
. Similarly Krabbenhoft and Damkilde [
8] citing for Siau [
13] show that
and
, although for higher humidity. Diffusion resistances obtained from the model presented in
Table 5,
Table 6,
Table 7 and
Table 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.
For comparison,
Figure 8 also shows the curves determined on the basis of the multi-Fickan model presented in [
10]. This model also describes the process quite well (especially for the radial and tangential directions), but the source term in this model has as many as seven parameters, which makes it very difficult to use it for the reverse method. The form of the source term proposed here has only two parameters, but the fit to the experimental data is even better. All this confirms the suitability of the proposed model for modelling moisture transport in wood.
8. 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 and for global . 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.