Abstract
At present, most high-performance cellulose matrix composites only use cellulose as reinforcement material, which is an obstacle to maximize the advantages of nanocellulose in structure and properties. The development of new functional nanocomposites with cellulose as the main component can better meet people’s needs for high-performance and degradable composites, which requires a comprehensive and thorough understanding of cellulose. Considering the limitations of physical experiments, we performed molecular dynamics simulation of the uniaxial tensile behavior of the cellulose system at three different strain rates (10−4/ps, 10−5/ps, and 10−6/ps), and the stress-strain responses of cellulose systems at different strain rates are obtained. The effect of the strain rate on the mechanical properties of amorphous cellulose system during the tensile processes is analyzed. The deformation mechanism of cellulose amorphous system during the tensile processes is characterized by the energy changes of the different terms including dihedral angle torsion term, bond tensile term, angle bending term, and nonbond term. Structural evolution of the cellulose crystal system during the tensile processes is used to explain the failure mechanism of cellulose. The kinetic simulation results show that the mechanical properties of the cellulose amorphous system increase with the increase of strain rate. Compared with the strain rate of 10−5/ps, the elastic modulus of the system increases by 6.73 GPa at the strain rate 10−4/ps. During the tensile processes, cellulose amorphous region adapts to the applied load mainly through the stretching of the cellulose macromolecular chains, i.e., the deformation of bond lengths and bond angles, without any breakage of the molecular chains. The main causes of chain lengthening at different strain rates are different. The failure of cellulose is caused by the slip and rearrangement of some molecular chains in the crystal structure.
1. Introduction
With the rapid development of industry, traditional materials have been unable to meet the increasing performance requirements of people and the extensive use of traditional materials will aggravate the damage to the environment. There is a growing demand for new, biodegradable, and high-performance materials. The ultralight, ultra-high-strength, and biodegradable composites based on cellulose meet this demand [1, 2]. As the most abundant renewable biomaterial resource on earth, natural cellulose has received a lot of attention from researchers due to its outstanding physicochemical properties such as high aspect ratio (usually 10–100 but some cellulose can reach 1000, e.g., bacterial cellulose), high strength, and high functional activity on the surface and can be used as a substrate to synthesise composite materials with excellent properties. In the past decades, the research on its physical and chemical properties has never stopped [3–5]. It has been shown that cellulose composites offer significant improvements in mechanical properties, thermal stability, and degradability compared to conventional materials [6–8]. Alemdar and Sain [9] mixed nanocellulose and starch in different proportions to obtain nanocomposites, and the results showed that the properties of the composites have been greatly improved compared with that of the pure starch materials. Saba et al. [10] prepared composite materials by mixing nanocellulose with epoxy resin to investigate the effect of nanocellulose content (0%-1%) on the properties of composite materials. The results showed that when the mass fraction of nanocellulose was 0.75%, the mechanical properties of composite materials are improved. However, most of the current development uses cellulose as the reinforcement material. To develop new functional nanocomposites with cellulose as the main component, maximize the superiority of nanocellulose structure and performance. This requires a comprehensive and thorough understanding of cellulose. The insolubility of cellulose makes it challenging to extract cellulose from the cell wall for experimental characterization. If chemical treatment (such as strong alkali and heating) is used to remove cell wall components to obtain cellulose, this process will change the properties of the extracted cellulose resulting in biased test results [11]. Coupled with the small size of cellulose, it has been a challenge to characterize its mechanical properties through experimental tests [12]. Computer-based molecular dynamics simulation not only provides a method to solve these problems but also can reveal the properties of cellulose from the molecular level, providing a theoretical basis for the development of high-performance cellulose-based materials. Previous work has shown that molecular dynamics simulation is an effective way to analyze the mechanical behavior of materials at the microscopic scale [13, 14]. Fu et al. [15, 16] used molecular dynamics method to construct the indenter to simulate the nanoindentation of Cu/Ni multilayer films. The results of this study have a profound impact on unveil potential mechanisms related to the excellent mechanical performance and expanding the application prospects of metallic multilayered films. This demonstrates the effectiveness of molecular dynamics simulation methods in studying the properties of materials. Poma et al. [17] measured the elastic modulus of cellulose crystals by molecular dynamics and compared the simulation results with experimental data to verify the feasibility and accuracy of the molecular dynamics approach in studying the properties of cellulose. Wu et al. [18] used the fully atomistic model of the cellulose to calculate its elastic properties. The effects of chain length, strain regime, and calculation method on axial modulus prediction using atomistic simulation are analyzed. They found that chain length, boundary conditions, and calculation method affect the magnitude of the predicted axial modulus and the uncertainty associated with that value. However, her study was only for a single cellulose chain. Tanaka and Iwata [19] measured the elastic modulus of natural cellulose crystals under different force fields using the molecular simulation technique. The measured value of elastic modulus under the CFF91 force field varies from 124 GPa to 155 GPa, and the measured value in the COMPASS force field is 124.6 GPa–130.3 GPa that was nearly equal to the observed value of 138 GPa. The analysis of the models of different sizes obtained the conclusion that the lateral (that is, intermolecular) interactions between cellulose chains were found to play an important role in the expression of the mechanical properties of cellulose crystal. Zhang et al. [20] used molecular dynamics under ReaxFF reaction field to simulate the tensile behavior of the amorphous structure of cellulose, and measured the elastic modulus of 9.37 GPa and the shear modulus of 3.94 GPa. The deformation mechanism of the structure is not involved in her research. Most of the current work focuses on the mechanical properties of cellulose. However, little research has been performed on the tensile behavior and subsequent failure mechanisms of cellulose which is an obstacle to thorough understanding of the potential properties of cellulose and the development of new cellulose-based materials.
In this study, the Materials Studio software [21] was used to construct the crystalline and amorphous models of cellulose. The LAMMPS [22] software package was employed to simulate the uniaxial tensile processes of the model, and it captures the stress response of the model to tensile strain in the tensile directions. The elastic modulus, Poisson’s ratio, yield stress and strain, and ultimate stress and strain were calculated for the cellulose amorphous system (it should be noted that the mechanical properties of cellulose crystal structures have not been analyzed in this study because they have been covered in existing studies). The energy changes of each component during the simulation of amorphous system are discussed to reveal the tensile mechanism. Finally, the stress-strain response and the tensile evolution of the crystal structure are discussed to derive the cellulose failure mechanism. The temperature of the tensile simulation is normal room temperature (300 K).
2. Molecular Dynamics Simulation
2.1. Construction of Amorphous Cellulose Model
Natural cellulose is a linear chain of ringed glucose molecules. The repeat unit is comprised of two anhydroglucose rings ((C6H10O5)n, n = 10,000 to 15,000) linked together through oxygen covalently bonded to C1 of one glucose ring and C4 of the adjoining ring (1–4 linkage) and so-called the β-1,4 glucosidic bond [3] (Figure 1). The cellulose chain includes crystalline and amorphous regions. The molecules in the crystalline region are arranged neatly and tightly, and can present a clear X-ray diffraction pattern, while the molecules in the amorphous region are irregularly arranged with large gaps. There is no clear division between the crystalline and amorphous regions in the cellulose chain, and it is a gradual transition state. It is difficult to discuss the crystalline and amorphous states separately during the experiments, whereas the molecular simulation technique allows to easily unveil the influence of two forms on the research problem.

Mazeau and Heux [23] performed molecular dynamics simulations for cellulose with different degrees of polymerisation (DP), and the results showed that the physicochemical properties obtained during the simulation do not differ significantly. Chen et al. [24, 25] used molecular dynamics method to construct cellulose chain with DP20 to investigate the mechanical properties. The results showed that the simulated values were basically consistent with the experimental values. Combining its natural state and the results of studies on the degree of polymerization of cellulose, the cellulose chains of DP15, DP30, and DP45 are established, respectively (Figure 2). Table 1 shows detailed information on cellulose chains with DP of 15, 30, and 45. The amorphous model was constructed with three cellulose chains of DP15, DP30, and DP45. The amorphous model was constructed using the method of constructing amorphous polymers, first proposed by Theodorou and Suter [26]. The target density of the models is set at 1.5 g/cm3, which is made to the density of pure cellulose. The model is shown in Figure 3. The size of the simulated box is 51.2 Å × 51.2 Å × 51.2 Å.

(a)

(b)

(c)

(d)

2.2. Construction of Cellulose Iβ Crystal Model
There are two main types of cellulose crystals in nature, Iα and Iβ. Cellulose Iα crystals are triclinic crystals with only one cellulose chain, while cellulose Iβ crystals are monoclinic crystals with two cellulose chains parallel to each other. The Iα crystal phase is high in bacteria cellulose and algae cellulose. In contrast, the Iβ crystalline phase is mainly found in some higher plants, like cotton, trees, and bushes [27]. It is clear that cellulose Iβ crystals are the focus of our study. The initial structure of the cellulose Iβ crystal lattice used in the present investigation was the one provided by Nishiyama et al. [28]. Table 2 displays the lattice parameter of cellulose Iβ models employed in this work. Then, we stacked the unit-cells into the size of 4 × 4 × 8 (Figure 4). The crystal system contains 10,752 atoms.

2.3. Simulation Details
The whole simulation process is divided into two stages. Before the tensile simulation, the equilibrium relaxation stage of the system should be carried out to minimize the energy of the high-energy configurations in the model. The force and energy stop tolerances are 10−6 (kcal/mol)/Å and 10−4 kcal/mol, respectively. The cellulose model was structurally optimized, and then, performed a dynamic simulation of 300 K in the NVT (constant number of atoms, volume, and temperature) ensemble for 100 ps with timestep 1 fs. After the equilibrium relaxation stage, the model structure tended the nature cellulose.
On the basis of the previous full relaxation, the tensile strain is applied in X-direction of the simulation cell. 107 times (the strain rate 10−4/ps) tensile dynamics simulate were performed on the model under the NPT ensemble with timestep 1 fs, the strain rate of the tensile deformation is controlled by changing the tensile times from 107 to 109 to achieve strain rate from 10−4/ps to 10−6/ps. The simulated system is equilibrated in the NPT ensemble for time 100 ps with 1 atmosphere of pressure applied in the two nondeformed directions to achieve the system stability. After the tensile equilibration, the stress-strain data for the simulated system is extracted. Make stress-strain curve and perform calculations of related mechanical properties. The lengths of the two nondeformed directions are used to calculate the Poisson’s ratio of the response.
ReaxFF force field is used in both structural optimization and tensile simulation of the model, which has been proved to be suitable for molecular simulation of cellulose [29–31]. Periodic boundary conditions are applied in all directions. The particles in the system will be copied by simple “projection” under periodic boundary conditions, but the calculation process is not involved. To a certain extent, the limitation of the atomic model scale can be eliminated using this boundary condition, and the computer simulation time can be reduced. The temperature and pressure are controlled in the simulation process by Andersen and Berendsen method, respectively [32, 33].
3. Molecular Simulation Results
3.1. Effect of Strain Rate on Mechanical Properties of the Amorphous Cellulose
The stress-strain behavior of the amorphous system subject X-direction strain at a rate of 10−4/ps is shown in Figure 5. Through the stress fitting curve, the yield stress 0.29 GPa can be obtained, which corresponds to the strain 0.033, the ultimate stress 0.39 GPa, and the strain 0.15. The strain on the left side of the yield limit is identified as a small deformation region, and its mechanical properties can be calculated by data. On the right, it is identified as a large deformation region, and the internal structure of the model will change permanently. In the elastic deformation stage, the slope of the fitting line is elastic modulus. The elastic modulus of the cellulose amorphous system was 8.78 GPa at the strain rate of 10−4/ps, which was consistent with the available literature data [34]. From Figure 5, it can be seen that when the tensile simulation is performed in the X-direction, the contraction is performed in the Y- and Z-directions, which can be found to be consistent according to the curve, and the Poisson’s ratio VY/X = VZ/X = 0.2 is calculated. The stress-strain data at three strain rates of 10−4/ps, 10−5/ps, and 10−6/ps are presented together with the deformations in the other two nondeformed directions (Figure 6). In order to more intuitively show the effect of strain rate on its mechanical properties, Tables 3 and 4 show the yield point, ultimate stress, elastic modulus, and Poisson’s ratio of the system at different strain rates.


According to Table 3, the strain rate increased the elastic modulus of the cellulose amorphous system becomes larger. When the strain rate was larger, the system quickly reaches the yield limit, and the internal molecular chains do not have enough time to adapt to the strain by changing the chain structure, bond length and bond angle, corresponding to a small strain, which leads to a higher elastic modulus. When the strain rate of the simulated system decreased by an order of magnitude, the corresponding elastic modulus values do not show such a pattern—order of magnitude decreased. The trends of yield stress and ultimate stress were the same as the elastic modulus, both decreased in value as the strain rate decreased. The strain value at which the yield point appears increased sequentially for the three strain rates. The small strain rate causes the simulated system to stretch slowly and the strain value at which the yield point appears was large. The law of strain value change for the ultimate stress is also the same. From the data in Table 3, it can be concluded that the mechanical properties of the cellulose amorphous system are better at the strain rate of 10−4/ps. At this strain rate, the elastic modulus of the system is 8.78 GPa, which is 4.28 times of the strain rate of 10−5/ps. Poisson’s ratio values in the Y-direction and Z-direction at the same strain rate were basically the same, indicating that the deformation of the cellulose amorphous system is isotropic. It is interesting to note that the Poisson’s ratio in the unstretched direction of the simulated system increased as the strain rate decreased. With the decreased strain rate, the corresponding strain at the yield point increased, and the deformation in the unstretched direction increased, leading to the increase of Poisson’s ratio.
3.2. Effect of Strain Rate on Tensile Deformation Mechanism of Amorphous Cellulose
According to Figure 6, the stress in the simulated system is along greater than 0 throughout the tensile processes (regardless of the magnitude of the strain rate). Throughout the simulation, the macromolecular chains in the amorphous region did not undergo failure behavior and adapted to the loaded stress by changing the bond lengths and bond angles in the macromolecular chains. Stretching is a dynamic process and it is difficult to monitor the structural changes in the molecular chains (stretching torsion of bond length, bond angle, and dihedral angle), and changes in the nonbonding interactions between molecules in real-time. However, it is can be described the process by calculating the energy of each part of the system and to explore the intrinsic stretching mechanism of amorphous system. Energy changes of each part of cellulose amorphous system are calculated at three strain rates (Figure 7).

(a)

(b)

(c)

(d)
As can be seen from the graph above, the cellulose amorphous system has very smooth energy changes at the strain rate 10−6/ps, with the data basically fluctuating above and below the initial value. This means that at the strain rate 10−6/ps, the cellulose amorphous system is subjected to less stress and has enough time to adapt to the tension by changing the structure of the molecular chain. The energy of each component changes all the time but with little fluctuation. Comparing the energy changes of the system at the strain rates 10−4/ps and 10−5/ps, it is found that the bond stretching term does not change much at the strain rates 10−5/ps, while the angle bending term, nonbonding term and dihedral angle torsion term change dramatically at the strain rates 10−4/ps and 10−5/ps. It shows that the bond was more stable than the bond angle and the dihedral angle in the cellulose amorphous system; in other words, it requires more energy to stretch than the bond angle and dihedral angle. When the strain rate is 10−5/ps and 10−6/ps, the small stress cannot cause the drastic change of the bond. Comparison of the energy changes at the three strain rates shows that the bond stretching term, angle bending term, dihedral angle torsion term, and nonbonding term change more at the strain rates 10−4/ps, suggesting that the change in bond length, bond angle, and dihedral angle within the molecule is the main reason for the stretching and lengthening of the system. At the strain rates 10−5/ps and 10−6/ps, the bond stretching term does not change much, and the stretching length of the system at this point is mainly caused by changes in the intramolecular bond angle and dihedral angle.
3.3. Effect of Strain Rate on the Tensile Deformation Mechanism of Crystalline Cellulose
The stress values during the stretching of the crystalline structure appear negative (Figure 8(a)), which implies a partial failure of the cellulose chains in this region during the tensile processes. To illustrate the failure mechanism in the cellulose crystal region, correlating key points within the stress-strain behavior to the evolution of the crystals’ atomic configuration during deformation it is possible to identify the failure mechanisms of the crystal structure at different strain rates.

(a)

(b)
As shown in Figure 8(a), there are specific points on the stress-strain plot identified at each strain rate: the origin (0), the ultimate point (1, 2), the inflection point (3, 4, 5), and 100% strain (6). The corresponding atomic structures viewed from the Y-direction at each of these points are shown in Figure 8(b) (the middle area is colored to make the motion trajectory more visible). Data from 50% to 95% were not shown because the stress remains around zero, and thus provides little additional information. It is worth noting that the height and width of each structural snapshot was different, which is caused by the selection of different strain points and the Poisson’s ratio of the system. For all strain rates, multiple peaks are observed within the first 50% of tensile deformation where the number, positions and magnitudes of the peaks vary with strain rate. It can be clearly seen that the number of peaks and valleys in the stress-strain plot decreased as the strain rate decreased. At large strain rate, the whole stretching time was relatively short and the strain was large, and the van der Waals force damage process between the layers of the crystal structure chain is sequential, with the middle layer being destroyed first before going to the sides, leading to the failure of some of the cellulose chains. The smaller the strain rate and the longer the simulation time, the milder the strain. The loading strain is uniformly distributed throughout the crystal structure, and the interlayer damage to the chains tends to act more simultaneously. This leads to a smaller number of peaks and valleys at a strain rate of 10−6/ps. When uniaxial stretching of the X-direction is performed, there are no covalent bonds in the macromolecular chain to accommodate the strain (the effect of covalent bonding will only be reflected in the stretching in the Z-direction). All strain is accommodated by breaking the interaction forces between molecular chains (hydrogens bond, van der Walls force) causing sliding, reorientation, and local displacement of cellulose chains, ultimately leading to failure of the crystalline region. At three strain rates, the stress increased from the origin 0 to the ultimate point 1. The atomic structure does not change significantly during this process but the lattice spacing was an increase. With the continuous increased of strain, the van der Waals forces between the chain layers and the hydrogen bonding between the molecular chains are disrupted, leading to the slip of the cellulose chains. These slips result in local relaxation of the stress (e.g., after point 1 in Figure 8(b)). The stress then increased with further strain as deformation occurs through localized displacement of the chains. However, the peak stress decreases a lot relative to point 1. This is the result of partial failure of the cellulose chain. As for point 6, the van der Waals forces between the molecular chain layers and the hydrogen bonding interactions between the molecular chains in the crystal structure have been completely destroyed at this point. From the corresponding snapshots, it was easy to see that the distances between cellulose chain layers are of different sizes, and the structure was loose. Comparing the atomic snapshots at all strain rates, it was found that the smaller the strain rate the looser the atomic snapshot shows the crystal structure. The smaller the strain rate, the longer the whole simulation is stretched, and the strain will have time to act on the atoms in the molecular chain, which will change the bond lengths and bond angles in the system to some extent. In other words, the van der Waals forces, hydrogen bonds, bond lengths, and bond angles in the crystal structure will all adapt to the strain at this point, eventually leading to the loosening of the structure.
4. Conclusions
(1)During the tensile processes of the cellulose amorphous system at three strain rates, there is a significant positive correlation between the mechanical properties of the system and the strain rate. The mechanical properties of the amorphous system were best at a strain rate of 10−4/ps, when the modulus of elasticity is increased by 4.28 times compared to that at a strain rate of 10−5/ps.(2)The deformation mechanism of cellulose amorphous system during the tensile processes was characterized by the energy changes of the different terms. The amorphous system exhibits plasticity characteristics, allowing elongation without breakage of the molecular chains. The reasons for the tensile elongation of the simulated system were different for different strain rates. At a strain rate of 10−4/ps, changes in intramolecular bond lengths, bond angles, and dihedral angles are the main contributors to the elongation of the system by stretching, and at smaller strain rates of 10−5/ps and 10−6/ps, bond length changes contribute less to the elongation of the system.(3)At the same stretching length, cellulose chain failure occurs in the crystalline region and the failure takes the form of slip and rearrangement of cellulose chains. Lower strain rates lead to a looser crystalline structure.Data Availability
The ReaxFF force field data used to support the findings of this study are included within the supplementary information files. The cellulose Iβ crystal structure (deposition number: 810597) used in this study was downloaded from the Cambridge Crystallographic Data Centre. Other relevant data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interests.
Acknowledgments
This work was financially supported by the Bintuan Science and Technology Program (2021CB036), Innovation Research Team Program of Tarim University President Funded (TDZKCX202202), and the President’s Fund Project of Tarim University (TDZKSS202111).
Supplementary Materials
The field resx file in the supplementary material is the force field file during the kinetic simulation, and the cellulose Iβ-815079 file is the crystal file used in the model construction. (Supplementary Materials)