Currently, cellulose exists in the largest amount of resources on Earth. In recent years, a technology for manufacturing this cellulose from wood to fiber has been developed, and it is expected that cellulose can be applied to various fields as a kind of sustainable materials . This type is called a cellulose nanofiber (CNF) and has a diameter of nanometer order . CNF has the highly hierarchical structure, and the minimum component of CNF is called cellulose microfibril (CMF). Multiple cellulose molecules are polymerized to form a linear molecular chain, and 30 to 40 chains are regularly formed into a bundle. By intermolecular hydrogen bonds between chains, a molecular sheet is formed in the direction perpendicular to the axis of the chain. CNF is considered to be used as a structural member for automobiles and is expected to construct a fiber-reinforced composite material . It will be an alternative material for carbon or glass fiber, having the advantage concerning resource amount and low environmental load. The report on the evaluation of experimental mechanical properties of CNF has been increasing trend in recent years . However, the mechanical properties of CMF are more required because CMF is dispersed in a matrix of composite and is nanoscale, and its experiment is not easy . Therefore, CMF has been researched by a numerical simulation. For example, there is a report that estimates Young’s modulus and tensile strength by tensile simulation in the axial or molecular sheet direction . It is also pointed out that twist around the axis exists in CMF  - , and it is reported that the cause of twisting is discussed in terms of hydrogen bonding . In considering that CNF will be integrated in actual structural and mechanical equipment, it will be also important to know the material’s behavior in which it is subject to various deformations such as bending and twisting as well as tensile loading. It is reported that CNF’s experimental shear modulus is between 1.8 and 3.8 GPa , and the simulated one by all-atom model is 1.6 GPa . Although knowledge on tensile properties of CMF is increasing, the material’s behavior in torsion or bending deformation is quite different  . In addition, the actual CMF is tangled complicatedly to form a hierarchical structure. In natural material, through a lot of hierarchical structures inside one material, torsional force or deformation (i.e. chirality) is successfully transmitted . Thus, it is worthwhile investigating hierarchical transmission mechanism in cellulose materials. Therefore, it is important not only to understand only one CMF but also to elucidate the transmission mechanism of force and deformation between many CMFs.
In this research, we focus on the torsional characteristics of CMF and the transmission mechanism of shear stress and strain during torsional deformation is investigated in a hierarchical structure. Torsional deformation simulation of CMF by molecular dynamics (MD) method is performed to calculate shear modulus of elasticity. Dependence on temperature and rotation direction is also investigated. Then, we assume a state in which two CMFs are ideally arranged in parallel, and create a computational model of hierarchical structure. For that model, we also evaluate the dependence on temperature and discuss the bond strength and toughness occurring in the hierarchical structures. Furthermore, to clarify the transmission mechanism between components of the hierarchical structures, only one of CMFs is twisted, and we observe how shear stress and strain are transmitted to the other CMF.
2. Computational Methods and Models
For the conformation of cellulose nanofibers, we use the potential functions (force field) including stretch, bending, torsion, van der Waals interaction (Buckingham) and electrostatic (Coulomb) interaction   . The details of potential function forms, parameters and variables are omitted here and the readers can refer to the former studies    and our methodology  as well. In this research, we use the hybrid molecular model developed based on the method of united atom model for the repeating unit of cellulose molecule. In order to accelerate the calculation, hydrogen (H) atoms in the six-membered ring are integrated into the carbon (C) atom therein. On the other hand, the H atoms in the hydroxy (OH) group which will form a hydrogen bond are not integrated and preserve enough accuracy. Hydrogen bonds are reproduced and calculation speed is improved by this method. Repetitive unit of cellulose molecule is modelled as shown in Figure 1, and its parameters are shown in Table 1.
The cellulose Iβ type, which is a main component of plant-derived natural cellulose, is the present objective. The crystal structure of cellulose Iβ is monoclinic. Therefore, 10 or 15 molecules, using the repeat unit model of Figure 1, are polymerized in the y axis direction to prepare a molecular chain. Subsequently, molecular chains are arranged in a crystal structure around the central molecular chain in the xz plane, and a CMF of five layers is made as shown in Figure 2. The conditions of those CMFs are shown in Table 2. CMFs with degrees of polymerization of 10 and 15 are named Model 1 and Model 2, respectively. Generally, the actual degree of polymerization of CMF is 200 - 300. But, in this study, we assume simple twisting of linear-shaped material with almost circular cross-section. The shorter model is preferable to exclude any complicated mode of deformation other than simple twist and it will make the discussion simple. We can assume, when simple twist is applied, the cross-sectional shape is still axisymmetric around the rotation axis, so that all points at the same radius are displaced by the same amount along the circumference path and no axial displacement occurs. Thus, CMF is modelled with relatively small degree of polymerization, like 10 (Model 1) or 15 (Model 2) as stated above, in this study.
Figure 1. Repeating unit of cellulose molecule reproduced by Hybrid model.
Figure 2. CMF model of cellulose Iβ. (a) Model 1: degree of polymerization = 10; (b) Model 2: degree of polymerization = 15.
Table 1. Parameters of repeating unit of cellulose molecule in the hybrid model.
Table 2. Parameters of cellulose micro fibril (CMF).
3. Calculation Condition
3.1. Torsional Simulation of Single CMF
At first, structural relaxation to obtain thermal equilibration is performed on the Model 2 of CMF by using MD method. This calculation is performed for 0.1 ns from the initial state with temperature control at 10, 100, 200 and 300 K. After that relaxation calculation, as shown in Figure 3 (left), an atomic group at one end (yellow-colored atoms in the figure), which includes just one repetition unit of cellulose molecule, is spatially fixed in all directions (in x, y, and z directions), while a constant angular velocity around the y axis is given to the other end (blue-colored atoms) and a torsion simulation is performed. For the rotation direction, it is done for both plus and minus angular velocities. Those directions are schematically shown by arrows in Figure 3 (right). The given constant angular velocity is ω = +2.0 or −2.0 rad/ns, and the calculation time is for 0.3 ns. The temperature condition is 10, 100, 200 and 300 K. In this research, to evaluate the shear modulus, it is necessary to calculate the atomic shear stress occurring in the rotation direction of each atom in the CMF. Therefore, we assume that the CMF is placed in a system of cylindrical coordinate. Each atomic stress tensor components in the Cartesian coordinates should be converted into those in the cylindrical coordinates and the average value calculated is obtained for all the atoms for each step.
3.2. Hierarchical Structure Simulation Using Two CMFs
Hierarchical structure simulation with multiple CMFs is performed using the Model 1 of CMF. First, a hierarchical structure is created by combining two equivalent CMF structures which has already relaxed in the same way as Model 2 stated in previous subsection. Then, both tensile and torsion simulations are performed on the hierarchical structure. The tensile simulation model is shown in Figure 4. Two groups of five molecular chains located at each ends in x direction of the combined structure (yellow-colored in the figure) are given an opposite constant speed in the x direction, and the structure is stretched. The given constant velocities are ν = ±15 m/s, and the calculation time is for 0.1 ns. The temperature condition is 10, 100, 200, 300, 400 and 500 K.
Configuration of the torsion simulation model for the hierarchical structure is shown in Figure 5. The CMF on the left is named CMF1 and that on the right CMF2. Both CMF1 and CMF2 are fixed in all the directions along x, y, and z axes at one end (yellow-colored atoms in the figure), while only CMF 1 is given a constant angular velocity around the y axis at the other end (blue-colored atoms). That is, the other end of CMF2 is free of boundary constraint. The given constant angular velocity is ω = +20 rad/ns (it is done only in anticlockwise rotation, viewing from +y), and the calculation time is for 0.03 ns. The temperature condition is 10, 100, 200 and 300 K.
In addition, in this research of the hierarchical structure, “transmission ratio” of stress (force) or strain (deformation) between components (i.e. between CMF1 and CMF2) will be discussed. Here, the transmission ratio is defined as the ratio of the value of CMF2 to the value of CMF1 (since CMF1 is driven first). Suppose, for example, that the averaged shear stress of CMF1 and CMF2 at the final step of the loading are τ1 and τ2, then the transmission ratio of the stress between them R is expressed by R = τ2/τ1.
Figure 3. MD simulation model for torsion test (an example of 5-layer model in 10 K).
Figure 4. MD simulation model for the tensile simulation of hierarchical structure.
Figure 5. MD simulation model for torsion simulation of hierarchical structure.
4. Result and Discussion
4.1. Torsional Simulation of Single CMF
Shear stress-shear strain diagrams under positive and negative angular velocities are shown in Figure 6 and Figure 7. Table 3 shows the shear moduli calculated within the strain range from 0 to 0.05. Even when the strain is zero, the stress is not zero. This is because twist angle already occurs during structural relaxation anticlockwisely in viewing from +y axis and the residual stress exists inside. It is observed that, regardless of the twisting direction, the higher the temperature is, the lower stress value is obtained at the same strain. This indicates that, when the temperature rises, the CMF becomes soft and easily deformed. For any angular velocity condition, the shear modulus at the temperature less than 300 K is always smaller than that at 10 K. The value of shear modulus obtained here is about 1.4 GPa at 10 K, and about 0.8 - 1.0 GPa between 100 and 300 K. They are somewhat smaller than the experimental value, 1.8 - 3.8 GPa . The calculated value reported for CMF using MD with the all atom model is 1.6 GPa . This slight discrepancy is caused by inhomogeneous deformation of the whole system due to the influence of the fixed and velocity-constrained ends.
4.2. Hierarchical Structure Simulation Using Two CMFs
4.2.1. Tensile Simulation
The stress-strain diagram obtained by tensile simulation is shown in Figure 8. The maximum stress obtained for each temperature is shown in Table 4. The relationship between toughness and temperature is shown in Figure 9. Here, the toughness value is defined as the strain energy until the tensile strain reaches 0.3. The toughness in Figure 8 is shown by the strain energy per the initial volume of CMFs.
With respect to Figure 8, the stress increases with increasing strain up to the maximum stress, but then declines. Even when tensile deformation is applied, the original crystal shape of individual CMF is almost preserved. From this fact,
Figure 6. Shear stress-shear strain diagram for ω = +2.0 rad/ns.
Figure 7. Shear stress-shear strain diagram for ω = −2.0 rad/ns.
Figure 8. Stress-strain diagram for tensile simulation of hierarchical structure.
Figure 9. Relationship between toughness and temperature for tensile simulation.
Table 3. Shear modulus G [GPa].
Table 4. Maximum stress for tensile simulation obtained for each temperature.
the hydrogen bond between CMF surfaces tends to be locally stretched and broken there, more possibly than bonds inside the CMF crystal. As a result, a large reduction in stress occurs. Under the temperature between 100 and 400 K, there is no big difference in the maximum stress values. The reason why the maximum stress at 10 K is particularly low is that the intermolecular hydrogen bond between the CMF surfaces is weakly formed there, compared to other temperatures. It means that the deformation resistance between the CMFs has weakened during combining process of the hierarchical structure. At 500 K, no stress reduction is observed after reaching the maximum values. In this temperature, the molecules in the CMF are close to the thermal decomposition, so that the entire hierarchical structure uniformly deforms, without concentration on the hydrogen bond between CMFs. In Table 4, the maximum stress value at 300 K is 0.68 GPa. Using the all-atom model, the tensile strength in the direction of the molecular sheet inside the CMF crystal was estimated from 0.40 to 0.90 GPa . It is recognized that the hybrid model used in this research maintains important hydrogen bonds and accurately reflects the limit value of hydrogen bonding in the cellulose.
Since the manner of degradation of stress after reaching the maximum value is quite different between temperatures as shown in Figure 9, we will mention about the temperature dependence of the toughness value. Generally speaking, in crystals of small molecules, a solid state in low temperature changes to a liquid state in high temperature by increasing diffusion of the molecules inside the structure. On the other hand, in polymer crystals, atomic vibration in a molecule becomes more active as temperature increase, and so each single molecular chain moves in more flexible manner. However, the whole molecular chain is not able to diffuse or move in large distance due to entanglement and so on.
In CMF as well, the toughness also increases as the temperature rises, by virtue of flexibility of molecules. However, when the temperature is 400 K or more, the toughness value saturates at the same limit because the thermal decomposition of CMF into separated molecules has been reached in such high temperature.
4.2.2. Torsion Simulation
Transmission ratio of shear stress and shear strain from CMF1 to CMF2 at temperature from 10 to 300 K are shown in Figure 10. At 10 K, the shear strain is not transmitted at all, and the shear stress induces the counterpart to produce a stress in the opposite rotation direction (from rotational and fibrillar axis, it will be like a pair of mechanical gears transmitting torque). In this case, for very large deformation, the bond between surfaces of two CMFs is broken due to large separation and the CMF2 return to the original shape quickly. The transmission ratio of shear strain increases as the temperature rises. On the other hand, the transmission ratio of the shear stress decreases over 200 K. This indicates that the shear stress is not greatly transmitted at high temperature. Therefore, it is inferred that in a high temperature shear strain or deformation can be transmitted via small stress. Figure 11 shows the atomic arrangement for each time at 300 K and a schematic of torsion deformation in hierarchical structure. As can be seen in Figure 11, it is confirmed that CMF1 shows torsional deformation, while CMF2 behaves like bending so that deflection occurs in the negative z direction. It shows possibility that the rotation of one CMF (torsional moment) is directly converted to the bending moment of adjacent CMF.
Figure 10. Difference in transmission ratio for each temperature.
Figure 11. Snapshots of hierarchical structure in torsion simulation at 300 K and a schematic of torsion deformation in hierarchical structure.
In this study, the mechanical properties of cellulose microfibrils (CMF), as the basic component of cellulose nanofibers (CNF), are investigated. In particular, torsional deformation and its transmission mechanism between components in hierarchical structure are discussed based on the results obtained by MD simulations.
· The shear modulus is about 0.7 - 1.0 GPa and there is no clear temperature dependence between 100 and 300 K;
· The strength of the hierachical structure in the sheet direction does not much depend on the temperature, but its toughness estimated from strain energy until it is broken shows strong temperature dependency;
· Transmission of torsional deformation between components in the hierarchical structure is well observed but is dependent on temperature. The transmission ratio of mechanical values (ratio between averaged stress or strain occurring in two adjacent CMFs) increases as temperature rises, but that of shear stress tends to decrease above a certain temperature (in this model, above 200 K).
In our opinion, it is possible to design a new composite material by using the shear modulus of CMF which is difficult to evaluate by experiments but is able to be estimated by MD as done in our study. This study will lead us to a better understanding of the dynamic behavior of CMF in hierarchical structure.
This study is supported by the Kansai University Grant-in-Aid for progress of research in graduate course, 2017 (Apr.)-2018 (Mar.).
 Klemm, D., Kramer, F., Moritz, S., Lindström, T., Ankerfors, M., Gray, D. and Dorris, A. (2011) Nanocelluloses: A New Family of Nature-Based Materials. Angewandte Chemie International Edition, 50, 5438-5466.
 Moon, R.J., Martini, A., Nairn, J. Simonsen, J. and Youngblood, J. (2011) Cellulose Nanomaterials Review: Structure, Properties and Nanocomposites. Chemical Society Reviews, 40, 3941-3994.
 Lee, K.-Y., Aitomäki, Y., Berglund, L.A., Oksman, K. and Bismarck, A. (2014) On the Use of Nanocellulose as Reinforcement in Polymer Matrix Composites. Composites Science and Technology, 105, 15-27.
 Milanez, D.H., Morato do Amaral, R., Lopes de Faria, L.I. and Gregolin, J.A.R. (2013) Technological Indicators of Nanocellulose Advances Obtained from Data and Text Mining Applied to Patent Documents. Materials Research, 16, 635-641.
 Guhados, G., Wan, W. and Hutter, J.L. (2005) Measurement of the Elastic Modulus of Single Bacterial Cellulose Fibers Using Atomic Force Microscopy. Langmuir, 21, 6642-6646.
 Heiner, A.P. and Teleman, O. (1996) Comparison of the Interface between Water and Four Surfaces of Native Crystalline Cellulose by Molecular Dynamics Simulations. Pure and Applied Chemistry, 68, 2187-2192.
 Yui, T., Nishimura, S., Akiba, S. and Hayashi, S. (2006) Swelling Behavior of the Cellulose Iβ Crystal Models by Molecular Dynamics. Carbohydrate Research, 341, 2521-2530.
 Hadden, J.A., French, A.D. and Woods, R.J. (2014) Effect of Microfibril Twisting on Theoretical Powder Diffraction Patterns of Cellulose Iβ. Cellulose, 21, 879-884.
 Shklyaev, O.E., Kubicki, J.D., Watts, H.D. and Crespi, V.H. (2014) Constraints on Iβ Cellulose Twist from DFT Calculations of 13C NMR Chemical Shifts. Cellulose, 21, 3979-3991.
 Uto, T., Mawatari, S. and Yui, T. (2014) Theoretical Study of the Structural Stability of Molecular Chain Sheet Models of Cellulose Crystal Allomorphs. The Journal of Physical Chemistry B, 118, 9313-9321.
 Paavilainen, S., Róg, T. and Vattulainen, I. (2011) Analysis of Twisting of Cellulose Nanofibrils in Atomistic Molecular Dynamics Simulations. The Journal of Physical Chemistry B, 115, 3747-3755.
 Northolt, M.G., Boerstoel, H., Maatman, H., Huisman, R., Veurink, J. and Elzerman, H. (2001) The Structure and Properties of Cellulose Fibres Spun from an Anisotropic Phosphoric Acid Solution. Polymer, 42, 8249-8264.
 Zhao, Z., Shklyaev, O.E., Nili, A., Mohamed, M.N.A., Kubicki, J.D., Crespi, V.H. and Zhong, L. (2013) Cellulose Microfibril Twist, Mechanics, and Implication for Cellulose Biosynthesis. The Journal of Physical Chemistry A, 117, 2580-2589.
 Kannam, S.K., Oehme, D.P., Dobblin, M.S., Gidley, M.J., Bacic, A. and Downton, M.T. (2017) Hydrogen Bonds and Twist in Cellulose Microfibrils. Carbohydrate Polymers, 175, 433-439.
 Wang, J.-S., Wang, G., Feng, X.-Q., Kitamura, T., Kang, Y.-L., Yu, S.-W. and Qin, Q.-H. (2013) Hierarchical Chirality Transfer in the Growth of Towel Gourd Tendrils. Scientific Reports, 3, Article No. 3102.
 Neyertz, S., Pizzi, A., Merlin, A., Maigret, B., Brown, D. and Deglise, X. (2000) A New All-Atom Force Field for Crystalline Cellulose I. Journal of Applied Polymer Science, 78, 1939-1946.
 Neyertz, S. and Brown, D. (1995) A Computer Simulation Study of the Chain Configurations in Poly (Ethylene Oxide)-Homolog Melts. The Journal of Chemical Physics, 102, 9725-9735.
 Koyama, A., Yamamoto, T., Fukao, K. and Miyamoto, Y. (2001) Molecular Dynamics Studies on Local Ordering in Amorphous Polyethylene. The Journal of Chemical Physics, 115, 560-566.
 Saitoh, K., Ohno, H. and Matsuo, S. (2013) Structure and Mechanical Behavior of Cellulose Nanofiber and Micro-Fibrils by Molecular Dynamics Simulation. Soft Nanoscience Letters, 3, 58-67.