Atomistic Study of Mechanical Behaviors of Carbon Honeycombs

Abstract

With an ultralarge surface-to-volume ratio, a recently synthesized three-dimensional graphene structure, namely, carbon honeycomb, promises important engineering applications. Herein, we have investigated, via molecular dynamics simulations, its mechanical properties, which are inevitable for its integrity and desirable for any feasible implementations. The uniaxial tension and nanoindentation behaviors are numerically examined. Stress–strain curves manifest a transformation of covalent bonds of hinge atoms when they are stretched in the channel direction. The load–displacement curve in nanoindentation simulation implies the hardness and Young’s modulus to be 50.9 GPa and 461±9 GPa, respectively. Our results might be useful for material and device design for carbon honeycomb-based systems.

1. Introduction

Graphene is known widely, due to its excellent mechanical nature, as a so-called “miracle material”, with many of its characteristics measured experimentally or theoretically exceeding those obtained in other materials: A Young’s modulus of 1 TPa and intrinsic strength of 130 GPa [1,2], high stiffness [3] and fracture strain [4,5], a normal-auxeticity mechanical phase transition [5], etc. Due to these predominant properties, graphene has promising potential to be employed in various applications: Paints and coatings of nanocomposites [4,6], flexible electronics [7], and bioapplications [8,9,10]. Nevertheless, considering the difficulty in engineering synthesis of large-area and high-quality graphene, so far, there is little practical application of industrially produced graphene. Consequently, it is worth focusing on allotropes of graphene [11,12,13], which are more stable and feasible to produce. Recently, a stable carbon allotrope, namely, carbon honeycomb (CHC), a three-dimensional graphene, was synthesized. There are a few reports on of its physical absorption [14,15], electronic band structure [16], thermoelectric performance [17], thermal transport properties [18,19], and phonon properties [15].

According to the cell patterns in the plane perpendicular to the cell axis, carbon honeycombs can be categorized into two sets: Armchair CHC (ac_n CHC) and zigzag CHC (zz_m CHC) [20]. When the number of armchair or zigzag lines (n in ac_n or m in zz_m) varies, the cell sizes of CHCs and, thus, their mechanical properties change correspondently. Krainyukova et al. [14] synthesized this carbon allotrope by deposition of vacuum-sublimated graphite and proposed periodic and random structures of CHCs. Regarding the mechanical properties of CHC, some analytical studies have been carried out [16,19,20,21,22,23,24]. Karfunkel et al. first proposed a family of hypothetical zz_m–sp²–sp³ CHC structures and revealed that the carbon modifications are as stable as diamond by solid-state semi-empirical SCF methods using the modified neglect of diatomic overlap (MNDO) Hamiltonian [23]. Park et al. studied the mechanical properties of zz_m–sp²–sp³ CHCs with different m (where m can be an integer or a half integer) using ab initio pseudopotential as well as the environment-dependent tight-binding method, making clear that the carbon allotrope is elastically stable and has a fairly high shear modulus [24]. Pang et al. made a systematical analysis of failure strength and strong anisotropic Poisson’s effect of ac_n–sp²–sp³ CHCs with different cell sizes via molecular dynamics simulation using optimized reactive empirical bond-order potential [19]. Gu et al. used molecular dynamics (MD) simulation with modified reactive empirical bond-order potential to give the stress–strain curves of symmetrical and asymmetrical ac₂–sp²–sp³ CHCs, ac₃–sp²–sp³ CHCs, and zz₂–sp²–sp³ and zz_2.5–sp²–sp³ CHCs [21]. Zhang et al. [16] built the in-plane compression and out-of-plane nanoindentation tests of ac_n–sp² and zz_m–sp²–sp³ CHCs via MD analysis with the Adaptive Intermolecular Reactive Empirical Bond Order potential (AIREBO) [25] potential. Meng et al. investigated the out-of-plane compression behaviors of both ac_n–sp²–sp³ and zz_m–sp²–sp³ CHCs using MD simulation with AIREBO potential [20]. Despite these studies, a full and clear understanding of the inner mechanism of the deformation of CHCs is still lacking but desirable due to its various promising applications as an engineering material. In this study, we used MD simulation with AIREBO potential to investigate the transformation from an ac₂–sp²–sp³ CHC to an ac₂–sp² CHC, which causes the change of tensile deformation behavior of ac₂ CHCs. In the meanwhile, we gave the stress–strain curves, Young’s modulus, and Poisson’s ratio of ac₂ CHCs, which show a great agreement with previous studies. We also built an out-of-plane nanoindentation test to investigate the plastic deformation behaviors of ac₂–sp²–sp³ CHCs, giving the hardness and Young’s modulus.

2. Crystal Structure of Carbon Honeycomb

Carbon honeycomb, just as its name implies, is a kind of tubular structure with a honeycomb-like pattern through a top–down perspective. The tube wall can be regarded as a graphene-like monolayer, thereby deemed to be a 3D version of graphene. One thing about its structure that remains to be discussed is whether three adjacent tube walls are bound up with each other via sp³ or sp² bonding (Figure 1b,c). Undoubtedly, atoms that build up the tube shells are combined with neighboring carbon atoms by an sp² bond, as is the case in graphene, whereas those who pull contiguous graphene-like walls together, which are called hinge atoms, give rise to more options, one kind of typical sp³ bonds with non-hinge atoms (tube atoms) or one variant of diamond-like sp³ bonds with both adjacent hinge and non-hinge atoms. In the case of sp² bonds, hinge atoms do not interact with each other, allowing for electron exchange with only adjacent tube atoms. In the case of sp³ bonds, it is not only nearby tube atoms but also nearest-neighbor hinge atoms that get bonded with the center hinge atom, resulting in a diamond-like atom-hinge that ties up three around the atomic walls. According to the carbon patterns in the plane perpendicular to the axis of cell channel, carbon honeycombs can be categorized into two sets: Armchair CHC (ac_n CHC) and zigzag CHC (zz_m CHC) [20], as shown in Figure 1. Furthermore, carbon honeycombs can be subdivided into sp² and sp²–sp³ CHCs. In zz_m CHCs, carbon atoms of graphene-like walls are bonded by sp² bonds, while hinge atoms can be bonded with adjacent atoms only by sp³ bonds. Therefore, zz_m CHCs always refer to zz_m–sp²–sp³ CHCs. However, the other case is different. In ac_n CHCs, hinge atoms can be bonded with atoms around them by sp² or sp³ bonds. Therefore, ac_n CHCs can be subdivided into ac_n–sp² CHCs and ac_n–sp²–sp³ CHCs. In the case of ac_n–sp²–sp³ CHCs, hinge atoms can be distributed symmetrically or asymmetrically in a crystal cell, according to which ac_n–sp²–sp³ CHCs can also be subdivided into a symmetrical one and asymmetrical one, as shown as Figure 1b,c,e,f.

Herein, we implemented a density functional theory (DFT) calculation for both primitive cells of sym-ac₃–sp² and sym-ac₃–sp²–sp³ CHCs, as shown in Figure 2b,c, using PWmat [26,27], to examine their stabilities. For both types, cell relaxation and atom relaxation were exerted in subsequence to determine their respective energy-minimized configurations. Subsequently, self-consistent calculation was carried out to optimize the electron interaction in both primitive cells, as visualized in Figure 1d,e. As for the sym-ac₃–sp²–sp³ CHC, stronger sp3 bonds between the hinge atom and nearest-neighbor tube atoms can be observed than those between two adjacent hinge atoms, the bond angle of which is measured as 103.35°. Meanwhile, exactly strong sp² bonds among graphene-like atoms can be distinguished from sp³ bonds, with denser electron density around tube atoms observed (red-yellow hexagons in Figure 1d). As for the sym-ac₃–sp² CHC, hinge atoms are not attached to each other, with weaker sp2 bonds between hinge and non-hinge atoms observed than those between only tube atoms. The DFT calculation results, i.e., that the cohesive energy is –154.19 eV/atom of the sym-ac₃–sp² CHC, –154.84 eV/atom of sym-ac₃–sp²–sp³ CHC, respectively, indicate that two types of configurations are both likely to exist stably. Additionally, the slightly stronger cohesive interaction implies the sym-ac₃–sp²–sp³ CHC structure is relatively more stable.

3. Molecular Dynamics Simulations

We used molecular dynamics simulations to do the investigations of mechanical behaviors of the sym-ac₃–sp² CHC and sym-ac₃–sp²–sp³ CHC. The reason we focused on them is that in the work of Krainyukova et al. [14], a periodic carbon structure identical to the sym-ac₃–sp² CHC according to their synthetic CHC samples, proving that ac₃ CHCs are more likely to exist in the nature. The method of MD simulations is well established to investigate the elastic and plastic behaviors of sym-ac₃–sp² and sym-ac₃–sp²–sp³ CHCs. Initially, we carried out uniaxial tension tests for both types of CHCs, simultaneously determining the Young’s moduli of the sym-ac₃–sp²–sp³ CHC to be 551 ± 4 GPa and of sym-ac₃–sp² CHC to be 542±4 GPa in the channel direction (Figure 1g,h), following which we carried out tension tests in the channel direction of two kinds of CHCs, giving stress–strain curves, respectively. Interestingly, a transformation of hinge atoms from sp³ to sp² was observed on the atomic scale, which can account for the yield stage that can be observed in the sym-ac₃–sp²–sp³ CHC, while cannot be observed in the sym-ac₃–sp² CHC. Furthermore, a nanoindentation simulation was devised and the load–displacement curve was plotted to determine the hardness and Young’s modulus of the CHC, which are 50.9 GPa of hardness and 461 ± 9 GPa of Young’s modulus, respectively. Meanwhile, we discussed the plastic behavior of CHC in the process of nanoindentation.

We used a molecular dynamics (MD) software, LAMMPS [28], which is an open-source code, to simulate two typical mechanical tests, uniaxial tension tests of the sym-ac₃–sp² CHC and sym-ac₃–sp²–sp³ CHC and nanoindentation of the ac₃–sp²–sp³ CHC. As for the tension simulation, given that we have set the boundary conditions as periodic in three directions, and thus, the mechanical response of the system is independent from the simulation scale, in order to accelerate the calculation, the simulation box was chosen as 3.408 nm × 2.951 nm × 4.181 nm, containing 2856 atoms. Both end-faces perpendicular to the direction of tension were freed from all kinds of forces. The whole system was minimized for energy and relaxed at specified temperatures (10 K and 300 K) and zero pressure in the NPT ensemble (where the number of particles, pressure, and temperature are constant) using AIREBO [25] potential (cutoff radius = 2.5 Å), to eliminate inner stress. The strain rate is set as 0.01 ps⁻¹. Considering that temperature effect will cause random distribution of atoms, which will prevent us from observing the bonding transformation of hinge atoms, in addition to at room temperature, uniaxial tension tests were also carried out at 10K, at which little thermal motion is allowed so that a perfect configuration might remain.

We also simulated the entire process of nanoindentation on an ac₃–sp²–sp³ CHC along the channel direction. A thick CHC plate measured as 17.82 nm × 15.43 nm × 11.60 nm was set up in a simulation box with a vacuum layer 9.00 nm thick in addition (Figure 2). The indenter is a cubic diamond sphere, 8.00 nm in diameter, with [001] direction aligned with the channel direction. The velocity of this indenter is set at 0.25 Å/ps. There is a 0.50-nm-thick layer on the bottom of the slab of CHC, all forces on which are zeroed. Nonperiodic conditions were selected in three directions in order to simulate an isolated system. Before the mechanical test, an energy minimization had been finished and the total ensemble had been relaxed efficiently without applied stress at 300 K, using an NVT ensemble (where the number of particles, volume and temperature are constant) and NPT ensemble with zero pressure successively, with AIREBO-morse [29] potential (cutoff radius = 2.5 Å).

4. Results and Discussion

4.1. Uniaxial Tension

Strain–Stress Curve

From uniaxial tension, we obtained the engineering strain–stress curves of both the sym-ac₃–sp² CHC and sym-ac₃–sp²–sp³ CHC (Figure 3), which can depict the elastic pattern directly and characterize the plastic properties indirectly, gaining an insight into the similarities and differences between the mechanical behaviors of two kinds of CHCs. Our results are in great agreement with the stress–strain curves from a previous study [21], in which the fracture strength and strain equal to our results are exhibited. Although the stress–strain curves of ac_n–sp² CHCs with different sizes have been studied systematically, there have been a limited number of studies of the yield process, i.e., the transformation from sp³ to sp², which is what we focus on. We also calculated the Young’s moduli of both CHCs, 542 ± 4 GPa for the sym-ac₃–sp² CHC and 551 ± 4 GPa for the sym-ac₃–sp²–sp³ CHC (300 K). From a previous study [16], the Young’s modulus of the CHC, which is assumed to be a nanoscale cell solid, can be calculated by the cell wall width via an analytic method. In our cases (wall width is determined as 7 Å), the analytical result is 560 GPa, in a great agreement with our simulation results.

When the tension starts, both CHCs are stretched quasilinearly under homogeneous and elastic deformation; herein, the initial curve of the sym-ac₃–sp²–sp³ CHC deviates a little from that of the sym-ac₃–sp² CHC. Then, in a certain strain range (from 0.072 to 0.132 for case at 10 K), a yield stage occurs in the deformation of the sym-ac₃–sp²–sp³ CHC, which indicates a plastic deformation has initiated. Later in this article, the mechanism therein will be discussed in detail. After the yield step, the sym-ac₃–sp²–sp³ curve behaves just like a sym-ac₃–sp² curve does, even breaking off at the same fracture strength, which is around 61 GPa (56 GPa for case at 300 K), with a retardation measured as around 0.06 in strain for the case at 10 K, equal to the length of the preceding yield stage. Thus, the mechanical behavior of the sym-ac₃–sp²–sp³ CHC can be distinguished from that of the sym-ac₃–sp² CHC by the fracture retardation caused by yield terrace mentioned above.

In order to illustrate the mechanism of this yield stage, a slice through some sheets of graphene-like walls has been exhibited in Figure 4. Both cases at 10 K and 300 K exhibit the same yield stage, so in order to avoid random distribution of atoms caused by thermal motion, which will prevent us from observing the detailed transfomation of hinge atoms, we explain the yield process based on the tension at 10 K. The color of atoms characterizes the channel-direction component of the stress on each atom, of which the colormap is linearly related to the value of stress. Through the variation of the channel-direction stress, we can justify whether the hinge atoms interact with other hinge atoms, so that we can determine the bond type of hinge atoms. For hinge atoms with an sp³ bond, the neighboring hinge atoms must be bonded with each other, thus exhibiting a high channel-direction stress; while for the case of sp², hinge atoms do not interact with each other, so there will be a low and even zero channel-direction stress on them. In addition, sp² bonding is in-plane (two-dimentional), while sp³ is three-dimentional; thus, we can also tell them apart by their geometry and nearest neighbors. Three snapshots at different states of strain, corresponding to three characteristic points in Figure 4, respectively, are displayed. At the beginning of tension (Figure 4a), hinge atoms share most of the distortion energy through the extension of the sp³ bond between two hinge atoms, while the sp³ configuration remains. However, when the engineering strain reaches around 0.072, some sp³ pairs of hinge atoms break up, the distortion energy of graphene-like atoms begins to increase gradually, and simultaneously, the yield stage starts. When the strain attains 0.132, almost all sp³ pairs are disconnected, most of the deformation is shared by the atomic cellular walls, and the yield stage terminates at the same time. By contrast, the sym-ac₃–sp² CHC, without sp³ interaction, does not experience such a transformation, and atomic tubes instead of hinge atoms bear most proportions of tensile force. Compared with graphene, although CHC has a lower fracture strength and modulus, it can obtain a higher fracture strain than graphene through this transformation, which may be applied in flexible coating.

4.2. Nanoindentation

At present, nanoindentation is commonly used for the research of mechanical properties of materials on the nanoscale [30,31,32]. There are two main reasons this methodology has been put into widespread usage [33]. As for nanoindentation, the stress applied by indenter is nonhomogeneous and can thus create an elastically physically isolated volume. Secondly, not only the location of elastic instability, but also the slip behaviors can be predicted via nanoindentation.

In a previous study [16], nanoindentation behaviors of different sizes of ac_n–sp² and zz_m–sp²–sp³ CHCs were reported. However, the indent depth was small (4 Å) and the unloading process was not displayed. Herein, we simulated the process of nanoindentation in the channel direction of the sym-ac₃–sp²–sp³ CHC, an indentation depth of about 40 Å, containing two steps, acting load and unloading, via LAMMPS [28], under the confinement mentioned in the section “Molecular Dynamics Simulations”. The result is displayed in Figure 5, which we used to determine the hardness and Young’s modulus of the CHC in the channel direction.

Nanoindentation hardness is the average pressure the material can support under an external load, defined as the ratio of indentation load P and projected contact area A_p [30]. From the load–displacement curve, the hardness H can be attained using the maximum load as:

$$\mathrm{H}=\frac{P_{max}}{A_p}.$$

(1)

In addition, Young’s modulus can also be calculated indirectly via reduced modulus E_r, which can be expressed as:

$$E_r=\frac{\sqrt{\pi }}{2\beta }\frac{S}{\sqrt{A_p}},$$

(2)

where S is termed as contact stiffness, which can be obtained as the slope at the tip of unloading curve, β is a constant that depends on the geometry of indenter, for a spherical indenter, β = 1. Furthermore, we can determine the substrate’s modulus using following relationship that reflects impacts of both indenter’s and substrate’s moduli:

$$\frac{1}{E_r}=\frac{1-{\nu }_i^2}{E_i}+\frac{1-{\nu }_s^2}{E_s},$$

(3)

where the subscript i refers to indenter, while s refers to substrate.

The nanoindentation hardness is determined as 50.9 GPa in

Table 1. Hardness and Young’s modulus calculated from nanoindentation. Calculated hardness on the nanoscale.

Projected Contact Area (Å2)	Maximum Load (nN/ Å)	Hardness (GPa)
3689	1878.1	50.9

and

Table 2. Hardness and Young’s modulus calculated from nanoindentation. Comparison between Young’s moduli from tension test and nanoindentation.

Contact Stiffness (nN/Å)	Reduced Modulus (GPa)	Young’s Moudulus (From Tension Simulation) (GPa)	Young’s Modulus (From Nanoindentation) (GPa)
228 ± 3	285 ± 4	551 ± 4	461 ± 9

, which is higher than quartz (9.25 GPa) and even some ion-beam-irradiation hardened metal (over 20 GPa) [34]. In order to analyze the plastic behavior of the ac₃–sp²–sp³ CHC in the process of nanoindentation, five snapshots are displayed in Figure 5. In the case of loading, the amounts of atoms which are identified as simple cube lattice are increasing with the dent produced by spherical diamond indenter deepening. This phenomenon indicates a transformation of covalent bonds from a previous normal type to a more compact one. Subsequently, although the indenter retreated and the dent recovering partially in an elastic way, the number of simple cubic atoms is maintained almost the same number as there is when the indenter penetrates the deepest part, which illustrates that a plastic deformation has occurred and cannot rebound.

5. Conclusions

We discussed the elastic and plastic behaviors of sym-ac₃–sp² CHCs and especially the sym-ac₃–sp²–sp³ CHC which is slightly more stable according to our DFT calculation and previous experimental study, using LAMMPS to carry out a typical uniaxial tension and nanoindentation simulation. In the case of uniaxial tension, we discussed the yield stage in the stress–strain curve of the sym-ac₃–sp²–sp³ CHC caused by a intriguing transformation of covalent bonds of hinge atoms from strong sp³ to comparatively weak sp², and we also determined the Young’s moduli at room temperature, 542 ± 4 GPa for the sym-ac₃–sp² CHC and 551 ± 4 GPa for the sym-ac₃–sp²–sp³ CHC, which are in great agreement with a previous analytical study. Then, in the discussion of nanoindentation, we determined the hardness (50.9 GPa) and Young’s modulus (461 ± 9 GPa) via the load–displacement curve and gave an insight into the plastic deformation caused by the indenter on the nanoscale, which can be explained by the permanent conversion of covalent bonds from sp³ type to sp² type. These findings provide an insight into the relationship between the covalent bond type of hinge atoms in carbon honeycomb and the plastic behavior of carbon honeycomb.