This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Exploration of mechanical, thermal conduction and electromechanical properties of graphene nanoribbon springs

Brahmanandam Javvajia Bohayra Mortazavi a Timon Rabczukb and Xiaoying Zhuanga,∗
aInstitute of Continuum Mechanics, Leibniz Universität Hannover, Appelstr. 11, 30167 Hannover, Germany
bInstitute of Structural Mechanics, Bauhaus University Weimar, Marienstrasse 15, 99423 Weimar, Germany.
Corresponding author E-mail: [email protected]

Abstract

Recent experimental advances [Liu et al., npj 2D Materials and Applications, 2019, 3, 23] propose the design of graphene nanoribbon spring (GNRS) to substantially enhance the stretchability of pristine graphene. GNRS is a periodic undulating graphene nanoribbon, where undulations are of sinus or half-circles or horseshoe shapes. Besides those, GNRS geometry depends on design parameters, like pitch’s length and amplitude, thickness and joining angle. Because of the fact that parametric influence on the resulting physical properties are expensive and complicated to be examined experimentally, we explore the mechanical, thermal and electromechanical properties of GNRS using molecular dynamics simulations. Our results demonstrate that horseshoe shape design of GNRS (GNRH) can distinctly outperform the graphene kirigami design concerning the stretchability. The thermal conductivity of GNRS were also examined by developing a multiscale modeling, which suggests that the thermal transport along these nanostructures can be effectively tuned. We found that however, the tensile stretching of GNRS and GNRH does not yield any piezoelectric polarization. The bending induced hybridization change results in a flexoelectric polarization, where the corresponding flexoelectric coefficient is 25%25\% higher than graphene. Our results provide a comprehensive vision to the critical physical properties of GNRS and may help to employ the outstanding physics of the graphene to design novel stretchable nanodevices.

1 Introduction

During the last decade, graphene [1, 2, 3] the two-dimensional (2D) form of sp2 carbon atoms has attracted astonishing attentions of scientific and industrial communities, stemmed from its extraordinarily high mechanical [4] and thermal [5] conduction properties along with exceptional electronic and optical characteristics. In particular, graphene can exhibit remarkably high Young’s modulus and tensile strength of 1000 GPa and 130 GPa [4] respectively, along with superior thermal conductivity of around 4000 W/mK [6] that outperforms diamond and other conventional materials. The exceptional physical and chemical properties of graphene, promoted the experimental and theoretical endeavours to fabricate novel graphene’s 2D counterparts, and subsequently explore their intrinsic properties and application prospects. As a results of scientific accomplishments, 2D materials family is commonly considered as the most vibrant class of materials, in which new members are continuously introduced, either theoretically predicted or experimentally fabricated. Worthy to remind that despite the outstanding properties of pristine graphene, it is not an ideal material and naturally shows few drawbacks, such as the lack of an electronic band-gap and brittle failure mechanism [4, 7]. In addition, the ultra-high thermal conductivity of graphene also prohibits its application for thermoelectric energy generation.
Nonetheless, it has been also practically proven that graphene can show largely/finely tunable electronic, mechanical, thermal, optical and electromechanical properties, with accurately controlled mechanical straining [8, 9, 10, 11, 12, 13], defect engineering [14, 15, 16, 17] or chemical doping [18, 19, 20, 21, 22]. We remind that for many centuries, springs have been playing pivotal roles in the design and fabrication of various kind of devices. The importance of springs originates from the fact that, while the mechanical properties of a material is considered as its inherent property and thus invariable, when it is shaped in the form of springs the subsequent structures can show superior stretchability and diverse mechanical responses. In particular, the design of spring like structures have been known as one of the most effective ways to design highly stretchable and flexible moving components.
For the employment of graphene in flexible nanoelectronics, its ductile and highly rigid mechanical properties are undesirable. Therefore, engineering of the graphene design in order to improve its stretchability is a critical issue [23, 24, 25, 26, 27]. To address this challenge, in a most recent exciting experimental advance by Liu et al. [28] the old idea of spring design has been applied for the case of graphene to enhance its stretchability and flexibility via a nanowire lithography technology. This experimental advance consequently raises questions concerning that how the design of graphene springs can be improved to reach higher degrees of stretchability. In addition, the electronic and heat transport properties of these novel nanostructures should be also examined, in order to provide comprehensive visions for the design of nanodevices. As a common challenge in the electronic apparatus, the thermal conductivity of employed components should be high to effectively dissipate the excessive heats. On the opposite side, low thermal conductivity is a key requirement for the enhancement in the efficiency of thermoelectric energy conversion. As a common barrier, it is well known that the evaluation of the mechanical and transport properties of graphene springs by the experimental tests are not only complicated but also time consuming and expensive as well. This study therefore aims to investigate the mechanical response and heat conduction properties of graphene springs via conducting extensive classical molecular dynamics simulations. Since the graphene is the frontier and symbolic member of 2D materials, commonly the experimental and theoretical methodologies that are applied for the graphene can be extended for the other members of 2D materials family. We are thus hopeful that the obtained results by this study may sever as valuable guides for the future theoretical and experimental studies on the design of 2D materials spring like structures.

2 Methods

We conducted classical molecular dynamics (MD) simulations to evaluate the mechanical properties and thermal conductivity of graphene springs, using the large scale atomic/molecular massively parallel simulatior package [29]. To this aim, we used the optimized Tersoff potential parameter set proposed by Lindsay and Broido [30] for introducing the atomic interaction between carbon. This version of Tersoff potential is not only highly computationally efficient, but also can yield accurate results for the mechanical and thermal properties of graphene. We analysed the mechanical response by conducting the uniaxial tensile simulations at room temperatures. For the evaluation of mechanical properties, we modified the cutoff of Tersoff potential from 0.18 nm to 0.20 nm to remove an unphysical stress pattern and moreover accurately reproduce the experimental results for the tensile strength of pristine graphene [31]. In this case, the time increment of MD simulations was set at 0.25 fs. Before applying the loading conditions, all structures were equilibrated using Nosé-Hoover thermostat method. For the loading condition, a constant engineering strain rate of 1×\times108 s-1 was applied, by increasing the periodic size of the simulation box along the loading direction in every time step. Virial stresses at every step were recorded and averaged over every 20 ps intervals to report the stress-strain relations. To evaluate the thermal conductivity, we used equilibrium molecular dynamics (EMD) method. The time increment of EMD simulations was set at 0.25 fs. In the EMD method, the heat flux vector was calculated via:

𝐉(t)=i(ei𝐯i+12i<j(𝐟ij(𝐯i+𝐯j))𝐫ij)\mathbf{J}(t)=\sum_{i}\left(e_{i}\mathbf{v}_{i}+\frac{1}{2}\sum_{i<j}\left(\mathbf{f}_{ij}\cdot\left(\mathbf{v}_{i}+\mathbf{v}_{j}\right)\right)\mathbf{r}_{ij}\right) (1)

where eie_{i} and 𝐯i\mathbf{v}_{i} are respectively the total energy and velocities of atom ii, 𝐟ij\mathbf{f}_{ij} and rijr_{ij} are the interatomic force and position vector between atoms ii and jj, respectively. In the approach, first the structures were equilibrated at a constant volume and room temperature using the Berendsen thermostat method. Before the evaluation of thermal conductivity, in order to remove the effects of previously applied thermostat, we used constant energy (NVE) simulations. For the evaluation of effective thermal conductivity, individual simulations were conducted for 1 ns under the NVE ensemble. The EMD method relies on relating the ensemble average of the heat current auto-correlation function (HCACF) to the thermal conductivity kk, via the Green-Kubo expression:

kαβ=1VkBT20Jα(0)Jβ(t)𝑑tk_{\alpha\beta}=\frac{1}{Vk_{B}T^{2}}\int_{0}^{\infty}\left\langle J_{\alpha}(0)J_{\beta}(t)\right\rangle dt (2)

where kBk_{B} is Boltzmann’s constant, TT is the simulation temperature, and VV is the total volume of the graphene spring. In order to reach a converged thermal conductivity, several independent simulations with uncorrelated initial velocities were carried out and the acquired HFACFs were averaged.

3 Results

Refer to caption
Figure 1: Unit cell representation and definition of various structural parameters for (a) sinus shape and (b) horseshoe shape graphene nanosprings.

Fig. 1(a) and (b) shows the atomic unit-cell representation for graphene nanoribbon (GNR) in sinus shape and horseshoe shape, respectively. Sinus shape of GNR (GNRS) obtained from cutting the infinite graphene sheet using two sine curves that are parallel to each other. The parameters defining the sine curve are pitch length (sp)(s_{p}) and amplitude (sa)(s_{a}). The locus of each point normals with a constant length (st)(s_{t}) creates a parallel sine curve. The variable sts_{t} defines the thickness of GNRS. The choice of (sp)(s_{p}) and (sa)(s_{a}) is arbitrary. However, (st)(s_{t}) should be less than the radius of the curvature of the sine curve, which is to avoid producing bigger arcs that cover the crests and troughs (cusps) of the sine curve. Using these variables for GNRS, we define a sinus shaped region on a pristine graphene sheet and removed atoms other than this region, which creates an initial atomic configuration for GNRS. The volume of GNRS is defined as the area under the parallel sine curve times the thickness of pristine graphene sheet 3.3Å3.3~{}\text{\AA} [32].
The horseshoe shape design of GNR (GNRH) composes by connecting two circular arcs of the inner radius (hr)(h_{r}) with the intersecting angle (hθ)(h_{\theta}). When hθ=0h_{\theta}=0^{\circ}, GNRH looks like a series of connected semi-circles and for hθ>0h_{\theta}>0^{\circ} and hθ45h_{\theta}\leq 45^{\circ}, the horseshoe design is obtained. For hθh_{\theta} more than 4545^{\circ}, the semi-circles merge each other, which is not desirable. These parameters define the shape of a single horseshoe curve. Another curve with radius hr+hth_{r}+h_{t} creates a parallel curve, where hth_{t} defines the thickness of GNRH. We construct the horseshoe-shaped region on a pristine graphene sheet and removed atoms other than this region, which creates an initial atomic configuration for the GNRH system. After the initial preparation of spring structures, we removed the carbon atoms bonded with a single carbon atom along the lateral edges. Since these atoms can lead to instability in simulations. It is worth noting that the carbon atoms belong to the curved edges of the GNRS and GNRH system are not terminated with hydrogen atoms.
With the defined geometrical parameters and initial atomic configurations for GNRS and GNRH, we consider the following parameter sets for estimating the mechanical and thermal properties. Parameter set 𝐬𝐩𝐬𝐭𝐬𝐚\mathbf{s_{p}-s_{t}-s_{a}}: The starting value for sps_{p} is 9 nm (where the size effects are absent in pristine graphene [13]), and sas_{a} is 2.5 nm, where at least 10 carbon rings accompanying in the lateral direction. The minimum possible value of sts_{t} is 1.6 nm for this combination of sps_{p} and sas_{a} for having a reasonable thickness for these spring systems. Further, we increase the values of spsp, stst and sasa by integral multiples from 1 to 5. Parameter 𝐡θ\mathbf{h}_{\theta}: We consider hθh_{\theta} as 0,15,300^{\circ},15^{\circ},30^{\circ} and 4545^{\circ} by keeping hrh_{r} at 2.5 nm and hth_{t} as 1.6 nm. Using this choice of parameters, we explore the effect of hθh_{\theta} on mechanical and thermal properties. Parameter 𝐡𝐫𝐡𝐭\mathbf{h_{r}-h_{t}}: In this set, hrh_{r} and hth_{t} parameters are scaled by integral multiples from 1 to 5 starting from 2.5 and 1.6 nm at fixed hθh_{\theta}. We fix the value of hθh_{\theta} from previous parameter set, which has shown exceptional mechanical properties. The spring structures from the different parameter sets are made by cutting from graphene sheets oriented along the zigzag direction. Our comparative results discussed in the supplementary information†document however confirm that the orientation of original graphene sheet show negligible effects in the estimated properties.

3.1 Mechanical properties

Refer to caption
Figure 2: Stress-strain curves for (a) parameter set spstsas_{p}-s_{t}-s_{a} in GNRS, (b) hθh_{\theta} and hrhth_{r}-h_{t} for GNRH. (c) Standard deviation of the zz-coordinates for selected GNRS and GNRH systems. (d) Variation of rupture strain and tensile strength (TS) for GNRS and GNRH systems with arc length. Solid lines indicate RS and dashed lines correspond to TS.
Refer to caption
Figure 3: Tensile deformed atomic configurations for (a-e) 183.2518-3.2-5 GNRS and (f-j) 53.2455-3.2-45^{\circ} GNRH. The strain levels for the atomic snapshots are as follows: (a and f) 0.0, (b) 0.18, (c) 0.30, (d) 0.42, (e) 0.54, (g) 1.42, (h) 2.0, (i) 2.17 and (j) 2.23. Color bar indicate the stress values.

Fig. 2(a) shows the stress-strain curves for the parameter set spstsas_{p}-s_{t}-s_{a}. Under stretching of GNRS system, the stress values remain low up to a strain of 0.3. This behavior is in contrast to that of pristine graphene. The given deformation stretches the bond lengths in pristine graphene and increases the stress levels. In GNRS systems, the given deformation deflects the atomic system instead of stretching the bonds, which is similar to the earlier observations for graphene kirigami [33]. Further, the recent study on nonlinear vibrations for helical graphene nanoribbon shows a transformation of softening type to hardening type response between the amplitude to frequency variation during the increase in mechanical loading[34]. In GNRS, the initial plateau in stress-strain curve up to a strain of 0.30.3 corresponds to the soft nonlinearity. Because of the transformation from softening to hardening type natural vibrations, GNRS system start stretching due to the increased loading. Fig. 3(a) visualizes the initial (strain is 0) atomic configuration for 183.2518-3.2-5 GNRS. Visual molecular dynamics package [35] has been used to generate the atomic snapshots. At a strain of 0.18, because of the deflections, a higher number of crests and troughs are observed in Fig. 3(b). Further straining to 0.30, deflects the GNRS without increasing the peaks (Fig. 3(c)). For strain less than 0.30, there is no significant increase in the stress distributions. When strain is more than 0.30, there is an increase in stress due to the bond stretching. Fig. 3(d) shows the locations of stress concentrations (near peaks) in the atomic configurations. Atomic configuration in Fig. 3(e) shows the multiple bond failures at strain 0.54. The failure process for other members in this parameter set is similar to the configuration 183.2518-3.2-5 GNRS. However, there are differences in the tensile strength (TS) and onsite of failure/rupture strain (RS) values due to the changes in the geometrical parameters. We use arc length as the variable to discuss the parametric dependence of RS and TS on sps_{p}, sts_{t} and sas_{a}. Fig. 2(d) shows that RS is nearly converged to 0.45 for arc lengths larger than 170 nm. Whereas, TS has a very large variation with arc length, from 33 to 7 GPa, which is due to the increase of parameter sts_{t} in GNRS design. As sts_{t} increases, the interaction between stress centers near crest and troughs of GNRS decreases, which decrease the total system stress. Further increase in the size parameters will converge to a constant value. Finally, for GNRS systems, the stress concentrates near the peak portions lead to global failure by breaking the bonds in carbon rings. The RS for the experimentally manufactured GNRS is about 0.350.35 for a system with thickness 5050 nm [28]. This observation is in close agreement with the converged RS value of 0.430.43 for spring design, which shows that our simulation predictions are accurate enough to explore the theoretical understanding for these novel design.
Fig. 2(b) shows the stress-strain response for parameter set hθh_{\theta}, where hrh_{r} and hth_{t} values are at 2.5 and 1.6 nm, respectively. When hθ=0h_{\theta}=0^{\circ}, there is no stress rise up to a strain of 0.3. For strain range 0.3 to 0.6, the deformation in atomic configuration rises the system stress followed by a failure. The strain range for the non-zero portion of the stress-strain curve shifts between 0.5 to 1.0 when hθh_{\theta} is 1515^{\circ}. For hθh_{\theta} equal to 3030^{\circ}, the non-zero portion of the stress-strain curve span the strain range of 0.94 to 1.5. For 45 GNRH, this strain range increased to 2.4. Figs. 3(f)-(j) shows the atomic configurations for a GNRH with a 45 connecting angle. The closeness between semi-circular segments in GNRH develops strong repulsion interactions compared to GNRS. Such repulsion largely deflects the atomic system. Further, mechanical stretching reduces atomic deflections by maintaining the stress levels via transforming the smooth circular GNRH segments to sharp peaks. Fig. 3(g) shows the atomic configuration with several peaks in GNRH. For strain greater than 1.5, Fig. 2(b) shows a linear stress-strain response due to the bond stretching. At strain 2, GNRH system looks like a combination of thread and knots, as seen in Fig. 2(h), where stress concentrates near the thread portions. Further increase in strain in GNRH, leads to bond failure. Figs. 2(i) and (j) at strain levels of 2.17 and 2.23 shows the complete failure of GNRH.
Interestingly, GNRH with various hθh_{\theta} maintained the stress levels when increasing the strain range. We consider varying the width parameter hth_{t} in GNRH by keeping the other parameters constant to check its influence on the mechanical properties. The range of hth_{t} is limited by the choice of other two parameters. For example, consider hrh_{r} equal to 55 nm and hθh_{\theta} is 4545^{\circ}. The maximum available value of hth_{t} is 44 nm. When hth_{t} is greater than 44 nm, the two circular cross-sections of GNRH unitcell overlap with each other, which is not desirable. We consider hth_{t} values as 1.61.6, 2.42.4 and 3.23.2 and 44 nm and the corresponding RS values are 2.642.64, 2.442.44, 2.172.17 and 1.91.9, respectively. TS values are noted as 53.8353.83, 38.0638.06, 24.8124.81 and 17.4917.49 GPa. (see supplementary information†) When varying hθh_{\theta}, the RS increased with very low effect on TS (see circle markers with dotted and solid line in Fig. 2(d)). Whereas, variation of hth_{t} effecting both RS and TS in GNRH. The increase in thickness, increases the separation between stress centers and decreases the TS, which lead to early failure and decrease in RS. This finding implies that changing of hth_{t} strongly influence the both RS and TS of GNRHs. The observations concerning the width effect on the mechanical response are in close agreement with the eralier report based on MD simulations[36].
When compared to GNRS, GNRH shows higher stress levels (dotted lines in Fig. 2(d)). From the structural point of view, GNRH differs from GNRS in two factors, one is the smoothness of undulations and the second is the closeness between the undulations. The smoothness of circular arcs in GNRH makes the stress to distribute across all the boundary atoms. The increased number of atoms with higher per atom stress values, increase the total stress in the atomic system. In the case of GNRS, the lower number of atoms with higher per atom stress near the peaks of sine curve makes the total stress lower. From Fig. 2(d), TS values for 2.51.602.5-1.6-0^{\circ} GNRH and 91.62.59-1.6-2.5 GNRS are 39.72 and 31.93 GPa, respectively, which represents that the smoothness factor accommodates more number of atoms with high stress levels in horseshoe shape design. The closeness between the undulations increase the deflections in the atomic system, which help to avoid the bond stretching and stress rise. These deflections in GNRS and GNRH systems are measured using the standard deviation of the zz-coordinates (Z¯)(\bar{Z}) [37], which is defined as Z¯=i=1N(ziz0)2/N\bar{Z}=\sqrt{\sum_{i=1}^{N}\left(z_{i}-z_{0}\right)^{2}/N}, where ziz_{i} is the zz-coordinate of the ithi^{\text{th}} atom and z0z_{0} is the averaged zz-coordinate over NN atoms. Fig. 2(c) plots the computed Z¯\bar{Z} with respect to strain for the selected GNRS (91.62.5)(9-1.6-2.5) and GNRH (parameter set hθh_{\theta} ) systems. Z¯\bar{Z} initially increases with strain, which represents that the energy of given mechanical straining used to increase the deflections in both GNRH and GNRS systems. After reaching a maximum deflection, given tensile loading starts stretching the atomic system and decreasing the deflections which decrease Z¯\bar{Z}. However, the magnitude of Z¯\bar{Z} for 00^{\circ} GNRH configuration is high compared to 91.62.59-1.6-2.5 GNRS, which supports that the repulsion interactions in GNRH are heavier compared to GNRS. Z¯\bar{Z} increasing with an increase in parameter hθh_{\theta}. Z¯\bar{Z} is highest for 4545^{\circ} GNRH.
The very strong repulsions exist between the semi-circular rings due to the minimum spacing. The very high deflections and smoothness of circular cross-sections helps to avoid the stress concentrations in GNRH, which helps to enhance the mechanical properties. As a total, very high value of RS is noted for GNRH. With increasing the system size (arc length of GNRH), RS tending to converge to a value of 22, which is about 17%17\% higher than the graphene kirigami design [38, 39], keeping the stress-levels identical.

3.2 Thermal properties

Refer to caption
Figure 4: Room temperature thermal conductivity for (a) parameter set spstsas_{p}-s_{t}-s_{a} in GNRS, (b) hθh_{\theta} and (c) hrhth_{r}-h_{t} for GNRH as a function of correlation time. (d) Thermal conductivity for graphene nanosprings at the room temperature as a function of arc lengths. The insets in (d) illustrate the finite element modeling results for temperature distribution of GNRS and GNRH in diffusive regime.

We next study the thermal transport through the GNRS and GNRH systems. Fig. 4 shows the EMD results of effective thermal conductivity of GNRS for few studied parameter sets, as a function of correlation time. In this case, we normalized the effective conductivity with respect to that of the pristine graphene to better illustrate the suppression rate. As it is clear, for the samples with lower thermal conductivity the convergence occurs at lower correlation times. For 91.62.59-1.6-2.5 GNRS, κ\kappa is about 0.01750.0175 times the thermal conductivity of pristine graphene κ0\kappa_{0} (from Fig. 4(a)). We employed a square sheet of pristine graphene with 1616 nm side length to estimate κ0\kappa_{0}. The estimated value of κ0\kappa_{0} is 1000±1001000\pm 100 W/(mK), which is an average over 88 independent simulations. This value is in close agreement with earlier reports using EMD method [40]. Note that due to the implementation of inaccurate heat-flux formula for many-body interactions in the LAMMPS tool, EMD method significantly underestimates the the thermal conductivity of graphene [41]. The suppressed thermal transport along the GNRS is due to the phonon-boundary scattering in these systems, which is in accordance with the previous reports concerning the graphene kirigami [33]. Also, there exist a very strong conversion for the phonon mode polarization from out-of-plane to in-plane and opposite phase for the left and right parts of the unitcell [42, 43]. Such conversion of phonon modes loose the transporting energy through the scattering with localized phonon modes near the edges. As a result the transport of heat flux is low and reduce the thermal conductivity. When sts_{t} increases, the ratio of atoms on the boundaries to the total number of atoms decreases, which results in lower phonon-boundary scattering rate and subsequently facilitate the heat transport.
As proposed in the previous work [33], we use a microscale continuum model of graphene spring to evaluate the effective thermal conductivity. This evaluation carried within the diffusive range, in which the phonon-boundary scattering vanishes. To this aim, a system is modeled within the finite element (FE) approach to establish connections between the effective thermal conductivity and nanoribbon’s arc length. We apply inward and outward heat-fluxes on the two opposite sides of GNRS as the boundary conditions. Using the measured temperature gradient along the heat transfer direction, the effective thermal conductivity was computed from the Fourier’s law. We then used a first order rational curve fitting to extrapolate the atomistic results (circular markers in Fig. 4(d) that correspond to the averaged κ/κ0\kappa/\kappa_{0} over several samples of spstsas_{p}-s_{t}-s_{a} and the standard deviation among them) dominated by the phonon-boundary scattering to the diffusive transport by the FE simulations. As shown in Fig. 4(d), this approach could provide a very accurate estimation of thermal transport at different arc lengths, and reveals that the phonon-boundary scattering starts to vanish at large arc lengths.
For GNRH systems with varying hθh_{\theta} is shown in Fig. 4(b). κ/κ0\kappa/\kappa_{0} for 2.51.602.5-1.6-0^{\circ} GNRH is about 0.01690.0169, which is nearly same as 91.62.59-1.6-2.5 GNRS. The constant thickness and similar scattering effects in these two spring systems produce the thermal conductivity nearly equal. Keeping thickness constant and increasing the joining angle hθh_{\theta} to 1515^{\circ}, κ/κ0\kappa/\kappa_{0} decreased from 0.01690.0169 to 0.01110.0111. As hθh_{\theta} increases, there is an increase in the radius of curvature of the junction region that connects the two semi-circular segments. The phonon transport through this increased curvature experiences significant scattering, which reduces the heat transfer and κ\kappa. The thermal conductivity for hθh_{\theta} 3030^{\circ} and 4545^{\circ} is nearly the same.
We examine thermal conductivity for the GNRH samples used in estimating the effect of width on mechanical properties in Section 3.1. The effective thermal conductivity for 2.51.6452.5-1.6-45^{\circ} and 51.6455-1.6-45^{\circ} are 0.00840.0084 and 0.00550.0055. This represent that increase in hrh_{r}, increases the radius of curvature and produce more edge localized phonon modes. These modes do not contribute for thermal transport, as a result the thermal conductivity decreases for 51.6455-1.6-45^{\circ} GNRH system (see supplementary information†) However, the increase in hth_{t} increase the number of phonon modes in GNRH keeping the density of edge localized modes same. This reduces the edge scattering and increase the phonon transport, thus increase in κ\kappa [44]. Fig. 4(c) plots κ\kappa for increasing values of hrh_{r} and hth_{t} keeping joining angle hθh_{\theta} as 4545^{\circ}. As the GNRS radius and thickness increasing, the available region for heat transport increases, which helps to lower the scattering and rise κ/κ0\kappa/\kappa_{0} from 0.00840.0084 to 0.02160.0216. This rise of κ\kappa is small compared to the GNRS systems due to the large curvature induced scattering. We repeat the FE modeling for GNRH with hθh_{\theta} as 4545^{\circ} similar to GNRS. The fitting between atomistic results and FE modeling is very good. However, GNRH fitting is converged at significantly larger cut lengths compared to GNRS. This proves that, curvature induced scattering reduces κ\kappa in GNRH.

3.3 Electromechanical properties

The nanoscale electromechanical properties (piezoelectricity and flexoelectricity) have been gaining intense attention due to their ability to sustain large deformations. This feature adds many different applications in the energy conversion process. These properties are limited in pristine graphene due to the crystallographic symmetry. The structural and chemical modifications break the symmetry and induce polarization under mechanical deformations [45, 46, 47, 48, 49, 17, 50]. The bending deformation of pristine graphene induces a polarization due to the change of hybridization state of the carbon atom, known as pyramidalization [51, 48, 52]. In GNRS and GNRH systems, cut patterns cancel the symmetry and promises for electromechanical coupling. In this section, we subject the GNRS and GNRH system to both tensile and bending deformations to obtain the respective polarization variations.
To calculate the polarization in the atomic system, we utilize the charge-dipole (CD) model along with the short-range bonded interactions (Tersoff potential). According to CD model, each atom ii is associated with charge qiq_{i} and dipole moment 𝐩i\mathbf{p}_{i} [53, 54]. The mathematical CD formulation involves the various contributions from charge-charge, charge-dipole and dipole-dipole interactions to the total system short-range interaction energy. The minimization of energy function gives the numerical values of qiq_{i} and 𝐩i\mathbf{p}_{i}. The complete details about the CD model and estimation of charges and dipole moments can be found in [17, 52] and references therein.
For applying deformation, we have added left and right rectangular regions to the GNRS and GNRH systems by discarding the periodic boundary condition used in Section 2. These regions have equal sts_{t} or hth_{t} to the spring systems with 11 nm length along the spring longitudinal direction. The left and right regions helps to hold the given displacement, particularly during bending test, and relax the remaining system. We define a load cycle by prescribing the displacement of atoms to left and right regions for 11 ps time period, followed by a relaxation of 22 ps time period. Because of the non periodic boundaries, we perform simulations at different repetitions of sinus and horseshoe shapes in spring systems. These simulations help us to study the size effect on electromechanical properties. For every load cycle, we note the evolution of atomic configuration, corresponding charge and dipole moments. The total polarization of the atomic system is the sum of all atomic dipole moments divided by the volume of atomic system.

Refer to caption
Figure 5: (a) Variation of polarization PxP_{x} with tensile strain. (b) Bending induced polarization PzP_{z} response with respect to strain gradient. (c) Dependence of flexoelectric coefficient μzxzx\mu_{zxzx} with arc length of GNRS and GNRH systems. Legends indicate the respective atomic configurations used.

For tensile deformation, we apply a displacement in longitudinal direction ux=ϵ˙tloadl0u_{x}=\dot{\epsilon}t_{\text{load}}l_{0} to the atomic system, where ϵ˙\dot{\epsilon} is strain rate equal to 1×\times108 s-1 as used in Section 2, tloadt_{\text{load}} is the loading time (1 ps) and l0l_{0} is the initial length of the atomic system in longitudinal direction. The load cycles continued to reach the strain ϵxx\epsilon_{xx} of 0.40.4 for 91.62.59-1.6-2.5 and 2.51.602.5-1.6-0^{\circ} systems. The strain limit 0.40.4 corresponds to the linear rise in stress-strain response for these systems as shown in Fig. 2(a). At each load cycle, the polarization is measured and the variation with strain is plotted in Fig. 5(a). The variation in polarization with strain is nearly negligible for both GNRS and GNRH systems. The coefficient of variation for the polarization response is nearly equal to 1 for both GNRS and GNRH, similar to non-piezoelectric pristine graphene [17]. The cancellation of polarization at the sinus and horseshoe cut patterns make these systems as non-piezoelectric materials.
For bending deformation, we supply the following out-of-plane displacement field to the atomic system

uz=Kx22,u_{z}=K\frac{x^{2}}{2}, (3)

where xx represent the atom coordinate in the longitudinal direction, KK represent the inverse of curvature (strain gradient) of the bending plane. After prescribing the bending deformation, the atoms belongs to the left and right termination regions are held fixed. Whereas, the interior atoms are allowed to relax to energy minimizing positions using the conjugate-gradient algorithm. For the energy minimized configuration, we note the charges and dipole moments for each atom. From these, the relationship between polarization to strain and strain gradient as

Pz=dzxzϵxz+μzxzxϵxzx,P_{z}=d_{zxz}\epsilon_{xz}+\mu_{zxzx}\frac{\partial\epsilon_{xz}}{\partial x}, (4)

where PzP_{z} is the out-of-plane polarization, μzxzx\mu_{zxzx} is the out-of-plane bending flexoelectric coefficient, dzxzd_{zxz} is the piezoelectric coefficient and ϵxz\epsilon_{xz} is the strain component. Here the piezoelectric contribution (dzxzϵxz)(d_{zxz}\epsilon_{xz}) is removed from the total polarization, because the given bending deformation leads to a linear variation in ϵxz\epsilon_{xz} along the xx direction. The linear variation demonstrates that the induced deformation is symmetric and the resulting polarization due to strain is canceled out [52]. From Eq. 3, the strain gradient is given as

ϵxzx=12x(uzx+uxz)=122uzx2=12K=Keff,\frac{\partial\epsilon_{xz}}{\partial x}=\frac{1}{2}\frac{\partial}{\partial x}\left(\frac{\partial u_{z}}{\partial x}+\frac{\partial u_{x}}{\partial z}\right)=\frac{1}{2}\frac{\partial^{2}u_{z}}{\partial x^{2}}=\frac{1}{2}K=K_{\text{eff}}, (5)

where uxu_{x} is zero because of fixing the atom positions belongs to the left and right boundary and KeffK_{\text{eff}} is the effective strain gradient, which is equal to half of the given value of KK. Substituting Eq. 5 along with zero piezo contribution in Eq. 4 leads to

Pz=μzxzxKeff.P_{z}=\mu_{zxzx}K_{\text{eff}}. (6)

The polarization PzP_{z} at various strain gradients KeffK_{\text{eff}} is plotted in Fig. 5(b). The dipole moment 𝐩i\mathbf{p}_{i} on atom rises due to the pyramidalization. The bending deformation transforms the atomic hybridization of carbon atom from sp2sp^{2} to sp3sp^{3}. The bending of the bond between carbon atoms forces the valence electron to develop a bonding interaction with neighboring atoms. This interaction allows a mixing between valence (π)(\pi) and bonding (σ)(\sigma) electrons which lead to πσ\pi-\sigma interactions. The πσ\pi-\sigma interaction modifies the charge state of carbon atoms and the local electric fields, which leads to the dipole moment of that atom. The increased deformation increases the πσ\pi-\sigma interaction, which rise the dipole moment. Fig. 6(a) and (b) show the distribution of dipole moment for the unitcells next to the center of the atomic systems. Both sides of the lateral edges have opposite dipole moments, and their magnitude is a function of bending displacement (from Eq. 3). The observed linear variation between polarization and strain gradient for spring systems is similar to the pristine graphene variation. The slope of this variation gives the flexoelectric coefficient μzxzx\mu_{zxzx}. Numerical value of μzxzx\mu_{zxzx} for 91.62.59-1.6-2.5 and 2.51.602.5-1.6-0^{\circ} systems are 0.00340.0034 and 0.00350.0035 nC/m, which is 25%25\% higher than pristine graphene [52]. The hybridization angle θσπ\theta_{\sigma\pi} is the angle between a fixed out-of-plane point to one of the bond between carbon atoms. For the initial or flat system, this angle is exactly 9090^{\circ}. The deformation of bond changes this angle. The change in angle for pristine and 2.51.602.5-1.6-0^{\circ} system are 4.624.62 and 6.386.38^{\circ}, respectively. The increment in θσπ\theta_{\sigma\pi}, increases the dipole moment of the system [51], which increase the flexoelectric coefficient. Note that we calculated the hybridization angle at the same location in both the systems.

Refer to caption
Figure 6: Atomic configurations colored with dipole moment pzp_{z} for a strain gradient of 0.0060.006 nm-1. (a) 91.62.59-1.6-2.5 GNRS design. GNRH design with inner radius 2.52.5 nm, thickness 1.61.6 nm and connecting angle hθh_{\theta} as (b) 00^{\circ}, (c) 1515^{\circ}, (d) 3030^{\circ} and (e) 4545^{\circ}. Dashed line indicate the middle portion of the atomic system.

From Fig. 5(b), the slope of the polarization to strain gradient curve is decreasing with the increase of hθh_{\theta} in GNRH systems. The ratio of change in pyramidalization angle in GNRH to pristine grpahene decreases as 1.271.27 , 1.141.14 and 1.051.05 and 0.850.85 with hθh_{\theta} as 00^{\circ}, 1515^{\circ}, 3030^{\circ} and 4545^{\circ}, respectively. The decrease in θσπ\theta_{\sigma\pi} decreases, the dipole moment distribution, as seen from Figs. 6(c) to (d), which decreases the flexoelectric coefficient. For hθ=45h_{\theta}=45^{\circ}, atoms on the lateral boundaries near the central line does not have dipole moments. The strong repulsions between these edges cancels the effect of pyramidalization, which decreases the flexoelectric coefficient.
In order to check the dependence of hth_{t} on the polarization, we consider 2.52.402.5-2.4-0^{\circ} and 2.53.202.5-3.2-0^{\circ} GNRH configurations. For these configurations, the polarization variation and flexoelectirc coefficient are nearly equal to the 2.51.602.5-1.6-0^{\circ} GNRH. The increase in thickness does not effecting the induced polarization and flexoelectiric coefficient (see supplementary information†). Since the increased thickness unable to change the pyramidalization angle further, makes the polarization comparable with smaller thickness system.
Further, the result of flexoelectric coefficient with length variation is given in Fig. 5(c) for GNRS and GNRH (hθ=0)(h_{\theta}=0^{\circ}). This result represent that there is a boundary effect when arc length is less than 100100 nm. For systems with arc length about 3030 nm, the local electric fields is strongly effected by the left and right region of atoms. Where the imposed boundary condition constrain the natural motion of the interior atoms, which restricts the process of pyramidalization and controls the dipole moment. When increasing the arc length this effect is slowly nullifying and the atomic configuration deforms to generate the dipole moments. For systems with lengths higher than 100100 nm, the boundary effect if completely negligible and the flexoelectric coefficients turn into a stable value. Finally, the flexoelectric coefficient of GNRS and GNRH-00^{\circ} system is 0.250.25 times higher than the pristine graphene.

4 Conclusion

Motivated by a latest experimental advance, we performed extensive classical molecular dynamics simulations to explore the mechanical, thermal conductivity and electromechanical properties of graphene nanoribbon springs (GNRS). In particular, we examined the effects of different GNRS design parameters on their physical properties. We found that by optimal design of GNRS systems, they can yield higher stretchability in comparison with kirigami counterparts. Horseshoe shape design of GNRS were found to show the better stretchability in comparison with other design strategies. In the aforementioned case, large deflections due to the strong repulsions between semi-circular rings help to keep the load bearing at larger strain levels. Our analysis of deformation process reveals that the stress concentrations occurring near the peak portions of GNRS induce local failure of carbon bonds and lead to final failure of structure. The thermal conductivity of GNRS were found to be substantially dependent on the width of nanoribbon’s width, due to the fact that phonon boundary scattering dominates the thermal transport. On this basis, we could establish the connection between the effective thermal conductivity of GNRS as a function of nanoribbon’s width size, by extrapolating the molecular dynamics results to the diffusive heat transfer model by the finite element. This approach can be used to effectively estimate the thermal conductivity, but also suggest that the thermal transport of GNRS can be effectively tuned by changing the design parameters. The negligible variation between polarization and strain proves that GNRH and GNRS systems are non-piezoelectric similar to graphene. A linear variation of polarization with strain gradient is observed in bending deformation test. The flexoelectric coefficient for GNRS and GNRH-00^{\circ} is 25%25\% higher than graphene. The decrease in μzxzx\mu_{zxzx} with increasing hθh_{\theta} is due to the decrease of pyramidalization angle. Our extensive theoretical results highlight the superior stretchability, finely tunable thermal transport and improved flexoelectric coefficient of GNRS, and suggest them as highly attractive components to design flexible nanodevices. The obtained results will hopefully guide future theoretical and experimental studies, to extend the idea of nanoribbon springs for graphene and other 2D materials.

Conflicts of interest

There are no conflicts to declare.

Acknowledgments

The authors gratefully acknowledge the sponsorship from the ERC Starting Grant COTOFLEXI (No. 802205). We acknowledge the support of the cluster system team at the Leibniz Universität of Hannover, Germany in the production of this work.

References

  • [1] K. S. Novoselov, “Electric Field Effect in Atomically Thin Carbon Films,” Science, vol. 306, pp. 666–669, oct 2004.
  • [2] A. K. Geim and K. S. Novoselov, “The rise of graphene,” Nature Materials, vol. 6, pp. 183–191, mar 2007.
  • [3] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Reviews of Modern Physics, vol. 81, pp. 109–162, jan 2009.
  • [4] C. Lee, X. Wei, J. W. Kysar, and J. Hone, “Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene,” Science, vol. 321, pp. 385–388, jul 2008.
  • [5] A. A. Balandin, “Thermal properties of graphene and nanostructured carbon materials,” Nature Materials, vol. 10, pp. 569–581, aug 2011.
  • [6] S. Ghosh, W. Bao, D. L. Nika, S. Subrina, E. P. Pokatilov, C. N. Lau, et al., “Dimensional crossover of thermal transport in few-layer graphene,” Nature Materials, vol. 9, pp. 555–558, jul 2010.
  • [7] P. Zhang, L. Ma, F. Fan, Z. Zeng, C. Peng, P. E. Loya, et al., “Fracture toughness of graphene,” Nature Communications, vol. 5, p. 3782, sep 2014.
  • [8] F. Guinea, “Strain engineering in graphene,” Solid State Communications, vol. 152, pp. 1437–1441, aug 2012.
  • [9] C. Metzger, S. Rémi, M. Liu, S. V. Kusminskiy, A. H. Castro Neto, et al., “Biaxial Strain in Graphene Adhered to Shallow Depressions,” Nano Letters, vol. 10, pp. 6–10, jan 2010.
  • [10] V. M. Pereira and A. H. Castro Neto, “Strain Engineering of Graphene’s Electronic Structure,” Physical Review Letters, vol. 103, p. 046801, jul 2009.
  • [11] S. Barraza-Lopez, A. A. Pacheco Sanjuan, Z. Wang, and M. Vanević, “Strain-engineering of graphene’s electronic structure beyond continuum elasticity,” Solid State Communications, vol. 166, pp. 70–75, jul 2013.
  • [12] F. Guinea, M. I. Katsnelson, and A. K. Geim, “Energy gaps and a zero-field quantum Hall effect in graphene by strain engineering,” Nature Physics, vol. 6, pp. 30–33, jan 2010.
  • [13] B. Javvaji, P. Budarapu, V. Sutrakar, D. R. Mahapatra, M. Paggi, G. Zi, et al., “Mechanical properties of Graphene: Molecular dynamics simulations correlated to continuum based scaling laws,” Computational Materials Science, vol. 125, pp. 319–327, dec 2016.
  • [14] A. Lherbier, S. M.-M. Dubois, X. Declerck, Y.-M. Niquet, S. Roche, and J.-C. Charlier, “Transport properties of graphene containing structural defects,” Physical Review B, vol. 86, p. 075402, aug 2012.
  • [15] A. W. Cummings, D. L. Duong, V. L. Nguyen, D. Van Tuan, J. Kotakoski, J. E. Barrios Vargas, et al., “Charge Transport in Polycrystalline Graphene: Challenges and Opportunities,” Advanced Materials, vol. 26, pp. 5079–5094, aug 2014.
  • [16] A. Cresti, N. Nemec, B. Biel, G. Niebler, F. Triozon, G. Cuniberti, et al., “Charge transport in disordered graphene-based low dimensional materials,” Nano Research, vol. 1, pp. 361–394, nov 2008.
  • [17] B. Javvaji, B. He, and X. Zhuang, “The generation of piezoelectricity and flexoelectricity in graphene by breaking the materials symmetries,” Nanotechnology, vol. 29, p. 225702, jun 2018.
  • [18] T. O. Wehling, K. S. Novoselov, S. V. Morozov, E. E. Vdovin, M. I. Katsnelson, A. K. Geim, et al., “Molecular Doping of Graphene,” Nano Letters, vol. 8, pp. 173–177, jan 2008.
  • [19] X. Miao, S. Tongay, M. K. Petterson, K. Berke, A. G. Rinzler, B. R. Appleton, et al., “High Efficiency Graphene Solar Cells by Chemical Doping,” Nano Letters, vol. 12, pp. 2745–2750, jun 2012.
  • [20] X. Wang, L. Zhi, and K. Müllen, “Transparent, Conductive Graphene Electrodes for Dye-Sensitized Solar Cells,” Nano Letters, vol. 8, pp. 323–327, jan 2008.
  • [21] F. Schedin, A. K. Geim, S. V. Morozov, E. W. Hill, P. Blake, M. I. Katsnelson, et al., “Detection of individual gas molecules adsorbed on graphene,” Nature Materials, vol. 6, pp. 652–655, sep 2007.
  • [22] D. Soriano, D. V. Tuan, S. M.-M. Dubois, M. Gmitra, A. W. Cummings, D. Kochan, et al., “Spin transport in hydrogenated graphene,” 2D Materials, vol. 2, p. 022002, may 2015.
  • [23] C. Feng, Z. Yi, L. F. Dumée, F. She, Z. Peng, W. Gao, and L. Kong, “Tuning micro-wrinkled graphene films for stretchable conductors of controllable electrical conductivity,” Carbon, vol. 139, pp. 672–679, nov 2018.
  • [24] K.-M. Hu, Z.-Y. Xue, Y.-Q. Liu, P.-H. Song, X.-H. Le, B. Peng, et al., “Probing built-in stress effect on the defect density of stretched monolayer graphene membranes,” Carbon, vol. 152, pp. 233–240, nov 2019.
  • [25] Z. Liu, Z. Qian, J. Song, and Y. Zhang, “Conducting and stretchable composites using sandwiched graphene-carbon nanotube hybrids and styrene-butadiene rubber,” Carbon, vol. 149, pp. 181–189, aug 2019.
  • [26] L.-C. Jia, L. Xu, F. Ren, P.-G. Ren, D.-X. Yan, and Z.-M. Li, “Stretchable and durable conductive fabric for ultrahigh performance electromagnetic interference shielding,” Carbon, vol. 144, pp. 101–108, apr 2019.
  • [27] Z. Yang, Y. Niu, W. Zhao, Y. Zhang, H. Zhang, C. Zhang, et al., “Nanoparticle intercalation-modulated stretchable conductive graphene fibers with combined photoelectric properties,” Carbon, vol. 141, pp. 218–225, jan 2019.
  • [28] C. Liu, B. Yao, T. Dong, H. Ma, S. Zhang, J. Wang, et al., “Highly stretchable graphene nanoribbon springs by programmable nanowire lithography,” npj 2D Materials and Applications, vol. 3, p. 23, dec 2019.
  • [29] S. Plimpton, “Fast Parallel Algorithms for Short-Range Molecular Dynamics,” Journal of Computational Physics, vol. 117, pp. 1–19, mar 1995.
  • [30] L. Lindsay and D. A. Broido, “Optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene,” Physical Review B, vol. 81, p. 205441, may 2010.
  • [31] B. Mortazavi, Z. Fan, L. F. C. Pereira, A. Harju, and T. Rabczuk, “Amorphized graphene: A stiff material with low thermal conductivity,” Carbon, vol. 103, pp. 318–326, jul 2016.
  • [32] M. Ishigami, J. H. Chen, W. G. Cullen, M. S. Fuhrer, and E. D. Williams, “Atomic Structure of Graphene on SiO 2,” Nano Letters, vol. 7, no. 6, pp. 1643–1648, 2007.
  • [33] B. Mortazavi, A. Lherbier, Z. Fan, A. Harju, T. Rabczuk, and J. C. Charlier, “Thermal and electronic transport characteristics of highly stretchable graphene kirigami,” Nanoscale, vol. 9, no. 42, pp. 16329–16341, 2017.
  • [34] R. Liu, J. Zhao, L. Wang, and N. Wei, “Nonlinear vibrations of helical graphene resonators in the dynamic nano-indentation testing,” Nanotechnology, vol. 31, no. 2, p. 025709, 2020.
  • [35] W. Humphrey, A. Dalke, and K. Schulten, “Vmd: Visual molecular dynamics,” Journal of Molecular Graphics, vol. 14, no. 1, pp. 33–38, 1996.
  • [36] Y. Chu, T. Ragab, and C. Basaran, “The size effect in mechanical properties of finite-sized graphene nanoribbon,” Computational Materials Science, vol. 81, pp. 269–274, 2014.
  • [37] J. Chen, J. H. Walther, and P. Koumoutsakos, “Strain engineering of kapitza resistance in few-layer graphene,” Nano Letters, vol. 14, no. 2, pp. 819–825, 2014.
  • [38] P. Z. Hanakata, E. D. Cubuk, D. K. Campbell, and H. S. Park, “Accelerated Search and Design of Stretchable Graphene Kirigami Using Machine Learning,” Physical Review Letters, vol. 121, no. 25, p. 255304, 2018.
  • [39] P. Z. Hanakata, E. D. Cubuk, D. K. Campbell, and H. S. Park, “Erratum: Accelerated Search and Design of Stretchable Graphene Kirigami Using Machine Learning [Phys. Rev. Lett. 121, 255304 (2018)],” Physical Review Letters, vol. 123, p. 69901 (E), 2019.
  • [40] L. F. C. Pereira and D. Donadio, “Divergence of the thermal conductivity in uniaxially strained graphene,” Physical Review B, vol. 87, no. 12, p. 125424, 2013.
  • [41] Z. Fan, L. F. C. Pereira, H.-Q. Wang, J.-C. Zheng, D. Donadio, and A. Harju, “Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations,” Physical Review B, vol. 92, p. 094301, sep 2015.
  • [42] N. Yang, X. Ni, J.-W. Jiang, and B. Li, “How does folding modulate thermal conductivity of graphene?,” Applied Physics Letters, vol. 100, no. 9, p. 093107, 2012.
  • [43] J. Zhao, J. Wu, J.-W. Jiang, L. Lu, Z. Zhang, and T. Rabczuk, “Thermal conductivity of carbon nanocoils,” Applied Physics Letters, vol. 103, no. 23, p. 233511, 2013.
  • [44] Z. Guo, D. Zhang, and X.-G. Gong, “Thermal conductivity of graphene nanoribbons,” Applied Physics Letters, vol. 95, no. 16, p. 163103, 2009.
  • [45] M. T. Ong and E. J. Reed, “Engineered Piezoelectricity in Graphene,” ACS Nano, vol. 6, no. 2, pp. 1387–1394, 2012.
  • [46] S. Chandratre and P. Sharma, “Coaxing graphene to be piezoelectric,” Applied Physics Letters, vol. 100, no. 2, pp. 12–15, 2012.
  • [47] Alamusi, J. Xue, L. Wu, N. Hu, J. Qiu, C. Chang, S. Atobe, H. Fukunaga, T. Watanabe, Y. Liu, H. Ning, J. Li, Y. Li, and Y. Zhao, “Evaluation of piezoelectric property of reduced graphene oxide (rGO)–poly(vinylidene fluoride) nanocomposites,” Nanoscale, vol. 4, no. 22, p. 7250, 2012.
  • [48] F. Ahmadpoor and P. Sharma, “Flexoelectricity in two-dimensional crystalline and biological membranes,” Nanoscale, vol. 7, no. 40, pp. 16555–16570, 2015.
  • [49] S. I. Kundalwal, S. A. Meguid, and G. J. Weng, “Strain gradient polarization in graphene,” Carbon, vol. 117, pp. 462–472, 2017.
  • [50] M. B. Ghasemian, T. Daeneke, Z. Shahrbabaki, J. Yang, and K. Kalantar-Zadeh, “Peculiar piezoelectricity of atomically thin planar structures,” Nanoscale, vol. 12, no. 5, pp. 2875–2901, 2020.
  • [51] T. Dumitrică, C. M. Landis, and B. I. Yakobson, “Curvature-induced polarization in carbon nanoshells,” Chemical Physics Letters, vol. 360, no. 1-2, pp. 182–188, 2002.
  • [52] X. Zhuang, B. He, B. Javvaji, and H. S. Park, “Intrinsic bending flexoelectric constants in two-dimensional materials,” Physical Review B, vol. 99, no. 5, p. 054105, 2019.
  • [53] A. Mayer, “Polarization of metallic carbon nanotubes from a model that includes both net charges and dipoles,” Physical Review B, vol. 71, no. 23, p. 235333, 2005.
  • [54] A. Mayer, “Formulation in terms of normalized propagators of a charge-dipole model enabling the calculation of the polarization properties of fullerenes and carbon nanotubes,” Physical Review B, vol. 75, no. 4, p. 045407, 2007.