Effect of Thermo-Hydro-Mechanical Treatment on Mechanical Properties of Wood Cellulose: A Molecular Dynamics Simulation

Abstract

Based on molecular dynamics, a water and cellulose model was constructed to provide more theoretical support for the behavior characteristics of cellulose properties in thermo-hydro-mechanical treatment. In this paper, dynamic simulations were carried out under the NPT ensemble at 4, 5.5, 8, and 12 MPa, respectively. Moreover, we analyze the effects on the mechanical properties of wood cellulose in terms of the hydrogen bond numbers, small molecule diffusion coefficients, end-to-end distances, and mechanical parameters of the water–cellulose model. The results indicate that the densification of the water–cellulose model gradually increases with increasing pressure. The effect of pressures on mechanical properties is mainly due to the formation of massive hydrogen bonds within the cellulose chain and between water and cellulose. This is reflected in the fact that water molecules are more difficult to diffuse in the cellulose, which therefore weakens the negative effect of large amounts of water on the cellulose. The increase in end-to-end distance represents the stiffness of the cellulose chains being strengthened. The mechanical parameters indicate an increase in wood stiffness to resist deformation better, while reducing tensile properties at the same time. The dynamic simulation results in this paper can well correspond to macroscopic experiments.

1. Introduction

As the only renewable green material, wood is used extensively because of its wide coverage of raw materials and strong maneuverability. The hygroscopic expansion of wood materials can form bonding interfaces, but at the same time, they are also vulnerable to corrosion cracking and dimensional deformation, which reduces the utilization rate of wood. Therefore, natural wood needs to be modified [1]. Heat treatment is the most common method used in wood modification [2]. Since the concept of heat treatment was first proposed by Stamm [3] in 1946, this method has been widely studied and continuously modified to better improve wood utilization, such as dimensional stability, moisture absorption, durability, etc. [4].

Hydrothermal treatment uses water vapor or liquid water as the heat transfer medium and modifies wood in the range of 160–240 °C [5]. In general, as the temperature of hydrothermal treatment increases, the dimensional stability and durability of wood increases significantly, but it will cause the loss of mechanical properties of wood [6]. This was confirmed by the study of Camiyani Black Pine, where the compression strength of Black Pine was reduced by 27% after 10 h treatment at 180 °C [7], while in the M. Gaff et al. study, the highest reduction in the impact bending strength of oak and spruce reached 32.2% and 39.8%, respectively [8].

To reduce the influence of hydrothermal treatment on the mechanical properties of wood, thermo-hydro-mechanical (THM) treatment can be used to modify wood [9]. Thermo-hydro-mechanical wood treatment is based on a combination of high temperatures, moisture, and mechanical loading, which leads to modified wood [10]. It has been demonstrated that THM treatment can better improve the mechanical strength of wood compared to dry-heating conditions [11]. Bao et al. applied 4–12 MPa of pressure while THM densified hybrid poplar wood at 160 °C [12]. The modulus of elasticity increased significantly compared to the wood without high pressure. Moreover, this study showed a positive correlation with the mechanical properties of the wood as the densification degree of wood enhanced.

By calculating the dynamic evolution of molecules, molecular dynamics can help us understand their motion. Cellulose is the primary component of wood. The study of cellulose movements by molecular dynamics assists us in both verifying the changes of wood in macroscopic experiments and predicting the impact of changes in conditions on the properties of wood. Wood cellulose belongs to higher plant cellulose, whose crystalline form is monoclinic crystal [13]. Cellulose is composed of crystalline and amorphous regions. The cellulose chains are disordered in the amorphous region. It is obviously observed that there are many pore structures in the system of cellulose [14], which causes a high frequency of interaction between water and cellulose in the amorphous region. Du et al. investigated the mechanical properties of amorphous cellulose by constructing an amorphous region from a cellulose chain with a degree of polymerization of 20 [15]. The mechanical properties obtained were similar to the experimental values. Therefore, this study focuses on the amorphous region of cellulose.

The accuracy of molecular dynamics simulation depends on the selection of force fields. Zhang et al. used the GLYCAM force field to study the effect of temperature on the thermal response of crystalline cellulose [16]. Chen et al. investigated the structural stability and thermal motion of natural cellulose in different potential hydrogen bonding patterns using the Gromos56Acarbo force field [17]. Wang et al. used the PCFF force field to simulate the cellulose at high temperatures, which investigated the effects of different temperatures on the mechanical properties of cellulose [18]. Meanwhile, it is also proved that the PCFF force field is suitable for the amorphous region of cellulose.

Available studies on molecular dynamics simulation mainly focused on the effects of temperature and transfer medium on wood heat treatment. The research on the effect of THM treatment on the properties of wood is limited to macroscopic properties such as wood size, hardness, and strength. However, the mechanism between heat and pressure and solid wood with high moisture has not been well explained. Exploring the effects of the combined mechanisms of pressure, temperature, and moisture on wood is crucial to improving wood utilization.

In this paper, a cellulose model was conducted under certain humidity in the con-trolled environment of THM treatment, which is used to analyze the micro-mechanism associated with the effect of THM treatment on the mechanical properties of the water–cellulose model that was studied. Moreover, we analyzed this model in terms of its hydrogen bond numbers, small molecule diffusion coefficients, and end-to-end distances. The main purpose of this paper is to study the changes in the macroscopic properties of wood in the THM treatment while providing stronger theoretical support for wood thermo-hydro-mechanical modification.

2. Materials and Methods

2.1. Model Establishment

Material Studios software (version 2020, Accelrys, San Diego, CA, USA) was used to simulate a cellulose chain with a degree of polymerization of 20. In order to investigate the effect of humidity on cellulose, a composite model with water contents of 25% was constructed [12]. The model was constructed using the method proposed by Theodorou et al. for the construction of amorphous polymers [19]. The Ewald method was adopted to charge addition. The Van der Waals addition method used the atom-based method. The water–cellulose composite model was established with a target density of 1.5 g/cm³ [20], which consists of a cellulose chain with 125 water molecules, as shown in Figure 1. The water molecules are displayed in ball-and-stick mode and the cellulose molecule is displayed in stick mode. The hydrogen and oxygen atoms are represented with the white and yellow spheres, respectively.

2.2. Dynamic Simulation

All simulations in this paper were based on the PCFF force field [18]. To minimize the energy of the system before the molecular dynamics simulation, the model must first completely relax and minimize its energy. First, the Smart algorithm was used to perform 5000 steps of operations, after which the model was geometrically optimized. Then we performed a dynamic relaxation at 433.15 K under the NVT ensemble for 1000 ps. One frame was output every 5000 steps. Huang et al. has proved that these parameters are reliable [21]. After all these steps, the system energy tended to be stable.

Based on the previous complete relaxation, we performed a 1000 ps under the NVT ensemble with a temperature of 433.15 K as pre-heating. After that, a 1000 ps dynamic simulation was performed under the NPT ensemble. The temperature was set to 433.15 K and the pressure was set to 4, 5.5, 8, and 12 MPa for simulation, respectively. The model with dynamic simulations at atmospheric pressure of 0.1 MPa as the untreated control group was prepared as well. Bao et al. tested the selection of temperature and pressure, demonstrating that the effect of the pressurized hydrothermal treatment on cellulose could be well explored [12]. The electronic effect was controlled by the Ewald method [22], the van der Waal force was calculated by the atom-based method [23], the temperature was controlled by the Anderson method [24], and the pressure was controlled by the Berendsen method [25].

3. Results and Discussion

3.1. Model Characterization

3.1.1. System Energy

In molecular dynamics, temperature and energy will change with time, which can be used as criteria for determining whether a system is in equilibrium. The variation of the water–cellulose system’s energy with time is shown in Figure 2a. The system equilibrium can be expressed by the energy convergence parameter δE [20], which is defined in Equation (1):

$$\delta E=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left|\frac{E_i-E_0}{E_0}\right|}$$

(1)

where N is the total number of steps of the simulation, $E_0$ represents the system’s initial energy, and $E_i$ represents the energy after the simulation i steps.

When $\delta E\le 0.001-0.003$, the system tends to be in equilibrium, while the simulation results are reliable. After analysis, the $\delta E$ for the water–cellulose model was 0.0015 at the last 200 ps.

The temperature–time variation of the water–cellulose composite model is shown in Figure 2b. The temperature, as shown in the figure, varied slightly over time while remaining within ±25 K, which indicated that the system was in equilibrium after energy relaxation.

As a result of geometric optimization and energy relaxation, the energy of the water–cellulose composite model is stable, and the subsequent dynamic simulation results are trustworthy.

3.1.2. Lattice Parameters and Density

The lattice cell of the water–cellulose model is a cube. The cell parameters show the length of the cell, which can be used to characterize the cell size. The variation in pressure affects the cell size as well as density, and the variation of the cell parameters under different pressures is shown in

Table 1. Cell parameters of water–cellulose model.

Pressure (MPa)	Cell Parameters (Å)
	The Length	The Width	The Height
Untreated	24.86	24.86	24.86
4	22.83	22.83	22.83
5.5	22.79	22.79	22.79
8	22.71	22.71	22.71
12	22.69	22.69	22.69

.

The density of the model reflects the densification degree of cellulose chains under different pressures, which can be used to analyze the mechanical properties of wood cellulose [26]. The density of the water–cellulose composite model after hydrothermal treatment under high pressure is reported in

Table 2. Density of water–cellulose model.

Pressure (MPa)	Density (g/cm3)
	Final	Average	Std. Dev.
Untreated	0.979	0.978	0.020
4	1.242	1.264	0.014
5.5	1.264	1.270	0.015
8	1.268	1.280	0.014
12	1.285	1.288	0.013

.

Where Final represents the density of the model at the end of the dynamic simulation and Average is the average density after dynamic simulation. Std. Dev. is the standard deviation, which indicates the dispersion degree of the density at different moments. The lower the standard deviation, the more reliable the density value.

Analyzing the data in Table 1, the cell parameters decreased from 24.86 to 22.69. The cell parameters decrease as the pressure increases, which indicates that the cell is compressed under high pressure. From Table 2, it can be inferred that the average density of the water–cellulose model increases with increasing pressure. The cellulose volume decreases and the average density increases significantly under high pressure, which demonstrates that the densification degree of wood can be greatly enhanced.

3.2. Hydrogen Bond Analysis

Glucose is the minimum component unit of cellulose, which contains multiple hydroxyl groups, hence the strong hydrogen bond in cellulose [27]. The hydrogen bond network is directly related to the mechanical properties of materials. Hydrogen bonds can be formed within the cellulose chain, between cellulose chains, as well as between cellulose and water molecules. The number of various hydrogen bonds under different pressures can be clearly seen in

Table 3. Number of hydrogen bonds in the water–cellulose model.

Pressure (MPa)	Number of Hydrogen Bonds
	Between Cellulose Chains	Between Water Molecular	Between Water–Cellulose	Total
Untreated	39	56	20	115
4	59	134	150	343
5.5	62	130	160	352
8	64	128	167	359
12	70	127	173	370

.

Water and cellulose form more hydrogen bonds than cellulose chains, which may contribute to moisture diffusion into the cellulose system to form hydrogen bonds. Significant quantities of moisture may weaken the original hydrogen bonding network, making the wood susceptible to fracture.

Table 3 demonstrates that the THM treatment significantly increased the number of hydrogen bonds of the THM cellulose in comparison with those of the untreated model. Table 3 also determines that hydrogen bonding between cellulose chains is higher at high pressure, which is probably owing to the increased cellulose densification caused by high pressure. The volume of cellulose decreases and intermolecular contacts boost, making it easier to form hydrogen bonds [28]. The formation of intermolecular hydrogen bonds not only enhances the structural stability, but also mitigates the negative effect of water intrusion on the mechanical properties of cellulose.

3.3. The Molecular Diffusion Coefficient

The diffusion of small molecules in polymers can be characterized by the diffusion coefficient. The size of the diffusion coefficient can be expressed by the motion range of small molecules in the system. The greater the diffusion coefficient, the more likely it is that water molecules will react with cellulose chains. The investigation of the diffusion coefficient of water molecules contributes to the movement of cellulose chains under high pressure [29]. The mean square displacement (MSD) is one of the important parameters to describe the motion of molecular chains, which is usually used to describe the centroid displacement of molecular chains within a section of time, and can be obtained from Equation (2):

$$\mathrm{MSD}={\displaystyle \sum _{i=1}^n\langle {\left|\overrightarrow{r_i}\left(t\right)-\overrightarrow{r_i}\left(0\right)\right|}^2\rangle }$$

(2)

where $\overrightarrow{r_i}\left(t\right)$ represents the displacement of the molecule i at time t; $\overrightarrow{r_i}\left(0\right)$ represents the initial displacement. The diffusion coefficient D of water molecules can be calculated through the Einstein equation [30], which is related to MSD, represented by Equation (3):

$$D=\frac{1}{6n}\underset{t\to \infty }{\lim }\frac{d}{dt}{\displaystyle \sum _{i=1}^n\langle {\left|\overrightarrow{r_i}\left(t\right)-\overrightarrow{r_i}\left(0\right)\right|}^2\rangle }$$

(3)

when t is sufficiently large, the diffusion coefficient can be solved as Equation (4),

$$D=\frac{k}{6}$$

(4)

where k represents the slope obtained by simulating the MSD curve with the least square method. Analyzing the water–cellulose model produced the MSD curves of water molecules under various pressures.

It can be seen from Figure 3 that the MSD curve of water molecules linearly correlated with time. The MSD curves were simulated to obtain the slope k and the correlation coefficient r², as well as the diffusion coefficients D at various pressures, which are shown in

Table 4. Diffusion coefficients of water–cellulose model under different pressures.

Pressure (MPa)	k	D	R-Square
Untreated	1.9112	0.3185	0.9993
4	1.6162	0.2694	0.9989
5.5	1.4387	0.2398	0.9952
8	0.7320	0.1220	0.9967
12	0.2752	0.0459	0.9983

.

From the analysis of Table 4, the diffusion coefficient of the pressurized model was smaller than those of the untreated model. Furthermore, the diffusion coefficient of water molecules decreases as the pressure increases, implying that water molecules move less flexibly in cellulose, as summarized in Table 4. As a strong polar molecule, water easily forms hydrogen bonds with the cellulose chain, thus inhibiting water movement in cellulose [31], as clearly reflected by the increased number of hydrogen bonds formed between water and cellulose.

The interaction between water and cellulose is increased as a result of the massive formation of intermolecular hydrogen bonds [32], while the diffusion coefficient of water molecules decreases as the adsorption capacity of cellulose to water increases. The rapid diffusion of water in cellulose accelerates the thermal movement of molecules, which is the primary reason for wood fracture. THM treatment reduces the diffusion of water molecules, which improves the stability of molecular chain movement and the deformation resistance of wood.

3.4. End-to-End Distance

End-to-end distance is the distance between the ends of a polymer chain, describing the degree of molecular curl. The curling of the cellulose polymer chain is owing to the non-bonding force in the molecular chain, such as hydrogen bond and van der Waals force, so that the energy of the whole system reaches the lowest value to ensure the stability of the system, which conforms to the entropy increase principle [33]. The end-to-end distance of the water–cellulose model is shown in Figure 4, and the cellulose presents a crimp state.

In this paper, end-to-end distance was used to characterize the stiffness of the cellulose chain. The higher stiffness of the cellulose chain at the microscopic level strengthens the deformation resistance of the wood, which means it is harder to deform due to external forces [34]. The Flory’s characteristic ratio [35] is defined as follows:

$$C=\frac{\langle h_0{}^2\rangle }{N{l}^2}$$

(5)

where $h_0$ is end-to-end distance, N denotes the number of cellulose chain molecules, and l denotes the length of a single bond. The larger value of C indicates the better stiffness of the cellulose chain. Given that N and l remain unchanged in this paper, the change in C is only related to $h_0$.

The end-to-end distance data at each pressure were extracted to obtain the graph shown in Figure 5.

Comparing the untreated model with the THM groups, end-to-end distance is demonstrated in Figure 5. With the enhancement of pressure, the amplitude of end-to-end distance keeps increasing. Hence, higher pressure lowered the curl degree of the cellulose. Large amounts of hydrogen bonds are formed, which makes a great contribution to the increasing end-to-end distance and increases the difficulty of rotation and displacement of cellulose chains. Above all, the equilibrium conformation of molecular chains depends more on the interaction energy rather than the conformation entropy. Non-bonding interaction of hydrogen bonds between water and cellulose is also the reason why cellulose chains are stretched [36]. As a result, the increase in end-to-end distance means that cellulose chain stiffness is strengthened, and its deformation resistance will also be enhanced.

3.5. Mechanical Properties

The mechanical parameters were obtained from the dynamic simulation results at different pressures, which were able to characterize the mechanical properties of the wood. Hooke’s law governs the stress–strain behavior of solid linear–elastic materials, so that the mechanical properties of cellulose can be calculated by Equation (6):

$$\left[C_{ij}\right]=\left[\begin{array}{cccccc}\lambda +2\mu & \lambda & \lambda & 0& 0& 0\\ \lambda & \lambda +2\mu & \lambda & 0& 0& 0\\ \lambda & \lambda & \lambda +2\mu & 0& 0& 0\\ 0& 0& 0& \mu & 0& 0\\ 0& 0& 0& 0& \mu & 0\\ 0& 0& 0& 0& 0& \mu \end{array}\right]$$

(6)

where λ and μ are known as Lamé constants and are often applied to compute mechanical parameters, such as Young’s modulus (E), shear modulus (G), Poisson’s ratio (γ), and so on. The formulas are as follows:

$$E=\frac{\mu \left(3\lambda +2\mu \right)}{\lambda +\mu }$$

(7)

$$G=\mu$$

(8)

$$K=\lambda +\frac{2}{3\mu }$$

(9)

$$\gamma =\frac{\lambda }{2\left(\lambda +\mu \right)}$$

(10)

As shown in Figure 6, each mechanical parameter of the water–cellulose model under different pressures to that of the untreated model was calculated to quantify the influence of the THM process.

Young’s modulus (E) is the ratio of stress to the corresponding strain of the material within the elastic limit, which is commonly applied to indicate whether the material is prone to deformation. A large value indicates that wood is more resistant to deformation. The shear modulus (G), as the ratio of shear stress to strain, also represents the stiffness of the material.

It can be seen from Figure 6a that the Young’s modulus and shear modulus increase with increasing pressure, which indicates that wood stiffness increases. The results showed that the E and G value of the THM model with 12 MPa increased by 12.0381 GPa and 4.4378 GPa compared with those of the untreated model, respectively. Both Young’s modulus and shear modulus increases that are caused by the enhancing densification degree of wood have been investigated [37]. The increase in hydrogen bonds between cellulose chains also increases the van der Waals force of the inter-chain contact [31]. The enhancement of the intermolecular force makes the entire cellulose system form a denser structure, resulting in the decrease in Young’s modulus and shear modulus. In consequence, the wood has a stronger stiffness to resist deformation. It is also consistent with the conclusion that pressurized hydrothermal treatment can improve the Young’s modulus and shear modulus of wood, which Bao et al. [12] studied.

Poisson’s ratio (γ) and the ratio of volume variable to shear variable (K/G) are used as the elastic constants to describe the transverse deformation of the material. The greater the value, the better the flexibility of the material. The results shown in Figure 6b excluded the increasing pressure between γ and K/G values, showing that the flexibility of wood is weakened. An updated review by Yu et al. confirms the reduction in wood toughness attributed to the degradation of cellulose at high temperatures [38].

Different application scenarios of wood have different requirements for its performance. Compared with other industries, wood has higher requirements of stiffness when used as a building material. Pressurized hydrothermal treatment can significantly improve the stiffness and deformation resistance of wood, allowing it to be used in applications that require high wood stiffness.

4. Conclusions

This paper studies the effect of pressurized hydrothermal treatment on the mechanical properties of wood cellulose. Model parameters, hydrogen bond numbers, molecule diffusion coefficients, end-to-end distances, and mechanical parameters of the water–cellulose model were analyzed. The conclusions obtained are as follows:

• High pressure improves the densification degree of wood. The total number of hydrogen bonds in the system, the number of hydrogen bonds between cellulose chains, and the combined hydrogen bonds of water and cellulose are on the rise, which improves the structural stability of cellulose. Reduced diffusion coefficient of water molecules in the cellulose indicates that high pressure weakens the cracking of the wood caused by the bulk diffusion of moisture.
• The stiffness of the polymer chain can reflect the macroscopic mechanical properties of the material from the microscopic level. The analysis found that the increase in pressure caused the continuous increase in the end-to-end distance of the cellulose chain and the decrease in the degree of curling, which was related to the formation of many hydrogen bonds between the cellulose chains. The increase in the end-to-end distance also means that the stiffness of the cellulose chain increases, indicating that the wood can have better deformation resistance at the macro level.
• The mechanical properties of the water–cellulose model were studied, and the Young’s modulus, shear modulus, Poisson’s ratio, and K/G value were calculated and analyzed. With the increase in the pressure, the Young’s modulus and the shear modulus were in an upward trend, and the increase was very large, which means that the deformation resistance and rigidity of the wood have been greatly improved. Moreover, Poisson’s ratio and K/G value decreased with the increase in pressure, and the toughness of the wood decreased. The pressurized hydrothermal treatment of wood can significantly improve its rigidity and deformation resistance, and what is more, its utilization rate in applications requiring high rigidity.