Mechanical Properties of Cellulose Nanocrystal (CNC) Bundles: Coarse-Grained Molecular Dynamic Simulation

: Cellulose nanocrystals (CNCs) is a promising biodegradable nanomaterial with outstanding physical, chemical, and mechanical properties for many applications. Although aligned CNCs can self-assemble into bundles, their mechanical performance is reduced by interfacial strength between CNCs and a twisted structure. In this paper, we employ developed coarse-grained (CG) molecular dynamic (MD) simulations to investigate the inﬂuence of twist and interface energy on the tensile performance of CNC bundles. CNC bundles of di ﬀ erent sizes (number of particles) are tested to also include the e ﬀ ect of size on mechanical performance. The e ﬀ ect of interfacial energy and twist on the mechanical performance shows that elastic modulus, strength, and toughness are more sensitive to twisted angle than interfacial energy. In addition, the e ﬀ ect of size on the bundle and twist on their mechanical performance revealed that both size and twist have a signiﬁcant e ﬀ ect on the results and can reduce the strength and elastic modulus by 75% as a results of covalent bond dissociation. In addition, a comparison of the broken regions for di ﬀ erent values of twist shows that by increasing the twist angle the crack propagates in multiple locations with a twisted shape.


Introduction
Cellulose is an abundant green polymer on earth mostly found in plant cell walls.Cellulose nanocrystals (CNCs) are crystalline structures that are formed by the stacking of many cellulose chains through hydrogen bonding [1].They have gained significant attention in recent years due to their outstanding physical, mechanical, and chemical properties such as high aspect ratio, high strength, high stiffness, and active surface for surface functionalization [1][2][3][4].Having a very active surface and high aspect ratio, aligned CNCs can easily assemble to CNC bundles [5] and be used in CNC reinforced composite materials.However, studies have shown that twisting of fibers can affect the interaction between fibers and reduce the overall performance of the bundle [6].This is particularly important for CNCs as many studies have shown that CNCs can easily twist [7][8][9].
In the past decade, substantial efforts have been placed on employing CNCs as neat cellulose-based materials [10] or as a reinforcement element in fiber nanocomposites [11][12][13].Many researchers employed experimental and theoretical approaches to characterize the mechanical properties of individual CNCs [14][15][16][17][18][19][20][21], bio-inspired CNCs [14][15][16][17][18][19][20][21], and CNC nanocomposite.For example, experimental and theoretical studies show that the Young modulus of a pristine CNC in the longitudinal direction is in the range of 110-200 GPa [14,21] and strength of a perfect CNC is in the range of 2-8 GPa [19,22].While in MD simulations the focus of the study is limited to the short length and small number of CNCs, some researchers have developed coarse-grained (CG) models to overcome this limitation and provide insights on the performance of long CNCs or the assembly of many CNCs [21].Many researchers developed CG models of CNCs with different coarsening details such as three beads [23][24][25][26] and one bead [27][28][29][30][31] per glucose ring for conformational [23,25,32], mechanical [24,26,31], and chemical properties [27,28,33].For example, a solvent free CG model of CNC was developed to understand the transition of crystalline to amorphous [28].Fan and Maranas, (2015), developed a unique CG model for conformational analysis of long fibrils [29].Wu et al. (2012), used the MARTINI CG model of CNC and studied the effect of the number of chains in a CNC on the twist angle and found that CNCs with more than 36 chains broke into two parts [25].A detailed CG model for mechanical and interfacial properties of CNC was developed to study the performance of bioinspired Bouligand and staggered structure of CNCs and it showed that the natural twist of CNCs improved the performance of materials [5].In another study, a spring-bead model for CNC was used to study the effect of length and surface energy on the performance and CNC nanopapers [34].
In this work, we employ the developed CG model for CNC [5,17,18,35,36] to understand the effect of twisting on the strength, elastic modulus, and toughness of CNC bundles.This study helps the researchers to understand the correlation between twist angle and bundle size, and the mechanical performance for optimum design of bundle size.Although to the best of our knowledge, this is the first study on the effect of twist on the CNC bundle, there have been many studies on the carbon nanotube (CNT) bundles and yarns [6,[37][38][39].For example, MD simulations of CNT bundle under axial tension showed that twisted CNT reduces the mechanical performance [37].A study on the twisted CNT yarn revealed that twisting could improve the interaction between tubes by increasing the interaction area [6].

The Coarse-Grained (CG) Model for Cellulose Nanocrystals (CNC)
The idea behind coarse-graining at the atomistic scale is to save computational time and energy by reducing the number of particles involved in the simulation while maintaining the critical physics of the original model (all-atom model).The CG model used in this study had been recently developed by Shishehbor and Zavattieri (2019), for mechanical and interfacial properties of CNCs [5].We mainly choose this model for the following reasons: (1) because it provides enough structural details of CNCs, (2) because of its tuneability of interfacial interaction between CNCs for surface functionalization, and (3) because it is available for the public through the nanoHUB toolkit [17,18].Here we explain the key features and equations of the model, however, for more details readers are referred to the original paper [5].
Figure 1a shows 3D and cross-section views of the molecular structure for a single CNC particle, formed by stacking of many cellulose chains (β-D glucopyranose) through hydrogen bonding (h-bonding) with reported unit cell parameters of a = 7.784 Å, b = 8.201 Å, c = 10.380Å, α = 90 • , β = 90 • , γ = 96.55• , at 293 K [1].While covalent bond (also called bonded interactions) is the dominant interaction in the chain direction, the interaction between chains is mainly through Van der Waals (VdW) and h-bonding (also called nonbonded interactions) [1,7].Previous studies suggest that the elastic modulus and strength of CNC in the chain direction ([001]) are mainly associated with covalent bonding, while VdW interactions and h-bonding have the main contributions for those properties in the lateral directions (110 and 1-10) [7,40].In addition, for the interaction between two or more CNC particles, VdW and h-bonding are the major governing interactions.In the CG model proposed for the mechanical and interfacial properties of CNCs [5], CG beads represent disaccharide molecules, and they are connected though different in order to model the mechanical properties of a single CNC (Figure 1b).Although adding more details (6 beads per cellobiose) would provide more information, the computational costs (number of particles and timestep) would increase dramatically.The current CG model would allow us to model big systems without sacrificing mechanical properties and interfacial properties as shown by [5].The interfacial properties, however, are modeled through nonlocal nonbonded potentials.Therefore, the total energy of the system (E) is defined through bonded (E b ) and nonbonded (E nb ) interactions (Equation ( 1)).
To accurately model the anisotropic mechanical properties of a single CNC in all directions, five different bonds were considered in the model and parametrized to fit the experimental and MD results.
For bonded interactions, Morse potential was considered to model dissociation of interactions, and therefore the total bonded energy terms for the model is as follow: where, D, r 0 , and α are the parameters of Morse potential.
To accurately model the anisotropic mechanical properties of a single CNC in all directions, five different bonds were considered in the model and parametrized to fit the experimental and MD results.For bonded interactions, Morse potential was considered to model dissociation of interactions, and therefore the total bonded energy terms for the model is as follow: where, D, r0, and α are the parameters of Morse potential.For nonbonded interactions Lennard-Jones (LJ) potential were considered as follow: where, σ is approximately the distance between two beads at equilibrium, ϵ is the energy term for LJ potential depth, and r is the distance between two beads.Figure 1c shows the assembly of 25 CNCs as a 5 × 5 CNC bundle where each particle is shown with a different color.From experimental observations, CNCs extracted from plants and trees, have a width of 2-5 nm (consisting of 16-36 chains) and a length of 50-500 nm [1].Therefore, the width and length of each particle in the model (2.5 nm width and 100 nm length) are in the range of the experimental observations.

Twisted Bundles of CNCs
The CG model for the bundle of CNCs studied here are bundles with a square cross-section of 2 × 2, 3 × 3, 4 × 4, and 5 × 5 CNC particles with 100 nm length.To study the effect of twist on the mechanical performance, bundles with different twisted angles (θ) of 0°, 90°, 180°, 270°, and 360° were generated.The twisted structure was generated by applying a transformation matrix to the original position of particles in the structure with θ = 0 as follow: For nonbonded interactions Lennard-Jones (LJ) potential were considered as follow: where, σ is approximately the distance between two beads at equilibrium, is the energy term for LJ potential depth, and r is the distance between two beads.Figure 1c shows the assembly of 25 CNCs as a 5 × 5 CNC bundle where each particle is shown with a different color.From experimental observations, CNCs extracted from plants and trees, have a width of 2-5 nm (consisting of 16-36 chains) and a length of 50-500 nm [1].Therefore, the width and length of each particle in the model (2.5 nm width and 100 nm length) are in the range of the experimental observations.

Twisted Bundles of CNCs
The CG model for the bundle of CNCs studied here are bundles with a square cross-section of 2 × 2, 3 × 3, 4 × 4, and 5 × 5 CNC particles with 100 nm length.To study the effect of twist on the mechanical performance, bundles with different twisted angles (θ) of 0 • , 90 • , 180 • , 270 • , and 360 • were generated.The twisted structure was generated by applying a transformation matrix to the original position of particles in the structure with θ = 0 as follow: where, x is the new position of particles, X is the original position of particles, and R θ is the transformation matrix around z axis by θ degree.Figure 2a shows a 5 × 5 bundle of CNCs with 0 • , 90 • , and 360 • twisted angle.It is worth mentioning that in this study we focused on the mechanical features of the CNC bundle rather than self-assembly and chemical aspects and we assumed that the structure had already twisted due to the shape and negative charge and then we studied the mechanical properties.In addition, the mechanical and interfacial properties for individual particles were in the range of experimental and MD results which shows that a negative charge does not have dramatic influence on the mechanical properties.For self-assembly simulations, where conformational analysis is more important than mechanical properties, CG beads could contain a negative charge (although this would increase the computational cost significantly), and therefore further analysis needs to be done for validation of the model.
mechanical properties.In addition, the mechanical and interfacial properties for individual particles were in the range of experimental and MD results which shows that a negative charge does not have dramatic influence on the mechanical properties.For self-assembly simulations, where conformational analysis is more important than mechanical properties, CG beads could contain a negative charge (although this would increase the computational cost significantly), and therefore further analysis needs to be done for validation of the model.

Simulation and Analysis Protocol
Initially, the bundles (shown in Figure 2a) are minimized by Hessian-free truncated Newton (HFTN) and conjugate gradient (CG) approaches and then equilibrated at 1 K temperature for three ns using a Nosé-Hoover thermostat and timestep of 5 fs.We used 1 K temperature and fixed the movement of the boundary atoms in the y-direction during equilibration to retain the initial twisted angle.We employed low temperature MD simulations to remove the effect of kinetic energy (thermal fluctuations) on the structure and mechanical properties.This approach had been used before by many authors [7,41,42] for carbon nanotube and CNC.After equilibration, the tensile displacement was at the speed of 0.00001 A/fs (1 m/s) applied on one side of the bundle (the moving boundary particles), while the particles on the opposite side (the fixed boundary particles) were being fixed and the thermostat particles were equilibrated as shown in Figure 2b.We examined different speeds in the range of 0.1-10 m/s and found a convergence of the results below 2 m/s with 10% error.To examine the effect of interfacial properties between CNC particles on the mechanical performance of CNC bundles, variations of ϵ (energy depth in LJ potential in Equation 3) from 1.25-10 kcal/mol were considered.In the original paper, the authors explained that the values of ϵ = 5.0 and ϵ = 10.0 showed natural twist and perfect interface between two CNCs, respectively [5].Therefore, here, the range of 1.25-10 kcal/mol was a suitable range to examine the effect of interfacial properties on the results.We employed a Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [43] MD package to carry out all the simulations and used OVITO package for all the visualization [44].

Simulation and Analysis Protocol
Initially, the bundles (shown in Figure 2a) are minimized by Hessian-free truncated Newton (HFTN) and conjugate gradient (CG) approaches and then equilibrated at 1 K temperature for three ns using a Nosé-Hoover thermostat and timestep of 5 fs.We used 1 K temperature and fixed the movement of the boundary atoms in the y-direction during equilibration to retain the initial twisted angle.We employed low temperature MD simulations to remove the effect of kinetic energy (thermal fluctuations) on the structure and mechanical properties.This approach had been used before by many authors [7,40,42] for carbon nanotube and CNC.After equilibration, the tensile displacement was at the speed of 0.00001 A/fs (1 m/s) applied on one side of the bundle (the moving boundary particles), while the particles on the opposite side (the fixed boundary particles) were being fixed and the thermostat particles were equilibrated as shown in Figure 2b.We examined different speeds in the range of 0.1-10 m/s and found a convergence of the results below 2 m/s with 10% error.To examine the effect of interfacial properties between CNC particles on the mechanical performance of CNC bundles, variations of (energy depth in LJ potential in Equation ( 3)) from 1.25-10 kcal/mol were considered.In the original paper, the authors explained that the values of = 5.0 and = 10.0 showed natural twist and perfect interface between two CNCs, respectively [5].Therefore, here, the range of 1.25-10 kcal/mol was a suitable range to examine the effect of interfacial properties on the results.We employed a Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [43] MD package to carry out all the simulations and used OVITO package for all the visualization [44].

The Effect of Interface Energy and Twist
After carrying out the simulations for each case study, the elastic modulus (fitted line to 0-0.01 strain in the stress-strain curves), strength (maximum force divided by the cross-section area), toughness (area below stress-strain curve), and failure strain (strain at the maximum force) are extracted from the stress-strain curves (Figure A1).The effect of interfacial energy ( ) and twist (θ) on the elastic modulus of a 2 × 2 bundle (Figure 3a), shows that elastic modulus is less sensitive to the , but significantly decreases as θ increases.For instance, for θ = 360, the elastic modulus (130 ± 10 GPa) drops almost 35% from its value for θ = 0 (200 ± 10 GPa).The results for the strength of 2 × 2 bundles shows more sensitivity to both and θ.While increasing the for θ = 0, 90, and 180 increases the strength, it decreases the strength for θ = 270 and 360.For example, for θ = 0, the strength increases from 4.5 GPa at = 25 kcal/mol to 6.0 GPa at = 10.0 kcal/mol (33% increase), while for θ = 360, the strength decreases from 3.5 GPa at = 1.25 kcal/mol to 2.5 GPa at = 10.0 kcal/mol (28% decrease).In addition, the result shows that θ has more effect than on the strength and, for the same , the strength decreases by the increase of θ.The values obtained here for elastic modulus (140-200 GPa) and strength (2.5-6 GPa) are in the range of the reported experimental and MD simulation results for elastic modulus (110-200 GPa) [14,21] and the strength of 2-8 GPa [19,22], respectively.The toughness results (Figure 3c) shows a similar trend observed in the strength.Toughness is sensitive to both and θ, and for θ = 0, 90, 180, and large values of (e.g., = 5.0 and 10.0) toughness increases by the increase of .However, for θ = 270 and 360, and large values of (e.g., = 5.0 and 10.0) toughness decreases.For small values of (e.g., = 5.0 and 10.0), the results for toughness fluctuates.In addition, as θ increases from 0 to 360, toughness decreases.For example, toughness decreases from 130 ± 10 Mj/m3 for θ = 0 to 60 ± 20 Mj/m3 for θ = 360.

The Effect of Interface Energy and Twist
After carrying out the simulations for each case study, the elastic modulus (fitted line to 0-0.01 strain in the stress-strain curves), strength (maximum force divided by the cross-section area), toughness (area below stress-strain curve), and failure strain (strain at the maximum force) are extracted from the stress-strain curves (Figure A1).The effect of interfacial energy (ϵ) and twist (θ) on the elastic modulus of a 2 × 2 bundle (Figure 3a), shows that elastic modulus is less sensitive to the ϵ, but significantly decreases as θ increases.For instance, for θ = 360, the elastic modulus (130 ± 10 GPa) drops almost 35% from its value for θ = 0 (200 ± 10 GPa).The results for the strength of 2 × 2 bundles shows more sensitivity to both ϵ and θ.While increasing the ϵ for θ = 0, 90, and 180 increases the strength, it decreases the strength for θ = 270 and 360.For example, for θ = 0, the strength increases from 4.5 GPa at ϵ = 25 kcal/mol to 6.0 GPa at ϵ = 10.0 kcal/mol (33% increase), while for θ = 360, the strength decreases from 3.5 GPa at ϵ = 1.25 kcal/mol to 2.5 GPa at ϵ = 10.0 kcal/mol (28% decrease).In addition, the result shows that θ has more effect than ϵ on the strength and, for the same ϵ, the strength decreases by the increase of θ.The values obtained here for elastic modulus (140-200 GPa) and strength (2.5-6 GPa) are in the range of the reported experimental and MD simulation results for elastic modulus (110-200 GPa) [14,21] and the strength of 2-8 GPa [19,22], respectively.The toughness results (Figure 3c) shows a similar trend observed in the strength.Toughness is sensitive to both ϵ and θ, and for θ = 0, 90, 180, and large values of ϵ (e.g., ϵ = 5.0 and 10.0) toughness increases by the increase of ϵ.However, for θ = 270 and 360, and large values of ϵ (e.g., ϵ = 5.0 and 10.0) toughness decreases.For small values of ϵ (e.g., ϵ = 5.0 and 10.0), the results for toughness fluctuates.In addition, as θ increases from 0 to 360, toughness decreases.For example, toughness decreases from 130 ± 10 Mj/m3 for θ = 0 to 60 ± 20 Mj/m3 for θ = 360.Figure 3d shows the results for failure strain.The results do not show a specific trend for failure strain by changing θ and , and the values are in the range of 0.03 ± 0.01 for all cases.The effect of θ and for a 5 × 5 bundle is also discussed here and compared with a 2 × 2 bundle.However, the effect of size and θ are discussed in more details in Section 3.2.Similar to 2 × 2 bundles, the effect of and θ on the elastic modulus, strength and toughness of 5 × 5 bundles (solid lines in Figure 4), shows less sensitivity to the .However, mechanical performance significantly decreases as the twist angle increases.
Figure 3d shows the results for failure strain.The results do not show a specific trend for failure strain by changing θ and ϵ, and the values are in the range of 0.03 ± 0.01 for all cases.The effect of θ and ϵ for a 5 × 5 bundle is also discussed here and compared with a 2 × 2 bundle.However, the effect of size and θ are discussed in more details in Section 3.2.Similar to 2 × 2 bundles, the effect of ϵ and θ on the elastic modulus, strength and toughness of 5 × 5 bundles (solid lines in Figure 4), shows less sensitivity to the ϵ.However, mechanical performance significantly decreases as the twist angle increases.Comparing the results of 2 × 2 (dashed lines) and 5 × 5 (solid lines) bundles shows that increasing the size of the bundle significantly decreases the mechanical performance.For example, for θ = 0, the elastic modulus decreases by 67% from 150 ± 10 GPa for the 2 × 2 bundles to 50 ± 20 GPa for the 5 × 5 bundles.In addition, the strength decreases by 75% from 3.0 ± 0.5 GPa for the 2 x 2 bundles to 0.75 ± 0.25 GPa for the 5 × 5 bundles.The toughness also decreases on average by 50% from 60 ± 15 Mj/m3 for the 2 × 2 bundles to 30 ± 10 Mj/m3 for the 5 × 5 bundles.There is less trend in the results for failure strain, however, for most θ and ϵ the failure strain decreases by 40%.
Figure 5 demonstrates the effect of θ on the failure mechanism of 5 × 5 bundles and for ϵ = 5.0.The comparison of the broken regions for different values of θ shows that by increasing θ the crack propagates in multiple locations with a twisted shape, while for θ = 0 the crack occurs in one plane.Comparing the results of 2 × 2 (dashed lines) and 5 × 5 (solid lines) bundles shows that increasing the size of the bundle significantly decreases the mechanical performance.For example, for θ = 0, the elastic modulus decreases by 67% from 150 ± 10 GPa for the 2 × 2 bundles to 50 ± 20 GPa for the 5 × 5 bundles.In addition, the strength decreases by 75% from 3.0 ± 0.5 GPa for the 2 × 2 bundles to 0.75 ± 0.25 GPa for the 5 × 5 bundles.The toughness also decreases on average by 50% from 60 ± 15 Mj/m 3 for the 2 × 2 bundles to 30 ± 10 Mj/m 3 for the 5 × 5 bundles.There is less trend in the results for failure strain, however, for most θ and the failure strain decreases by 40%.
Figure 5 demonstrates the effect of θ on the failure mechanism of 5 × 5 bundles and for = 5.0.The comparison of the broken regions for different values of θ shows that by increasing θ the crack propagates in multiple locations with a twisted shape, while for θ = 0 the crack occurs in one plane.The reduction in mechanical performance with respect to the increase of θ can be explained by bond dissociation as a result of twist.When a bundle is twisted, the CNC particles experience tension (length increases) and reduction in stiffness and load carrying capacity due to bond dissociation.Consequently, an increase in θ decreases mechanical performance.Although increasing the ϵ in some cases compensate the reduction in strength and stiffness due to twist, for large values of ϵ, CNCs break at the lateral direction (h-bond direction) as opposed to chain direction (covalent bond), and therefore increasing ϵ cannot fully recover the strength and stiffness of the bundle.

The Effect of Bundle Size and Twist
In this section, the influence of bundle size and twist on the mechanical performance of CNC bundles are evaluated.As discussed in Section 2.2, bundles with a square cross-section of 2 × 2, 3 × 3, 4 × 4, and 5 × 5 CNCs under a range of different values for ϵ and θ are examined here.Since we have examined the effect of ϵ and θ in Section 3.1.for bundles of 2 × 2 and 5 × 5, here, the effect of size and θ are evaluated for ϵ = 1.25 and 5.0 kcal/mol.The effect of size and θ on the elastic modulus for ϵ = 1.25 kcal/mol (Figure 6a), indicates that by increasing both size and θ, the elastic modulus significantly decreases, and the sensitivity of the results to the size of the bundle is more apparent as θ increases.For example, for θ = 0, the elastic modulus is 200 ± 10 GPa and 180 ± 15 GPa for 2 × 2 and 5 × 5 bundles, respectively (10% decrease); while for θ = 360, the elastic modulus is 140 ± 15 GPa and 30 ± 10 GPa, respectively (75% decrease).The reduction in mechanical performance with respect to the increase of θ can be explained by bond dissociation as a result of twist.When a bundle is twisted, the CNC particles experience tension (length increases) and reduction in stiffness and load carrying capacity due to bond dissociation.Consequently, an increase in θ decreases mechanical performance.Although increasing the in some cases compensate the reduction in strength and stiffness due to twist, for large values of , CNCs break at the lateral direction (h-bond direction) as opposed to chain direction (covalent bond), and therefore increasing cannot fully recover the strength and stiffness of the bundle.

The Effect of Bundle Size and Twist
In this section, the influence of bundle size and twist on the mechanical performance of CNC bundles are evaluated.As discussed in Section 2.2, bundles with a square cross-section of 2 × 2, 3 × 3, 4 × 4, and 5 × 5 CNCs under a range of different values for and θ are examined here.Since we have examined the effect of and θ in Section 3.1.for bundles of 2 × 2 and 5 × 5, here, the effect of size and θ are evaluated for = 1.25 and 5.0 kcal/mol.The effect of size and θ on the elastic modulus for = 1.25 kcal/mol (Figure 6a), indicates that by increasing both size and θ, the elastic modulus significantly decreases, and the sensitivity of the results to the size of the bundle is more apparent as θ increases.For example, for θ = 0, the elastic modulus is 200 ± 10 GPa and 180 ± 15 GPa for 2 × 2 and 5 × 5 bundles, respectively (10% decrease); while for θ = 360, the elastic modulus is 140 ± 15 GPa and 30 ± 10 GPa, respectively (75% decrease).The results for the strength (Figure 6b) also show sensitivity to both size and θ, and the strength decreases more rapidly for bigger bundles and larger θ.For example, for θ = 0, the strength is 4.75 ± 0.25 GPa and 3.5 ± 0.25 GPa for 2 × 2 and 5 × 5 bundles, respectively (25% decrease), while for θ = 360, the strength is 3.5 ± 0.5 GPa and 0.75 ± 0.25 GPa, respectively (75% decrease).The values obtained here for elastic modulus (175-200 GPa) and strength (3.5-5.0GPa) for θ = 0 are in the range of the reported experimental and MD simulation results for elastic modulus (110-200 GPa) [14,21] and the strength of 2-8 GPa [19,22], respectively.However, for 4 × 4 and 5 × 5 bundles and θ > 270, the effect of size and twist leads to a decrease of elastic modulus and strength below the minimum margins for pristine CNCs (110 GPa and 2 GPa for elastic modulus and strength, respectively).Toughness is also sensitive to both size and θ (Figure 6c), and it decreases with an increase of size and θ.However, the size dependency is clearer for larger θ (θ > 180) and for θ < 180 there are regions with higher toughness for a 5 × 5 bundle than 3 × 3 and 4 × 4 bundles (for example θ = 90).For θ = 360 (large θ), toughness drops approximately 60% from 65 ± 5 (Mj/m3) for 2 × 2 bundles to 20 ± 5 (Mj/m3) for 5 × 5 bundles.The results for failure strain (Figure 6d) indicate that the failure strain varies in the range of 0.025-0.03for a different size and θ.However, the sensitivity of failure strain to θ is higher than the sensitivity to the size.The results for the strength (Figure 6b) also show sensitivity to both size and θ, and the strength decreases more rapidly for bigger bundles and larger θ.For example, for θ = 0, the strength is 4.75 ± 0.25 GPa and 3.5 ± 0.25 GPa for 2 × 2 and 5 × 5 bundles, respectively (25% decrease), while for θ = 360, the strength is 3.5 ± 0.5 GPa and 0.75 ± 0.25 GPa, respectively (75% decrease).The values obtained here for elastic modulus (175-200 GPa) and strength (3.5-5.0GPa) for θ = 0 are in the range of the reported experimental and MD simulation results for elastic modulus (110-200 GPa) [14,21] and the strength of 2-8 GPa [19,22], respectively.However, for 4 × 4 and 5 × 5 bundles and θ > 270, the effect of size and twist leads to a decrease of elastic modulus and strength below the minimum margins for pristine CNCs (110 GPa and 2 GPa for elastic modulus and strength, respectively).Toughness is also sensitive to both size and θ (Figure 6c), and it decreases with an increase of size and θ.However, the size dependency is clearer for larger θ (θ > 180) and for θ < 180 there are regions with higher toughness for a 5 × 5 bundle than 3 × 3 and 4 × 4 bundles (for example θ = 90).For θ = 360 (large θ), toughness drops approximately 60% from 65 ± 5 (Mj/m3) for 2 × 2 bundles to 20 ± 5 (Mj/m3) for 5 × 5 bundles.The results for failure strain (Figure 6d) indicate that the failure strain varies in the range of 0.025-0.03for a different size and θ.However, the sensitivity of failure strain to θ is higher than the sensitivity to the size.
The effect of size and θ on the mechanical performance for = 5.0 kcal/mol are also demonstrated in Figure 7 (solid lines) and compared with the results for = 1.25 kcal/mol (dashed lines in Figure 7 and solid lines in Figure 6).The results show the same trend for both = 1.25 and 5.0 kcal/mol, indicating that bundle size and θ have more influence than .By increasing from 1.25 to 5.0 kcal/mol, the most affected results are observed in the strength, toughness, and failure strain of 2 × 2 bundles, and for other bundle sizes the change is not significant (less than 20%).Figure 8a shows the effect of size and θ on the failure mechanism for = 5.0.The comparison of the broken regions for both 2 × 2 and 5 × 5 bundles shows that crack regions move from a plane crack at the boundary to a twisted crack in the middle of the bundle.The effect of size and θ on the mechanical performance for ϵ = 5.0 kcal/mol are also demonstrated in Figure 7 (solid lines) and compared with the results for ϵ = 1.25 kcal/mol (dashed lines in Figure 7 and solid lines in Figure 6).The results show the same trend for both ϵ = 1.25 and 5.0 kcal/mol, indicating that bundle size and θ have more influence than ϵ.By increasing ϵ from 1.25 to 5.0 kcal/mol, the most affected results are observed in the strength, toughness, and failure strain of 2 × 2 bundles, and for other bundle sizes the change is not significant (less than 20%).Figure 8a shows the effect of size and θ on the failure mechanism for ϵ = 5.0.The comparison of the broken regions for both 2 × 2 and 5 × 5 bundles shows that crack regions move from a plane crack at the boundary to a twisted crack in the middle of the bundle.The reduction in mechanical performance with respect to the increase of size can be explained by bond dissociation as a result of twist.For a constant twisted angle, increasing the bundle size puts CNC particles under more tension (length of particles increases as a result of twist, especially for boundary particles) and there is a reduction in stiffness and load carrying capacity due to bond dissociation.

Figure 1 .
Figure 1.(a) 3D and cross-section view of the atomistic structure of cellulose nanocrystals (CNC) [1,40].Pink spheres denote oxygen atoms, blue spheres represent carbon atoms, and white spheres are hydrogen atoms.(b) The coarse-grained (CG) model of CNC, where each bead represents a disaccharide and is connected with 5 different internal bonds is shown with different colors.(c) A 3D and cross-section view of a 5 × 5 CNC bundle, where different particles are shown with a different color.

Figure 1 .
Figure 1.(a) 3D and cross-section view of the atomistic structure of cellulose nanocrystals (CNC) [1,41].Pink spheres denote oxygen atoms, blue spheres represent carbon atoms, and white spheres are hydrogen atoms.(b) The coarse-grained (CG) model of CNC, where each bead represents a disaccharide and is connected with 5 different internal bonds is shown with different colors.(c) A 3D and cross-section view of a 5 × 5 CNC bundle, where different particles are shown with a different color.

Figure 2 .
Figure 2. 5 × 5 bundles of CNC with different twisted angle (θ).(a) θ = 0°, 90°, and 360°; (b) side view of a 360° rotated bundle.Tensile displacement is applied on the boundary particles of one side while the boundary particles on the other side are fixed.

Figure 2 .
Figure 2. 5 × 5 bundles of CNC with different twisted angle (θ).(a) θ = 0 • , 90 • , and 360 • ; (b) side view of a 360 • rotated bundle.Tensile displacement is applied on the boundary particles of one side while the boundary particles on the other side are fixed.

Figure 7 .
Figure 7. Mechanical performance of twisted bundle of CNCs for different bundle size shown for  = 5.0 kcal/mol and  = 1.25 kcal/mol with solid and dashed lines, respectively.(a) elastic Modulus, (b) strength, (c) toughness, and (d) failure strain.