Theoretical Evaluation of Impact Characteristics of Wavy Graphene Sheets with Disclinations Formed by Origami and Kirigami

Evaluation of impact characteristics of carbon nanomaterials is very important and helpful for their application in nanoelectromechanical systems (NEMS). Furthermore, disclination lattice defects can generate out-of-plane deformation to control the mechanical behavior of carbon nanomaterials. In this study, we design novel stable wavy graphene sheets (GSs) using a technique based on origami and kirigami to control the exchange of carbon atoms and generate appropriate disclinations. The impact characteristics of these GSs are evaluated using molecular dynamics (MD) simulation, and the accuracy of the simulation results is verified via a theoretical analysis based on continuum mechanics. In the impact tests, the C60 fullerene is employed as an impactor, and the effects of the different shapes of wavy GSs with different disclinations, different impact sites on the curved surface, and different impact velocities are examined to investigate the impact characteristics of the wavy GSs. We find that the newly designed wavy GSs increasingly resist the kinetic energy (KE) of the impactor as the disclination density is increased, and the estimated KE propagation patterns are significantly different from those of the ideal GS. Based on their enhanced performance in the impact tests, the wavy GSs possess excellent impact behavior, which should facilitate their potential application as high-impact-resistant components in advanced NEMS.


Introduction
Impact resistance is an essential evaluation criterion in material design, especially for the application of materials developed for use in aircraft, automobiles, sports products, and machine elements. Since they were first reported [1], graphene sheets (GSs) have attracted enormous interest among scholars globally. Ideal GSs are one-carbon-atom-thick, twodimensional carbon nanomaterials with hexagonal lattice structures that exhibit excellent mechanical, electrical, chemical, and thermal properties [2][3][4]. In particular, the high-strainrate behavior [5] of GSs can be exploited for potential applications as high-impact-resistant components in nanoelectromechanical systems (NEMS). The difficulty and complexity of nano-scale impact simulations has meant that impact testing of carbon nanomaterials, such as GSs [6][7][8][9] and carbon nanotubes (CNTs) [10], has been performed based on molecular dynamics (MD) simulation, revealing the high impact resistances of both GSs and CNTs. In addition, a numerical spring-mass model has been adopted to describe the impact behavior

Simulation Models
Origami and kirigami, traditional Japanese artistic techniques, have been utilized in a wide range of natural science and engineering fields [40]. Miura-ori, conceived by Miura in 1970 [41], involves nonlinear creasing of materials to allow expansion, and folding in the vertical and horizontal directions to improve their strength. We adopted origami and kirigami in this study to form complicated paper models of GSs with disclinations, instead of the complex phase field crystal (PFC) method used for designing curved GS with lattice defects [20,23]. Positive and negative disclinations were introduced by removing and inserting wedge shapes from/into an ideal planar GS, respectively. The details of the proposed origami-and kirigami-based technique for forming wavy GSs with disclinations are presented in Figure 1a. Based on the paper models constructed using the origami-kirigami technique, atomic models of the disclination-incorporated wavy GSs were created, with different types of disclinations generating different out-of-plane deformations and yielding differently shaped GSs. A four-membered ring in the GSs was defined as an isolated positive disclination that forms a conical curved surface, whereas an eight-membered ring in the GSs was defined as an isolated negative disclination that forms a saddle-shaped curved surface, as shown in Figure 1b. Thus, each repeating unit of the wavy GSs was constructed as a disclination quadrupole combining two pairs of the four-and eight-membered rings, i.e., four isolated disclinations. Observing the structure along the z-axis from the positive direction, the two eight-membered rings are arranged in the upper right and lower left regions, and the two four-membered rings are arranged in the lower right and upper left regions. The four isolated disclinations connected in this manner were regarded as one basic wavy unit, and multiple units were connected to create a unique wavy GS model constructed as a network of disclinations. Although an ideal GS has a stable structure with a planar shape, the introduction of disclinations formed a GS with an inherently curved surface, as shown in Figure 2. Furthermore, the wavy GS with an assembly of disclinations exhibited greater stability than the ideal GS, with improvements in the impact characteristics of the GS that might possibly be attributed to the corrugated structure.
As shown in Figure 2, the impact tests involved impinging the fullerene C 60 (mass, 1.197 zg; diameter, 6 Å) on a GS with dimensions of 120 Å × 120 Å. Moreover, all four edges of the target GSs were fixed during the impact test. Impact tests were performed on an ideal GS and four types of wavy GSs (named as wavy1, wavy2, wavy3, and wavy4) with different disclination densities, as shown in Table 1. The central position in the z direction of each GS model was set to 0 Å. Because C 60 is initially positioned above the GS, it was necessary to place it at a sufficiently large distance to avoid the influence of the van der Waals (vdW) interaction forces between the carbon atoms of C 60 and the GS (20-30 Å) [42]; therefore, in this work, the height of C 60 was set to 30 Å. As detailed in Table 1, l was the distance between adjacent upward-facing four-membered rings (Figure 2), and the disclination density was defined as the number of disclinations per unit area. After performing the impact simulation with different types of GSs, the wavy4 model with l = 31.79 Å was selected for further analysis using different positions of impact and impact velocities. With respect to the points of impact, the lowest and highest points (z positions with the greatest magnitudes) nearest to the center of the wavy GS model were defined as A0 and A4, respectively; moreover, three evenly spaced points between A0 and A4 were defined as A1, A2, and A3 (black circles in Figure 2).  As shown in Figure 2, the impact tests involved impinging the fullerene C60 (mass, 1.197 zg; diameter, 6 Å) on a GS with dimensions of 120 Å × 120 Å. Moreover, all four edges of the target GSs were fixed during the impact test. Impact tests were performed on an ideal GS and four types of wavy GSs (named as wavy1, wavy2, wavy3, and wavy4) with different disclination densities, as shown in Table 1. The central position in the z

MD Analysis Conditions
The impact tests were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [43], and atomic-scale visualization was performed using the Open Visualization Tool (OVITO). The adaptive intermolecular reactive empirical bond order (AIREBO) [44] force field (Equation (1)), which has been extensively used for probing the mechanical properties of carbon nanomaterials in general, was employed for all the intra-and inter-molecular interactions.
where is the REBO potential of a covalent bond, is the Lennard-Jones (LJ) potential of a non-covalent bond, and is tortional potential that depends on the dihedral angle. Because the AIREBO potential has been applied for calculating dihedral angles, it was a suitable choice for use in evaluating the curved surface model characteristics of the constructed GSs. In addition, because the vdW interaction forces cannot be ignored, the interactions between the carbon atoms in C60 and the GSs were considered. The cut-off distance in the potential was set to 1.7 Å to maintain the smoothness of the wavy GS structures and avoid the scattering of atoms during the impact tests. A time step of 1.0 fs and a temperature of 5 K were employed in all the simulations. Prior to the impact tests with C60, the target GSs and C60 were fully relaxed and equilibrated at 5 K and a pressure of 0 atm in an isothermal-isobaric ensemble (NPT) for Model of wavy1 GS, l is 11.42 Å. Disclination density is 31.30 × 10 −3 /Å 2 .

MD Analysis Conditions
The impact tests were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [43], and atomic-scale visualization was performed using the Open Visualization Tool (OVITO). The adaptive intermolecular reactive empirical bond order (AIREBO) [44] force field (Equation (1)), which has been extensively used for probing the mechanical properties of carbon nanomaterials in general, was employed for all the intra-and inter-molecular interactions.
where is the REBO potential of a covalent bond, is the Lennard-Jones (LJ) potential of a non-covalent bond, and is tortional potential that depends on the dihedral angle. Because the AIREBO potential has been applied for calculating dihedral angles, it was a suitable choice for use in evaluating the curved surface model characteristics of the constructed GSs. In addition, because the vdW interaction forces cannot be ignored, the interactions between the carbon atoms in C60 and the GSs were considered. The cut-off distance in the potential was set to 1.7 Å to maintain the smoothness of the wavy GS structures and avoid the scattering of atoms during the impact tests. A time step of 1.0 fs and a temperature of 5 K were employed in all the simulations. Prior to the impact tests with C60, the target GSs and C60 were fully relaxed and equilibrated at 5 K and a pressure of 0 atm in an isothermal-isobaric ensemble (NPT) for Model of wavy2 GS, l is 18.06 Å. Disclination density is 11.40 × 10 −3 /Å 2 .

MD Analysis Conditions
The impact tests were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [43], and atomic-scale visualization was performed using the Open Visualization Tool (OVITO). The adaptive intermolecular reactive empirical bond order (AIREBO) [44] force field (Equation (1)), which has been extensively used for probing the mechanical properties of carbon nanomaterials in general, was employed for all the intra-and inter-molecular interactions.
where is the REBO potential of a covalent bond, is the Lennard-Jones (LJ) potential of a non-covalent bond, and is tortional potential that depends on the dihedral angle. Because the AIREBO potential has been applied for calculating dihedral angles, it was a suitable choice for use in evaluating the curved surface model characteristics of the constructed GSs. In addition, because the vdW interaction forces cannot be ignored, the interactions between the carbon atoms in C60 and the GSs were considered. The cut-off distance in the potential was set to 1.7 Å to maintain the smoothness of the wavy GS structures and avoid the scattering of atoms during the impact tests. A time step of 1.0 fs and a temperature of 5 K were employed in all the simulations. Prior to the impact tests with C60, the target GSs and C60 were fully relaxed and equilibrated at 5 K and a pressure of 0 atm in an isothermal-isobaric ensemble (NPT) for Model of wavy3 GS, l is 25.38 Å. Disclination density is 6.29 × 10 −3 /Å 2 .

MD Analysis Conditions
The impact tests were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [43], and atomic-scale visualization was performed using the Open Visualization Tool (OVITO). The adaptive intermolecular reactive empirical bond order (AIREBO) [44] force field (Equation (1)), which has been extensively used for probing the mechanical properties of carbon nanomaterials in general, was employed for all the intra-and inter-molecular interactions.
where is the REBO potential of a covalent bond, is the Lennard-Jones (LJ) potential of a non-covalent bond, and is tortional potential that depends on the dihedral angle. Because the AIREBO potential has been applied for calculating dihedral angles, it was a suitable choice for use in evaluating the curved surface model characteristics of the constructed GSs. In addition, because the vdW interaction forces cannot be ignored, the interactions between the carbon atoms in C60 and the GSs were considered. The cut-off distance in the potential was set to 1.7 Å to maintain the smoothness of the wavy GS structures and avoid the scattering of atoms during the impact tests. A time step of 1.0 fs and a temperature of 5 K were employed in all the simulations. Prior to the impact tests with C60, the target GSs and C60 were fully relaxed and equilibrated at 5 K and a pressure of 0 atm in an isothermal-isobaric ensemble (NPT) for Model of wavy4 GS, l is 31.79 Å. Disclination density is 3.96 × 10 −3 /Å 2 .

MD Analysis Conditions
The impact tests were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [43], and atomic-scale visualization was performed using the Open Visualization Tool (OVITO). The adaptive intermolecular reactive empirical bond order (AIREBO) [44] force field (Equation (1)), which has been extensively used for probing the mechanical properties of carbon nanomaterials in general, was employed for all the intra-and inter-molecular interactions.
is the REBO potential of a covalent bond, E LJ ij is the Lennard-Jones (LJ) potential of a non-covalent bond, and E tors kijl is tortional potential that depends on the dihedral angle. Because the AIREBO potential has been applied for calculating dihedral angles, it was a suitable choice for use in evaluating the curved surface model characteristics of the constructed GSs. In addition, because the vdW interaction forces cannot be ignored, the interactions between the carbon atoms in C 60 and the GSs were considered. The cutoff distance in the potential was set to 1.7 Å to maintain the smoothness of the wavy GS structures and avoid the scattering of atoms during the impact tests. A time step of 1.0 fs and a temperature of 5 K were employed in all the simulations. Prior to the impact tests with C 60 , the target GSs and C 60 were fully relaxed and equilibrated at 5 K and a pressure of 0 atm in an isothermal-isobaric ensemble (NPT) for 5 ps. The relaxation was conducted using the conjugate gradient (CG) method, and the accuracy (∆E) was set to 10 −20 . Moreover, all the impact tests were conducted using the NVE ensemble, i.e., Newton's equation.

Continuum Mechanics Method
A theoretical analysis of an ideal GS based on continuum mechanics was performed with the aim of verifying the accuracy of our MD simulations. According to the continuum mechanics method, an ideal GS can be simulated using a continuum plate model; hence, the governing equation of motion for a rectangular GS can be written as follows: where z(x, y, t) is the flexural deflection of the rectangular GS, t is the time, q(x, y, t) is the transverse load acting on the GS plate, and ρ is the density of the GS. D is the flexural rigidity of the GS plate, which is expressed as follows: where E is the elastic modulus of the GS, h is the thickness of the GS, and ν is Poisson's ratio for the GS. All edges of the rectangular GS are simply supported; therefore, the displacement and load function could be represented as follows: where Z mn (t) is the time-dependent coefficient, and Q mn (t) is a term corresponding to the Fourier series expansion. L a and L b (both equal to 120 Å) are the side lengths of the GSs. A concentrated load, F(t), is assumed to be located at the (ξ, η) position. Therefore, Substituting Equations (3)-(6) into Equation (2) resulted in the following expression: ..
where j = 1,2, and The solution to Equation (7) can be expressed as The deflection at an impact position (ξ, η) located in the middle of the GS plate can be expressed as Therefore, the contact deformation of the rectangular GS plate can be written where w(t), as expressed in Equation (12), gives the displacement of the impact mass.
where m 0 and v 0 F are the mass and initial velocity of the impactor C 60 , respectively. The Hertz law of contact [45] was employed to express the relationship between the contact force, F(t), and deformation, α(t), as follows: where K, the Hertzian contact stiffness, was determined using the radius of the impactor and the stiffness values of the impactor and GS plate. The energy Ψ(t) absorbed by the GS plate within time t can be expressed in term of the impact force F acting on the GS plate and the displacement w of the point of impact:

Comparison of Ideal GS Impact Tests via MD Simulation and Continuum Mechanics Method
Impact tests of an ideal GS using MD simulation and the continuum mechanics method were first compared. Figure 3 shows the calculated energies for the impact of an ideal GS with C 60 at an initial impact velocity of 10 Å/ps, along with the total kinetic energies (KEs) and potential energies (PEs) values of GS and C 60 , respectively, which suggest that the KE of the impactor can be absorbed as the KE and PE of GS. For GS, the impact-induced deformations were not completely preserved, and it returned to its initial shape, resulting in as absorption of the KE of the impactor; the KE of GS is a little higher than the KE of the impactor due to the newly generated negative PE of GS. For C 60 , because it was considered as a rigid body, the KE of C 60 after impact became completely 0.0 eV. In addition, the theoretical value of the initial KE of fullerene, m F v 0 F 2 /2 = 3.75 eV, was almost equal to total energy, which was depicted in Figure 3 (purple line) to facilitate this comparison.
In the MD simulation, the initial distance between C 60 and the ideal GS was set as 30 Å to eliminate the influence of the vdW interaction forces; therefore, a period of descent of the C 60 prior to its impact on the GS was considered. The KE of the GS during this period was a nonzero constant value prior to the start of the impact. The atoms in the GS were presumed to sensitively react to slight changes prior to impact and continue to experience some movement after the impact. When C 60 impinged on the GS, the KE of the GS considerably increased with increasing impact time. However, the theoretical results obtained using the continuum mechanics method showed that the KE of the GS increased from zero and converged to a value a little higher than the initial KE of the impinging C 60 . The KE of the GS prior to impact was zero in the theoretical calculations because the minimal motions of the atoms in the GS were neglected. A comparison of the post-impact KEs of the GS revealed reasonable agreement between the KEs obtained using the two methods.
In particular, the results satisfy the law of conservation of energy, indicating the suitability of the employed MD simulation for investigating the impact tests. Similar MD simulation impact tests were performed on the wavy GSs with disclinations, and the effects of disclination density, KE distribution, site of impact, and impact velocity were examined.

Effects of Disclination Density
Impact simulations on the ideal and four wavy GSs with different disclination densities (Table 1) were conducted using = 10 Å/ps to analyze the impact characteristics, as shown in Figure 4. For the ideal and wavy GSs, the points of impact were at the center and the A4 site near the center, respectively ( Figure 2). Figure 4a shows the KEs of each GS obtained during the impact tests. The A4 location corresponds to the tip of a protruding portion of the wavy GSs, and represents the region with the greatest z-coordinate value. Therefore, A4 is guaranteed to be in contact with any impactor for at least some time during the collision, even if the impactor is larger than C60 or has a flat tip. In addition, the pre-impact KEs of the ideal and wavy GSs are almost the same, presumably because the numbers of carbon atoms in the models were composed of approximately the same 10,000 carbon atoms; the Kes per carbon atom of the ideal and wavy GSs were approximately equal to each other. A difference can be observed in the impact time at which the KE of each GS began to increase because of the varying heights of the A4 sites of the different wavy models. Therefore, wavy4, which has the highest A4 site, was impacted the earliest, followed by wavy3; the ideal GS exhibited the longest impact time. At the end of each impact test, the KE of all the GS converged to approximately the same value, i.e., the initial KE of the C60 impactor.

Effects of Disclination Density
Impact simulations on the ideal and four wavy GSs with different disclination densities (Table 1) were conducted using v 0 F = 10 Å/ps to analyze the impact characteristics, as shown in Figure 4. For the ideal and wavy GSs, the points of impact were at the center and the A4 site near the center, respectively ( Figure 2). Figure 4a shows the KEs of each GS obtained during the impact tests. The A4 location corresponds to the tip of a protruding portion of the wavy GSs, and represents the region with the greatest z-coordinate value. Therefore, A4 is guaranteed to be in contact with any impactor for at least some time during the collision, even if the impactor is larger than C 60 or has a flat tip. In addition, the pre-impact KEs of the ideal and wavy GSs are almost the same, presumably because the numbers of carbon atoms in the models were composed of approximately the same 10,000 carbon atoms; the Kes per carbon atom of the ideal and wavy GSs were approximately equal to each other. A difference can be observed in the impact time at which the KE of each GS began to increase because of the varying heights of the A4 sites of the different wavy models. Therefore, wavy4, which has the highest A4 site, was impacted the earliest, followed by wavy3; the ideal GS exhibited the longest impact time. At the end of each impact test, the KE of all the GS converged to approximately the same value, i.e., the initial KE of the C 60 impactor.   During the impact simulation, the wavy1 model exhibited higher KE values than the other models. However, in terms of the change in KE immediately before and after the impact, that is, the impact-induced change in KE, the ideal and wavy1 GSs models displayed similar behavior, with the other models showing smaller changes. In essence, the ideal GS exhibited a greater capacity to absorb the KE of the impactor, whereas the wavy4 GS absorbed the least KE from the impactor, as shown in Figure 4b. ΔE is the amount of KE change of GSs before and after being impacted by C60, as shown in Figure 4a. The wavy GSs have high resistance to impact. This behavior can also be observed by comparing the KEs of C60 during the collisions, as shown in Figure 4c. The C60 that collided with the ideal GS lost almost all of its KE, whereas the C60 units that collided with the wavy2, wavy3, and wavy4 GSs retained KEs of approximately half of their original value after impact. Compared to the ideal GS, the wavy GSs were better able to resist the KE of the impactor.

Distribution of KE
The impact-related changes to the various GS model structures were elucidated using the KE distributions of the GSs at representative durations of 0.0, 0.5, 1.0, and 1.5 ps (Figure 5). Although the initial sizes of the five models are similar, different KE distributions are exhibited because of the different out-of-plane deformations produced by the introduced disclinations. For the ideal, disclination-free GS, the KE distribution expanded as a concentric circle from the central point of impact and was eventually distributed throughout the GS. The wavy1 model with the highest disclination density showed the most extended KE distribution among the four wavy GSs. However, the KE distribution is circular at 0.5 ps and spread in a cross-shaped pattern along the x and y directions, with negligible propagation in the diagonal direction. The wavy2, wavy3, and wavy4 models almost did not absorb KE from the fullerenes, which is consistent with the analysis presented in the previous section; therefore, no noticeable trends were observed in their KE During the impact simulation, the wavy1 model exhibited higher KE values than the other models. However, in terms of the change in KE immediately before and after the impact, that is, the impact-induced change in KE, the ideal and wavy1 GSs models displayed similar behavior, with the other models showing smaller changes. In essence, the ideal GS exhibited a greater capacity to absorb the KE of the impactor, whereas the wavy4 GS absorbed the least KE from the impactor, as shown in Figure 4b. ∆E is the amount of KE change of GSs before and after being impacted by C 60 , as shown in Figure 4a. The wavy GSs have high resistance to impact. This behavior can also be observed by comparing the KEs of C 60 during the collisions, as shown in Figure 4c. The C 60 that collided with the ideal GS lost almost all of its KE, whereas the C 60 units that collided with the wavy2, wavy3, and wavy4 GSs retained KEs of approximately half of their original value after impact. Compared to the ideal GS, the wavy GSs were better able to resist the KE of the impactor.

Distribution of KE
The impact-related changes to the various GS model structures were elucidated using the KE distributions of the GSs at representative durations of 0.0, 0.5, 1.0, and 1.5 ps ( Figure 5). Although the initial sizes of the five models are similar, different KE distributions are exhibited because of the different out-of-plane deformations produced by the introduced disclinations. For the ideal, disclination-free GS, the KE distribution expanded as a concentric circle from the central point of impact and was eventually distributed throughout the GS. The wavy1 model with the highest disclination density showed the most extended KE distribution among the four wavy GSs. However, the KE distribution is circular at 0.5 ps and spread in a cross-shaped pattern along the x and y directions, with negligible propagation in the diagonal direction. The wavy2, wavy3, and wavy4 models almost did not absorb KE from the fullerenes, which is consistent with the analysis presented in the previous section; therefore, no noticeable trends were observed in their KE distributions. Videos of the KE distribution during the impact tests obtained via MD simulation, provided as Supplementary Videos S1-S5, clearly illustrate the spread of the KE distribution. Therefore, in the ideal GS, the impact-induced KE spread concentrically and was ultimately dispersed throughout the GS, whereas in the wavy GSs, which absorbed similar amounts of KEs as the ideal GS, the KE propagation was limited in terms of direction and range. Nanomaterials 2022, 12, x FOR PEER REVIEW 13 distributions. Videos of the KE distribution during the impact tests obtained via MD ulation, provided as Supplementary Videos S1-S5, clearly illustrate the spread of th distribution. Therefore, in the ideal GS, the impact-induced KE spread concentrically was ultimately dispersed throughout the GS, whereas in the wavy GSs, which abso similar amounts of KEs as the ideal GS, the KE propagation was limited in terms of d tion and range.

Effects of Impact Location
Different impact characteristics, such as energy conversion and GS deforma were obtained for the collisions between C60 and different locations within the wavy The wavy4 model, which had the lowest disclination density, was selected as a represe tive GS to reveal the effects of impact location. The wavy4 model was impacted wit at = 10 Å/ps at five locations (A0-A4; Figure 2), and the results are shown in Figu The moments of impact at the different sites occurred at different times because o different distances between each impact location and C60. After the impact, the KEs o GS impacted at the A0, A1, A2, and A3 sites converged to higher KE values. The A0 corresponds to the lowest point on the surface (most negative z-coordinate value) o inverted cone in the wavy GSs; therefore, C60 collides with the A0 site in a manner sim to the dropping of a mass into a saucer, i.e., with minimal bouncing. After the first im of C60 on the GS at the A3 site, it bounces and impacts the GS again at a different loca which results in two energy changes in the curve corresponding to A3 prior to con gence. In the scenarios involving the A0, A1, and A2 sites, the post-impact KEs incr

Effects of Impact Location
Different impact characteristics, such as energy conversion and GS deformation, were obtained for the collisions between C 60 and different locations within the wavy GSs. The wavy4 model, which had the lowest disclination density, was selected as a representative GS to reveal the effects of impact location. The wavy4 model was impacted with C 60 at v 0 F = 10 Å/ps at five locations (A0-A4; Figure 2), and the results are shown in Figure 6. The moments of impact at the different sites occurred at different times because of the different distances between each impact location and C 60 . After the impact, the KEs of the GS impacted at the A0, A1, A2, and A3 sites converged to higher KE values. The A0 site corresponds to the lowest point on the surface (most negative z-coordinate value) of the inverted cone in the wavy GSs; therefore, C 60 collides with the A0 site in a manner similar to the dropping of a mass into a saucer, i.e., with minimal bouncing. After the first impact of C 60 on the GS at the A3 site, it bounces and impacts the GS again at a different location, which results in two energy changes in the curve corresponding to A3 prior to convergence. In the scenarios involving the A0, A1, and A2 sites, the post-impact KEs increase only slightly over time and nearly converged to the converged KE of location A3; however, the KE corresponding to location A4 converged to a value approximately half that of the converged KE of the other impact sites. This is because the collision of C 60 with site A4, which is the highest point on the wavy GS surface (most positive z-coordinate value), resulted in C 60 bouncing directly off the GS and not bouncing back onto it, unlike the behavior at the other impact sites. However, because location A4 is the highest point on the wavy GS surface, it can be presumed that this is the most probable location for impact, especially when the diameter of the impactor is considerably larger than that of C 60 . The KE absorption at A4 is lower than that at the other locations. Therefore, when a collision occurs at A4, the wavy GS can presumably retain the impact-induced shape change, that is, it can absorb more KE. By contrast, the protruding portion of the surface cannot easily resist the deformation, which means that there is little KE absorption. Therefore, KE absorption by the wavy GSs is strongly dependent on the site of impact. The wavy structure can be approximately classified into protruding areas and other areas, and we established that KE absorption by the wavy GSs was greater for collisions with surface depressions. only slightly over time and nearly converged to the converged KE of location A3; however, the KE corresponding to location A4 converged to a value approximately half that of the converged KE of the other impact sites. This is because the collision of C60 with site A4, which is the highest point on the wavy GS surface (most positive z-coordinate value), resulted in C60 bouncing directly off the GS and not bouncing back onto it, unlike the behavior at the other impact sites. However, because location A4 is the highest point on the wavy GS surface, it can be presumed that this is the most probable location for impact, especially when the diameter of the impactor is considerably larger than that of C60. The KE absorption at A4 is lower than that at the other locations. Therefore, when a collision occurs at A4, the wavy GS can presumably retain the impact-induced shape change, that is, it can absorb more KE. By contrast, the protruding portion of the surface cannot easily resist the deformation, which means that there is little KE absorption. Therefore, KE absorption by the wavy GSs is strongly dependent on the site of impact. The wavy structure can be approximately classified into protruding areas and other areas, and we established that KE absorption by the wavy GSs was greater for collisions with surface depressions.

Effects of Impact Velocity
The wavy GSs exhibited reduced Young's moduli and bending moments compared to those exhibited by the ideal GS. The wavy GSs can be considered as consisting of deformable shells. To investigate the effects of impactor velocity on the KEs of the GSs, the wavy4 model was impacted at site A4 by C60 at velocities, , of 10, 20, 30, 40, 50, and 60 Å/ps. As shown in Figure 7, the pre-impact KEs of all the GSs are similar. However, all the KEs of the GSs converged to different constant values. For comparison, short lines colored to match the corresponding KE curve have been included on the right in Figure 7; these lines indicate the values of 1 (2 ( ) ) ⁄ , which increase with impactor velocity. However, the converging value of the KE is close to 1 (2 ( ) ) ⁄ only at = 10 Å/ps, which can be explained by considering the conservation of energy. With increasing impact velocity, the difference between the value of 1 (2 ( ) ) ⁄ and the corresponding Figure 6. Comparison of KEs of the wavy4 GS model impacted by C 60 at v 0 F = 10 Å/ps at the A0, A1, A2, A3, and A4 sites.

Effects of Impact Velocity
The wavy GSs exhibited reduced Young's moduli and bending moments compared to those exhibited by the ideal GS. The wavy GSs can be considered as consisting of deformable shells. To investigate the effects of impactor velocity on the KEs of the GSs, the wavy1 model was impacted at site A4 by C 60 at velocities, v 0 F , of 10, 20, 30, 40, 50, and 60 Å/ps. As shown in Figure 7, the pre-impact KEs of all the GSs are similar. However, all the KEs of the GSs converged to different constant values. For comparison, short lines colored to match the corresponding KE curve have been included on the right in Figure 7; these lines indicate the values of m F v 0 F 2 /2, which increase with impactor velocity. However, the converging value of the KE is close to m F v 0 F 2 /2 only at v 0 F = 10 Å/ps, which can be explained by considering the conservation of energy. With increasing impact velocity, the difference between the value of m F v 0 F 2 /2 and the corresponding converging KE increased. This is because of the microvibrations of the C 60 colliding with the wavy GS at low impact velocities, and the rebounding of the C 60 from the GS at higher velocities; the wavy1 GS was ruptured when C 60 collided with the GS at v 0 F = 60 Å/ps. In particular, at higher impact velocities, the deformation of the GS was considerably greater than that achieved at lower impact velocities. The deformation energies of the GSs with different impact velocities are shown in the inset of Figure 7, confirming that the KEs of the GSs impacted by C 60 at higher velocities are retained at higher velocities.
converging KE increased. This is because of the microvibrations of the C60 colliding with the wavy GS at low impact velocities, and the rebounding of the C60 from the GS at higher velocities; the wavy4 GS was ruptured when C60 collided with the GS at = 60 Å/ps. Therefore, the rate of energy absorption presumably increased as the impact velocity decreased. In particular, at higher impact velocities, the deformation of the GS was considerably greater than that achieved at lower impact velocities. The deformation energies of the GSs with different impact velocities are shown in the inset of Figure 7, confirming that the KEs of the GSs impacted by C60 at higher velocities are retained at higher velocities.

Conclusions
In this work, new stable wavy GSs were generated by appropriately introducing disclinations into an ideal GS using a technique based on origami and kirigami. The impact characteristics of wavy GSs with incorporated disclinations were investigated using MD simulations, and the accuracy of these was verified using a continuum mechanics method. The disclination density significantly affected the impact characteristics of the wavy GSs and even determined their energy absorption and dispersion rates. Because disclinations can produce prominent out-of-plane deformations in GSs, they can reduce the spread of post-impact effects. Consequently, the KE distribution in the ideal GS was an expanding circle, whereas that in the wavy GS exhibited cross-shaped (local) propagation, indicating that the wavy GS structures offer control over the range of KE propagation. By comparing the impact results obtained at different points of impact, the C60 impactor was found to undergo various numbers of bounces and possess varying KE absorption, owing to the nonuniformity of the disclination-induced out-of-plane deformations in GSs.

Conclusions
In this work, new stable wavy GSs were generated by appropriately introducing disclinations into an ideal GS using a technique based on origami and kirigami. The impact characteristics of wavy GSs with incorporated disclinations were investigated using MD simulations, and the accuracy of these was verified using a continuum mechanics method. The disclination density significantly affected the impact characteristics of the wavy GSs and even determined their energy absorption and dispersion rates. Because disclinations can produce prominent out-of-plane deformations in GSs, they can reduce the spread of post-impact effects. Consequently, the KE distribution in the ideal GS was an expanding circle, whereas that in the wavy GS exhibited cross-shaped (local) propagation, indicating that the wavy GS structures offer control over the range of KE propagation. By comparing the impact results obtained at different points of impact, the C 60 impactor was found to undergo various numbers of bounces and possess varying KE absorption, owing to the nonuniformity of the disclination-induced out-of-plane deformations in GSs. Nonetheless, with respect to the impact velocity, the absorption ratio of the target GS decreased at velocities higher than v 0 F = 10 Å/ps. Therefore, the results of these impact tests indicate that these new wavy GS structures possess improved impact characteristics with respect to ideal GSs, which should promote their application as high-impact-resistant components in advanced NEMS.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this paper: NEMS nanoelectromechanical systems GSs graphene sheets CNTs carbon nanotubes C 60 buckminsterfullerene MD molecular dynamics KEs kinetic energies PEs potential energies