A Novel Space-Time Finite Element Algorithm to Investigate the Hygro-Mechanical Behaviours of Wood Fiber-Polymer 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, moisture-induced expansion, i.e., hygroexpansion, is one of the major drawbacks of wood-based composites. Unlike wood fibers, matrix components in wood-based composites, such as polymers, are much less susceptible to water absorption than wood [1-3]. The drastic difference in hygroelastic properties can even lead to serious problems in the case of wood fiber-polymer composites (WPC), such as buckling and warping of the installed WPC panels [4]. Therefore, studying the hygro-mechanical behaviours of WPC is an essential part of improving the durability of the materials.

It has been shown experimentally [5-8] 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, hygro-mechanical modeling is becoming increasingly important. Numerous hygro-mechanical modeling studies are conducted for wood-based materials [9-12] and for polymeric materials [13-15], but there is very little research on predicting the time-dependent hygro-mechanical behaviours of WPC. There is particularly little research that attempts to accurately predict the three-dimensional anisotropic responses of WPC subjected to water absorption.

Mbacké et al. [16] have used a numerical method to study the hygro-mechanical behaviours of homogenised composites, determining the time-dependent 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 multi-scale analytical model for stress-dependent moisture diffusion in fiber-reinforced 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]. Time-dependent 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 semi-discretisation. Like heat transfer problems, moisture diffusion in materials can be considered as a first-order parabolic problem, and perhaps the most commonly used methods for time integration are classified as single-step/single-solver methods (SS/SS) [20]. Such methods are based on solving a system of $n_{e q}$ linear equations leading to second-order 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 space-time finite element method, which was developed by Hughes and Hulbert [21]. The method is based on the time-discontinuous Galerkin finite element method (TDG-FEM), 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 space-time linear systems of equations and the relatively high computational cost of solving non-symmetric TDG matrices directly, several techniques [22-24] have been developed. As a result, there are now several applications of TDG-FEM [25-28] However, to our knowledge, there are no studies using TDG-FEM to investigate the hygro-mechanical behaviours of WPC.

The current study has two objectives: 1) to develop a complete and efficient TDG-FEM algorithm written in MATLAB for hygro-mechanical 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 TDG-FEM matrix equations in space and time for the hygro-mechanical behaviours of WPC is presented. Then, the predictor/multi-corrector is described using a multi-pass algorithm from [24] and applied to the first-order parabolic moisture diffusion problem. A suitable adaptive time-stepping 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 time-stepping scheme using the embedded solution approach is discussed. It is followed by model validation and quantitative analysis of model reliability.

There is no or very little commercial software capable of modeling such a time-dependent moisture diffusion problem using high-accuracy adaptive time integration methods. Therefore, to simulate the hygro-mechanical 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.

1.jpg

Figure 1. Overall FE-based algorithm. ← and → denote input and output respectively

3.1 3D FE formulation in spatial domain

For each element in space, i.e. Ω^e, an eight-node 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 semi-discretisation 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)

2.png

Figure 2. An eight-node 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, moisture-induced 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^{e-1} \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 weak-form 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, body-force is neglected (f_i=0), and there are no mechanical tractions $\left(\bar{t}_{i}=0\right)$ on the boundaries. The weak-form 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 non-trivial solutions arise from the moisture-induced strains from Eq. (5) and by substituting Eq. (5) into Eq. (23), the general 3D FE-discretised 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}\left|J^{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}\left|J^{e}\right| d \xi d \eta d \zeta$.    (27)

|J^e| is the determinant of J^e. $B_{m}^{e}$ is the strain-displacement 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 first-order 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}\left|J^{e}\right| d \xi d \eta d \zeta$    (31)

$D^{e}=\int_{\Omega^{e}}{ }^{T} B_{h}^{e} D_{h}^{e} B_{h}^{e}\left|J^{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 hygro-mechanical 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 time-discontinuous Galerkin finite element method (TDG-FEM). This method allows unconditionally stable, highly accurate time integration via using high-order p polynomial approximation with the accuracy order of 2p+1 at the end of each time step.

The formulation of TDG-FEM 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}^{N-1}\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.

3.png

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, …, N-1. The principle of TDG-FEM [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₀)=H₀. The L₂-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 high-order accuracy, the time shape functions are obtained thereby mapping of monomial basis {1, t, t², …, 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}^{i-k}$     (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 first-order 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^-1B 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 off-diagonal 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 predictor-corrector algorithm, the left-hand side of Eq. (51) is re-written 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 multi-pass algorithm, a Gauss-Seidel 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 re-written as:

$\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}}^{(l-1)}$    (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/multi-corrector (v-form) multi-pass algorithm is summarised in Algorithm 1.

In order to investigate and ensure the stability and accuracy of the proposed algorithm, a homogeneous first-order 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. A-stable, if $|\mathcal{A}|<1$ for all positive λΔt_n. The graphs of amplification factor obtained by TDG-FEM and by second-order accurate methods as a function of normalized frequency, λΔt_n are plotted in Figure 4. It shows that the solutions obtained by TDG-P2 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 time-independent 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

The iterative algorithms of TDG-P1 and TDG-P2 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.

4.jpg

Figure 4. Amplification factors of different methods

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 Guggenheim-Andersen-de Boer (GAB) three-parameter 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}{(1-k \times H / 100)(1+(c-1) 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

5.png

Figure 5. RVE model of WPC; fiber(brown), matrix(yellow). (a) full model with annotated multi-view 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.