OPEN ACCESS
Moistureinduced swelling, or hygroexpansion, has been known to greatly deteriorate the durability of wood fiberpolymer composites (WPC). It is generally very expensive to perform experiments to completely obtain the diffusion kinetics as the process occurs in a very extended period of time. For the first time, we have developed a spacetime finite element algorithm that employs time discontinuous Galerkin (TDG) method for timedependent 3D hygromechanical behaviours of WPC. The formulation of matrix equations in spatial and temporal domains are explained in detail. A block GaussSeidel iterative method is used in the predictor/multicorrector multipass algorithm, which efficiently yields unconditionally stable and highorder accurate solutions. The model is validated by comparing the predicted timedependent hygroexpansion with that obtained in a previous experimental study. The quantitative analysis ensures the reliability of model, based on a Fickian diffusion process. With our adaptive timestepping scheme that bases on the embedded solution from the multipass iterations, the model efficiently progresses the kinetics with relatively large time steps. A runtime of a few hours compared to about three months of actual laboratory experimentation confirms the novelty and robustness of our model.
spacetime finite element method, timediscontinuous Galerkin, hygromechanical modeling, wood fiberpolymer composites
As is well known, wood fiber reinforced composites have been used in the building materials industry for several decades due to their low cost, light weight, biodegradability, abrasion resistance, and good specific properties. However, moistureinduced expansion, i.e., hygroexpansion, is one of the major drawbacks of woodbased composites. Unlike wood fibers, matrix components in woodbased composites, such as polymers, are much less susceptible to water absorption than wood [13]. The drastic difference in hygroelastic properties can even lead to serious problems in the case of wood fiberpolymer composites (WPC), such as buckling and warping of the installed WPC panels [4]. Therefore, studying the hygromechanical behaviours of WPC is an essential part of improving the durability of the materials.
It has been shown experimentally [58] that the diffusion kinetics of water in WPC takes place over a very long period of time  from weeks to a few months, depending on the diffusivity of the materials. Therefore, laboratory experiments are sometimes too expensive to obtain fully meaningful results. Consequently, hygromechanical modeling is becoming increasingly important. Numerous hygromechanical modeling studies are conducted for woodbased materials [912] and for polymeric materials [1315], but there is very little research on predicting the timedependent hygromechanical behaviours of WPC. There is particularly little research that attempts to accurately predict the threedimensional anisotropic responses of WPC subjected to water absorption.
Mbacké et al. [16] have used a numerical method to study the hygromechanical behaviours of homogenised composites, determining the timedependent mechanical properties as a function of the moisture content of the materials. However, the model ignores the complexity of the anisotropy of materials and uses simple linear interpolation to obtain the resulting properties, leading to inaccurate predictions in the case of WPC. In another study [17], a multiscale analytical model for stressdependent moisture diffusion in fiberreinforced polymer matrix composites was proposed under the assumption that the fibers do not absorb moisture. Thus, it is not applicable to the case of WPC where both constituents respond to moisture transport to different extents.
The finite element method (FEM) has been used extensively in the modeling and analysis of a composite structure with randomly distributed short fiber reinforced composites [18, 19]. Timedependent problems are usually solved using FEM by first discretising the spatial domain and solving the resulting ordinary differential equation in discrete time steps. This method is commonly known as semidiscretisation. Like heat transfer problems, moisture diffusion in materials can be considered as a firstorder parabolic problem, and perhaps the most commonly used methods for time integration are classified as singlestep/singlesolver methods (SS/SS) [20]. Such methods are based on solving a system of $n_{e q}$ linear equations leading to secondorder accuracy by unconditionally stable implicit SS/SS methods such as the trapezoidal method.
However, in applications where time is greatly extended, the SS/SS methods can lead to an accumulation of errors when using large time steps. The most efficient way to solve the above problems is to discretise the temporal domain with an accuracy no less than that of the spatial domain. This is also known as spacetime finite element method, which was developed by Hughes and Hulbert [21]. The method is based on the timediscontinuous Galerkin finite element method (TDGFEM), which is possibly considered one of the best methods for time integration in a discretised temporal domain. It provides relatively high accuracy with good stability. Due to the large linear spacetime linear systems of equations and the relatively high computational cost of solving nonsymmetric TDG matrices directly, several techniques [2224] have been developed. As a result, there are now several applications of TDGFEM [2528] However, to our knowledge, there are no studies using TDGFEM to investigate the hygromechanical behaviours of WPC.
The current study has two objectives: 1) to develop a complete and efficient TDGFEM algorithm written in MATLAB for hygromechanical behaviours of WPC and 2) to validate the model by comparing the solutions with corresponding experiments, investigating and discussing the reliability with respect to Fick’s theory.
The present paper is organised as follows. First, the formulation of the TDGFEM matrix equations in space and time for the hygromechanical behaviours of WPC is presented. Then, the predictor/multicorrector is described using a multipass algorithm from [24] and applied to the firstorder parabolic moisture diffusion problem. A suitable adaptive timestepping scheme is briefly explained. Then, the overall algorithm with the necessary defined material parameters is used to calculate the moisture evolution and the resulting hygroexpansion of the composites. The convergence of the mesh is considered. The efficiency of the timestepping scheme using the embedded solution approach is discussed. It is followed by model validation and quantitative analysis of model reliability.
The general diffusion equation of moisture in a given material according to Fick’s second law reads:
$\frac{\partial W}{\partial H} \frac{\partial H}{\partial t}=\nabla \cdot\left(D_{h} \nabla H\right)$ (1)
where, ∂W/∂H is the moisture capacity that can be obtained through the equilibrium adsorption isotherm of the material, H represents the relative humidity ranging from zero to one, ∇H denotes the gradient vector of humidity. Also, D_{h} is the diffusion coefficient of moisture in the material or moisture diffusivity in m^{2}/s. In a 3D Cartesian system or X space:
$D_{h}=\left[\begin{array}{lll}D_{h, x} & 0 & 0 \\ 0 & D_{h, y} & 0 \\ 0 & 0 & D_{h, z}\end{array}\right]$ (2)
where, D_{h,x}, D_{h,y} and D_{h,z} are the diffusion coefficients in the respective directions, which are are assumed to be constant. Therefore, Eq. (1) can be rewritten for the 3D diffusion equation as:
$\frac{\partial H}{\partial t}=\frac{1}{\partial W / \partial H}\left(D_{h, x} \frac{\partial^{2} H}{\partial x^{2}}+D_{h, y} \frac{\partial^{2} H}{\partial y^{2}}+D_{h, z} \frac{\partial^{2} H}{\partial z^{2}}\right)$ (3)
The moisture absorption gives rise to hygroexpansion of material via:
$\boldsymbol{\varepsilon}_{h}=\left\{\begin{array}{l}\varepsilon_{h, x} \\ \varepsilon_{h, y} \\ \varepsilon_{h, z} \\ 0 \\ 0 \\ 0\end{array}\right\}=\Delta H\left\{\begin{array}{l}\beta_{x} \\ \beta_{y} \\ \beta_{z} \\ 0 \\ 0 \\ 0\end{array}\right\}$ (4)
where, ΔH is the change in humidity content, and β_{x}, β_{y} and β_{z} are the directional hygroexpansion coefficients of material. The hygroexpansion strains do not contribute to shear deformation, leading to zero shear strains in Eq. (4). Consequently, the hygroelastic stressstrain relationships in X of the total strain, ε, elastic strains, ε^{m} can, therefore, be written as:
$\sigma=D_{m} \varepsilon_{m}=D_{m}\left(\varepsilon\varepsilon_{h}\right)$ (5)
where,
$\sigma={ }^{T}\left\{\begin{array}{llllll}\sigma_{x} & \sigma_{y} & \sigma_{z} & \tau_{x y} & \tau_{y z} & \left.\tau_{z x}\right\}\end{array}\right.$ （6）
and
$\varepsilon_{m}={ }^{T}\left\{\frac{\partial u}{\partial x}, \frac{\partial v}{\partial y}, \frac{\partial w}{\partial z}, \frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}, \frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}, \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right\}$ (7)
The symbol ^{T} denotes transpose operation, σ denotes normal stresses, τ denotes shear stresses, u, v and w are the components of strain in x, y and z respectively, D^{m} is defined as:
$D_{m}$
$=\left[\begin{array}{llllll}1 / E_{x} & v_{x y} / E_{x} & v_{x z} / E_{x} & 0 & 0 & 0 \\ v_{y x} / E_{y} & 1 / E_{y} & v_{y z} / E_{y} & 0 & 0 & 0 \\ v_{z x} / E_{z} & v_{z y} / E_{z} & 1 / E_{z} & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 / G_{y z} & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 / G_{z x} & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 / G_{x y},\end{array}\right]^{1}$ (8)
where, E’s and G’s are elastic and shear moduli. v’s are Poisson’s ratios.
The equilibrium equations of a differential volume element that are satisfied over the physical domain are given in a tensor notation as:
$\sigma_{i j, j}+f_{i}=0$ (9)
where, σ_{ij} are the stress components in the Cartesian stress tensor, f_{i} is the component of the body force vector along the i direction, and j represents partial differentiation with respect to the j coordinate. Respectively, the essential and natural boundary conditions are:
$u_{i}(X)=\bar{u}_{i}$; on $\Gamma_{u}$ (10)
$\sigma_{i j} \hat{n}_{j}=\bar{t}_{i} ;$ on $\Gamma_{t}$ (11)
where, $u_{i}(X)$ is the displacement along the i direction of the point at coordinates X, $\bar{u}_{i}$ is the corresponding specified boundary displacement at the point, $\hat{n}_{j}$ is the component of unit normal vector in the j direction, and $\bar{t}_{i}$ is the specified boundary stress along the i direction.
Therefore, the governing equations of threedimensional hygromechanical behaviours in a matrix form can be obtained by substituting the corresponding stress components from Eq. (6) into that of Eq. (9).
There is no or very little commercial software capable of modeling such a timedependent moisture diffusion problem using highaccuracy adaptive time integration methods. Therefore, to simulate the hygromechanical behaviours of WPC, we have developed a complete algorithm written in MATLAB [29]. With the following series of FE formulations in space and time, the overall algorithmic structure is summarised in Figure 1.
Figure 1. Overall FEbased algorithm. ← and → denote input and output respectively
3.1 3D FE formulation in spatial domain
For each element in space, i.e. Ω^{e}, an eightnode hexahedron with natural coordinates Ξ as shown in Figure 2 is used. Trilinear Lagrange interpolation functions (ξ_{j}, η_{j}, ζ_{j}) of local node j are applied:
$\Psi_{j}^{e}(\Xi)=\frac{1}{8}\left(1+\xi \xi_{j}\right)\left(1+\eta \eta_{j}\right)\left(1+\zeta \zeta_{j}\right)$ (12)
The element displacements, i.e., U^{e}=^{T}{u^{e} v^{e} w^{e}}, are interpolated in terms of the shape functions in Eq. (12) and nodal displacements $u_{j}^{e}, v_{j}^{e}$ and $w_{j}^{e}$ as:
$u^{e}=\sum_{j=1}^{8} \Psi_{j}^{e} u_{j}^{e}, v^{e}=\sum_{j=1}^{8} \Psi_{j}^{e} v_{j}^{e}, w^{e}=$$\sum_{j=1}^{8} \Psi_{j}^{e} w_{j}^{e}$ (13)
The values of H are also interpolated from the nodes, in which semidiscretisation is used. The approximated humidity over Ω^{e} is as follows:
$H^{e}(X, t)=\Psi^{e}(X) H^{e}(t)$ (14)
where, H^{e}(t) contains the nodal humidity values at time t. The first order time derivative of Eq. (14) is, therefore, obtained as:
$\frac{\partial}{\partial t} H^{e}(X, t)=\Psi^{e}(X) \frac{\partial H^{e}(t)}{\partial t}=\Psi^{e}(X) \dot{H}^{e}(t)$ (15)
Figure 2. An eightnode hexahedron (brick) element, Ω^{e}, with the coordinate systems. Cartesian (x, y, z) or X (left), and natural (ξ, η, ζ) or Ξ (right) with nodal coordinates
With x, y, z physical coordinates of node j in X, the Jacobian matrix, J^{e} can be calculated as:
$J^{e}=\left[\begin{array}{lll}\Sigma_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \xi} x_{i} & \Sigma_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \xi} y_{i} & \Sigma_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \xi} z_{i} \\ \Sigma_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \eta} x_{i} & \Sigma_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \eta} y_{i} & \sum_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \eta} z_{i} \\ \Sigma_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \zeta} x_{i} & \Sigma_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \zeta} y_{i} & \Sigma_{i=1}^{8} \frac{\partial \Psi_{j}^{e}}{\partial \zeta} z_{i}\end{array}\right]$ (16)
The invertibility of Jacobian matrix in Eq. (16) is ensured by setting the initial shapes of all elements as close to a square as possible with all internal angles being 90 degrees. Additionally, moistureinduced strains do not contribute to shear deformation regarding to Eq. (4), meaning that anomalous distortion or twisting of an element is not possible.
So, the derivatives in X can be found via:
$\frac{\partial}{\partial X}=J^{e1} \frac{\partial}{\partial \Xi}$ (17)
The essential boundary conditions of the nodes on surface boundaries Γ are:
$\mathbf{U}(X, t)=\overline{\mathbf{U}} ;$ on $\Gamma_{u}$ (18)
$\mathbf{H}(X, t)=\overline{\mathbf{H}} ;$ on $\Gamma_{h}$ （19）
where, $\overline{\mathbf{U}}$ and $\overline{\mathbf{H}}$ are specified displacements and humidities imposed on the nodes of corresponding boundaries. Also, the initial conditions are:
$\mathbf{U}(X, t=0)=0 ; \mathbf{H}(X, t=0)=0 ;$ over $\Omega$ (20)
By Galerkin Method of Weighted Residuals, the weakform of Eq. (9) can be obtained through:
$\int_{\Omega} \delta U_{i}\left(\sigma_{i j, j}+f_{i}\right) d \Omega=0$ (21)
where, δU_{i} denotes the weight functions on the i direction. Integration by parts yields:
$\int_{\Omega} \delta U_{i, j} \sigma_{i j} d \Omega=\int_{\Omega} \delta U_{i} f_{i} d \Omega+\int_{\Gamma_{\mathrm{t}}} d U_{i} \overline{\mathrm{t}}_{\mathrm{i}} \mathrm{d} \Gamma$ (22)
In this work, bodyforce is neglected (f_{i}=0), and there are no mechanical tractions $\left(\bar{t}_{i}=0\right)$ on the boundaries. The weakform of the original problem in Eq. (9) is now simplified as:
$\int_{\Omega} \delta U_{i, j} \sigma_{i j} d \Omega=0$. (23)
Note that the nontrivial solutions arise from the moistureinduced strains from Eq. (5) and by substituting Eq. (5) into Eq. (23), the general 3D FEdiscretised linear elasticity equation can be written as:
$K_{m}^{e} d^{e}=F_{m}^{e}$ （24）
where,
$K_{m}^{e}=\int_{\Omega^{e}}{ }^{T} B_{m}^{e} D_{m}^{e} B_{m}^{e}\leftJ^{e}\right d \xi d \eta d \zeta$ (25)
$d^{e}=T\left\{u_{1}^{e} v_{1}^{e} w_{1}^{e} u_{2}^{e} v_{2}^{e} w_{2}^{e} \ldots u_{8}^{e} v_{8}^{e} w_{8}^{e}\right\}$, (26)
$F_{m}^{e}=\int_{\Omega e}{ }^{T} B_{m}^{e} D_{m}^{e} \varepsilon_{h}^{e}\leftJ^{e}\right d \xi d \eta d \zeta$. (27)
J^{e} is the determinant of J^{e}. $B_{m}^{e}$ is the straindisplacement transformation matrix given as:
$\underset{6 \times 24}{B_{m}^{e}}=\left[\begin{array}{ccccc}{\left[B_{m, 1}^{e}\right]} & {\left[B_{m, 2}^{e}\right]} & {\left[B_{m, 3}^{e}\right]} & \ldots & {\left[B_{m, 8}^{e}\right]} \\ 6 \times 3 & 6 \times 3 & 6 \times 3 & & 6 \times 3\end{array}\right]$ (28)
From $\varepsilon^{e}=B_{m}^{e} d^{e}$ and with Eqns. (7), (14), (16), (17) and (26), each $B_{m, i}^{e}$ in Eq. (28) is given as:
$\boldsymbol{B}_{m, i}=\left[\begin{array}{lll}\frac{\partial \Psi_{i}}{\partial \xi} \frac{\partial \xi}{\partial x} & 0 & 0 \\ 0 & \frac{\partial \Psi_{i}}{\partial \eta} \frac{\partial \eta}{\partial y} & 0 \\ 0 & 0 & \frac{\partial \Psi_{i}}{\partial \zeta} \frac{\partial \zeta}{\partial z} \\ 0 & \frac{\partial \Psi_{i}}{\partial \zeta} \frac{\partial \zeta}{\partial z} & \frac{\partial \Psi_{i}}{\partial \eta} \frac{\partial \eta}{\partial y} \\ \frac{\partial \Psi_{i}}{\partial \zeta} \frac{\partial \zeta}{\partial z} & 0 & \frac{\partial \Psi_{i}}{\partial \xi} \frac{\partial \xi}{\partial x} \\ \frac{\partial \Psi_{i}}{\partial \eta} \frac{\partial \eta}{\partial y} & \frac{\partial \Psi_{i}}{\partial \xi} \frac{\partial \xi}{\partial x} & 0\end{array}\right]$ (29)
Similarly, by using the above approach, the FE matrix equations of moisture diffusion can be obtained. With δH being the weight functions of H, and by substitution of Eq. (14) and Eq. (15) into Eq. (3), and no moisture flux being posed to any boundary, the following firstorder parabolic moisture diffusion equation can be obtained:
$M^{e} \dot{H}^{e}(t)+D^{e} H^{e}(t)=0$ (30)
where,
$M^{e}=\int_{\Omega^{e}}{ }^{T} \Psi^{e} \Psi^{e}\leftJ^{e}\right d \xi d \eta d \zeta$ (31)
$D^{e}=\int_{\Omega^{e}}{ }^{T} B_{h}^{e} D_{h}^{e} B_{h}^{e}\leftJ^{e}\right d \xi d \eta d \zeta$ (32)
$H^{e}={ }^{T}\left\{\begin{array}{llll}H_{1}^{e} & H_{1}^{e} & H_{3}^{e} \ldots & H_{8}^{e}\end{array}\right\}$ (33)
$\dot{H}^{e}=^{T}\left\{\dot{H}_{1}^{e} \quad \dot{H}_{1}^{e} \quad \dot{H}_{3}^{e} \ldots \quad \dot{H}_{8}^{e}\right\}$ (34)
$B_{h}^{e}$ is the transformation matrix used to obtain moisture gradients, which can be defined as:
$\begin{aligned}&B_{h}^{e} \\&3 \times 8\end{aligned}=\left\{\begin{array}{l}\frac{\partial}{\partial \xi} \frac{\partial \xi}{\partial x} \\ \frac{\partial}{\partial \eta} \frac{\partial \eta}{\partial y} \\ \frac{\partial}{\partial \zeta} \frac{\partial \zeta}{\partial z}\end{array}\right\}\left\{\begin{array}{lllll}\Psi_{1} & \Psi_{2} & \Psi_{3} & \ldots & \left.\Psi_{8}\right\}\end{array}\right.$ (35)
Upon the assembly operation for Eqns. (24) and (30) over Ω, the global matrix equations of hygromechanical behaviours can be given as:
$\mathbf{K} \mathbf{U}=\mathbf{F}$ (36)
and
$\mathbf{M} \dot{\mathbf{H}}+\mathbf{D} \mathbf{H}=\mathbf{0}$ (37)
where, K, F, M and D respectively denote the global stiffness matrix, global force matrix, global mass matrix and global diffusivity matrix, while vectors $\dot{\mathbf{H}}$, H of size n_{eq} and vector U of size 3 n_{eq} contain solutions of humidity rate, humidity and displacements respectively.
3.2 TDG formulation in temporal domain
To solve Eq. (37), we employ timediscontinuous Galerkin finite element method (TDGFEM). This method allows unconditionally stable, highly accurate time integration via using highorder p polynomial approximation with the accuracy order of 2p+1 at the end of each time step.
The formulation of TDGFEM is developed by partitioning the time domain for each time step, $t \in$=(0, t_{f}] into N intervals, I_{n}=(t_{n}, t_{n+1}), with time step size of Δt_{n}=t_{n+1}t_{n}. The trial solutions are members of the set of polynomial functions of order p:
$\mathcal{V}^{h}=\left\{\mathbf{H}^{h} \in \bigcup_{n=0}^{N1}\left(P^{p}\left(I_{n}\right)\right)^{n_{e q}}\right\}$ (38)
The method allows unknowns H^{h}(t) to be discontinuous across the interface of time intervals, and the approximations at time t_{n} are written as:
$\mathbf{H}^{h}\left(t_{n}^{+}\right)=\lim _{\delta t \rightarrow 0} \mathbf{H}^{h}\left(t_{n}+\delta t\right) ; \delta t>0$ (39.1)
$\mathbf{H}^{h}\left(t_{n}^{}\right)=\lim _{\delta t \rightarrow 0} \mathbf{H}^{h}\left(t_{n}\delta t\right) ; \delta t>0$ (39.2)
An illustration of time discontinuous approximation is provided in Figure 3.
Figure 3. Time discontinuous approximation of solutions
In such discontinuous solutions, the jump operator across each time interface t_{n} is defined as:
$\mathbf{H}^{h}\left(t_{n}\right)=\mathbf{H}^{h}\left(t_{n}^{+}\right)\mathbf{H}^{h}\left(t_{n}^{}\right)$ (40)
where, n=0, 1, 2, …, N1. The principle of TDGFEM [21] is to find $\mathbf{H}^{h}(t) \in \mathcal{V}^{h}$ such that for all $\mathbf{W}^{h}(t) \in \mathcal{V}^{h}$, such that:
$a\left(\mathbf{W}^{h}, \mathbf{H}^{h}\right)_{n}=l\left(\mathbf{W}_{n}^{h}\right)$ (41)
where,
$a\left(\mathbf{W}^{h}, \mathbf{H}^{h}\right)_{n}:=\left(\mathbf{W}^{h}, \mathbf{M} \dot{\mathbf{H}}^{h}\right)_{I_{n}}+\left(\mathbf{W}^{h}, \mathbf{D} \mathbf{H}^{h}\right)_{I_{n}}$$+\mathbf{W}^{h}\left(t_{n}^{+}\right) \cdot \mathbf{M H}^{h}\left(t_{n}^{+}\right)$ (42)
and
$l\left(\mathbf{W}_{n}^{h}\right):=\mathbf{W}^{h}\left(t_{n}^{+}\right) \cdot \mathbf{M} \mathbf{H}^{h}\left(t_{n}^{}\right)$ (43)
in which H^{h}(t_{0})=H_{0}. The L_{2}inner product for the time interval can be obtained as:
$\begin{aligned}(\mathbf{W}, \mathbf{H})_{I_{n}}=\int_{I_{n}} \mathbf{W H} d t &=\lim _{\delta t \rightarrow 0} \int_{t_{n}+\delta t}^{t_{n}+1\delta t} \mathbf{W H} d t \\ \delta t &>0 \end{aligned}$ (44)
From Eq. (43), it can be seen that the solutions of $\mathbf{H}^{h}\left(t_{n}^{}\right)$ from the end of the previous time step are to be used as initial conditions for the current time step. The approximations to H(t) are given as:
$\mathbf{H}^{h}(t)=\mathbf{H}_{n}^{}+\left[\boldsymbol{\varphi}(t) \otimes \mathbf{I}_{n_{e q}}\right] \overline{\mathbf{H}}, t \in\left(t_{n}, t_{n+1}\right)$ （45）
where, $\mathbf{H}_{n}^{}$ is a column vector of size n_{eq}, containing n_{eq} humidity values at the end of the previous time step, $\boldsymbol{\varphi}(t)=\left[\varphi_{0}(t), \varphi_{1}(t), \ldots, \varphi_{p}(t)\right]$ a row vector of size p+1, containing the time shape functions of each polynomial order, from 0 to p. $\overline{\mathbf{H}}=T\left[\overline{\mathbf{H}}_{0}, \overline{\mathbf{H}}_{1}, \ldots, \overline{\mathbf{H}}_{p}\right]$  a column vector of size (p+1) n_{eq}, in which, each $\overline{\mathbf{H}}_{i}$ denotes the vector of undetermined coefficients corresponding to humidity, $\otimes$ represents the Kronecker (tensor) product, $\mathbf{I}_{n_{e q}}$ is the identity matrix of order n_{eq}.
As to achieve an optimally efficient iterative algorithm with highorder accuracy, the time shape functions are obtained thereby mapping of monomial basis {1, t, t^{2}, …, t^{p}}. The time interval, t is normalized to an interval $\hat{t} \in(0,1)$, using linear scaling, i.e. $t=t_{n}^{+}\left(\Delta t_{n}\right) \hat{t}$. For choices of time shape functions, we adopt those defined in Eq. (56) in a previous study [24], which reads:
$\varphi_{i}(\hat{t})=\sum_{k=0}^{i}(1)^{k} k ! \hat{\gamma}^{k}\left(\begin{array}{l}i \\ k\end{array}\right)^{2} \hat{t}^{ik}$ (46)
for i=0, 1, 2, …, p and $\left(\begin{array}{l}i \\ k\end{array}\right)$ denotes the binomial coefficient, $\hat{\gamma}=1 /(2 p+1)$.
By the substitution of the trial functions from Eq. (45) into Eq. (41), the coupled matrix equations for each time interval are obtained as:
$(A \otimes \mathbf{M}+B \otimes \mathbf{D}) \overline{\mathbf{H}}=b \otimes \mathbf{M} \mathbf{H}_{n}^{}$ (47)
where, A and B, of size (p+1)×(p+1), are defined as follows:
$A=\int_{0}^{1} T_{\boldsymbol{\varphi}}(\hat{t}) \dot{\boldsymbol{\varphi}}(\hat{t}) d \hat{t}+{ }^{T} \boldsymbol{\varphi}(0) \boldsymbol{\varphi}(0)$ （48）
$B={ }^{T} B=\Delta t_{n} \int_{0}^{1}{ }^{T} \boldsymbol{\varphi}(\hat{t}) \boldsymbol{\varphi}(\hat{t}) d \hat{t}$ (49)
and
$b=\Delta t_{n} \int_{0}^{1}{ }^{T} \boldsymbol{\varphi}(\hat{t}) d \hat{t}$ (50)
From above, the number of linear equations in Eq. (47) is (p+1) times larger than the system of the second order accurate SS/SS algorithm such as that of generalized trapezoidal method. Plus, direct solutions of the large system of coupled equations $\overline{\mathbf{H}}_{i}$ require high computational costs. In this work, we have adapted the numerical implementation proposed in [24] to solve firstorder parabolic moisture diffusion equations. Taking Eq. (47) and multiplying both sides by ($A^{} \otimes \mathbf{I}_{n_{e q}}$) yields:
$\left(I_{p+1} \otimes \mathbf{M}+W \otimes \mathbf{D}\right) \overline{\mathbf{H}}=d \otimes \mathbf{M H}_{n}^{}$ (51)
where, I_{p}_{+1} is the identity matrix of order p+1, W=A^{1}B and
$d=A^{1} b=\Delta t_{n} \int_{0}^{1}={ }^{T}\left\{\hat{\gamma} \Delta t_{n}, \Delta t_{n}, 0,0, \ldots, 0\right\}$ (52)
Let D and E be the diagonal and offdiagonal parts of W; W=D+E, where $D=\gamma I_{p+1}=\widehat{\gamma} \Delta t I_{p+1}$, and
$E=\left[\begin{array}{llllll}0 & E_{01} & E_{02} & E_{03} & \ldots & E_{0 p} \\ E_{10} & 0 & E_{12} & E_{13} & \ldots & E_{1 p} \\ E_{20} & E_{21} & 0 & E_{23} & \ldots & E_{2 p} \\ E_{30} & E_{31} & E_{32} & 0 & \ldots & E_{3 p} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ E_{p 0} & E_{p 1} & E_{p 2} & E_{p 3} & \ldots & 0\end{array}\right]$ (53)
To develop a predictorcorrector algorithm, the lefthand side of Eq. (51) is rewritten as:
$\left(I_{p+1} \otimes[\mathbf{M}+\gamma \mathbf{D}]+E \otimes \mathbf{D}\right) \overline{\mathbf{H}}=d \otimes \mathbf{M H}_{n}^{}$ (54)
where, $\gamma=\widehat{\gamma} \Delta t$. To achieve for an efficient multipass algorithm, a GaussSeidel block matrix iterative scheme is employed. This is done by splitting matrix E into upper and lower triangular matrices, i.e., E=E^{U}+E^{L}. So, Eq. (54) can be rewritten as:
Algorithm 1. Predictor/multicorrector multipass algorithm at time t. 
Input: $\mathbf{M}^{*}, \mathbf{D}, \boldsymbol{\varphi}(\hat{t}), \gamma, E, p, l_{\max }$ Output: $\mathbf{H}_{t}^{h}$ $\mathbf{M}^{*} \leftarrow \mathbf{M}+\gamma \mathbf{D}$ Form effective mass matrix $\overline{\mathbf{H}}_{i}^{(0)} \leftarrow \mathbf{0}$ i=0, 1, 2, …, p for $l \leftarrow 1$ to l_{max} do for $i \leftarrow 1$ to p do $\widetilde{\mathbf{H}}_{i}^{(l)} \leftarrow \mathbf{H}_{n}^{}+E_{i i1}^{L} \overline{\mathbf{H}}_{i1}^{(l)}+\sum_{j=i+1}^{p} E_{i j}^{U} \overline{\mathbf{H}}_{j}^{(l1)}$ Predictor phase $\mathbf{M}^{*} \overline{\mathbf{H}}_{i}^{(l)} \leftarrow\mathbf{D} \widetilde{\mathbf{H}}_{i}^{(l)}$ Corrector phase end for end for $\mathbf{H}_{n+1}^{} \leftarrow \mathbf{H}_{n}^{}+\sum_{i=0}^{p} \varphi_{i}(1) \overline{\mathbf{H}}_{i}$ Update solutions 
$\left(I_{p+1} \otimes[\mathbf{M}+\gamma \mathbf{D}]+E^{L} \otimes \mathbf{D}\right) \overline{\mathbf{H}}^{(l)}=d \otimes \mathbf{M H}_{n}^{}E^{U} \otimes \mathbf{D} \overline{\mathbf{H}}^{(l1)}$ (55)
where, l=1, 2, …, l_{max}, and $\overline{\mathbf{H}}^{(l)}$ denotes the solutions obtained from the l^{th} iterate. With a given set of initial values, $\overline{\mathbf{H}}^{(0)}=0$ , Eq. (55) is repeated to determine $\overline{\mathbf{H}}^{\left(l_{\max }\right)}$. The solutions H^{h}(t) can then be calculated via Equation (45). The above iterative scheme is written in the form of the predictor/multicorrector (vform) multipass algorithm is summarised in Algorithm 1.
In order to investigate and ensure the stability and accuracy of the proposed algorithm, a homogeneous firstorder problem with one degree of freedom is taken:
$\dot{H}(t)+\lambda H(t)=0 ; \quad t_{n}^{+} \leq t \leq t_{n+1}^{}$ (56)
$H\left(t_{n}^{}\right)=H_{n}^{}$ (57)
For the current time step, the solutions at the end of the interval is written as:
$H^{h}\left(t_{n+}^{}\right)=\mathcal{A}^{(l)} H^{h}\left(t_{n}^{}\right)$ (58)
where, $\mathcal{A}^{(l)}$ is the amplification matrix at the l^{th} iterate. The algorithm is considered unconditionally stable, i.e. Astable, if $\mathcal{A}<1$ for all positive λΔt_{n}. The graphs of amplification factor obtained by TDGFEM and by secondorder accurate methods as a function of normalized frequency, λΔt_{n} are plotted in Figure 4. It shows that the solutions obtained by TDGP2 are the most accurate, with reference to the exact solutions.
The accuracy of the algorithm is measured by replacing H^{h} in Eq. (58) by exact values H. Then, the local truncation error of solution in each iteration pass l at the end of time step t_{n} is given as:
$\varepsilon^{(l)}=H\left(t_{n+1}^{}\right){ }_{\mathcal{A}}^{(l)} H\left(t_{n}^{}\right)$ (59)
and by using the Taylor’s series at about Δt=0 to obtain the expansion of Eq. (59), the local order accuracy is defined as:
$\tau^{(l)}\left(t_{n}\right)=C_{k} \Delta t^{k}+C_{k+1} \Delta t^{k+1}+\ldots$ (60)
where, C are timeindependent coefficients. The consistency of algorithm is measured as the order of k where k>0 and $\tau(t) \leq C_{k} \Delta t^{k}$ for $t \in\left[0, t_{f}\right]$. The orders of accuracy achieved by the proposed iterative algorithms of linear, when p=1, and quadratic, when p=2, are illustrated in Table 1.
Table 1. Accuracy and stability of iterative algorithms
Algorithm 
l_{max} 
Error Rate 
C_{k} 
Stability 
TDGP1 
1 
Δt 
0.05556 
Astable 
2 
Δt^{3} 
0.01698 
Astable 

3 
Δt^{3} 
0.01389 
Astable 

Direct Solution 
Δt^{3} 
0.01389 
Astable 

TDGP2 
1 
Δt^{2} 
0.02133 
Astable 
2 
Δt^{4} 
0.00064 
Astable 

3 
Δt^{5} 
0.00014 
Astable 

Direct Solution 
Δt^{5} 
0.00014 
Astable 
The iterative algorithms of TDGP1 and TDGP2 methods require three iteration passes (l_{max}=3) to achieve the same orders of accuracy and coefficients as that calculated by direct solving. Therefore, a high order of accuracy obtained by only a few iterations passes means that the algorithm is efficient and stable.
Figure 4. Ampliﬁcation factors of different methods
3.3 Adaptive timestepping method
Usually, moisture diffusion occurs in a very extended period of time in normal environmental conditions. This means that the smaller the time step, the more expensive the calculations. To minimise the simulation time with reliable accuracy, we employ the embedded pair approach. This method has been successfully implemented in the RungeKutta methods [3032], where the resulting parameters determined by the higherorder (q) accurate algorithm are used in the algorithm with a lower order (r) of accuracy.
From Table 1, since our TDGP2 iterative method with three passes, i.e., TDGP2(3) yields the optimal accuracy, the embedded solutions are, hence, of a lowerorder accuracy with two passes, i.e., TDGP2(2) such that r=q1. The difference:
$\mathbf{e}=\mathbf{H}_{n+1}^{}\widehat{\mathbf{H}}_{n+1}^{}$ (61)
is used as a local error estimate, where $\mathbf{H}_{n+1}^{}$ and $\widehat{\mathbf{H}}_{n+1}^{}$ are the solutions by TDGP2(3) and TDGP2(2) respectively. Using the standard rule to predict the new step size [3335] leads to:
$\Delta t_{n+1}=\alpha \Delta t_{n}\left(\frac{\eta_{\text {tol }}}{\eta}\right)^{1 / p}$ (62)
where, α is a safety factor, and p=2 for TDGP2. η_{tol} is the specified tolerance being equal to 0.1% in this work, and η is the relative local error estimate:
$\eta=\frac{\\mathbf{e}\_{E}}{\\mathbf{U}\_{E}} \times 100 \%$ （63）
Using the error measured in total energy norm, the estimated local error becomes:
$\\mathbf{e}\_{E}=\left({ }^{T} \mathbf{e} \mathbf{D} \mathbf{e}\right)^{1 / 2}$ (64)
and, thus, $\\mathbf{U}\_{E}$ is the maximum total energy norm recorded during the simulation:
$\\mathbf{U}\_{E}=\max _{i=0,1, \ldots, n+1}\left({ }^{T} \mathbf{H}_{i}^{} \mathbf{D} \mathbf{H}_{i}^{}\right)^{1 / 2}$ (65)
To make the time stepping method efficient and reliable, the new step size Δt_{n+1} from Eq. (62) is accepted only when η<η_{tol}. Also, to avoid drastic changes in step size, the limit of time step change is imposed, i.e., Δt_{n+1}≤βΔt_{n}, in which β is the control factor, which is taken as 4 in this work. In the case where the error does not satisfy the criterion, i.e. η>η_{tol}, the updated solutions from the current step are rejected, and the determination of $\mathbf{H}_{n+1}^{}$ is attempted with the new smaller step size Δt_{n+1}.
3.4 Material parameters and conditions
The orthotropic elastic constants of wood fibers are determined, using the model proposed by Charupeng and Kunthong [36]. The model requires bulk density and lignocellulosic composition of wood fibers, and completely estimates all the essential elastic constants of wood fibers. Then, the formulation of FE equations in 3D space for the initial configuration (t=0) are constructed, using representative volume element method (RVE), as shown in Figure 5(a). This method has been proven to be reliable with modeling the mechanical properties of wood fiber composites [3740]. With a given value of ϕ and densities of matrix and fibers, the physical parameters a, b, c, T and L can be determined. However, the full RVE model is rather expensive to simulate and unnecessary. We, instead, use a reduced model which is split about the three planes of symmetry in Cartesian system, as depicted in Figure 5(b). This reduces the total number of elements in calculations by a factor of 8.
The geometrical parameters shown in Figure 5(a) can be determined as follows. The volume fraction of fiber is determined from basic materials properties as:
$f=\frac{\frac{w_{f}}{\rho_{f}}}{\frac{w_{f}}{\rho_{f}}+\frac{1w_{f}}{\rho_{m}}}$ (66)
where, w_{f} and ρ_{f} are the weight fraction and bulk density of wood fibers, ρ_{m} is the density of matrix material.
The volume fraction is then used to approximate the initial geometrical parameters a and b of the model  for a constant thickness of fiber wall, a=b=t. The parameter c is determined by averaging the element sizes along the x and y direction. Given the volume of model V, fiber width ϕ and widthtothickness ration; k=ϕ/T, the average element size along the fiber direction c=L/N_{z}; N_{z} is the number of elements along the fiber direction, and $L=\frac{V}{(\phi+2 a)(\phi / k+2 b)}\quad.$ Finally, the corrected value of t is determined via:
$f=\frac{\phi^{2}(L2 c)}{k(\phi+2 t)(\phi / k+2 t)}$ (67)
where, according to literature [41], k=2.
In each complete iteration, the development of moisture, i.e., ΔH calculated by the aforementioned algorithm in Algorithm 1, gives rise to nonzero components in matrix F of Eq. (36). With the predefined global stiffness matrix K, the global displacements U are determined, and hence used to update the physical coordinates of the domain. This process is repeated until the simulation time is reached. The final geometry is postprocessed to obtain meaningful results.
All the linear systems of equations are solved by preconditioned conjugate gradient (PCG) method [42]. This method offers a relatively short period of time to converge and good reliability in solving symmetric positive definite matrices. The convergence of PCG calculation is ensured by setting the tolerance to 1×10^{7} and the maximum iteration number at 10,000. which are then used to update the physical coordinates for the next time step.
In order to illustrate a numerical example as well as validate the algorithm, the experimental results in Ref. [5] are taken for comparison. There are two matrix polymers: polylactic acid (PLA) and polypropylene (PP). The fibers are bleached birch kraft fibers. The essential material parameters are summarised in Table 2. The properties of polymers are assumed to be isotropic. The diffusion coefficients and elastic constants are assumed to be independent upon moisture absorption.
Table 2. Parameter inputs. Referenced studies from which values are taken are included. %C, %H and %L denote the weight fractions of cellulose, hemicellulose and lignin respectively
Constituent 
Parameter 
Value 
Ref. 
Birch 
ρ(kgm^{3}) 
690 
[43] 
%C 
71.7 
[43] 

%H 
26.1 
[43] 

%L 
2.0 
[43] 

ϕ(μm) 
21.3 
[44] 

D_{h,x}(10^{12}m^{2}s^{1}) 
2.26 
[45] 

D_{h,y}(10^{12}m^{2}s^{1}) 
2.26 
[45] 

D_{h,z}(10^{12}m^{2}s^{1}) 
5.65 
[45] 

PLA 
ρ(kgm^{3}) 
1230 
[46] 
E(GPa) 
1.93 
[46] 

ν 
0.35 
[47] 

D(10^{13}m^{2}s^{1}) 
1.77 
[48] 

∂W/∂H 
1×10^{2} 
[49] 

PP 
ρ(kgm^{3}) 
925 
[50] 
E(GPa) 
1.6 
[51] 

ν 
0.45 
[52] 

D(10^{13}m^{2}s^{1}) 
1.90 
[53] 

∂W/∂H 
8×10^{4} 
[54] 
The hygroexpansion coefficient is generally very small in polymeric materials [3, 55] compared to that of wood, and hence assumed to be zero. As the moisture capacity of fiber usually varies with relative humidity, we employ the GuggenheimAndersende Boer (GAB) threeparameter sorption Eq. [56], in which the equilibrium moisture content in cellulose fibers can be determined via:
$W=\frac{V_{m} \times c \times k \times H / 100}{(1k \times H / 100)(1+(c1) k \times H / 100)}$ (68)
where, H is the relative humidity, and for adsorption parameters of birch wood, V_{m}=0.058, c=8.9, k=0.78 [57]. So, ∂W/∂H, can be determined by differentiating Eq. (68).
Also, by using the FE algorithm [36] with inputs from Table 2 ρ, %C, %H and %L, the elastic constants of wood fibers are estimated and summarised in Table 3. The above properties of the homogenised fiber model are used to simulate the birch pulp in the experiments [5].
Table 3. Elastic constants of homogenised birch fibers
$E_{x}$ (GPa) 
$E_{y}$ (GPa) 
$E_{z}$ (GPa) 
$G_{y z}$ (GPa) 
$G_{zx}$ (GPa) 
$G_{xy}$ (GPa) 
$\nu_{z y}$ 
$\nu_{z x}$ 
$\nu_{x y}$ 
2.62 
2.62 
42.8 
4.23 
4.23 
1.19 
0.447 
0.447 
0.104 
Figure 5. RVE model of WPC; fiber(brown), matrix(yellow). (a) full model with annotated multiview projection. (b) actual physical domain used for simulation. ϕ, T and L are the width, thickness and length of the homogenised fiber respectively
With reference to Figure 5(b), the boundary conditions imposed at all time steps are:
$\mathbf{H}(\infty, y, z)=\mathbf{H}(x, \infty, z)=\mathbf{H}(x, y, \infty)=R H$ (69)
$\mathbf{u}(0, y, z)=\mathbf{v}(x, 0, z)=\mathbf{w}(x, y, 0)=0$ (70)
where, RH is the relative humidity, u, v and w are nodal displacements along the x, y and z directions, and ∞ indicates the boundary on a specific dimension. In other words, Eq. (69) are posed on to the nodes with coordinates $\left(\frac{\phi}{2}+a, y, z\right),\left(x, \frac{T}{2}+b, z\right)$, and $\left(x, y, \frac{L}{2}+c\right)$ where ϕ, T, L, a, b, c are parameters according to Figure 5(a). And from Eq. (70), u=0 are posed on to the nodes with coordinates (0, y, z), v=0 are posed on to the nodes with coordinates (x, 0, z), and w=0 are posed on to the nodes with coordinates (x, y, 0).
The initial conditions from Eq. (20) are used on all nodes except that with the imposed essential conditions above.
As for the FE discretisation, the dimensions of the RVE model are specified such that the total volume of the model equals the total volume of the specimens in the experiments. The uniformity of mesh is ensured by making all element as close to a square as possible. Using the concepts of RVE, the average outofplane hygroexpansion, is plotted against time for convergence tests as shown in Figure 6.
Figure 6. Effects of mesh geometry on the solutions obtained for 30% birch and 70% PP by weight, at RH=0.97. 3×2×10 indicates that the FE model has three elements, two elements and 10 elements along the x, y and z respectively
The choices of crosssectional dimensions xyplanes of elements are rather limited as to conserve the volume fraction of fiber and the fiber widthtothickness ratio W/T. Refining the mesh from 2×2×10 to 3×2×10 and increasing the number of elements along the fiber direction yield converged solutions to less than 1% as the dimensions of each element become more uniform. However, the 4×3×13 mesh has the greatest number of elements but yields a rather different trend of solutions. This is because, in order to conserve the volume fraction of fibers, such mesh leads to uneven dimensions of elements on the crosssections, and hence poor aspect ratios. Therefore, the minimum refinement of mesh is 3 by 2 by 12 elements.
The preferred mesh refinement is taken for the evaluation of the efficiency of timestepping algorithm. The initial time step is set to 100 seconds with β=4 to avoid large increases of step size, α=0.9 and η_{tol}=0.1%. The size of time step, Δt_{n}, are plotted against iteration number, n, as shown in Figure 7.
Figure 7. Time step progression vs step number. Labels of approximate day number in the simulated period are included for certain iterations. Runtime is approximately 6.5 hours
The results show that the step size rises rather rapidly for the first few steps to a certain range  around 10^{4} seconds or three hours. It then increases from hours to days and weeks, as implied by the labelled day numbers. Such increasing intervals are also used commonly by numerous similar experimental studies [57]. This suggests that the algorithm quickly determines the appropriate step size from the solutions and progresses according to the hygromechanical behaviours. The fact that our algorithm requires 6.5 hours of runtime to simulate the actual period of 90 days indicates good robustness and its potential for industrial applications  predictions of dimensional stability and hence durability of materials are essential but expensive for laboratory experiments.
Figure 8 shows the predicted timedependent outofplane hygroexpansion as a function of time, for a series of different fiber volume fractions and material constituents. All FE calculations are performed with RH=0.97, Δt_{1}=100s, α=0.9, η_{tol}=0.1% and β=4. The results are compared to corresponding experimental results.
Figure 8. Timedependent outofplane hygroexpansion predicted by our FEA model (solid lines) compared to the experiments [5] (dotted lines)
There is good agreement between our FEA calculations and the results of the experiments [5]. The model accurately predicts accurately the hygroexpansion towards equilibrium moisture content with slight differences of up to 10% in the first 1020 days. Within the aforementioned period, the reliability of the prediction is slightly lower, which is due to the fact that the magnitude of the diffusion coefficient in polymers can vary from 10^{15} to 10^{11} m^{2}/s [58, 59]. Moreover, water absorption, although very small, can lead to rearrangement of the chemical chains and thus to a change in diffusivity [2].
While the numerical model assumes perfect interfaces between the fibers and the matrix, this is usually not the case in practice. The interfacial compatibility of wood fibers is poor in a PP matrix but strong in a PLA matrix, as evidenced by SEM micrographs according to the literature [5, 60]. Thus, as hygroexpansion progresses, microcracks and interfacial debonding develop, providing more channels for a greater extent of moisture diffusion. This has also been shown by Suri and Perreux [61] that damage and cracks lead to increased water diffusion rate in fiber composites. As a result, FEA calculations give higher values compared to the experiments in the case of 30 wt% fibers and 70 wt% PP after 20 days.
It should also be noted that the simulation environments are assumed such that all the sides (boundary surfaces) of composite specimens are exposed to uniform relative humidity and that the specimen are kept at constant and standard temperature and pressure during the diffusion process. In terms of actual experimentation, the consistency with of such conditions is ensured by complying to the relevant standard ASTM D 57081 [5].
To further evaluate the reliability and accuracy of the modeling methods, the kinetics of diffusion is analysed. The analytical equation describing the diffusion mechanism and kinetics can be given as follows:
$\log \left(\frac{M_{t}}{M_{\infty}}\right)=\log (k)+n \log (t)$ (71)
where, M_{t} and M_{∞} denote moisture content at time t and that at the equilibrium respectively, k and n are constants [62]. For relative moistures of lower than 0.5 (M_{t}/M_{∞}≤0.5), the n value indicates the type of diffusion mechanism; n=0.5 for general Fickian diffusion processes [63]. Here, the value of M_{t}/M_{∞} is assumed to be equal to H_{t}/RH. The value of H_{t} is given by the average value of the nodal solutions from Eq. (37) at time t. Since the equilibrium value cannot exceed the relative humidity RH (see Eq. (19)), this value is set as the equilibrium value.
Figure 9. Log of relative moisture against log time in seconds including trend lines
Figure 9 shows log(H_{t}/RH) as a function of log(t) for the three composites. For all composites, the slopes of the trend line, i.e. n, are approximately equal to 0.5  with a deviation of at most 6%. This is in agreement with the experimental study [6], according to which numerous natural fibers/ PP composites have n value of about 0.5, thus obeying the kinetics theory of a Fickian diffusion process. There is a slight difference between the log(k) values from the trend line equations of birch fiber/PLA composites because the RVE models are dimensioned differently due to the difference in fiber weight fractions. Moreover, the n value for birch fiber/PP composites is quite different due to the low interfacial compatibility mentioned earlier. The above analysis confirms that our model can effectively capture the physics of the hygromechanical behaviours of WPC as well as the kinetics according to Fick’s theory.
We have successfully developed a validated finite element model based on semidiscretization of space and time to study the threedimensional hygromechanical behaviours of wood fiberpolymer composites. The solutions are continuous in space but can be discontinuous in time, using TDGFEM for the latter to achieve optimal accuracy.
A predictor/multicorrector using a GaussSeidel block matrix iterative scheme algorithm for solving the firstorder parabolic problem in temporal domain has also been presented. It requires only up to three iterative passes to achieve the same accuracy as the direct solution. Moreover, the adaptive timestepping method using embedded solution approach with the maximum energy norm criterion has been shown to effectively determine the appropriate size of time step and to progress efficiently throughout the simulation period. A runtime of a few hours compared to about three months of actual laboratory experiments confirms the efficiency and robustness in predicting the complex timedependent hygromechanical behaviours.
As for the validation of the model, the FEA results are compared with the corresponding cases in previous experimental studies and show good agreement. The reliability of the prediction can increase with the interfacial compatibility between fiber and matrix. Diffusion at early stages may be more difficult to predict due to the variability of material parameters in polymeric materials. Data analysis of the timedependent moisture content confirms that Fick’s theory is very good at predicting the behaviours, with the characteristic n value in agreement with experiments.
[1] Bledzki, A.K., Sperber, V.E., Faruk, O. (2002). Natural and wood fibre reinforcement in polymers. iSmithers Rapra Publishing.
[2] Weisenberger, L.A., Koenig, J.L. (1989). NMR imaging of solvent diffusion in polymers. Applied Spectroscopy, 43(7): 11171126.
[3] Yew, G.H., Yusof, A.M., Ishak, Z.M., Ishiaku, U.S. (2005). Water absorption and enzymatic degradation of poly (lactic acid)/rice starch composites. Polymer Degradation and Stability, 90(3): 488500. https://doi.org/10.1016/j.polymdegradstab.2005.04.006
[4] Klyosov, A.A. (2007). WoodPlastic Composites. John Wiley & Sons, 1698. https://doi.org/10.1002/9780470165935
[5] Almgren, K.M., Gamstedt, E.K., Berthold, F., Lindström, M. (2009). Moisture uptake and hygroexpansion of wood fiber composite materials with polylactide and polypropylene matrix materials. Polymer Composites, 30(12): 18091816. https://doi.org/10.1002/pc.20753
[6] Espert, A., Vilaplana, F., Karlsson, S. (2004). Comparison of water absorption in natural cellulosic fibres from wood and oneyear crops in polypropylene composites and its influence on their mechanical properties. Composites Part A: Applied Science and Manufacturing, 35(11): 12671276. https://doi.org/10.1016/j.compositesa.2004.04.004
[7] Leu, Y.Y., Chow, W.S. (2011). Kinetics of water absorption and thermal properties of poly (lactic acid)/organomontmorillonite/poly (ethylene glycol) nanocomposites. Journal of Vinyl and Additive Technology, 17(1): 4047. https://doi.org/10.1002/vnl.20259
[8] Duanmu, J., Gamstedt, E.K., Rosling, A. (2007). Hygromechanical properties of composites of crosslinked allylglycidylether modified starch reinforced by wood fibres. Composites Science and Technology, 67(1516): 30903097. https://doi.org/10.1016/j.compscitech.2007.04.027
[9] Šejnoha, M., Sýkora, J., Vorel, J., Kucíková, L., Antoš, J., Pokorný, J., Pavlík, Z. (2019). Moisture induced strains in spruce from homogenization and transient moisture transport analysis. Computers & Structures, 220: 114130. https://doi.org/10.1016/j.compstruc.2019.04.008
[10] Autengruber, M., Lukacevic, M., Füssl, J. (2020). Finiteelementbased moisture transport model for wood including free water above the fiber saturation point. International Journal of Heat and Mass Transfer, 161: 120228. https://doi.org/10.1016/j.ijheatmasstransfer.2020.120228
[11] Wang, N., Liu, W., Lai, J. (2014). An attempt to model the influence of gradual transition between cell wall layers on cell wall hygroelastic properties. Journal of Materials Science, 49(5): 19841993. https://doi.org/10.1007/s1085301378855
[12] Joffre, T., Isaksson, P., Dumont, P.J., Roscoat, S., Sticko, S., Orgéas, L., Gamstedt, E.K. (2016). A method to measure moisture induced swelling properties of a single wood cell. Experimental Mechanics, 56(5): 723733. https://doi.org/10.1007/s1134001501199
[13] Obeid, H., Clément, A., Fréour, S., Jacquemin, F., Casari, P. (2018). On the identification of the coefficient of moisture expansion of polyamide6: Accounting differential swelling strains and plasticization. Mechanics of Materials, 118: 110. https://doi.org/10.1016/j.mechmat.2017.12.002
[14] Sar, B.E., Fréour, S., Célino, A., Jacquemin, F. (2015). Accounting for differential swelling in the multiphysics modeling of the diffusive behavior of a tubular polymer structure. Journal of Composite Materials, 49(19): 23752387. https://doi.org/10.1177%2F0021998314546704
[15] Sar, B.E., Fréour, S., Davies, P., Jacquemin, F. (2014). Accounting for differential swelling in the multiphysics modelling of the diffusive behaviour of polymers. ZAMMJournal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 94(6): 452460. https://doi.org/10.1002/zamm.201200272
[16] Mbacké, M.A., Khazaie, S., Fréour, S., Jacquemin, F. (2020). Numerical modeling of the hygromechanical behavior of the 3D composite materials. Journal of Composite Materials, 54(28): 44574471. https://doi.org/10.1177%2F0021998320932270
[17] Youssef, G., Fréour, S., Jacquemin, F. (2009). Stressdependent moisture diffusion in composite materials. Journal of Composite Materials, 43(15): 16211637. https://doi.org/10.1177%2F0021998309339222
[18] Stahl, D.C., Cramer, S.M. (1998). A threedimensional network model for a low density fibrous composite. J. Eng. Mater. Technol., 120(2): 126130. https://doi.org/10.1115/1.2807000
[19] Termonia, Y. (1994). Structureproperty relationships in shortfiberreinforced composites. Journal of Polymer Science Part B: Polymer Physics, 32(6): 969979. https://doi.org/10.1002/polb.1994.090320601
[20] Hughes, T.J. (2012). The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publication, Inc.
[21] Hughes, T.J., Hulbert, G.M. (1988). Spacetime finite element methods for elastodynamics: Formulations and error estimates. Computer Methods in Applied Mechanics and Engineering, 66(3): 339363. https://doi.org/10.1016/00457825(88)900060
[22] Chien, C.C., Wu, T.Y. (2000). An improved predictor/multicorrector algorithm for a timediscontinuous Galerkin finite element method in structural dynamics. Computational Mechanics, 25(5): 430437. https://doi.org/10.1007/s004660050490
[23] Zhang, R., Wen, L., Xiao, J., Qian, D. (2019). An efficient solution algorithm for space–time finite element method. Computational Mechanics, 63(3): 455470. https://doi.org/10.1007/s0046601816038
[24] Kunthong, P., Thompson, L.L. (2005). An efficient solver for the highorder accurate timediscontinuous Galerkin (TDG) method for secondorder hyperbolic systems. Finite Elements in Analysis and Design, 41(78): 729762. https://doi.org/10.1016/j.finel.2004.09.003
[25] Sharma, V., Fujisawa, K., Murakami, A. (2021). Spacetime finite element method for transient and unconfined seepage flow analysis. Finite Elements in Analysis and Design, 197: 103632. https://doi.org/10.1016/j.finel.2021.103632
[26] Luo, H., Absillis, G., Nourgaliev, R. (2021). A moving discontinuous Galerkin finite element method with interface condition enforcement for compressible flows. Journal of Computational Physics, 445: 110618. https://doi.org/10.1016/j.jcp.2021.110618
[27] Baccouch, M., Temimi, H. (2021). A highorder space–time ultraweak discontinuous Galerkin method for the secondorder wave equation in one space dimension. Journal of Computational and Applied Mathematics, 389: 113331. https://doi.org/10.1016/j.cam.2020.113331
[28] Dong, Z., Li, H. (2021). A spacetime finite element method based on local projection stabilization in space and discontinuous Galerkin method in time for convectiondiffusionreaction equations. Applied Mathematics and Computation, 397: 125937. https://doi.org/10.1016/j.amc.2020.125937
[29] MATLAB, V9.10.0.1739362 (R2021a). 2021, Massachusetts: The MathWorks, Inc.
[30] Dormand, J.R., ElMikkawy, M.E.A., Prince, P.J. (1987). Families of RungeKuttaNystrom formulae. IMA Journal of Numerical Analysis, 7(2): 235250. https://doi.org/10.1093/imanum/7.2.235
[31] Higham, D.J. (1997). Regular RungeKutta pairs. Applied Numerical Mathematics, 25(23): 229241. https://doi.org/10.1016/S01689274(97)000627
[32] de Swart, J.J., Söderlind, G. (1997). On the construction of error estimators for implicit RungeKutta methods. Journal of Computational and Applied Mathematics, 86(2): 347358. https://doi.org/10.1016/S03770427(97)001660
[33] Gear, C.W. (1966). The numerical integration of ordinary differential equations of various orders (No. ANL7126). Argonne National Lab., Ill. https://doi.org/10.1090/S00255718196702254945
[34] Tsitouras, C., Papakostas, S.N. (1999). Cheap error estimation for RungeKutta methods. SIAM Journal on Scientific Computing, 20(6): 20672088. https://doi.org/10.1137/S1064827596302230
[35] Hull, T.E., Enright, W.H., Fellen, B.M., Sedgwick, A.E. (1972). Comparing numerical methods for ordinary differential equations. SIAM Journal on Numerical Analysis, 9(4): 603637. https://doi.org/10.1137/0709052
[36] Charupeng, N., Kunthong, P. (2021). A novel finite element algorithm for predicting the elastic properties of wood fibers. International Journal of Computational Materials Science and Engineering, 2150027. https://doi.org/10.1142/S2047684121500275
[37] Shokrieh, M.M., MoshrefzadehSani, H. (2016). An optimized representative volume element to predict the stiffness of aligned short fiber composites. Journal of Composite Materials, 50(23): 33013310. https://doi.org/10.1177%2F0021998315618456
[38] Ogierman, W., Kokot, G. (2018). Generation of the representative volume elements of composite materials with misaligned inclusions. Composite Structures, 201: 636646. https://doi.org/10.1016/j.compstruct.2018.06.086
[39] Naderi, M., Apetre, N., Iyyer, N. (2017). Effect of interface properties on transverse tensile response of fiberreinforced composites: threedimensional micromechanical modeling. Journal of Composite Materials, 51(21): 29632977. https://doi.org/10.1177%2F0021998316681189
[40] Schneider, M. (2017). The sequential addition and migration method to generate representative volume elements for the homogenization of short fiber reinforced plastics. Computational Mechanics, 59(2): 247263. https://doi.org/10.1007/s0046601613507
[41] Mansfield, S.D., Kibblewhite, R.P., Riddell, M.J. (2004). Characterization of the reinforcement potential of different softwood kraft fibers in softwood/hardwood pulp mixtures. Wood and Fiber Science, 36(3): 344358.
[42] Barrett, R., Berry, M., Chan, T.F., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo, R., Romine, C., van der Vorst H. (1996). Templates for the solution of linear systems: Building blocks for iterative methods. Math. Comput., 64(211): 13491352. https://doi.org/10.2307/2153507
[43] Asikainen, S.A., Fuhrmann, A., Ranua, M., Robertsén, L. (2010). Effect of birch kraft pulp primary fines on bleaching and sheet properties. BioResources, 5(4): 21732183.
[44] Danielewicz, D., SurmaŚlusarska, B. (2017). Properties and fibre characterisation of bleached hemp, birch and pine pulps: A comparison. Cellulose, 24(11): 51735186. https://doi.org/10.1007/s1057001714766
[45] Šejnoha, M., KUCíKOVÁ, L., Vorel, J., Sýkora, J., De Wilde, W.P. (2019). Effective material properties of wood based on homogenization. International Journal of Computational Methods and Experimental Measurements, 7(2): 167180. http://doi.org/10.2495/CMEMV7N2167180
[46] Van de Velde, K., Kiekens, P. (2002). Biopolymers: overview of several properties and consequences on their applications. Polymer Testing, 21(4): 433442. https://doi.org/10.1016/S01429418(01)001076
[47] Fung, Y.C. (1965). Foundations of Solid Mechanics. New Jersey: PrenticeHall, Inc.
[48] Yang, C., Xing, X., Li, Z., Zhang, S. (2020). A comprehensive review on water diffusion in polymers focusing on the polymer–metal interface combination. Polymers, 12(1): 138. https://doi.org/10.3390/polym12010138
[49] Cairncross, R.A., Becker, J.G., Ramaswamy, S., O’Connor, R. (2006). Moisture sorption, transport, and hydrolytic degradation in polylactide. In TwentySeventh Symposium on Biotechnology for Fuels and Chemicals, pp. 774785. https://doi.org/10.1385/ABAB:131:1:774
[50] Maddah, H. A. (2016). Polypropylene as a promising plastic: A review. American Journal of Polymer Science, 6(1): 111. https://doi.org/10.5923/j.ajps.20160601.01
[51] Oksman, K., Skrifvars, M., Selin, J.F. (2003). Natural fibres as reinforcement in polylactic acid (PLA) composites. Composites science and technology, 63(9): 13171324. https://doi.org/10.1016/S02663538(03)001039
[52] Kumar, M., Gaur, K.K., Shakher, C. (2015). Measurement of material constants (Young's modulus and Poisson's ratio) of polypropylene using digital speckle pattern interferometry (DSPI). Journal of the Japanese Society for Experimental Mechanics, 15: s87s91. https://doi.org/10.11395/jjsem.15.s87
[53] Yi, X., Pellegrino, J. (2002). Diffusion measurements with Fourier transform infrared attenuated total reflectance spectroscopy: Water diffusion in polypropylene. Journal of Polymer Science Part B: Polymer Physics, 40(10): 980991. https://doi.org/10.1002/polb.10161
[54] Suherman, S., Peglow, M., Tsotsas, E. (2010). Measurement and modelling of sorption equilibrium curve of water on PA6, PP, HDPE and PVC by using FloryHuggins model. Reaktor, 13(2): 8994. https://doi.org/10.14710/reaktor.13.2.8994
[55] Weitsman, Y.J. (2011). Fluid Effects in Polymers and Polymeric Composites. Springer Science & Business Media. http://doi.org/10.1007/9781461410591
[56] Timmermann, E.O. (2003). Multilayer sorption parameters: BET or GAB values? Colloids and Surfaces A: Physicochemical and Engineering Aspects, 220(13): 235260. https://doi.org/10.1016/S09277757(03)000591
[57] Bratasz, Ł., Kozłowska, A., Kozłowski, R. (2012). Analysis of water adsorption by wood using the GuggenheimAndersonde Boer equation. European Journal of Wood and Wood Products, 70(4): 445451. https://doi.org/10.1007/s001070110571x
[58] Duval, S., Camberlin, Y., Glotin, M., Keddam, M., Ropital, F., Takenouti, H. (2000). Characterisation of organic coatings in sour media and influence of polymer structure on corrosion performance. Progress in Organic Coatings, 39(1): 1522. https://doi.org/10.1016/S03009440(00)000941
[59] De Rosa, L., Monetta, T., Bellucci, F. (1998). Moisture uptake in organic coatings monitored with EIS. Materials Science Forum, 289: 315326. https://doi.org/10.4028/www.scientific.net/MSF.289292.315
[60] Huda, M.S., Drzal, L.T., Misra, M., Mohanty, A.K. (2006). Woodfiberreinforced poly (lactic acid) composites: Evaluation of the physicomechanical and morphological properties. Journal of Applied Polymer Science, 102(5): 48564869. https://doi.org/10.1002/app.24829
[61] Suri, C., Perreux, D. (1995). The effects of mechanical damage in a glass fibre/epoxy composite on the absorption rate. Composites Engineering, 5(4): 415424. https://doi.org/10.1016/09619526(94)00014Z
[62] Comyn, J. (Ed.). (1985). Polymer Permeability. Springer Science & Business Media. http://doi.org/10.1007/9789400948587
[63] Crank, J. (1956). The Mathematics of Diffusion. Oxford, UK: Cambridge University Press.