Ultrahigh Ballistic Resistance of Twisted Bilayer Graphene

: Graphene is a good candidate for protective material owing to its extremely high stiffness and high strength-to-weight ratio. However, the impact performance of twisted bilayer graphene is still obscure. Herein we have investigated the ballistic resistance capacity of twisted bilayer graphene compared to that of AA-stacked bilayer graphene using molecular dynamic simulations. The energy propagation processes are identical, while the ballistic resistance capacity of the twisted bilayer graphene is almost two times larger than the AA-bilayer graphene. The enhanced capacity of the twisted bilayer graphene is assumed to be caused by the mismatch between the two sheets of graphene, which results in earlier fracture of the ﬁrst graphene layer and reduces the possibility of penetration.


Introduction
Materials for ballistic protection have been an important research topic throughout human history [1][2][3][4][5][6]. In order to stop high-speed projectiles, the material should be able to absorb large amounts of kinetic energy. Most current body armor is mostly made of lightweight fibrous materials exhibiting high stiffness and a high strength-to-weight ratio [2]. The fibrous materials can disperse the impact energy rapidly to the area far from the impact center through conic deformation and elastic wave. Considering that graphene is lightweight and has extremely high stiffness and strength [7][8][9][10][11][12][13][14][15], it seems to be a natural protective material.
Recently, Lee et al. demonstrated experimentally that the specific penetration energy of multilayer graphene is around 10 times larger than literature values for macroscopic steel sheets, which means graphene can be an extraordinary candidate of armor material [16][17][18][19]. The reason is attributed to the excellent impact-energy-delocalization property of graphene arising from its superior in-plane speed of sound, high strength, stiffness, and structural anisotropy [17].
Bilayer graphene (BLG) also illustrates advanced physical properties [20][21][22][23][24][25][26]. When the two sheets of bilayer graphene are twisted relative to each other by a certain angle, which forms the called twisted bilayer graphene (tBLG) [27][28][29], the intrinsic strength along transverse direction can be apparently increased with a slight change in the in-plane properties [30]. Consequently, tBLG may be harder to be penetrated than bilayer graphene without losing the superiority in the energy-delocalization property. This indicates that the tBLG should be able to absorb more energy during impact and thus possesses better ballistic resistance capacity [31].
A few efforts have addressed the ballistic properties of graphene membranes. Haque et al. analyzed the penetrating impact of 1-3 layers of graphene by molecular dynamics (MD) simulations [32]. Similar petal-shaped crack orifices were observed in the LIPIT. Meguid et al. is concerned with the ballistic behavior of multilayer graphene polymer composite [33]. They revealed that the ballistic impact resistance of polyethylene by Crystals 2021, 11, 206 2 of 9 over eight-fold with a single graphene layer membrane. Azevedo et al. analyzed the ballistic performance of a carbon-based material named penta-graphene [34]. Their results show that the fracture pattern is more spherical. Zhang et al. developed the transverse impact responses of three kinds of bilayers (h-BN/h-BN, graphene/graphene, and hybrid graphene/h-BN) twisted by molecular dynamics simulations [35]. Despite these efforts, the ballistic resistance performance of twisted bilayer graphene is still not fully understood. Motivated by the outstanding electronic properties of twisted bilayer graphene, we aim to study their impact resistance behaviors. By means of molecular dynamics simulations, we found that tBLG with a certain twist angle is a better ballistic protection material than BLG. The responses of tBLG and the referenced BLG under transverse impact loading and the mechanism for the enhanced capacity are investigated using molecular dynamic simulations [36].

Model and Molecular Dynamics Simulations
Molecular dynamics simulations were conducted using the large-scale atomic/molecular massively parallel simulator (LAMMPS) [37]. Through time and ensemble average of motion trajectory, the dynamic properties of the system are obtained. It can model extensive systems using a variety of interatomic potentials and boundary conditions [38,39]. A unit cell of tBLG was prepared according to Uchida's method [40]. Then it was further extended to a scale of around 80 nm to avoid the influence of the reflected elastic waves. The projectile was a diamond ball with a radius r = 3.57 nm. A periodic boundary condition was applied to the in-plane (x and y) directions, and the out-of-plane directions (z) were fixed at 6 nm. All simulations of impacts were conducted under micro-canonical (NVE) ensemble to maintain energy conservation in the entire system. Before the simulation, the system was fully relaxed in an isothermal-isobaric ensemble (NPT) for 40 ps at a temperature of 300 K and pressure of 1 atmosphere for both atomistic positions and simulation box dimensions in the in-plane directions to minimize the residue stresses. The time-step was set at 0.001 ps. For both BLG and tBLG samples, the same impacting processed were simulated twenty times, respectively, with different relaxation times.
The interaction of C-C bonds in a diamond ball and graphene layers is described by the second-generation reactive empirical bond order (REBO) potential, which can precisely describe the in-plane elastic properties of graphene and diamond [41][42][43]. The Lennard-Jones (LJ) potential was used for the rest of the van der Waals interactions, with parameters of 0.023049 eV and 2.852 Å for the interaction between diamond and graphene and 0.00284 eV and 3.4 Å for that between the two graphene layers [44]. This is a widely used potential for describing the interaction between two graphene layers [44]. Then the projectile was fired toward the center of samples transversely with an initial speed of 2.2 km/s. The simulation temperature was maintained at 300 K using a Nose-Hoover thermostat [45]. Other parameters were set the same for all the simulations. The penetration analyses were performed using OVITO [46].

Calculation of Perimeter of Crack Orifice
Common damage caused by a pass-through projectile during an impact test was an irregular left-over crack orifice. We introduced the quantity of perimeter of such an orifice as a descriptor of the damage caused by the applied projectile. We assumed the perimeter of one single atom was 1 (as a unit). For graphene lattice, there were three nearest atoms for each atom. Whenever an atom lost the nearest atom because of fracture or other reasons, it contributed 1/3 to the perimeter. The atoms with three nearest neighbor atoms contributed zero to the perimeter. First, we counted the nearest atoms for each atom within the impact center. Then we added up their corresponding contributions to the perimeter to had the quantity.

Atomistic Structures
The atomic configurations of AA-stacked bilayer graphene and twisted bilayer graphene tBLG in a unit cell are illustrated in Figure 1a,b, respectively, after full relaxation. Our results suggest that AA-stacked bilayer graphene is thermodynamically stable in our NPT ensemble at 300 K and 1 atmosphere pressure under the applied boundary conditions. Crystals 2019, 9, x FOR PEER REVIEW 3 of 9

Atomistic Structures
The atomic configurations of AA-stacked bilayer graphene and twisted bilayer graphene tBLG in a unit cell are illustrated in Figure 1a,b, respectively, after full relaxation. Our results suggest that AA-stacked bilayer graphene is thermodynamically stable in our NPT ensemble at 300 K and 1 atmosphere pressure under the applied boundary conditions.

Angle Effect on the Impact
We examined the angle effect on the impact performance of twisted bilayer graphene. The impact processes were simulated on the samples with 9 different twist angles , including 0.0°, 3.89°, 7.34°, 9.43°, 13.17°, 18.73°, 21.79°, 25.04° and 30.0°. The bullet was simulated by a diamond ball with a radius of 3.57 nm. We invested the impact using several impact speeds. The ballistic limit fell between 2.0 km/s and 2.5 km/s, so the initial velocities of the bullet were chosen within this range.
When the bullet velocity was 2.0 km/s, none of the samples were penetrated. When it came to 2.2 km/s, the tBLG samples with a twist angle of 7.34°, 18.73°, 21.79° and 25.04° were able to stop the bullet. However, the BLG ( = 0.0°) and other tBLG samples ( = 3.89°,9.43°, 13.17°, and 30.0°) failed to stop the bullet. When the bullet velocity reached 2.5 km/s, all of the samples were penetrated.

Penetration Comparison
To investigate the different performance on the ballistic resistance between BLG and tBLG, the bullet velocity of 2.2 km/s was fixed in all the following simulations. Among those non-penetrated samples in the bullet velocity of 2.2 km/s, tBLG with a twist angle 21.79° was selected to compare with the referenced BLG. A movie of the impacting process can be seen in the supplementary information.
After the ball impacted on the sample, an elastic wave radially propagated at highspeed, and a conic deformation (cone wave) followed afterward in good agreement with the experiment and previous simulation results for graphene [17,19]. To eliminate accidental errors caused by certain configurations of atoms, the simulations of impacting processed for BLG and tBLG were repeated twenty times, respectively. In each repeated simulation, the system was equilibrated with a different number of time-steps to obtain a different atomic configuration with a slight change in the position and velocity. Statistically, the ballistic resistance capacity of tBLG was 187.5% better than the referenced BLGthe possibility for the resistance of penetration was 40% for BLG and was 75% for tBLG.
According to the penetration situation, all the simulations were classified into four cases for further comparison, namely penetrated cases for BLG, non-penetrated cases for BLG, penetrated cases for tBLG and non-penetrated cases for tBLG. In all the cases, the first graphene layer was fractured. Setting the time zero point as when the bullet was 10 Å away from the first graphene layer, the average fracture of the first graphene layer

Angle Effect on the Impact
We examined the angle effect on the impact performance of twisted bilayer graphene. The impact processes were simulated on the samples with 9 different twist angles θ, 25.04 • and 30.0 • . The bullet was simulated by a diamond ball with a radius of 3.57 nm. We invested the impact using several impact speeds. The ballistic limit fell between 2.0 km/s and 2.5 km/s, so the initial velocities of the bullet were chosen within this range.
When the bullet velocity was 2.0 km/s, none of the samples were penetrated. When it came to 2.2 km/s, the tBLG samples with a twist angle of 7.

Penetration Comparison
To investigate the different performance on the ballistic resistance between BLG and tBLG, the bullet velocity of 2.2 km/s was fixed in all the following simulations. Among those non-penetrated samples in the bullet velocity of 2.2 km/s, tBLG with a twist angle θ = 21.79 • was selected to compare with the referenced BLG. A movie of the impacting process can be seen in the supplementary information.
After the ball impacted on the sample, an elastic wave radially propagated at highspeed, and a conic deformation (cone wave) followed afterward in good agreement with the experiment and previous simulation results for graphene [17,19]. To eliminate accidental errors caused by certain configurations of atoms, the simulations of impacting processed for BLG and tBLG were repeated twenty times, respectively. In each repeated simulation, the system was equilibrated with a different number of time-steps to obtain a different atomic configuration with a slight change in the position and velocity. Statistically, the ballistic resistance capacity of tBLG was 187.5% better than the referenced BLG-the possibility for the resistance of penetration was 40% for BLG and was 75% for tBLG.
According to the penetration situation, all the simulations were classified into four cases for further comparison, namely penetrated cases for BLG, non-penetrated cases for BLG, penetrated cases for tBLG and non-penetrated cases for tBLG. In all the cases, the first graphene layer was fractured. Setting the time zero point as when the bullet was 10 Å away from the first graphene layer, the average fracture of the first graphene layer started Crystals 2021, 11, 206 4 of 9 at 2.08 ps for tBLG and 2.39 ps for BLG, showing little difference between penetrated and non-penetrated cases. It appeared that the first graphene layer of tBLG could actually provide less protection to the second layer than BLG. Correspondingly, in these penetrated cases, the average fracture of the second graphene layer in tBLG started earlier than BLG, which was 3.28 ps for tBLG and 3.51 ps for BLG. However, this seemed contradictory with the earlier fracture of the first graphene layer and the lower penetration probability for tBLG, which needs further investigation.

Energy Propagation Process
The advantage of graphene serving as armor material is its excellent performance in energy delocalization [17,19]. Therefore, we were particularly interested in the energy propagation process in tBLG and BLG. Figure 2a is the energy absorbed by the first graphene layer in both penetrated and non-penetrated cases for BLG and tBLG, respectively. In the end, the first graphene layer in tBLG absorbed less energy, owing to its earlier fracture. However, there was little difference between BLG and tBLG before the fracture of the first graphene layer as well as between the penetrated and non-penetrated cases. It seemed strange that the tBLG performed worse in terms of energy absorption yet possessed a lower penetration probability. Figure 2b is the energy and out-of-plane displacement profile of the second graphene layer at t = 3.0 ps. The profile is illustrated along the line that goes through the impact center parallel to the armchair direction of the second graphene layer. The comparison of the energy curves and displacement curves in Figure 2b indicates that the energy was mostly transferred through a cone wave rather than an elastic wave. The propagation speeds of cone wave, which also represents the dissipation speeds of energy along a longitudinal direction, for BLG and tBLG, in both penetrated and nonpenetrated cases, were nearly identical. The energy profile along the zigzag direction can be seen in Figure S2 in the supplementary material. Moreover, the phonon density of state also showed a significant similarity between BLG and tBLG, which means an equivalent energy dissipation speed in terms of atomic vibration between BLG and tBLG. Above all, tBLG showed neither advantage nor a disadvantage in energy delocalization. Namely, by twisting a certain angle between two sheets of graphene, it did not lose its superiority in the energy-delocalization property.
Crystals 2019, 9, x FOR PEER REVIEW 4 of 9 started at 2.08 ps for tBLG and 2.39 ps for BLG, showing little difference between penetrated and non-penetrated cases. It appeared that the first graphene layer of tBLG could actually provide less protection to the second layer than BLG. Correspondingly, in these penetrated cases, the average fracture of the second graphene layer in tBLG started earlier than BLG, which was 3.28 ps for tBLG and 3.51 ps for BLG. However, this seemed contradictory with the earlier fracture of the first graphene layer and the lower penetration probability for tBLG, which needs further investigation.

Energy Propagation Process
The advantage of graphene serving as armor material is its excellent performance in energy delocalization [17,19]. Therefore, we were particularly interested in the energy propagation process in tBLG and BLG. Figure 2a is the energy absorbed by the first graphene layer in both penetrated and non-penetrated cases for BLG and tBLG, respectively. In the end, the first graphene layer in tBLG absorbed less energy, owing to its earlier fracture. However, there was little difference between BLG and tBLG before the fracture of the first graphene layer as well as between the penetrated and non-penetrated cases. It seemed strange that the tBLG performed worse in terms of energy absorption yet possessed a lower penetration probability. Figure 2b is the energy and out-of-plane displacement profile of the second graphene layer at t = 3.0 ps. The profile is illustrated along the line that goes through the impact center parallel to the armchair direction of the second graphene layer. The comparison of the energy curves and displacement curves in Figure  2b indicates that the energy was mostly transferred through a cone wave rather than an elastic wave. The propagation speeds of cone wave, which also represents the dissipation speeds of energy along a longitudinal direction, for BLG and tBLG, in both penetrated and non-penetrated cases, were nearly identical. The energy profile along the zigzag direction can be seen in Figure S2 in the supplementary material. Moreover, the phonon density of state also showed a significant similarity between BLG and tBLG, which means an equivalent energy dissipation speed in terms of atomic vibration between BLG and tBLG. Above all, tBLG showed neither advantage nor a disadvantage in energy delocalization. Namely, by twisting a certain angle between two sheets of graphene, it did not lose its superiority in the energy-delocalization property.

Crack Orifice Patterns
In order to get insight into the fracture mechanism, the evolution of the crack-orifice patterns during the impact process is demonstrated in Figure 3. Normally three to six cracks were initiated near the impact center and propagate outwards along the zigzag direction of graphene layers, creating the same number of triangular-shaped petals. The

Crack Orifice Patterns
In order to get insight into the fracture mechanism, the evolution of the crack-orifice patterns during the impact process is demonstrated in Figure 3. Normally three to six cracks were initiated near the impact center and propagate outwards along the zigzag direction of graphene layers, creating the same number of triangular-shaped petals. The damaged area spread quickly and finally became much larger than the impact area (outlined in dashed line in Figure 3). Such behavior agrees well with the experimental results [17]. Crystals 2019, 9, x FOR PEER REVIEW 5 of 9 damaged area spread quickly and finally became much larger than the impact area (outlined in dashed line in Figure 3). Such behavior agrees well with the experimental results [17]. Due to the high-speed impact of the bullet, three to six cracks were initiated near the impact center. The huge kinetic energy destroyed the graphene bond in the center of the impact, then due to the damaged area spread quickly, some atoms near the edge of petals were divorced from the petals, forming some fragments in the damaged area. Those fragments trapped in the impact region owing to the high velocity of the bullet were further pushed by the bullet and result in high-stress concentration points in the second layer, as shown in Figure S4a. It should be noticed that the cracks in the second layer always initiated from the position where some fragments of the fractured first layer existed. The crack-orifice patterns of all the cases simulated in this work are presented in Figure S4b.

Chaos Degree
In penetrated cases, there were more fragments trapped in the impact region on the first layer than those in non-penetrated cases. The fragments were assumed to play an important role in the fracture of the second layer. Here we use the term "chaos degree", defined as the total in-plane perimeter of atoms within the impact region (See Methods section), to evaluate the influence of fragments.
Considering that both cracks and fragments can be expressed by the in-plane perimeter of atoms, chaos degree could be helpful to explore the failure mechanism. The results are shown in Figure 4, and the data in each time point were averaged over all the simulations that belong to the same case. The dramatic increase of the chaos degree at around 2 ps was caused by the fracture of the first layer. The subsequent decrease was caused by the petals and fragments moving out from the impact region. Because of the earlier fracture of the first layer, the increase and decrease of the chaos degree for tBLG were both earlier than that for BLG.
For both tBLG and BLG, the peak chaos degree in penetrated cases was larger than in non-penetrated cases. This means that the penetration of the second layer was relevant to the worse chaos degree. The decreasing speed of chaos degree showed a slight difference between penetrated and non-penetrated cases for tBLG. Whereas for BLG, the decreasing speed of chaos degree in penetrated cases was much slower than that in nonpenetrated cases. This indicated that the slower escape speed of fragments and petals in BLG was one of the reasons for the fracture of the second layer. Due to the high-speed impact of the bullet, three to six cracks were initiated near the impact center. The huge kinetic energy destroyed the graphene bond in the center of the impact, then due to the damaged area spread quickly, some atoms near the edge of petals were divorced from the petals, forming some fragments in the damaged area. Those fragments trapped in the impact region owing to the high velocity of the bullet were further pushed by the bullet and result in high-stress concentration points in the second layer, as shown in Figure S4a. It should be noticed that the cracks in the second layer always initiated from the position where some fragments of the fractured first layer existed. The crack-orifice patterns of all the cases simulated in this work are presented in Figure S4b.

Chaos Degree
In penetrated cases, there were more fragments trapped in the impact region on the first layer than those in non-penetrated cases. The fragments were assumed to play an important role in the fracture of the second layer. Here we use the term "chaos degree", defined as the total in-plane perimeter of atoms within the impact region (See Methods section), to evaluate the influence of fragments.
Considering that both cracks and fragments can be expressed by the in-plane perimeter of atoms, chaos degree could be helpful to explore the failure mechanism. The results are shown in Figure 4, and the data in each time point were averaged over all the simulations that belong to the same case. The dramatic increase of the chaos degree at around 2 ps was caused by the fracture of the first layer. The subsequent decrease was caused by the petals and fragments moving out from the impact region. Because of the earlier fracture of the first layer, the increase and decrease of the chaos degree for tBLG were both earlier than that for BLG.

Kinetic Energy Dissipation
The total kinetic energy of atoms within the impact region is shown in Figure 5. The data were averaged over all the simulations that belong to the same case. All the curves exhibited two main peaks. Before the first peak, the atoms were sped up by the bullet consecutively until they reach the same speed as the bullet. Then the atoms started to slow down owing to the reduction of the speed of the bullet. When the first graphene layer fractured, the stored strain energy was released and transformed into the kinetic energy of the petals and fragments, resulting in the formation of the second peak in the curves.
Owing to the earlier fracture of the first layer, the fragments of tBLG had slightly larger kinetic energy than that of BLG. As a result, the fragments of tBLG possessed higher mobility and were more likely to move out from the impact region, making tBLG less likely to be penetrated. Above all, the lower probability of tBLG to create fragments with high chaos degree in the impact region and the higher mobility of the fragments should be responsible for the enhanced capacity of tBLG.  For both tBLG and BLG, the peak chaos degree in penetrated cases was larger than in non-penetrated cases. This means that the penetration of the second layer was relevant to the worse chaos degree. The decreasing speed of chaos degree showed a slight difference between penetrated and non-penetrated cases for tBLG. Whereas for BLG, the decreasing speed of chaos degree in penetrated cases was much slower than that in non-penetrated cases. This indicated that the slower escape speed of fragments and petals in BLG was one of the reasons for the fracture of the second layer.

Kinetic Energy Dissipation
The total kinetic energy of atoms within the impact region is shown in Figure 5. The data were averaged over all the simulations that belong to the same case. All the curves exhibited two main peaks. Before the first peak, the atoms were sped up by the bullet consecutively until they reach the same speed as the bullet. Then the atoms started to slow down owing to the reduction of the speed of the bullet. When the first graphene layer fractured, the stored strain energy was released and transformed into the kinetic energy of the petals and fragments, resulting in the formation of the second peak in the curves.
Owing to the earlier fracture of the first layer, the fragments of tBLG had slightly larger kinetic energy than that of BLG. As a result, the fragments of tBLG possessed higher mobility and were more likely to move out from the impact region, making tBLG less likely to be penetrated. Above all, the lower probability of tBLG to create fragments with high chaos degree in the impact region and the higher mobility of the fragments should be responsible for the enhanced capacity of tBLG.
The earlier fracture of the first layer was assumed to be caused by the mismatch between two sheets of graphene. During the transverse impact, the distance between two layers of graphene was reduced greatly, leading to intensive interaction between atoms in two layers. Because of the perfect overlap of the two layers of graphene and the symmetry of lattice for BLG, each atom on the first layer suffered identical force generated by atoms in the second layer. All the forces were perpendicular to the graphene plane. For tBLG, the symmetry was broken down. The forces generated by the second layer were different for each atom in the first layer and no longer perpendicular to the graphene plane, whose component force in the plane helped to break the bonds between atoms and caused the earlier fracture of the first layer.
Owing to the earlier fracture of the first layer, the fragments of tBLG had slightly larger kinetic energy than that of BLG. As a result, the fragments of tBLG possessed higher mobility and were more likely to move out from the impact region, making tBLG less likely to be penetrated. Above all, the lower probability of tBLG to create fragments with high chaos degree in the impact region and the higher mobility of the fragments should be responsible for the enhanced capacity of tBLG.  Besides the symmetry, the mismatch of lattice also resulted in a decrease in the crack number, which may have reduced the formation of fragments. As shown in Figure 3 and Figure S4, in most cases, six main cracks formed in the first graphene layer in BLG, but only three to four in tBLG. The cracks always propagated along the six zigzag directions, as illustrated by the red arrows in Figure 1. For tBLG, when cracks in the first layer propagated along different zigzag directions, the atoms near the cracks suffered different resistance forces generated from the second graphene layer. Therefore, the cracks propagated preferentially along certain zigzag directions and resulted in a decrease in the crack number.

Conclusions
We assessed the ballistic impact resistance performance of twisted bilayer graphene using molecular dynamics simulations under various conditions. We revealed that the second layer of twisted bilayer graphene with the twist angle of 21.79 • was harder to penetrate than the second layer of the referenced bilayer graphene by simulating their responses under the high-speed transverse impact, which was selected to always penetrate the first layer. It was observed that the penetration of the second layer was related to the fragments of the first layer within the impact region. The fragments in the twisted bilayer graphene were less likely to be trapped in the impact region due to their higher mobility resulting from the earlier fracture of the first layer than in the referenced bilayer graphene. The statistical analysis showed that the ballistic resistance capacity of the twisted bilayer graphene was almost double that of AA bilayer graphene. Our atomistic results and insights may be helpful in the design of graphene-based armor materials.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4 352/11/2/206/s1, Figure S1: Configurations and monitor points, Figure S2: Energy propagation process along the x-direction of tBLG samples, Figure S3: Energy propagation process at 2 km/s, Figure S4: Stress distribution of lower graphene layer for tBLGs with different twist angle, Figure  S5 Author Contributions: Q.P. and Q.C. conceived the idea of the paper; Q.P. and Q.C. designed and performed the simulations and wrote the original draft; S.P. and Q.P. discussed the results and modified the manuscript. All authors have read and agreed to the published version of the manuscript.