Effects of size and surface on the auxetic behaviour of monolayer graphene kirigami

Graphene is an active element used in the design of nano-electro-mechanical systems (NEMS) owing to its excellent in-plane physical properties on mechanical, electric and thermal aspects. Considering a component requiring negative Poisson’s ratio in NEMS, a graphene kirigami (GK) containing periodic re-entrant honeycombs is a natural option. This study demonstrates that a GK with specific auxetic property can be obtained by adjusting the sizes of its honeycombs. Using molecular dynamics experiments, the size effects on the auxetic behaviour of GK are investigated. In some cases, the auxetic difference between the hydrogenated GK and continuum kirigami (CK) is negligible, in which the results from macro CK can be used to predict auxetic behaviour of nano kirigami. Surface effect of GK is demonstrated from two aspects. One is to identify the difference of mechanical responses between the pure carbon GK and the hydrogenated GK at same geometry and loading condition. Another is from the difference of mechanical responses between the GK model and the CK model under same loading condition and geometric configuration. Generally, surface energy makes the GK possess higher variation of auxetic behaviour. It also results in higher modulus for the GK as comparing with that of the CK.

Scientific RepoRts | 6:35157 | DOI: 10.1038/srep35157 Single/few-layer graphene is a typical two dimensional (2D) isotropic material and has attracted much attention due to its excellent physical properties since it was obtained in laboratory 25 . For an ideal graphene, it is not an in-plane auxetics 5 . In recent years, graphene kirigami (GK) has been investigated due to its peculiar properties on deformation, electric conductivity and thermal conductivity [19][20][21][22][23] . As we tailor a kirigami of the 2D material into a re-entrant honeycomb ( Fig. 1a) 6,[26][27][28] , auxetic behaviour will emerge. The question is, will the size effect still be significant on mechanical properties when the re-entrant cellular structure is miniaturized to nanoscale? Besides, will surface energy plays an important role in mechanical behaviour of the tailored GK like it in 2D material with caves? In present study, the deformation analysis of a continuum kirigami (CK) with its geometry matching that of GK is carried out using finite element method (FEM) [29][30][31][32] aligning with molecular dynamics(MD) experiments of the GKs, aiming at revealing the influences of the sizes of microstructure and surface energy [33][34][35] on the mechanical properties of such GK.

Models and Methods
Geometric model of kirigami. The mechanical behaviour of the GK shown in Fig. 1a or continuum model shown in Fig. 1c 6,26-28 depends on surface energy distribution and five independent variables, i.e., θ (angle between oblique bar and vertical bar), w θ (width of oblique bar), w (width of vertical bar), l G (gap of adjacent vertical bars, it reflects the length of vertical bar) and l θ (it reflects length of oblique bar). The deformation of the atomic system (in Fig. 1a) will be calculated using the result from MD simulation. The deformation of the CK (in Fig. 1c) will be solved via FEM. The engineering strain along y-direction is controlled within 5%. Hence, the equivalent modulus of kirigami can be calculated using the principle of minimum potential energy (P e ). The engineering strain along x-direction, i.e., ε x , is calculated by the relative expansion of the blue points in Fig. 1a. The engineering strain along y-direction, i.e., ε y , is calculated by the time integration of the velocity of the upper side of specimen. The Poisson's ratio in x-y plane is the ratio between − ε x and ε y . To reflect the effect of surface energy on the equivalent modulus (in Eq. (4)), modulus ratio, i.e., the ratio between the equivalent modulus and the in-plane modulus of the two-dimensional matrix (i.e., the ideal graphene with modulus of ~995.8 GPa at 8 K), is defined as R m in Eq. (5). The Poisson's ratio and the modulus ratio of the atomic system and the continuum system will be compared to show the effect of surface energy. Geometric size effect is considered by varying the five independent geometric variables. The Poisson's ratio and modulus ratio can be obtained using the following equations. Where L x0 and L xt are the initial and current horizontal distance of the red and blue points in Fig. 1a. t is the time of loading. V 0 is the volume occupied by the atoms in statistic area for finding the potential energy of the atoms (P e ). The thickness of graphene is set to be 0.335 nm in this study.   36 , is adopted. In simulation, the AIREBO 37 potential function is used to describe the interaction among the hydrogen and carbon atoms. The time step for integration is 0.0005 ps. The procedure of MD experiment contains following steps: Step 1) Build the kirigami model (within the green frame in Fig. 1a) with specified geometry. Once the hydrogenation model is adopted to the simulation, the red bond of edge in Fig. 1b will be bonded with hydrogen atoms; Step 2) Set the two vertical sides (left and right sides) of the nanostructure as free boundaries, and the two horizontal sides as periodic boundaries (Let l a = l b + l c ); Step 3) Reshape the structure by minimizing the potential energy of the nanostructure. The conjugated-gradient algorithm is adopted for update the positions of the atoms; Step 4) Relax the nanostructure in a Nosé-Hoover thermostat for 200 ps with temperature of 8 K; Step 5) Provide a strain ratio of e r = 0.005%/ps along y-direction on specimen. The total motion costs 1000 ps; Step 6) Calculate the engineering strains along x-and y-directions according to the positions of the blue points in Fig. 1a. Calculate the Poisson's ratio according to Eq. (3). Calculate the equivalent elastic modulus of the kirigami using Eq. (4). And find the ratio between the equivalent modulus and that of the ideal graphene.
To show the effect of surface energy on the mechanical properties of GK, hydrogenation schemes are added. If there is no carbon atom on the graphene bonded with hydrogen atom, the model is labelled as "+ 0H"; if each carbon atom on the edge (i.e., unsaturated carbon atom) is bonded with a hydrogen atom, the model is labelled as "+ 1H". The potential energy and volume of the system with respect to the "+ 1H" model are the summation of those of carbon atoms only. All the results will be compared with those of a macrostructure with geometric similarity to the GK. The Poisson's ratio and the modulus ratio of the structure of CK are analysed by finite element method. The analysis procedure of a deformed continuum using FEM includes the following steps: Step a) Build the CK structure (shown in Fig. 1c) which has the same geometric configuration with that of the GK; Step b) Mesh the CK structure with finite elements (plane stress elements); Step c) Apply displacement constraints: the lower side is fixed along y-direction. The node located at the center of the lower side is fully fixed. Provide the nodes on the upper side of the structure with specified displacement, which reflects the same engineering strain along y-direction, i.e., no more than 5%; Step d) Find the deformation of the CK structure by solving its equilibrium equation; Step e) Calculate the Poisson's ratio using the displacement of the labelled points in Fig. 1c and equation (3).
Calculate the modulus ratio using equations (5).

Schemes for size effect analysis.
To show the size effect of the GK, the following schemes are considered: 120° and 150°. The microstructures are shown in Fig. 2a; The microstructures are shown in Fig. 2c;

Theoretical analysis
To estimate the mechanical behaviour of GK, a theoretical analysis is introduced here. In Fig. 3, the right part with yellow and grey areas (including adjacent bars) is a quarter of the unit cell of kirigami. The equivalent mechanical properties, e.g., equivalent modulus and Poisson's ratio are determined by the mechanical and geometric parameters of the oblique and vertical bars. We suppose the total strain is along y direction, i.e., ε y , is specified and label the equivalent moduli of yellow and grey areas as E 1 and E 2 , respectively. The two moduli can be obtained using MD simulation. The length weights of the two areas are w 1 and w 2 (= 1 − w 1 ), respectively. Hence, the equivalent modulus of the kirigami can be calculated using Let ε y1 and ε y2 represent the line strains of the yellow and grey areas along y-direction, respectively. They can be expressed as y y y y The line strain of the half part along x-direction, x-strain, can be calculated using: Where ε x1 (< 0) is the x-strain of the vertical bar adjacent to the yellow area. The value can be obtained using MD simulation. ε x2 is the x-strain of the grey area. For the structure with small deformation, it can be obtained using the following formulation.
Hence, the equivalent Poisson's ratio can be expressed as

Numerical results
Results of Scheme 1_ θ varies. In this scheme, θ is chosen from the angles, i.e., 30, 60, 90, 120 and 150 degrees, respectively. The other parameters are given in Section "Schemes for size effect analysis". After simulation, the Poisson's ratio and the modulus ratio are calculated and shown in Table 1. Table 1 lists the Poisson's ratios and modulus ratios of the GK with different values of θ. For each Poisson's ratio of GK in the three models, the value increases with the monotonous increase of θ. When θ is less than 90°, the Poisson's ratio is negative, i.e., GK shows auxetic behaviour. At θ = 90°, the Poisson's ratios are positive, rather than zero. The reason is that the oblique bar, which currently is horizontal, has a rotary angle which leads to θ > 90°. The table also demonstrates the value of R p is sensitive to the surface energy. For example, for "+ 0H" model, R p is 0.318, which is over twice of 0.149 of "+ 1H" model, and is far greater than 0.035 of CK model. Similar differences among the values of R p of the three models can also be found at any other values of θ. It is concluded that the GK without hydrogenation (+ 0H model) has the highest absolute value of R p among the three models. The hydrogenated GK (+ 1H model) has lower absolute value of R p than the pure carbon kirigami does, but the value is higher than that of the CK (analysed using FEM) if θ is not 60°.
For the modulus ratios, R m , of kirigami in the three models, increase of them is seen with the increase of θ. The reason is that the deformation of the structure is caused by the rotation of the oblique bars when θ is lower, while, at higher value of θ, the deformation of GK is mainly caused by the axial linear strain of bars. Hence, increment of potential energy is higher at higher value of θ. According to Eq. (4), the modulus is higher when θ is greater. Comparing the three models' modulus ratios at the same value of θ, the GK with "+ 1H" model has the highest modulus. It is about 10~20% higher than that of GK with "+ 0H" model. The reason can be found from the variation comparison of potential energy of the two models, i.e., in the "+ 1H" model, the carbon atoms on the edge have more steep increases as they are bonded with hydrogenation atoms. The modulus ratios of the GK with or without hydrogenation are far greater than that of the CK model under the same geometric parameter settings. It demonstrates that the surface influences the modulus of a 2D material significantly.
Results of Scheme 2_ width of oblique bar varies. In this scheme, only N(w θ ) changes from 2 to 7, indicating the width of oblique bars increases monotonously (Fig. 2b) while the other parameters are fixed. After simulation, the Poisson's ratio and the modulus ratio are calculated and shown in Table 2.
When comparing R p and R m of the three models with the same geometric configurations, the surface effect can also be found. Similar conclusions as those in Scheme 1 can be found. Interestingly, the difference between R p of GK with "+ 1H" model and that of the CK model is very small whatever the width of oblique bar is. It evidences that the auxetic behaviour of a 2D material shown in Fig. 2b can be estimated using experiment on a macro-continuum model.  Now, we examine the size effect on the mechanical properties of kirigami via adjusting the width of the oblique bars. Due to θ = 60° (less than 90°), the kirigamis of all the three models show auxetic behaviour when the width of oblique bar varies. It is also found that the absolute value of R p decreases with the increase of the width of oblique bar. The reason can be explained from Eq. (10). In the model, the geometric parameters, e.g., w, l θ , θ, l G , and E 1 and ε y are constants. Both w 2 /w 1 and E 2 increase with the increase of the width of oblique bar. Negative scalar ε x1 decreases with the increase of the width of oblique bar. Hence, the absolute value of the R p is lower when the width of oblique bar is wider. The modulus ratios of the three models are also different. The difference between the moduli of the GKs with and without hydrogenation is smaller than that between the GK models and CK model.

Results of Scheme 3_ width of vertical bar varies.
In this scheme, N(w) is chosen from 1, 3, 5, 7, i.e., the width of vertical bar increases monotonously. The rest parameters are fixed. After simulation, the Poisson's ratio and the modulus ratio are calculated and shown in Table 3.
In the scheme, the absolute values of R p do not vary monotonously, and the value peaks at N(w) = 3. The reason can be identified via the analysis with Eq. (10). For example, when N(w) = 1, the width of vertical bar is less than that of oblique bar, the vertical bar plays a major role in the deformation of GK. The rotary angle of the oblique bar is low because of small tension on the vertical bar (σ y ). As N(w) = 3, the vertical bar has stronger tensile force which leads to higher rotary angle of oblique bar. Therefore, the R p at N(w) = 3 is larger than that at N(w) = 1. When N(w) ≥ 5, the rotary angle of oblique bar is slightly higher than that at N(w) = 3. Simultaneously, w/l θ is over 1.6 times of that at N(w) = 3. From the second item in the left part of Eq. (10), we can find that the wider width of the vertical bar (N(w) ≥ 5) leads to smaller Poisson's ratio of GK. On the other hand, the difference between the Poisson's ratios of GK with "+ 1H" model and CK model is obvious no matter if θ = 60° or not in the models, which implies that the surface effect cannot be neglected when using the continuum model to estimate the auxetic behaviour of GK with different widths of vertical bar.
For the modulus ratios of the three models, they decrease with the increase of width of vertical bar. The reason is that the deformation of GK mainly depends on the rotation of the oblique bar during loading along y-direction. Hence, the variation of potential energy of the vertical bar decreases with the increase of its width. Therefore, the variation of potential energy of system decreases, too, which results in decreasing of the modulus according to Eq. (4).

Results of Scheme 4_ length of vertical bar varies.
In this scheme, the length of vertical bar N(l G ) is chosen from 6, 12, 18, 24, 30 and 36. The other parameters are fixed. After simulation, the Poisson's ratio and the modulus ratio are calculated and shown in Table 4.
When the length of vertical bar increases or N(l G ) changes from 6 to 36, the absolute values of Poisson's ratio of the three models increase monotonously. The reason can be revealed by the analysis of either Eq. (10) or the deformation of Fig. 3b. For example, for the second item of the left part of Eq. (10), when w, l θ , E 1 and E 2 (< E 1 ), are constants, the increase of w 1 leads to the increase of R p . It suggests that the 2D material can have higher auxetic effect by providing with longer vertical bars. In this scheme, the difference between the values of R p of the hydrogenated GK and CK model is very small. One can estimate the axuetic behaviour of the GK using the experiments on the continuum model with the same geometric configurations according to this principle.
The modulus ratios of the three models increase with the increase of the length of vertical bar. But the maximal value is reached when w 1 tends to be 1.0. The maximal value is E 1 /E matrix (= (w/(w+ l θ )) × (E vert_bar /E matrix ) ≈ 20%) according to Eqs (5 and 6) and Fig. 3b. The hydrogenated GK still has the highest modulus among the three models. The CK model has the lowest modulus, which is no more than 20% of that of the hydrogenated GK.  Table 3. Comparisons of Poisson's ratio and the modulus ratio for the three models in Scheme 3. Results of Scheme 5_ length of oblique bar varies. In this scheme, the length of oblique bar varies, i.e., N(l θ ) is chosen from 12, 15, 18, 24 and 27. The other parameters are fixed. After simulation, the Poisson's ratio and the modulus ratio are calculated and shown in Table 5.
The results of R p of CK with different geometric parameters show that the absolute value of R p decreases slightly with the increase of the length of oblique bar. Although the similar surface effects still exist, demonstrated by the absolute value of R p of hydrogenated GK is lower than that of the pure carbon GK, but greater than that of CK model, the value of R p does not change monotonously with respect to the variation of the length of oblique bar. For example, the absolute value of pure carbon GK peaks at N(l θ ) = 24. While, the maximal absolute value of hydrogenated GK appears at N(l θ ) = 21. The reason can be revealed from the configurations of the models shown in Fig. 4. For instance, at N(l θ ) = 21, the two vertical bars on the same column are not attracted to be bonded together. When N(l θ ) = 24 or 27, the two vertical bars get close to each other, with the distance of only 0.298 nm for N(l θ ) = 24 (Fig. 5a) or 0.265 nm for N(l θ ) = 27. The two adjacent ends of vertical bars are not bonded. However, the van der Waals interaction between the two adjacent ends becomes very strong, which leads to the decrease of the value of θ. Hence, the Poisson's ratio at N(l θ ) = 24 is greater than that at N(l θ ) = 21. When l G keeps unchanged, larger l θ leads to smaller variation of θ after relaxation. That is why the Poisson's ratio at N(l θ ) = 24 is higher than that at N(l θ ) = 27. The electron density distribution nearby the adjacent ends of vertical bars shown in Fig. 5a reveals that no new carbon-carbon covalent bond is generated. Hence, after deformation, the two ends are separated. However, the variation of potential energy is greater than that of the model with N(l θ ) = 21. It results in   Table 5. Comparisons of Poisson's ratio and the modulus ratio for the three models in Scheme 5. higher modulus according to Eq. (4). This is verified by the values of the R m listed in Table 4. On the other hand, if the value of N(l θ ) is far greater than 24, the oblique bars are softer and the two ends of the vertical bars can get closer, which provide a chance to generate new covalent bonds among the unsaturated carbon atoms on the ends (Fig. 5b). In such condition, the topology of the kirigami changes and the new GK may alter mechanical properties, such as becoming stiffer or even non-auxetic. When the GK is hydrogenated ("+ 1H" model), the repulsion exists between the adjacent vertical bar and oblique bar at the joint whilst attraction exists between the two adjacent vertical bars. A model with longer oblique bar will have higher deformation of the oblique bar due to repulsion between the oblique and vertical bars and stronger attraction at adjacent ends of vertical bars. Hence, the peak value of Poisson's ratio appears at N(l θ ) = 21. On the other hand, the two adjacent ends of vertical bars are attracted due to van der Waals interaction, rather than covalent bonds, the potential energy of the system changes slightly. The modulus ratio decreases with the increase of N(l θ ). The modulus ratio of CK model has similar variation. But the modulus ratio of hydrogenated GK is still much higher than that of CK.

Conclusions
After analysis the Poisson's ratio and modulus ratio of GK with re-entrant honeycomb microstructure, the dependence of the two mechanical properties on the sizes of the microstructure is revealed. And the mechanical response of the CK with geometric similarity to the GK is also calculated for finding the surface effect of GK through comparisons of the responses. The results show that the specified Poisson's ratio and modulus of GK can be obtained by adjusting the sizes of microstructure. Some conclusions can be drawn as follows (1) Of all the schemes, pure carbon GK, being with the highest surface energy, has the highest absolute value of Poisson's ratio among the three models. The absolute value of Poisson's ratio of the hydrogenated GK is higher than that of the CK due to the higher surface energy of the GK than that of the CK. Considering the effect of θ in scheme 1, when θ is less than 90°, GK shows auxetic. (2) In general, the modulus ratio of hydrogenated GK is about 10~20% higher than that of the pure carbon GK.
The modulus ratio of CK is far less than that of the GK when they have the same geometric configurations. It demonstrates the surface influence on the modulus of 2D nano materials. (3) In Scheme 2, the difference between the Poisson's ratios of hydrogenated GK and CK models is very small, meaning that the Poisson's ratios depend on the width of oblique bar, slightly. It indicates that the auxetic behaviour of a 2D nanomaterial can be estimated using the experiment on a macro-continuum model with respect to the width of oblique bar. For the modulus ratios of GK and CK models, when the width of oblique bar increases, the absolute value of R p decreases, but the modulus ratio increases, monotonically. (4) When increasing the width of vertical bar like in scheme 3, the peak value of Poisson's ration appears at N(w) = 3. The difference between the Poisson's ratios of GK models and CK model is obvious. As for the modulus ratios of the three models, they decrease with the increase of width of vertical bar. (5) When the length of vertical bar N(l G ) increases from 6 to 36, the absolute values of Poisson's ratio of the three models increases monotonously. It suggests that the 2D material can have higher auxetic effect with longer vertical bars as in Scheme 4. In this condition, the difference between the values of Poisson's ratio of the hydrogenated GK and CK model is very small. The axuetic behaviour of GK can be estimated using the experiments on the CK model with the same geometric configurations. The modulus ratios of the three models increase with the increase of the length of vertical bar. (6) In Scheme 5, the absolute value of Poisson's ratio varies slightly with the increase of the length of oblique bar (or N(l θ )). If the value of l G is small (i.e., vertical bars are short), the two adjacent vertical bars in pure carbon GK model may bonded together. It results in the sharp increase of both Poisson's ratio and modulus ratio. But it does not happen in either the hydrogenated GK model or the CK model.