# Finite Element Modelling and Mechanical Characterization of Graphyne

Authors: Ricardo Couto; Nuno Silvestre
Journal: Journal of Nanomaterials
 (2016)

---

## Abstract

Graphyne is an allotrope of carbon with excellent mechanical, electrical, and optical properties. The scientific community has been increasingly interested in its characterization and computational simulation, using molecular dynamics (MD) simulations and density functional theory (DFT). The present work presents, for the first time (to the authors’ knowledge), a finite element (FE) model to evaluate the elastic properties of graphyne. After presenting a brief literature review on the latest developments of graphyne and its mechanical characterization through computational methods, the FE model of graphyne sheet is presented in detail and the calculation of its elastic properties described. The linear elastic properties (Young’s modulus, Poisson’s ratio, bulk modulus, and shear modulus) obtained from the proposed FE models are in general agreement with those previously obtained by other authors using more complex computational models (MD and DFT). The influence of van der Waals (vdW) interatomic forces on the linear elastic properties of planar graphyne is negligible and can be disregarded if small strain hypothesis is adopted. The FE models also show that graphyne exhibits marginal orthotropic behavior, that is, “quasi-isotropic” behavior, a fact that agrees with the conclusions reported by other researchers.


---

## Body

## 1. Introduction

Providing the basis for life in our planet, carbon is found in many known life forms. Its numerous hybridizations (sp, sp2, and sp3) exhibit strong chemical bonds and display a rich variety of arrangements. Besides the well-known crystalline forms, such as diamond and graphite, carbon also appears in numerous allotropes, namely, graphene, fullerenes, and nanotubes [1]. These allotropes have distinct mechanical [2, 3], electronic [4, 5], optical [6], and chemical [7] properties. 2D graphene is a monolayer of carbon atoms organized in a honeycomb lattice and is one of the strongest materials ever tested [8]. With the advent of nanotechnologies and following the discovery/synthesis of graphene [9], plenty of 2D carbon-based materials still are pursued, tested, and simulated, either experimentally or computationally. Graphene/boron nitride [10], graphene oxide [11], aluminum nitride monolayers [12], and hydrogenated graphene are some of these carbon-based materials. A recent form of 2D carbon-based material is known as graphyne, which is still hard to synthesize (the scientific community has conceived of graphene for many years before its isolation and characterization by Novoselov et al. [9]; up to 2004, it has likely been unknowingly produced in small quantities, through the use of pencils and other similar applications of graphite; before 2004, graphene was mainly investigated within the framework of computational tools; therefore, it is foreseen that graphyne will also be synthetized in the near future, as graphene was in 2004). This material has several distinct atomic structures. It is a one-atom thick material (like graphene), but it has acetylenic linkages (C-C triple bond) between 33% and 100% of all bonds. Among the most known configurations of graphyne are α-, β-, and γ-graphyne. The aforementioned linkages have great influence on its mechanical properties and heat conduction. Regarding electrical properties, graphyne displays a band gap in the semiconductor range and its band gap can be changed through mechanical strain [13]. Extensive studies have been proposed in order to better explore its properties and future applications, by either experimental test or computational simulations.Graphyne is not limited to single atomic bonds and hexagonal lattice exhibited by graphene. Because graphyne has acetylenic linkages, it gets weaker than graphene. And this scenario is even the worst in case of graphdiyne and graphyne-3 because they have a higher number of acetylenic linkages. Because of the scarce (limited) production of graphyne [14], researchers must rely mostly on computational methods to assess the graphyne mechanical properties (stiffness and strength). Recently, some authors applied density functional theory (DFT) calculations and performed molecular dynamics (MD) simulations [15–18] to study the mechanical behavior of several types of graphyne. To the authors’ best knowledge, there is not a single application of finite element (FE) method to characterize the mechanical behavior of graphyne. As stated by Peng et al. [16], FE models could be very useful in evaluating the mechanical behavior of graphyne. Therefore, the main goals of this paper are (i) to propose an original FE model of graphyne and (ii) to calculate its elastic properties. In summary, the following three objectives are assumed:(i)
To develop a consistent FE model to simulate the mechanical behavior of graphyne.
(ii)
To compare the results of the FE model with others available in literature and obtained from either MD or DFT calculations.
(iii)
To assess the influence of nonbonded forces (vdW forces) on the elastic properties of graphyne.Initially, the literature review about the applications and mechanical characterization of graphyne is presented in Section 2. Then, Section 3 describes the proposed FE model of a graphyne sheet (regarding the fact that many possible structures of graphyne could be assumed, γ-graphyne will be adopted for the present computational study; in Sections 3 and 4, we shall designate γ-graphyne simply as graphyne; hence, this work seeks to provide a new path to analyze this promising carbon allotrope; to the authors’ knowledge, this is the first work on the application of FE method to study graphyne mechanical properties). Several particular issues, such as the FE type, boundary conditions, and interatomic forces, are addressed in detail. This reference model enables the computation of graphyne’s elastic properties without the effect of vdW forces. Then, these weak forces are included and their influence is assessed. After that, graphyne sheet is subjected to several tests (uniaxial tests, shear test, and biaxial test) and the results of these calculations are presented in Section 4 and compared with MD and DFT values found in literature. This comparison enables the validation of the FE model. Finally, some conclusions regarding the proposed FE model are presented.


## 2. Literature Review on Graphyne

Graphyne, graphdiyne, graphyne-3, graphyne-4, and graphyne-5 are allotropes of carbon, all belonging to a family called graphyne-N. These structures are planar sheets of sp and sp2 bonds of carbon atoms arranged in a crystal lattice [14]. The main difference between these types of graphyne is the number of acetylenic links existing between the hexagonal lattices of carbon atoms (see Figure 1).Figure 1
Types of acetylenic links.Graphyne-N has 33% of the C-C bonds of graphene replaced by one acetylenic unit [19]. Although having a fixed amount of C-C bonds, there are several structures of graphyne. The most known are α-, β-, and γ-graphyne and 6,6,12-graphyne (see Figure 2). Through the performance of both MD and DFT analyses, this allotrope has shown good mechanical properties in terms of Young modulus, fracture, and strength [15, 17, 18]. According to Peng et al. [14], only small fragments of graphyne have been synthesized up to now. Due to their limited amount and dimension, such samples were not studied in detail to find their properties (mechanical, electrical, optical, and thermal) experimentally. Thus one must rely mostly on theoretical and computational data. Graphyne’s mechanical properties suggest it can be used as nanofiller in certain composite materials, allowing not only the reduction of strength but also the increase of ductility.Figure 2
Atomic structures of (a) graphene, (b)α-graphyne, (c) β-graphyne, (d) γ-graphyne, and (e) 6,6,12-graphyne [20].
(a)
(b)
(c)
(d)
(e)Regarding its electronic properties, it has been shown through first principle calculations that graphyne ribbons with finite width have band gaps in the semiconductor range [21]. The band gap can be altered if mechanical strain is applied [13]. However, structural modifications of graphyne through the adsorption on Cu, Ni, and Co surfaces help modulate the energy gap, as was found by Lazić and Crljen [5] applying a DFT method with the novel nonlocal vdW-density functional correlation (vdW-DF). These authors also found that the binding of graphyne to the abovementioned metal surfaces is stronger than that of graphene, which makes graphyne more suitable for the building of future nanoelectronic components. This stronger binding to substrates leads to a lower sensitivity of graphyne to wrinkling phenomenon when the substrates are deformed. Given the performed computational studies that foresee a tunable band gap, the future looks bright in terms of electronic applications of graphyne. Graphyne semiconductor properties will be determined by the sheets size upon fabrication. Regarding its large elastic strain region coupled with its electronic properties, we may predict the future implementing of graphyne on sensors that can be useful in a wide range of possibilities, including temperature sensing. From a MD perspective, Hu et al. [22] found that graphyne nanotubes (GNTs) have an unprecedentedly low thermal conductivity, even lower than any values displayed by other carbon nanotubes (CNTs), being defected, doped, or chemically functionalized. This low thermal conductivity persists even if the tubes are infinitely long and/or have large diameters, up to 18 nm. Such discovery may be valuable for future applications, such as thermoelectrical devices for energy conversion. On the branch of optical applications, Bhattacharya et al. [6] recently found that, through the systematical substitution of boron (B) and nitrogen (N) into the structure of the graphyne family, this allotrope displayed very particular optical properties, compared to other carbon-based structures. Due to the presence of BN, the optical band gap is tuned, and BN doped graphyne exhibits a strong absorption peak in a wide UV-region for perpendicular polarization. Hence, graphyne shall also find applications in future UV light protection.After the development of reactive empirical bond order (AIREBO) potential for hydrocarbons [24], the mechanical, electrical, and vibrational properties of both graphene and graphyne started to be able to be analyzed in depth. Shao et al. [25] performed DFT calculations to compare the mechanical behavior of graphene and graphyne and to find a resembling way of reacting to the change of temperature on both materials. This research also found that the closer to 0 K the carbon allotrope got, the higher its mechanical properties become. Such discovery, even if only through computational methods, may in the future be a great contribute to aeronautics, as long-range aircrafts spend most of their flight time at negative temperatures (in Celsius degrees range). This field seems to have a great potential for new developments. Much like graphene, graphyne can be used to desalinize water, as Xue et al. [26] predicted through MD simulations and first principles modeling. Xue et al. [26] also found a 100% rejection of ions in sea water, using pristine graphyne. Therefore, graphyne shall provide a new solution to alleviate the shortage of fresh water around the globe.Focusing now on graphdiyne (see Figure3), which was first predicted by Haley et al. [27], it is the carbon allotrope with poor mechanical properties, with far less interesting elastic properties than those of graphyne or graphene [17]. Its stiffness is nearly 30% that of graphyne with a noticeably reduced strength according to Cranford and Buehler [28] and Pei [29]. The less stiff mechanical behavior of graphdiyne is due to the larger acetylenic linkage that results in weaker bonding. However, given its other interesting properties, graphdiyne is considered separately from the rest of the graphyne family. As for electronic properties, graphdiyne can have its band gap varied proportionally to externally applied strain, similarly to graphyne, though the band gap is quite different [29]. Graphdiyne is a direct gap semiconductor with its band gap equal to 0.47 eV. It was found that band structure of graphdiyne is tunable by a uniform strain, even if the graphdiyne remains in its own plane and does not undergo out-of-plane deformations. The band gap increases with increasing strain value, which results from the decreased orbital overlap between C atoms when strain increases. Although such results were obtained through computational simulations, graphdiyne has already been synthetized successfully in very limited samples [30, 31].Figure 3
Chemical structure of graphdiyne [27].After it was synthesized on copper substrates, the resulting graphdiyne-based composite material was analyzed through full atomistic MD in order to find its mechanical properties [2]. These authors concluded that graphdiyne can enhance Cu performance in a synergetic way, by combining Cu with layers of graphdiyne. Due to the confinement of the Cu crystal during deformation, the failure strain of graphdiyne-based composite also increased substantially. Another significant difference between graphyne and graphdiyne was due to Yue et al. [32]. They found out that graphdiyne is impermeable to water and ions due to its small pore size. Thus, graphdiyne should not be used either as a desalinator or as a mean to purify water. Since graphdiyne shows evidence of low effective electron mass [29], the electrons should move fast, even if a small potential field is applied. Therefore, one of the limitations of silicon based transistors is outclassed by graphdiyne. Such limitation lies in the fact that the maximum switching frequency on current transistors is limited by the velocity of the electrons across the transistor channel relative to the length of the channel. In terms of sensors, and knowing that graphdiyne has a desired electromechanical coupling, we are not far from glimpsing its application on several mechanically and electronically driven systems, for example, thermal control systems. Regarding the abovementioned doping of graphdiyne using boron and nitrogen, patterned sheets with different properties can be fabricated and be used to control large amounts of current flow by harnessing the electrons flow through the sheets.


## 3. Mechanical Characterization of Graphyne

In the previous section, the general behavior and possible applications of graphyne were highlighted. The present section is directed towards the characterization of mechanical behavior of graphyne. Although being weaker than graphene, the mechanical properties of graphyne always attracted the interest of many researchers [14]. The acetylenic links present in graphyne are mostly responsible for its mechanical, thermal, electrical, and optical properties. Although the mechanical properties of graphyne might be dependent on the temperature, this thermal dependency is yet to be fully understood. Shao et al. [25] studied the temperature-dependent elastic properties and ultimate strength of graphyne using first principles calculations combined with quasiharmonic approximations. They used a carbon lattice containing 24 atoms and the range of temperatures employed was in the range of 0 K–1000 K. Unlike many materials, the thermal expansion coefficient (TEC) of graphyne at low temperatures is negative, being low in the range of 0 K–350 K. The TEC values increase within the remaining range of temperatures. Regarding the elastic properties, Shao et al. [25] obtained values for Young’s modulus (E = 250.9 N/m) and ultimate tensile strength (σ
U
T
S = 81.2 GPa) at room temperature (graphyne is a material of discrete (atomic) nature; hence its thickness cannot be considered conceptually; this is the reason why the elastic modulus is measured in N/m; however, a thickness can be given to the graphyne sheet, thus turning the unit of the elastic modulus into GPa). Regarding the influence of temperature on these properties, Young’s modulus suffered a loss of 10.79% (from 0 K to 1000 K) while the ultimate tensile strength decreased 12.1% (Shao et al. [25] obtained E = 350.0 N/m and σ
U
T
S = 119.2 GPa for graphene at room temperature; varying the temperature from 0 to 1000 K, graphene was capable of keeping its superior mechanical properties: E and σ
U
T
S decreased by 2.2% and 4.03%, resp.) (from 0 K to 1000 K).Zhang et al. [33] performed a somewhat similar study to obtain the influence of temperature on the strain-rate and fracture strength for α-, β-, γ-, and 6,6,12-graphynes. For comparison purposes, a graphene sheet was also simulated (see Figure 2). With such a variety of graphyne structures, they were able to understand the overall influence of the amount of acetylenic links on the mechanical properties of each graphyne. Square sheets with 20 nm side and periodic boundary conditions were adopted and AIREBO potential [24] was used in MD simulations. The results, though using a different method, were in agreement with those found by Shao et al. [25]: they found that both fracture strength and Young’s modulus decrease with increasing temperature. They also concluded that the erosion of graphyne’s mechanical properties with increasing temperature was higher than those of graphene. Regarding the comparison between different graphynes, it can be mentioned that the increasing amount of acetylenic linkage causes a faster decay in the mechanical properties of graphyne. It means that the mechanical properties of α-graphyne decrease the most for increasing temperature. In their simulations, Zhang et al. [33] adopted a range of temperatures between 1 K and 1200 K. In this range, Young’s moduli of graphene decreased only by 10% while the graphyne’s moduli dropped from 37.4% (γ-graphyne) to 76.7% (α-graphyne). Similar trends were found for fracture stress and strain. Therefore, α-graphyne is the weaker structure of graphyne family [33].Very recently, Hou et al. [34] developed a molecular mechanics model to study the elastic properties of graphyne-N. The ultimate objective of this work was to compare their model results with others, either using MD [15, 18, 28] or DFT [10, 17, 32]. Hou et al. [34] studied the in-plane stiffness of 10 distinct structures of graphyne and intended to capture the degradation of this property with the increasing number of acetylenic linkages. They also assessed the shear stiffness and Poisson’s ratio. Having acquired their results from the developed equations, which they claim to be capable of linking macroscopic properties with graphyne-N atomic structure in a direct manner, they concluded that their predicted results matched the ones available in literature (MD and FDT simulations). As previously shown by Yang and Xu [17] and later confirmed by Yue et al. [32], the in-plane stiffness decreases with increasing number of acetylenic links present in graphyne. Hou et al. [34] obtained a similar trend of results but further extended the calculations up to graphyne-10. Regarding the shear stiffness, these researchers achieved a trend similar to that obtained for in-plane stiffness. It should be highlighted that the shear stiffness and Poisson’s ratio of graphynes are scarcely studied properties but, by no means, less important. Bearing this in mind, the present work will provide some new results related to both properties. Hou et al. [34] also obtained results for Poisson’s ratio of graphyne-N. Although the trend is realistic (Poisson’s ratio increases with the number of acetylenic linkages), the resultsper si require further studies.Also very recently, Asadpour et al. [23] studied the mechanical properties of γ-graphyne and two analogous systems of boron nitride (BN) sheets through DFT analysis. The three analyzed structures were (i) sheet of γ-graphyne, (ii) sheet of BN, called BN-yne, and (iii) sheet composed by hexagonal rings of BN linked by acetylenic links (graphyne-BN sheet). An image of such structures can be seen in Figure 4, where the gray spheres represent the carbon atoms while the blue and pink spheres represent the boron and nitride atoms, respectively. Through their simulations, they found that graphyne has the highest in-plane stiffness (E = 190.7 N/m) among the three analyzed structures (Young’s modulus of BN-yne and graphyne-BN were E
=
162.4 N/m and E
=
163.0 N/m, resp.). Regarding Poisson’s ratio, graphyne-BN had the highest value (nearly 0.50) while graphyne had the smallest one (0.41). Regarding the shear modulus and bulk modulus, graphyne displayed the highest values: the bulk modulus was K
=
122.7 N/m and the shear modulus was G
=
77.0 N/m. For these properties, the BN-yne structure displayed the smallest values. Finally, it was found by Asadpour et al. [23] that graphyne undergoes plastic deformation for high tensile loads.Figure 4
(a) Graphyne structure, (b) BN-yne structure, and (c) graphyne-BN structure [23].
(a)
(b)
(c)Still being limited, the number of studies on the mechanical properties of graphyne has been growing in recent years. Although a limited range of properties have been studied, it is well known that computational methods (MD, FDT, and molecular mechanics) are capable of filling this gap. Besides its mechanical properties, recent years have shown a renewed interest of the scientific community in graphyne’s electrical [21, 31, 35–38] and optical [6, 13, 14] properties. Graphyne-BN is also being increasingly studied nowadays. While graphyne is still hard to synthesize, much can be expected from this allotrope.


## 4. Finite Element Modeling of Graphyne

Computational modelling of carbon nanomaterials is very useful because the assessment of their mechanical behavior needs rigorous experiments at nanoscale, which still are complex and too expensive. Among the computational methods, MD and DFT proved to be the most effective. They have been applied to evaluate the mechanical behavior of graphene, fullerenes, and graphyne. However, these methods require intensive computations, as they depend not only on the electronic structure of atoms (ab initio, DFT) but also on several thermodynamic variables (MD), such as pressure, temperature, and energy. In order to avoid such computational efforts and time consuming analyses, the finite element (FE) method has been applied to assess the mechanical behavior of carbon structures at nanoscale. The application of FEs to investigate the behavior of fullerenes (tubes and buckyballs) and graphene has been very successful in the last decade [39]. The FE models are divided into two main groups: (i) beam FE models, in which the atomic interactions are simulated through beam (1D) elements distributed within the hexagonal lattice of the 2D molecular structure, and (ii) shell FE models, in which the atomic interactions and hexagonal lattice are taken into account directly in the constitutive laws of shell (2D) elements that form the 2D molecular structure. Because the present work falls within the scope of the first group, we will limit our literature review to beam FE models. Since the publication of the seminal works by Odegard et al. [40, 41], much has been done in the development of beam FE models for nanotubes and graphene. Regarding graphene, excellent models have been proposed by Scarpa et al. [42], Georgantzinos et al. [43], Chandra et al. [44], Rouhi and Ansari [45], Theodosiou and Saravanos [46], and Zhang et al. [47]. However, FE method has not yet been applied to study the mechanical behavior of graphyne, even if for its simple linear elastic behavior. The main objective of the present paper is to present a FE model of graphyne.To achieve this goal, a graphyne sheet is modelled to study its linear elastic properties. The sheet resembles that investigated by Cranford and Buehler [15] using MD: it has a square shape with each side being approximately 10 nm wide. The FE software chosen was Ansys [48]. As is known, graphyne has three different types of covalent bonding: aromatic, single, and triple. The first type (aromatic bonds) occurs inside the hexagonal molecules of graphene, while the second (single) and third (triple) types occur on the acetylenic links forming the γ-graphyne. These three types of C-C bonds can be viewed in Figure 5.Figure 5
Graphyne cells with (a) different interatomic bonds and (b) bond-angles.
(a)
(b)To simulate such bonds, finite element Beam4 [48] was chosen. This 2-node element was proven useful as it is a uniaxial element with tension, compression, torsion, and bending deformations, having 6 degrees of freedom per node (3 translations and 3 rotations in X, Y, and Z directions). In order to enable accurate modelling of graphyne sheet, the lengths of the atomic bonds were assumed to be the same as those adopted by Cranford and Buehler [15]. The length of each bond (equilibrium distance between atoms) is as follows: (i) aromatic (L = 1.49 Å), (ii) single (L = 1.48 Å), and (iii) triple (L = 1.19 Å).To obtain the covalent forces being exerted by the interatomic bonds, we follow the formulation proposed Odegard et al. [40, 41], which was successfully implemented to evaluate the stiffness of CNTs and graphene. The total molecular potential energy of the molecular model is given by(1)
V
=
∑
K
ρ
ρ
-
P
2
+
∑
K
θ
θ
-
Θ
2
+
∑
D
I
J
1
2
ρ
I
J
ρ
12
-
ρ
I
J
ρ
6
,where (i) P and Θ are the undeformed interatomic length of each bond and the undeformed bond-angle, respectively, (ii) ρ and θ are the distance and bond-angle after stretching and angle variance, respectively, (iii) K
ρ and K
θ represent the force constants associated with the stretching and angle variance of each bond and bond-angle, respectively, and (iv) D
I
J and ρ
I
J are the well depth and natural vdW distance for interaction between atoms I
J. As was stated by Odegard et al. [41], only the bond stretching, bond-angle variation, and vdW parameters were considered in (1) because other energy terms have a negligible contribution to the total molecular potential energy. In order to study the influence of vdW forces in the elastic properties of graphyne, we consider two models: (i) reference model (RM), which considers covalent bonding forces and disregards vdW forces, and (ii) van der Waals model (vdWM), which considers both covalent bonding forces and vdW forces. Therefore, RM considers the first two terms of (1) while vdWM considers the three terms. Odegard et al. [40] proved that the axial stiffness representing primary bonds and the bond-angle variance interactions may be determined as a function of the force constants. In order to calculate Young’s modulus for bond-stretching interaction, the following formula is used:(2)
E
α
A
=
2
K
ρ
Lwhich depends on bond length (L), cross-sectional area (A), and the stretching force constant (K
ρ). According to Odegard et al. [40], this constant has a value of K
ρ
=
469 Kcal/mol/Å2. Besides the C-C bond-stretching interaction, the bond-angle variation interaction also should be taken into account when calculating bond-angle interaction Young’s Modulus, given by(3)
E
β
A
=
32
K
θ
L
sin
⁡
Θ
2
2
.This equation is applied to calculate bond-angle interaction Young’s Modulus for each one of the three covalent bonds. In this equation, and besides the bond length (L) and cross-sectional area (A), it also depends on the undeformed bond-angle (Θ) and the angle variation constant, which is K
θ
=
63 Kcal/mol/rad2, according to Odegard et al. [41]. It should be noted that Θ is purely dependent on the angle between two adjacent covalent bonds, acting on the same carbon atom (see Figure 5(b)).It is now important to mention that both (2) and (3) give the axial stiffness (EA) of a given C-C bond. Therefore, Young’s modulus value depends on the value of cross section area while the opposite is also true (the value of A depends on the value of E). In conclusion, what really matters is the axial stiffness value (EA) and not the individual values of E and A. Therefore, the cross section area A = 1 Å2 is assumed and the 2nd moment of area of the cross section (circular) is I = 1/4π Å4. The cross section properties for aromatic, single, and triple covalent bonds implemented in FE model are given in Table 1.Table 1
Geometrical properties of beam cross section and Young’s modulus for aromatic, single, and triple bonds.
Bound type
Aromatic
Single
Triple
A [Å2]
1.0
1.0
3.0
I [Å4]
0.07958
0.07958
0.71620
E
α [nN/Å2]
97.109
96.457
25.852
E
β [nN/Å2]
70.507
70.507
0With the geometrical properties calculated, we can now determine Young’s moduli to simulate the covalent forces being exerted by the interatomic bonds. This is done by using (2) and (3). The values of (i) bond-stretching interaction Young’s Modulus E
α and (ii) bond-angle interaction Young’s Modulus E
β are given in Table 1. Note that, because of their alignment (orientation with θ = 180°—see Figure 5(b)) in acetylenic links leads, the triple bonds have E
β = 0.Having finished the inputs of geometrical and mechanical properties of covalent bonds, we now proceed to the spatial modeling of the graphyne sheet. Figure6(a) shows the beam element mesh in the RM. Note that (i) each aromatic bond in a hexagon was modelled with a single beam FE while (ii) each acetylenic link between two adjacent hexagons comprises three beam FEs, one corresponding to the triple bond and two corresponding to the simple bonds (see Figure 5(a)). In the FE mesh, the global x-axis (horizontal) is aligned with thearmchair direction of the graphyne sheet while the global y-axis (vertical) is aligned with thezigzag direction. All the displacements in the z-axis direction (out-of-plane) are set null.Figure 6
(a) Reference model (RM) without vdW forces and (b) model with vdW forces (vdWM).
(a)
(b)Regarding the vdWM (model with vdW forces), other FEs should be added to the RM (reference model) to simulate the effect of these forces, that is, to model the interaction of a carbon atom with its neighbor (nonbonded) atoms. The FE chosen is designated by Link180 [48]. This is a tridimensional truss element used to simulate uniaxial tension-compression forces. It has two nodes and three degrees of freedom per node (displacements in x-, y-, and z-axes).Regarding the calculation of mechanical properties, the bar axial stiffness due to vdW interaction is also obtained from the formulation by Odegard et al. [41], following the well-known 12-6 potential of Lennard-Jones: (4)
E
w
A
=
2
ρ
I
J
D
I
J
ρ
-
ρ
I
J
2
1
2
ρ
I
J
ρ
12
-
ρ
I
J
ρ
6in which ρ is the deformed interatomic distance, D
I
J is the well depth of interaction between carbon atoms I and J, and ρ
I
J is the vdW equilibrium distance for interaction between carbon atoms I and J. In case of carbon atoms, Odegard et al. [41] stipulate that D
I
J
=
0.07 Kcal/mol and ρ
I
J
=
3.55 Å.Because the vdW forces vary nonlinearly withρ (interatomic distance) and the hypothesis of small deformations is considered, we also assumed the value of ρ as being equal to the interatomic distance before the loading is applied to the sheet. Figure 7 shows the truss bars simulating vdW forces acting on carbons hexagons (Figure 7(a)) and vdW forces acting on carbon atoms of acetylenic linkages (Figure 7(b)). Special attention should be taken to the value of ρ as it is the maximum distance a bond of this kind can take; once surpassed, the vdW interaction becomes negligible. This explains why there are pairs of atoms that are not connected by bars (marginal effect of vdW forces). A closer look at Figure 7 shows that there are three relevant interatomic distances (bar lengths) in terms of calculation of the vdW interaction Young’s modulus E
w: ρ
=
2.5807 Å is the interatomic distance between neighbor aromatic bonds (hexagon atoms; see Figure 7(a)), ρ
1
=
2.572 Å is the interatomic distance between neighbor simple bonds (Figure 7(b)), and ρ
2
=
2.970 Å is the interatomic distance between neighbor triple bonds (Figure 7(b)). Again assuming A = 1 Å2, the values of E
w (Young’s modulus for vdW interactions) calculated from (4) are E
w
=
0
.
5945 nN/Å2 (between aromatic bond atoms), E
w
·
p
=
0
.
1372 nN/Å2 (between simple bond atoms), and E
w
·
c
=
0
.
6132 nN/Å2 (between triple bond atoms). The FE mesh of vdWM is depicted in Figure 6(b). Both RM and vdWM comprise 2830 carbon atoms. The graphyne sheet displays a total length L
i = 10.28 nm (x-axis direction) and height H
i = 10.14 nm (y-axis direction). This was the graphyne shape closest to the square shape of the MD model studied by Cranford and Buehler [15].Figure 7
Bars simulating vdW forces between (a) aromatic bond atoms and (b) acetylenic linkage atoms.
(a)
(b)


## 5. FE Analyses

The main aim of this work is to extract the elastic properties of graphyne using FE models. Because a planar sheet of graphyne is used, it is possible to determine the six elastic properties: Young’s moduli (E
x, E
x), Poisson’s ratio (ν
y
x, ν
x
y), shear modulus (G), and bulk modulus (K). In order to achieve this goal, the graphyne sheet should be submitted to four tests shown. However, in each test there are different procedures regarding the imposition of boundary conditions and applied loading, so they will be described separately. The boundary conditions and loadings applied in each test are schematically shown in Figure 8 for the case of RM (equal to vdWM). The four tests are the following:(i)
Uniaxial tensile test in the armchair direction (x-axis) to determine the elastic modulus E
x and Poisson’s ratio ν
y
x (see Figure 8(a)): unit horizontal forces (1 nN) were applied to all the 82 atoms on the right edge of the sheet (F
x = 82 nN), while the atoms on the left edge were restrained from displacement in x-axis but allowed to be displaced in y-axis direction. To avoid any undesirable rigid-body vertical displacement, the atoms at mid-height of the right and left edges of the sheet were prevented from being displaced along y-axis. In this test, the following displacements were extracted from the analyses: d, the displacement of the right edge of the sheet in the x-axis direction, and a, the displacement of top and bottom edges of the sheet in y-axis direction.
(ii)
Uniaxial tensile test in the zigzag direction (y-axis) to determine the elastic modulus E
y and Poisson’s ratio ν
x
y (see Figure 8(b)): unit vertical forces (1 nN) were applied to all the 83 atoms on the top edge of the sheet (F
Y = 83 nN) and the atoms on the bottom edge were prevented from being displaced in y-axis but allowed to be displaced in x-axis direction. The atoms at midwidth of the top and bottom edges of the sheet were fixed against displacements along x-axis, in order to avoid any undesirable rigid-body movement of the sheet. From this analysis, the following displacements are given: c, the displacement of the sheet top edge in y-axis direction, and b, the displacement of right and left edges of the sheet in x-axis direction.
(iii)
Biaxial tensile test (X
Y-plane) to calculate the bulk modulus K (see Figure 8(c)): unit horizontal forces (1 nN) were applied to all the 82 atoms on the right and left edges of the sheet (F
x = ±82 nN), and unit vertical forces (1 nN) were applied to all the 83 atoms on the top and bottom edges of the sheet (F
Y = ±83 nN). Despite the static equilibrium of the sheet, both x- and y-axis displacements of the carbon atom in the sheet center were prevented in order to exclude any possibility of rigid-body movement during the simulation. From this analysis, one must get L
f, length of the deformed sheet L
f in the x-axis direction, and H
f, height of the deformed sheet in the y-axis direction.
(iv)
Shear test (X
Y-plane) to calculate the shear modulus G (see Figure 8(d)): unit vertical forces (1 nN) were applied to all the 82 atoms on the sheet right (downward, F
Y = −82 nN) and left (downward, F
Y = +82 nN) edges, and unit horizontal forces (1 nN) were also applied to all the 83 atoms on the sheet top (rightward, F
X = +83 nN) and bottom (leftward, F
X = −82 nN) edges (this system of forces results in a clockwise moment of 83 nN × 10.14 nm = 841.62 nN·nm and an anticlockwise moment of 82 nN × 10.28 nm = 842.96 nN·nm, which gives a marginal value for the moment 1.34 nN·nm). The atoms in the center of each edge were prevented from being displaced perpendicularly to its edge. From this analysis, one must get α
X, the angle of rotation of horizontal (top and bottom) edges of the deformed sheet, and α
Y, the angle of rotation of vertical (right and left) edges of the deformed sheet.After running the four tests (analyses) for each model (RM and vdWM), the deformed shape of the graphyne sheet in each test is shown in Figure 9 for RM (vdWM is similar) and the numerical results are presented in Table 2 for both RM and vdWM: (i) a, b, c, and d (uniaxial tests), (ii) L
f and H
f (biaxial test), and (iii) α
X and α
Y (shear test).Table 2
FE test results of RM and vdWM.
FE test
Results
Reference model (RM)
vdW model (vdWM)
Uniaxial tensile (x-axis)
a [Å]
0.7562
0.74988
d [Å]
3.6405
3.6165
Uniaxial tensile (y-axis)
b [Å]
0.7789
0.77314
c [Å]
3.8781
3.8539
Biaxial tensile
H
f [Å]
104.008
103.992
L
f [Å]
105.094
105.082
Shear
α
X [rad]
0.11471
0.11390
α
y [rad]
0.11452
0.11372Figure 8
Boundary conditions and applied loading in (a) uniaxial tensile test in the armchair direction (x-axis), (b) uniaxial tensile test in the zigzag direction (y-axis), (c) biaxial tensile test (X
Y-plane), and (d) shear test (X
Y-plane).
(a)
(b)
(c)
(d)Figure 9
Initial and deformed shapes of graphyne sheet (RM): (a) uniaxial tensile test in armchair direction, (b) uniaxial tensile test in zigzag direction, (c) biaxial tensile test, and (d) shear test.
(a)
(b)
(c)
(d)


## 6. Elastic Properties of Graphyne

In this section, the elastic properties of graphyne are obtained and discussed. Using the test results and a set of basic elasticity equations, the elastic properties can be obtained:(i)
Elastic modulusE
x and Poisson’s ratio ν
y
x (uniaxial tensile test in the armchair direction, x-axis) are(5)
E
x
=
σ
x
ε
x
=
F
x
/
H
i
d
/
L
i
=
F
x
L
i
d
H
i
ν
y
x
=
ε
y
ε
x
=
-
2
a
/
H
i
d
/
L
i
=
-
2
a
L
i
d
H
i
.
(ii)
Elastic modulusE
y and Poisson’s ratio ν
x
y (uniaxial tensile test in the zigzag direction, y-axis) are(6)
E
y
=
σ
y
ε
y
=
F
y
/
L
i
c
/
H
i
=
F
y
H
i
c
L
i
ν
x
y
=
ε
x
ε
y
=
-
2
b
/
L
i
c
/
H
i
=
-
2
b
H
i
c
L
i
.
(iii)
Bulk modulusK (biaxial tensile test, X
Y-plane) is(7)
K
=
σ
-
Δ
=
σ
x
+
σ
y
/
2
A
f
-
A
/
A
=
1
/
2
F
x
/
H
i
+
F
y
/
L
i
L
f
H
f
-
H
i
F
i
/
H
i
F
i
=
L
i
F
x
+
H
i
F
y
2
L
f
H
f
-
H
i
F
i
.
(iv)
Shear modulusG (shear test, X
Y-plane) is(8)
G
=
τ
x
y
γ
=
1
/
2
F
x
/
L
i
+
F
y
/
H
i
α
x
+
α
y
.In these formulas, (i) σ
x and σ
y are the normal stresses in x- and y-axis directions, respectively, (ii) ε
x and ε
y are the extensions in x- and y-axis directions, respectively, (iii) σ
- is the bulk stress (average in-plane stress), (iv) Δ is superficial extension of the sheet, (v) A and A
f are the initial and final sheet area, respectively, and (vi) τ
x
y and γ are the shear stress and strain (sheet distortion), respectively. It should be noted that the stresses σ
x, σ
y, σ
-, and τ
x
y, are given in nN/nm (or N/m) simply because the thickness of graphyne sheet has no physical meaning. Graphyne (like graphene) is a discrete atomic structure, not a continuum plate; thus the familiar meaning of “thickness” does not hold true. Therefore, the values of elastic modulus are given in N/m (Yang and Xu [17]), not in N/m2. If a similarity between 2D nanomaterials, like graphyne and graphene, and continuum plates is assumed, then a thickness value of 3.2 Å is usually adopted (equilibrium distance between sheets)—this value comes from adhesion energy results, presented in Cranford and Buehler [15]. With this “approximation,” the modulus values can be calculated in N/m2 (or Pa) by distributing (i.e., dividing) the forces F also by the thickness t = 3.2 Å.Finally, using the data presented in Table2 and (5)–(8), the elastic properties of graphyne were calculated. All the results are displayed in Table 3 for both RM and vdWM, as well as those elastic properties values obtained by other researchers using either MD or DFT simulations. Cranford and Buehler [15], Zhang et al. [18], Yang and Xu [17], and Asadpour et al. [23] obtained values of elastic moduli and Poisson’s ratio. To the authors’ knowledge, the study by Asadpour et al. [23] is the only one available for the bulk modulus and shear modulus assessment.Table 3
Elastic properties of graphyne.
Elastic property
RM
vdWM
Cranford and Buehler [15]
Zhang et al. [18]
Yang and Xu [17]
Asadpour et al. [23]
Units
E
x
228.4 (713.7)
229.9 (718.5)
201.4 (629.4)
(508.0)
~155.0
190.69
N/m (GPa)
ν
y
x
0.421
0.421
0.410
—
E
y
208.5 (651.6)
209.8 (655.7)
247.0 (772.0)
(505.0)
~150.0
—
N/m (GPa)
ν
x
y
0.396
0.396
—
—
—
—
—
K
164.5
165.5
—
—
—
122.7
N/m
G
70.5
71.0
—
—
—
77.0
N/mAfter having obtained the elastic properties of graphyne using FE models, some discussions regarding these values should be undertaken on the following issues:(i)
Influence of vdW forces on the elastic properties of graphyne.
(ii)
Orthotropic behavior of graphyne.
(iii)
Validation of FE models of graphyne.The first main discussion regards the influence of nonbonded interatomic forces on the elastic properties of graphyne. From the comparison between the elastic properties (E
x, E
y, ν
y
x, ν
x
y, K, G
x
y) obtained from RM (reference model only considers bonded interatomic forces) and vdWM (both bonded interatomic forces and vdW interatomic forces are considered), it may be easily concluded that the influence of vdW interatomic forces on the elastic behavior of graphyne is negligible: the maximum difference in results is about 0.6%. Regardless the direction we may choose, the effect of vdW forces is very weak and accounts only for 1 N/m (difference in elastic moduli). As for Poisson’s ratio, bulk modulus, and shear modulus, their influence is even smaller. However, it should be recalled that our results assumed that the bar axial stiffness did not depend on its extension (vdW forces vary linearly with the interatomic distances) because we assumed small strain hypothesis. If we wanted to calculate the nonlinear elastic properties of graphene, we would have to relax the small strain hypothesis. Additionally, when studying the behavior of highly bended graphyne sheets (out-of-plane behavior) or even the behavior of multilayer graphyne, the consideration of vdW interatomic forces is mandatory because neighbor (nonbonded) atoms stay in equilibrium only due to vdW forces. In summary, it can be concluded that the consideration of vdW interatomic forces can be disregarded in the modelling of graphyne in-plane behavior (small strain). Because of this evidence, we will now proceed by describing the results obtained for RM only.The observation of the results shown in Table3 shows that we are in the presence of an orthotropic material: the values of Young’s modulus in armchair (x-axis) and zigzag (y-axis) directions are not equal. The idea that the armchair direction of graphyne is stiffer than its zigzag one (E
x
>
E
y) should also be highlighted. Their difference is roughly 20 N/m, which corresponds to nearly 10%. Comparing graphyne Young’s moduli with those of graphene (according to Yang and Xu [17], 325 N/m in armchair direction and 315 N/m in zigzag direction), we conclude that graphyne is much more flexible (less stiff) than graphene. This well-known characteristic is a result of the flexibility of acetylenic links, which exist in graphyne and are absent in graphene. However, it can also be stated that graphyne is “more orthotropic” than graphene as the ratio E
x
/
E
y is higher for graphyne than for graphene. Regarding Poisson’s ratio, it can be concluded that graphyne exhibits a high transversal contraction in both directions: Poisson’s ratio varies around 0.40. Because Poisson’s ratio becomes smaller for the zigzag direction of the graphyne sheet than for its armchair direction (ν
y
x
>
ν
x
y), it can be concluded that graphyne is also orthotropic regarding the transversal contraction. In this case, the percentage difference (about 7%) is smaller than that for Young’s moduli (about 10%). Therefore, graphyne is a slightly orthotropic material but, given the small (7–10%) difference, it can also be assumed as a “quasi-isotropic” material.Before going into detailed explanations of the shear modulus (G) and bulk modulus (K), it is useful to define the stress-strain relation for two-dimensional (2D) media, that is, without a 3rd dimension. Let us consider a 2D homogeneous isotropic sheet, for which the relation between the stress tensor and strain tensor reads(9)
σ
i
j
=
2
G
ε
i
j
+
K
-
G
ε
k
k
δ
i
j
i
,
j
,
k
=
1,2
,where index k satisfies the sum convention. This tensor relation depends on the 2D bulk modulus K and shear modulus G. Similarly, the strain-stress relation is given by(10)
ε
i
j
=
1
+
ν
E
σ
i
j
δ
i
j
-
ν
E
σ
k
k
δ
i
j
i
,
j
,
k
=
1,2
,where ν and E are 2D Poisson’s ratio and Young’s modulus, respectively. Clearly, there are only two independent moduli and the comparison between (11) and (12) yields the following formulas:(11)
G
=
E
2
1
+
ν
,
(12)
ν
=
K
-
G
K
+
G
.Replacing G (given by (13)) in (14) and solving this with respect to K, the following definition is obtained for the 2D bulk modulus:(13)
K
=
E
2
1
-
ν
.This definition, obtained from 2D elasticity, differs from the one derived from 3D elasticity and given by(14)
K
=
E
3
1
-
2
ν
.Using these 2D elasticity relations ((11) and (13)), we should bear in mind that we are not dealing with either plane-strain or plane-stress elasticity. Such specifications only arise if we deal with 3D elasticity. Eischen and Torquato [49] also derived relations between 2D and 3D moduli for a 3D isotropic material. They proved, among other results, that the 2D shear modulus G (either in plane-strain or plane-stress) is equal to the 3D shear modulus. However, the bulk moduli relations are not so simple. The 2D bulk modulus K is related to either the bulk modulus K
P
ε (for 3D plane-strain) or bulk modulus K
P
σ (for 3D plane-stress) by the relations(15)
K
P
ε
=
K
+
G
3
,
K
P
σ
=
9
K
G
3
K
+
4
G
.Concerning the shear modulus obtained for the RM, it can be compared with that obtained from the assumption of “quasi-isotropic” behavior of the graphyne sheet, in which G depends on Young’s modulus and Poisson’s ratio through (11). If this formula was applied to the armchair direction (E
x = 228.4 N/m, ν
y
x = 0.421), G = 80.4 N/m would be obtained. On the other hand, the application of (11) to the zigzag direction (E
y = 208.5 N/m, ν
y
x = 0.396) led to G = 74.7 N/m. The shear modulus obtained from RM (G = 70.5 N/m) is lower than both “quasi-isotropic” predictions (14% and 6%). These results also show the orthotropic nature of graphyne. However, and given the small (6–14%) difference in G, the shear behavior of graphyne is also “quasi-isotropic.”Likewise, the 2D bulk modulus obtained for the RM can also be compared with that calculated from the assumption of “quasi-isotropy,” in whichK depends on Young’s modulus and Poisson’s ratio through (13). If this formula is applied to the armchair direction (E
x = 228.4 N/m, ν
y
x = 0.421), it is obtained as K
=
197.2 N/m. The application of (13) to the zigzag direction (E
y = 208.5 N/m, ν
y
x = 0.396) leads to K = 172.6 N/m. The bulk modulus obtained from the RM (K = 164.5 N/m) is much lower than the former “quasi-isotropic” prediction (20%) and slightly lower than the latter “quasi-isotropic” prediction (4.9%). This fact also puts in evidence the orthotropic behavior of graphyne.In order to validate the proposed FE models, the results are compared with those obtained by other researchers. These results comprise MD results using AIREBO potentials (Cranford and Buehler [15], Zhang et al. [18], Yang and Xu [17], and Hou et al. [34]), or using generalized gradient approximation (Asadpour et al. [23]). The proposed FE model can be viewed as an approximation to these more complex models, which take into account many other parameters at atomic scale and MD simulations (e.g., displacement increment, time step, number of time steps, force field, equilibration phase, type of thermostat, MD ensemble, and cutoff function). Despite the difference of complexity between the FE models and MD models, it can be concluded from the observation of Table 3 that the elastic properties of graphyne derived from our FE model compare fairly well with those obtained by more rigorous approaches. Similarly to the results determined by Zhang et al. [18] and Yang and Xu [17], we also attained Young’s modulus for the zigzag direction slightly smaller than that for armchair direction. Yet, the results by Cranford and Buehler [15] are in disagreement with those by Zhang et al. [18] and Yang and Xu [17] and ours: Young’s modulus they obtained for the zigzag direction is higher than the value for the armchair direction. On the other hand, the absolute difference between both Young’s moduli (armchair and zigzag) determined by Cranford and Buehler [15] is higher than ours, but with opposite trend. The absolute difference between both Young’s moduli (armchair and zigzag directions) determined by Zhang et al. [18] and Yang and Xu [17] is lower than ours, but with a similar trend.Regarding Poisson’s ratio, the obtained values (0.40–0.42) match perfectly with that (0.41) found by Asadpour et al. [23], using DFT. Comparing the bulk modulus of graphyne with other available values, we can observe that although being somewhat different, our value is of the same order of magnitude of that presented by Asadpour et al. [23]. Finally, the value of the shear modulus obtained in our RM is remarkably close to the one obtained by Asadpour et al. [23]. Furthermore, Hou et al. [34] also analyzed the shear behavior of several graphyne sheet configurations and got a value of G around 75 N/m, which also agrees very well with the one obtained from our RM.


## 7. Conclusion

This paper presented, for the first time (to the authors’ knowledge), a finite element model to study the elastic properties of graphyne. After a brief introduction, the literature review on the latest developments of graphyne and the recent works dedicated to its mechanical characterization were presented, mostly focused on computational methods. Then, a FE model for graphyne sheets was presented in detail and the calculation of the graphyne’s elastic properties—Young’s moduli (E
x, E
x), Poisson’s ratio (ν
y
x, ν
x
y), shear modulus (G), and bulk modulus (K)—was conducted and further explained. Finally, the results of FE models were compared with those obtained from other methods (MD and DFT) and available in the literature. In brief, the main conclusions are the following:(i)
The FE models of graphyne proved to be simple but accurate. The linear elastic properties (Young’s modulus, Poisson’s ratio, bulk modulus, and shear modulus) obtained from the proposed FE models were in general agreement with those previously obtained by other authors using more complex computational models (MD and DFT).
(ii)
The influence of vdW interatomic forces on the linear elastic properties of planar graphyne is negligible and can be disregarded if small strain hypothesis is adopted. In case of either in-plane large deformations (non-linear elasticity) or out-of-plane bending of the graphyne sheet, the effect of vdW forces may be relevant.
(iii)
The FE models indicated that graphyne exhibits orthotropic behavior, a fact that agrees with the conclusions reported by other researchers. However, the differences between the orthotropic constants and the “quasi-isotropic” predictions were not meaningful, showing that graphyne may be modelled with “quasi-isotropic” behavior.



---

*Source: 7487049-2016-05-17.xml*