Molecular dynamics simulation of mechanical properties of carbon nanotube reinforced cellulose

Cellulose, hemicellulose, and lignin are the major chemical components in wood paper. Various types of wet and dry strength additives are used to enhance the optical and mechanical properties of recycled paper. One of the possible materials is the carbon nanotube. In order to explore the probability of the use of carbon nanotubes as reinforcing materials and to understand how carbon nanotubes affect the mechanical properties of paper, a single-walled carbon nanotube is inserted into a Iβ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{\beta }$$\end{document} cellulose nanocrystal, and its mechanical properties are studied by using energy minimization and molecular dynamics (MD) simulations. During simulations, the crystals are stretched in the axial direction at a constant speed, and stress and strain are computed and recorded at the atomic level. Our results indicate that carbon nanotube can significantly enhance the mechanical properties of paper.


Introduction
In the paper industry, millions of tons of paper are recycled and reused in a recycling process to save energy and natural wood resources. One of the problems that the paper industry faces today is that the mechanical strength of recycled fibers is greatly weakened by contaminants and by previously intensive drying processes [1].
To improve the mechanical strength of recycled cellulose fibers, various types of drying or wet strength additives, at the nanoscale, are mixed with cellulose fibers in the papermaking process. These strength additives include polyacrylamide, polyvinyl, polyamide, and polyamide-epichlorohydrin, as well as different types of fillers. One of the possible fillers or reinforce materials is carbon nanotubes (CNTs), which possess extremely high mechanical strength. It is possible that a small amount of CNTs could be mixed with cellulose to achieve novel properties by combining properties from the two different constituents. In fact, it was found that the mechanical, chemical, electrical, optical, magnetic, and electro-and magneto-optical properties of the CNTs have extensive applications in many different areas. For example, theoretical and experimental results reveal that the Young's modulus of CNTs could be as high as 1-5 TPa 2 and the tensile strength could be as high as 400 GPa [3,4]. Moreover, their aspect ratio is very large, and their diameter scale may be comparable with various cellulose chains. Thus, CNTs naturally suggest themselves for use as one component of composites to enhance the mechanical properties [5,6]. Indeed, people have experimentally studied the mechanical properties of polyimide/ single-wall CNT composites [7]. They found that the storage modulus increased by 65% upon the addition of 1% CNTs by volume [8]. Some authors simulated the stress-strain behavior of polyimide CNT composites by using molecular dynamics (MD) simulations and found that long carbon CNTs could significantly reinforce mechanical properties [9]. It suggests that cellulose mixed with carbon nanotubes may result in high mechanical properties.
Several attempts to use nanotubes in the cellulose community have been made recently. [10,11] used cellulose to build nanotubes and provided evidence that cellulose nanotubes have an advantage, as compared with cellulose nanofibers, in hydrophobic property by using density functions and molecular dynamics (MD) calculations [5,6,10,[12][13][14][15][16][17][18][19][20]. In order to employ the high mechanical strength of CNT, Shishehbor and Pouranian [21] used cellulose to wrap around CNT and obtained tunable and strong mechanical properties in their MD simulations.
In this work, we propose to directly and simply replace the cellulose chain with a CNT to build a CNT-reinforced cellulose crystal. The advantage of replacing cellulose chains with CNTs is that the cellulose crystal structure and properties are kept unchanged. It is expected that new materials made of CNT-cellulose crystals can be created. The present work is not done in labs and is limited to computer simulations. Atomistic modeling is used to investigate the stress-strain behavior of cellulose nanocrystals and their composites with a single-walled CNT. We will focus on numerical simulations of the performance of composites of CNT-reinforced cellulose and on the evaluation of their improvement regarding mechanical properties. A series of MD simulations is conducted to study how CNT influences the mechanical properties of composites of CNT and cellulose.
Our results indicate that CNT can significantly enhance the mechanical properties of cellulosic materials, which is potentially important for recycled paper.
The next section will introduce how to construct cellulose crystal, CNT crystal, and CNT-reinforced cellulose. The method of molecular dynamics simulations is also briefly introduced. "Results" section presents simulation results of densities, energies, and Young's moduli of the three crystals. "Discussion and conclusion" section draws a conclusion.

Cellulose nanocrystal
Cellulose can be presented in different crystalline phases as I , I , II , and III I . Native plant cellulose is a mixture of I .and I . Their ratio of I to I depends on original fiber sources. In wood fibers, I is a major component and adapted in this simulation work.
The simulated cellulose nanocrystal has fifty chains, each chain with nineteen repeated β-(1 → 4)-linked D-glucose units [24] as shown in Fig. 2. The total atom number in cellulose crystal is 19050, which are located in a parallelepiped box with a triclinic symmetry. The basis three vectors of the simulation box are A = (x, y, z) = (38.92, 0, 0),Å B = (− 4.642, 40.741, 0) Å, and C = (0, 0, 93.42) Å, respectively. In other words, the width of the built nanocrystal is five times of the unit cell in the a and b crystallographic directions and nine times of unit cell in the c crystallographic direction (see Fig. 2). The cellulose chains are along the c-direction of the unit cell or the z-direction of the simulation box.

Carbon nanotube
To verify the liability of the method, a carbon nanotube crystal is built, and its mechanical strength is tested and compared with experimental results.
First, a single-walled carbon nanotube is built by mapping the atoms of a 2D graphene sheet to the surface of a cylinder. Different types of CNT can be denoted by two integers ( n 1 , n 2 ) , called chirality, where the n 1 and n 2 refer to the chiral vector in the primitive Bravais lattice of the 2D graphene sheet. The diameter of the carbon nanotube depends on the integers of n 1 and n 2 . More details can be found in [25] work. Following their work, carbon nanotube builder software was developed at Western Michigan University [9] and is used here again. Two different To test the CNT model, a single CNT is simulated, and its mechanical properties are tested and compared with other authors' results. In MD simulations, the CNT axis is initially set along the z-direction of a simulation box, as shown in Fig. 3, and its stress as a function of stain in the z-direction is obtained from simulations.

Composite of CNT-reinforced cellulose
In order to build composites of CNT and cellulose, one of the cellulose chains is removed from the cellulose nanocrystal as shown in Fig. 4a and replaced by a singlewalled CNT of (4, 4) as shown in Fig. 4b. The single CNT is shown in Fig. 3. Its axis is placed along the same axis as the original cellulose chain to keep the cellulose structure unchanged. In addition, diameter of single walled CNT with ( 5,5) is also suitable for replacing cellulose chain. There are 19,277 in the (4,4) CNT cellulose and 19,425 atoms in the (5,5) CNT cellulose.

MD simulations
An open software of molecular dynamics simulation package called LAMMPS 26 is employed to simulate various physical properties, such as potential energy, total energy, density, and hydrogen bond energy.
In simulations, the popular CHARMM carbohydrate force fields, all36_carb and all36_cgenff [27,28], are adopted for cellulose and carbon CNTs. In the Charmm force fields, total energy E t in the cellulose system can be presented by The first term in the right side of the equation represents the bonding energy, the second term is the energy contribution due to the vibration of bonded two atoms (only the Charmm force field has this term), the third term is the angular energy between two bonds, the fourth term stands for the dihedral energy, the fifth term is improper torsion energy, the sixth term is van der Waals non-bonded Lennard-Jones (LJ) potential energy, and the last term is electric Coulombic interaction energy.k b ,k UB ,k a ,k d , and k i are the spring coefficients of corresponding terms; ij is the depth of energy well of Lennard-Jones potential; R min,ij is its minimum energy positions; q i is the charge of atom i; 0 is dielectric constant;r 0 ,s 0 , 0 ,ψ 0 , and ψ 0 are corresponding equilibrium positions or angles. Both the LJ and electric Coulombic forces in the sixth and last term of the right side of Eq. 1 are cut off at 12 Å. From the gradient of total energy, the forces acting on each atom can be calculated. Then, based on the forces, a leap-frog algorithm is utilized to update the velocities and positions of all atoms.
Hydrogen bonds (H-bonds) are critically important for the stability of cellulose structure. The Deriding model with three body forces [29,30], is selected to track the hydrogen bonds.
In MD simulation, first, an energy minimization is imposed on all atoms, then they run in a NVT ensemble for 40,000 timesteps at a temperature of 298 K. Next, the system runs at a constant pressure of 1 atm and at the same temperature for the additional 260,000 timesteps, allowing the system to reach an equilibrium state. During the constant pressure run, 1 atm of pressure is kept in the x-, y-, and z-directions, and a zero pressure is kept in the xy-, yz-, and xz-directions. Volume or density and different types of energies can be measured in the NPT time integration via the Nose/Hoover algorithm [31,32]. One step time interval is 0.5 fs. During simulations, periodic boundary conditions are imposed on both cellulose and CNTs in all three directions. It is noted that a correct periodic condition for CNT presents a model of an infinitely long CNT in simulations.

Stress and strain
The purpose of the present work is to use CNT to reinforce the mechanical properties of cellulose. Since the CNT lies along the axial or z-direction of cellulose, the stress zz and Young's moduli in the reinforced z-direction are our focus for cellulose crystal, CNT crystal, and CNT-cellulose composite.
At the atomic level, the stress tensor zz on atom I can be calculated from its two-, three-, and four-body forces by the first term of the right hand of the equation above is kinetic energy contribution (where v z is the velocity in the z-direction) and the second term is the viral contribution, which can be computed by The first term in the right hand of the equation is a pairwise energy contribution where the summation is over the n p neighbors of atom 1, r 1 , and r 2 are the position of the 2 atoms, respectively, and F 1 and F 2 are the force on the 2 atoms from the pairwise interaction. The second term is a bonding contribution of a similar form for the n b bonds. The third term is the n a angular bond contribution; the fourth term is the n d dihedral energy contribution, and the fifth term is the n i improper contribution.
During crystal stretching in the z-or axial direction, the stress in this direction is calculated and averaged over all atoms at each time step.
In simulations, after the system arrives at its equilibrium state, the cellulose nanocrystal is stretched at its both ends . 4 a One cellulose chain is removed from a cellulose crystal (left); b a single-walled carbon nanotube is inserted into the crystal (right). The CNT is shown in Fig. 3 along the z-direction at room temperature while 1 atm pressure is kept in the x-and y-directions. During stretching, the one-time step is in 0.25 fs, and the stretching strain rate is at 0.001 ps . The total time steps for stretching are 400,000, and the total displacement is 9.342 Å. The results of stress and strain are computed and recorded. Young's modulus is evaluated from the slope of the initial linear portion up to 2% of the small strain. The same applies to single CNT and CNT-reinforced cellulose. It is noted that like water molecules in ices at a low temperature, there are no covalent bonds between each cellulose chain in its crystal. There are hydrogen bonds (H-bonds) only between each cellulose chain. The cohesion between cellulose chains mainly comes from hydrogen bonds. There are no covalent bonds between cellulose chains and CNT either. When we stretch the CNT cellulose crystal, the hydrogen bonds are broken first without breaking the covalent bonds of CNT since its covalent bonds are much stronger than H-bonds. Therefore, the present model is correct for H-bonddominated materials such as the ices and cellulose.

Density
In the NPT simulations, the total number of atoms, temperature, and pressure are kept constants; volume is varied and computed; and then the density of the cellulose nanocrystal is predicted. The results of density as a function of the time step are plotted along with temperature in Fig. 5. It shows that after 200,000 timesteps, the values of density and temperature are stable with only small fluctuations, indicating that an equilibrium state has been reached. In the present work, the predicted average density of the cellulose nanocrystal is 1.68 g cm 3 . Experimental results for the density of I cellulose nanocrystals are 1.676 g cm 3 reported by Nishiyama et. al. (2002) and 1.625 g cm 3 reported [33]. [34] used a quantum mechanical density function to theoretically predict the density of 1.61 g cm 3 . It is clearly demonstrated that our simulation results are in an excellent agreement with others' experimental and numerical results. Figure 5 also shows that the average temperature is well controlled around 298 K as expected and that hydrogen bond energy is stable when the system reaches an equilibrium state.

Young's modulus
We first produce the relationship between stress and strain for cellulose crystals by stretching the simulation box in the z-direction. During stretching, the stress in this direction is calculated and averaged over all atoms at each time step. The results for stress, averaged over every 100 steps, against strain are depicted for the cellulose crystal in Fig. 6.
The slopes of stress-strain curves at their initial portion represent the Young's moduli. Our simulation result for Young's modulus in the axial or z-direction for I cellulose nanocrystals is 144.6 GPa. In literature, experimental values of Young's modulus have been determined by using different experimental techniques. For instance, x-ray diffraction measurements [35,36,37,38,39], of the axial elastic modulus measured along the longitudinal To explore the function of H-bonds, the same simulation is conducted, except that all H-bonds are removed. Figure 6 also compares the stress-stain curve of cellulose crystals with that of cellulose without H-bonds. The slope of the curve is much larger for cellulose than for cellulose with removed H-bonds. The Young's modulus for later is 52.82 GPa, which is 63% lower than the normal cellulose. It is also seen that tensile strength is 7.45 GPa for cellulose crystal and only 0.85 GPa for that without H-bonds, which is 89% lower. The total H-bond energy and number of H-bonds against time during stretching are recorded and plotted in Fig. 6 (right) too. It shows that during stretching, H-bonds are broken and the number of H-bonds reduces from 10,414 to 9840 while H-bond energy increases due to H-bonds breaking. This dynamic information during the tensile strength numerical test demonstrates that H-bonds play a critical role in the contribution of mechanical strength. This dynamic information was not reported by others although many authors [40], Eichhorn and Davies [41] have revealed the important role of H-bonds in mechanical properties.

Young's modulus of single CNT
To verify the model of CNT, following the same procedure, MD simulations are carried out for a single carbon nanotube. During stretching, the results of stress against strain are reported in Fig. 7. From their initial slopes, we obtain that Young's module is 1.31 TPa for (4,4) CNT and 1.20 TPa for (5,5) CNT.
The experimental measurement from a bending test on Young's modulus was an average value of 1.28 ± 0.59 TPa. A tight-binding calculation produced a value of the order of 1 TPa (Henandez and Subbaswamy 1998). In addition, the predicted density in our simulations is 1.38 g∕cm 3 , as compared with the experimental data range of 1.3-1.4 g∕cm 3 .
Our results agree well with those of other authors, indicating that the present model of CNTs is trustworthy.

CNT-cellulose composite
Finally, the result of the stress-strain curve is presented for the composite of CNT-reinforced cellulose and compared with that of cellulose crystal in Fig. 8.
The slopes of the initial stress-strain curves are larger for the composites of CNT-reinforced cellulose than for the pure cellulose crystal. The Young's modulus increases from E = 144.6 GPa for the pure cellulose crystal to E = 171.97 GPa for the composite of (4,4) CNT-reinforced cellulose. This is a 19% increase in Young's modulus by replacing a cellulose chain with a CNT although only 2% of the cellulose chain is replaced by CNT. The tensile strength increases from 7.40 to 7.98 GPa. This is an 8% increase in tensile strength. For (5,5) CNT cellulose, the Young's module increases to 165.91 GPa, and the tensile strength increases to 8.67GPa. This is a 17% increase in tensile strength and a 15% increase in Young's modulus. It is evident that carbon nanotube can significantly reinforce mechanical strength.

Discussion and conclusion
In this work, the molecular dynamics simulation package LAMMPS is employed to simulate the mechanical properties of I crystalline cellulose in NPT ensembles, where CHARMM force fields are adopted. The predicted densities and Young's modulus of cellulose crystals and CNT crystals from the simulations are in excellent agreement with the experimental results. Subsequently, one cellulose chain is replaced by a single-walled carbon nanotube in the cellulose crystal, the Young's modulus increases by 19%, and the tensile strength also largely increases. It is demonstrated that CNT can dramatically enhance the mechanical properties of the cellulose system. Also, we have shown that the number of H-bonds decreases and the H-bond energy increases due to H-bonds continuously breaking during the numerical breaking test.
It is noted that if CNTs are randomly distributed in space and mixed with cellulose, mechanical strength should be smaller than the present case where the CNT enhances the strength in the axis direction.
It is pointed out that the periodic boundary condition is imposed on both cellulose chains and CNT in the composite, implying that the CNT is infinitely long in the present model. The infinitely long CNT is realized by bonding the atoms in the simulation box to the ghost atoms on the outside of the simulation box. The interaction between the cellulose and CNT atoms is the nonbonded LJ forces. The large mechanical strength of CNT itself contributes to the improvement of the cellulose composite with the infinitely long CNT model.
It is recognized that adding polar function groups, such as OH or COOH groups, on the surface of CNT will further improve the adhesive interaction between cellulose and CNT atoms. We have finished this type of work, and the results will be reported in a different article.