Compound influence of surface and flexoelectric effects on static bending response of hybrid composite nanorod

Abstract

Nanoscale beams and rods are extensively used in several nano-electro-mechanical systems (NEMS) and their applications such as sensors and actuators. The surface and flexoelectricity phenomena have an extensive effect on nanosized structures and are related to their scale-dependent characteristics. This article presents the effect of different surface parameters and flexoelectricity on the electrostatic response of graphene-reinforced hybrid composite (GRHC) nanorods (NRs) using the theory of linear piezoelectricity, Euler-Bernoulli (EB), and Galerkin residual method. Based on these theories, the theoretical and finite element (FE) model is produced to investigate the static bending deflection of GRHC NRs when subjected to point and uniformly distributed load (UDL) considering different boundary conditions: cantilever (FC), fixed-fixed (FF), and simply supported (SS). This proposed FE model provides a useful tool for analyzing and investigating the outcomes of analytical models, which are found to be in good agreement. Our results presented in this article reveal that the effect of surface and flexoelectricity on the static bending response of GRHC NRs is noteworthy. These effects diminish with increased thickness/diameter of NR, and hence, these effects can be neglected for large-sized structures. The results presented here would help to identify the desired electrostatic response of GRHC NRs in terms of static bending response for a range of NEMS using different loading and boundary conditions as well as graphene volume fraction. This current study offers pathways for developing new proficient novel GRHC materials with enhanced control authority and present models can be exploited for numerous other materials as well as line-type structural systems such as beams, wires, rods, column/piers, and piles to study their global response.

Introduction

Over the last few decades, there is increasing research attention in flexible piezoelectric structures because of their inherent electro-mechanical coupling aspects that offer various opportunities to develop multifunctional next-generation nano-electro-mechanical systems (NEMS) with energy harvesting applications.^1,2 The smart piezoelectric composite nanostructures show size-dependent properties such as flexoelectric and surface effects, which classical continuum mechanics cannot account for as compared to conventional bulk fiber-reinforced composite structures. In addition, due to the coupling of piezoelectricity and flexoelectricity at nano-level, the continuum and finite element (FE) modeling of such novel composite nanostructures get complicated. Hence, to make the optimum usage of the piezoelectric composite materials, it is important to acquire a well prior knowledge of these scale-dependent properties. Piezoelectric composite structures have widespread usage in mechanical, civil, marine, aerospace, and other industrial applications. A piezoelectric micro/nanoscale structure such as beams, plates, rods, and panels is found to be having massive capabilities for numerous applications including energy harvesters which contain sensors, actuators, parallel plate capacitors, and nanogenerators.^3–10 Nanoscale beams and rods are extensively used in several NEMS applications. Precisely from the past couple of years, piezoelectric/non-piezoelectric materials have fascinated great interest for producing energy harvesters to generate the electric response on the application of mechanical load to the structure. In contrast to the phenomenon of the piezoelectric effect, an impulsive electric response in terms of polarization is exploited by a strain gradient known as the flexoelectric effect and it is more significant to accomplish the purpose of energy harvesting. The flexoelectric effect occurs in insulating as well as dielectric materials^11,12 whereas the inversion symmetry plays a vital role. By breaking this inversion symmetry, one can generate the polarization even in the centrosymmetric crystals (i.e. in non-piezoelectric materials) because of the strain gradient. Such a gradient term shows substantial results at the nanolevel and hence the flexoelectric effect is acknowledged as a scale-dependent phenomenon. Firstly, the flexoelectric effect was studied in 1950s ¹³ but acknowledged very less in research interest for a long period because of its weak effect. Very recently from 2000s, the flexoelectric effect captivated much more attention from the scientific and nanoscale application point of view. Tagantsev ¹¹ investigated the flexoelectric coefficients of dielectrics, ferroelectric, and materials such as relaxors and he observed that these coefficients depend on dielectric coefficient. These results confirmed the estimations made by Ma and Cross ¹⁴ through experimental characterization. Maranganti et al. ¹⁵ developed a comprehensive framework based on variational principle for dielectrics considering the flexoelectricity and offered solutions for the governing relations of a centrosymmetric continuum medium through Green’s function. Meanwhile, the advantage of flexoelectricity over piezoelectricity is, the piezoelectric effect which only occurs in 20 non-centrosymmetric point groups while the flexoelectric effect occurs in all dielectric with 32 crystalline point groups, interesting fact about the flexoelectricity is that one can generate piezoelectricity even in the non-piezoelectric materials. ¹⁶ For instance, Majdoub et al. ¹⁷ observed the flexoelectric effect for energy harvesters such as piezoelectric beam/ribbon due to substantial enrichment of the piezoelectric coefficient. In addition to flexoelectricity, Shen and Hu ¹⁸ recognized a general formulation for dielectrics with surface effects taken into account at nanoscale. In case of unimorph piezoelectric energy harvesters, Wang and Wang ¹⁹ proposed an analytical model with consideration of the flexoelectricity at nanolevel. Their outcomes revealed that flexoelectricity shows a noteworthy effect for piezoelectric cantilever beams in the energy harvesting application.

On the other hand, in addition to flexoelectricity, the influence of surface parameters is extensively renowned to substantially affect the physico-mechanical characteristics of nanostructures. The effects of different parameters of surface such as residual stresses, surface modulus, lame parameters, and surface piezoelectricity engrossed great research attention from fundamental applications. For example, Gurtin and Ian Murdoch ²⁰ proposed the continuum model considering the linear surface elasticity theory. They assumed the deformable surface having zero or negligible thickness, which is adhered to bulk material considering the perfect bond amongst the surface and bulk material. The surface effect can also be known as a size-dependent phenomenon. Due to reduced geometrical dimensions to the nanoscale, the surface effects are primarily liable to produce the enhanced scale-dependent electro-mechanical response from the base material. To distinguish the crystal and amorphous structure, Izumi et al. ²¹ carried out the molecular dynamic (MD) simulation for evaluating the surface elastic coefficient as well as surface stresses. As the values/magnitudes of surface elastic moduli and surface stresses are different for different elements such as nickel (Ni), aluminum (Al), palladium (Pd), silver (Ag), platinum (Pt), copper (Cu), gold (Au), etc., including positive and negative sign. It depends upon the surface orientations of elements and crystal structures (crystallographic direction). For example, Shenoy ²² studied the different elements for values of surface residual stresses and elastic moduli. They showed the negative value of surface stresses for nickel in case of specific crystallographic direction and crystal surface while all other elements (Al, Ag, Cu, Au, Pd, Pt) possess positive value. The influence of surface stress on static/dynamic characteristics of elastic and piezoelectric materials is investigated by several researchers,^23–26 and they observed that the surface effect will influence the performance of energy-harvesting elements as its size is minimized to nanolevel. Several authors^21,25,27,28 investigated the approximate approach to compute the surface elastic constant and Ru ²⁹ has predicted the same approach by considering the surface layer thickness and thermoelastic dissipation factor into account for the analysis of nanowire. Afterward, for demonstrating the influence of surface parameters including surface lame constants and residual surface stresses on the bending/modal analysis of circular and rectangular plates, Liu and Rajapakse ³⁰ presented the FE and analytical methods. Chen ³¹ considered the consequence of a plane boundary of a piezoelectric body molded as a thin layer means surface layer with quantified bulk properties and derived the state-space formulation for a transfer relation concerning the state vectors at the top/bottom surface. The author also stated the relationship amongst the surface layer thickness and the surface piezoelectric coefficient. Yan and Jiang ³² have adopted the Euler-Bernoulli (EB) model for static bending analysis of the piezoelectric cantilever beam by taking the residual surface effects into account, but they have not considered the effect of flexoelectricity. Pan et al. ³³ and Nan and Wang ³⁴ introduced the characteristic length to specify the surface constant over the bulk material properties. From the study of extensive literature related to classical continuum mechanics, it is observed that these theories are not able to consider the small-scale effect of nano-scaled structures due to the absence of different scale parameters. Considering the insufficiencies of classical continuum theories, the higher-order non-classical continuum theories such as non-local elasticity theory (NET) of Eringen, shear deformation theory (SDT), nonlocal strain gradient theory (NSGT), and modified coupled stress theory (MCST) which give more precise outcomes by considering size effects, have been recommended to investigate the different responses of nanostructures. For instance, Ebrahimi and Dabbagh^35–37 investigated the hygro-thermal and viscoelastic wave propagation properties of single- and double-layered graphene sheets (SLGSs and DLGSs) employing NET and NSGT. In addition, the value of wave frequency, phase velocity, and escape frequency of DLGSs were also obtained. Moreover, they performed parametric studies to examine the effect of different parameters such as length scale and nonlocal parameters, wave number, moisture concentration, temperature gradient, structural damping, Winkler and Pasternak coefficients, as well as axial load on the wave propagation response of SLGSs and DLGSs. Based on NET, Ebrahimi et al.^38–41 presented a new two-step porosity dependent homogenization technique to investigate the wave propagation responses of functionally graded (FG) porous nanobeams in the presence of axial pre-load. They reported that if the porosity volume fraction is increased then the dispersion responses of FG nanobeams were decreased. Afterward, they investigated the wave propagation response of smart magnetostrictive sandwich nanoplates based on NSGT. In this, they reported that for low wave numbers the magnetostriction significantly affect the dispersion responses of smart nanoplate. Using same NSGT framework, they investigated the viscoelastic behavior of nanostructures such as FG nanobeams and nanoplates which are affected by the relationship between nonlocal time and space. These studies reveal that an increase in the nonlocal parameter can reduce the loss factor in a significant way due to the coupling between the nonlocalities in the spatial and temporal domains. Most recently, Naskar et al. ⁴² introduced a very powerful semi-analytical “Extended Kantorovich method” in conjunction with MCST to analyze the static and dynamic response of flexoelectric FG plate at the nanoscale in the presence of surface and piezoelectric effect.

Novoselov et al. ⁴³ performed the ground-breaking experimental analysis of a single layer graphene sheet which fascinated enormous response from both academia and industry. It is also called mankind of the century because of its distinctive thermo-electro-mechanical properties. As graphene is made from allotrope of carbon, each atom is contributed to a chemical reaction from both sides, which is attributed to its two-dimensional (2-D) structure. In the present day, graphene becomes broadly acknowledged as the utmost outstanding material to form the graphene-based composite due to its high mechanical stiffness. To improve and economize the efficiency of composite/matrix, graphene nanosheets (GNS), graphene/graphite oxide (GO), and graphene platelets (GPLs) which are derivatives of graphene were surrounded into the different types of matrix.^44–49 Whereas Rafiee et al. ⁴⁸ accomplished the systematic experimentation to analyze the buckling behavior of graphene-based composite beam and they conveyed there is ∼52% enrichment of buckling capacity of the composite beam by adding graphene nanofillers with 0.1% weight fraction into the epoxy. Using algebraic polynomials and the Ritz method, the thermal post-buckling for examining the nonlinear thermal stability of graphene-based composite beams under constant temperature raise was performed by Kiani and Mirzaei ⁵⁰ and Zhang et al. ⁵¹ Zhang et al. ⁵¹ presented the mechanical analysis including vibration, bending, and buckling of GO-based FG beams using the theory of first-order shear deformation. In this, using the modified Halpin-Tsai (HT) method, they estimated the effective mechanical properties of the graphene-based composite. Several authors^52–54 studied the static and dynamic analysis of advanced composite, multi-material lattices, and FG materials by using different theories as well as a powerful machine learning tool for stochastic analysis. Wang et al. ⁵⁵ established a 2-D elastic model with consideration of the constant distribution of graphene in each layer and also shear strain as well as the thickness strain, that is, plane stress state in each layer for free vibration and bending analysis of layered graphene composite beams. In conjunction with these investigations, interphase is a multifaceted area that can be used amongst the matrix and reinforced fiber for finding additional upgraded effective properties of the proposed composite. For finding more accurate as well as precise properties of the composite, no slippage assumption is generally taken into consideration amongst the reinforcement, matrix, and interphase. In that manner, Chen et al. ⁵⁶ investigated the composite surrounded by GO to study gradient interphase not only to improve the distribution of carbon fiber/epoxy interface but also stress transfer features. By using different micromechanical and numerical models, Shingare and Naskar ⁵⁷ predicted the overall piezoelastic properties of hybrid graphene-reinforced composite. They also reported lead zirconate titanate (PZT-5H)-based epoxy composite with and without considering graphene interphase. They showed a significant enhancement because of graphene nanofillers in traditional composites. Thus, graphene can be utilized as interphase and nanofillers to enhance the effective properties of composites. Nowadays, apart from conventional composite and FG structures, hybrid nanocomposite attracted great research interest due to the incorporation of micro- and nano-scaled reinforcement in the matrix phase. For instance, Ebrahimi et al.^58–62 performed extensive research investigation on static stability, free vibration, and postbuckling analysis of multi-scale hybrid (MSH) nanocomposites and their different structural elements such as beams, plates, and shells composed of both macro- and nano-scale reinforcements such as carbon fiber (CF), glass fiber (GF), GO powder, and carbon nanotube (CNT) in a polymeric matrix. Most recently, Ebrahimi and Dabbagh ⁶³ and Ebrahimi et al.^64,65 presented an analytical solution for investigating the free vibration and buckling analysis of MSH nanocomposites plates reinforced with CNT nanoparticles. They considered the different micromechanical models such as rules-of-mixture, modified HT, and Eshelby-Mori-Tanaka model to compute the effective properties. In these above studies, they employed different analytical and FE theories such as EB beam theory, classical plate theory, Navier’s solution, Galerkin’s and Rayleigh-Ritz method based new refined higher-order SDT to envisage the buckling and vibration response of MSH nanoplates subjected to clamped-clamped and simply supported edge support conditions. Moreover, they also considered the effect of different parameters such as porosity, agglomeration of CNTs, straight and wavy CNTs as well as effect of viscoelastic properties of polymer and wavy shape of the CNTs to acquire desired response of MSH nanocomposite beams, plates, and shells. Shingare and Naskar ⁶⁶ also investigated the piezoelectric and surface effects on a hybrid graphene-reinforced composite plate to study its static/dynamic responses. The above all studies related to FG, MSH, and advanced composites were mainly focused on static and dynamic behavior, but they didn’t consider the electromechanical response considering size-dependent properties such as flexoelectric and surface effects.

Taking inspiration from the concept of MSH nanocomposites, authors considered three-phase hybrid nanocomposite, namely, “Graphene-Reinforced Hybrid Composite (GRHC)” composed of graphene as nanofillers or interphase and PZT-5H as active piezoelectric fiber incorporated in the epoxy matrix. On careful examining the available scientific literature, it is noticed that no single study is available on such piezoelectric GRHC nanorod considering size-dependent properties such as flexoelectric and surface effects to improve its electromechanical properties which is an undeniable inspiration behind the present research. The novelty of the current study proposes the analysis of static bending analysis of flexoelectric GRHC nanorod (NR) (hereinafter the “hybrid NR”) considering different boundary conditions (BCs): cantilever that is, free-clamped (FC), fixed-fixed (FF), and simply supported (SS) and loading conditions: point (P) and uniformly distributed load (UDL). According to the practical point of view, these loading and BCs are important in automotive, aerospace, and spacecraft engineering structures. These structures are mainly subjected to the penetration of external load and UDL with all mentioned BCs. The objective of the present work is to extend the base knowledge of the gradient of strain/electric response to promote the importance of surface and flexoelectric effects in the nanostructure. The surface effects comprising the surface stress, modulus, and piezoelectricity are integrated into the model. From the literature, it is also observed that the analytical models are restricted to simple geometries and BCs while these models are also far more complicated for systems with somewhat more complex geometry. It is also well-known that nanoparticles appear with irregular surfaces when their dimensions are in the nanometer range. Therefore, to further elucidate the surface effect and better characterization of the electrostatic behavior of nanostructures, the FE method (FEM) should be a better option. In case of intricate geometries of beam and BCs, usually accounted for NEMS applications, a multipurpose in-house FE model needs to be developed. However, the traditional FEM can only provide numerical solutions without considering the surface and flexoelectricity effect. This thing keeping in mind, we developed a new FE formulation in the present article by incorporating the surface parameters as well as piezoelectric and flexoelectric effects with the traditional one.^30,67 The Galerkin residual method is the most widely used technique to solve the problem with suitable approximation in analytical relations without using commercial FE software. Therefore, to avoid the complexity due to consideration of all dimensions of structures, we limited our FE analysis to one dimension (1-D). In the context of dimensional reduction, 1-D modeling always referred to line-type members such as beams, wires, rods, column, and piles, 2-D modeling referred to plate-type members such as walls, slabs, etc. (assuming suitable problem type: plane stress, plane strain, or axisymmetric) while 3-D modeling referred for the structural member by assuming appropriate dimensions to analyze prototype. For this reason, the authors developed an in-house FE code by assuming 1-D case based on the Galerkin residual method to validate the results obtained from an analytical formulation based on EB theory and compared both sets of results.

The results obtained from the present research work would offer new insights to engineer the domain configurations for tailoring the desired static electromechanical responses of the novel GRHC considering surface and flexoelectric effects. This has been demonstrated by comparison of different sets of results such as: (i) conventional composite nanorod (without surface and flexoelectric effects), (ii) flexo composite nanorod (considering only flexoelectric effect), and (iii) flexo-surface composite nanorod (considering surface and flexoelectric effects). From the computational development point of view, the impact of the above-mentioned aspects is quite significant. Thus, most importantly, the current results reveal that the incorporation of surface effect dominantly influences the electromechanical response of nanorods compared to the flexoelectric effect and both the effects should be considered to study the structural response of any nanostructure. Figure 1 illustrates the flowchart explaining the layout of the work carried out in the research article.

Figure 1.

Flowchart of electromechanical analysis of hybrid NRs.

Electromechanical response of GRHC nanorod

Effective properties of GRHC

The determination of effective elastic and piezoelectric properties of graphene-reinforced hybrid composite (GRHC) is required priori and therefore, its effective properties are predicted first. We considered graphene sheets as nanofillers or interphase and PZT-5H as active piezoelectric fiber incorporated in the epoxy matrix to form GRHC. Figure 2 demonstrates a single GRHC representative volume element (RVE) which is assumed to be extracted from composite laminate. Recently, Shingare and Naskar^6,57 studied the electrostatic response of GRHC and its structures accounting for the piezoelectric and flexoelectric effects. In this, they envisaged the overall effective properties of GRHC using analytical and numerical models such as two- and three-phase mechanics of materials (MOM) models, Halpin-Tsai (HT), rules-of-mixture (ROM), modified rules-of-mixture (mROM) as well as FE modeling. The effective properties of three-phase GRHC are evaluated by changing volume fractions of graphene interphase ( ${\mathrm{v}}_{\mathrm{g}}$ ) and PZT fibers ( ${\mathrm{v}}_{\mathrm{p}}$ ), that is, by taking the value of ${\mathrm{v}}_{\mathrm{g}}$ in terms of ${\mathrm{v}}_{\mathrm{p}}$ .

Figure 2.

Schematics of three-phase GRHC RVEs for estimating effective properties.

The detailed analytical and FE micromechanical model for the development of three-phase GRHC is not shown here and for more details regarding MOM and FE modeling, readers are referred to Shingare and Naskar. ⁵⁷ From this study, authors reported that these three-phase GRHCs show overall enhanced effective properties as compared to two-phase PZT-based composites, that is, without consideration of graphene interphase. Therefore, this three-phase hybrid composite is utilized for studying the static bending response of GRHC nanorod and their essential effective properties are mentioned in Table 2.

Governing equations

Considering the effects of flexoelectricity and surface parameters, the electric Gibb’s energy density function is characterized for materials with different effects: bulk and surface effect. In case of bulk material, the following internal energy density function is written by neglecting the higher-order strain gradient ⁶⁸ :

$${\mathrm{U}}_{\mathrm{b}}=\frac{1}{2}{\sigma }_{\mathrm{ij}}{\epsilon }_{\mathrm{ij}}+\frac{1}{2}{\tau }_{\mathrm{ijk}}{\eta }_{\mathrm{ij},\mathrm{k}}-\frac{1}{2}{\mathrm{D}}_{\mathrm{i}}{\mathrm{E}}_{\mathrm{i}}$$

(1)

wherein ${\epsilon }_{\mathrm{ij}}$ and ${\sigma }_{\mathrm{ij}}$ signify the traditional strain and stress tensor. ${\eta }_{\mathrm{ij},\mathrm{k}}$ and ${\tau }_{\mathrm{ijk}}$ signify the higher-order strain and stress gradient, ${\mathrm{D}}_{\mathrm{i}}$ and ${\mathrm{E}}_{\mathrm{i}}$ denote the electric flux density and field vector. These components are expressed as:

$${\epsilon }_{\mathrm{ij}}=\frac{1}{2}({\mathrm{u}}_{\mathrm{j},\mathrm{i}}+{\mathrm{u}}_{\mathrm{i},\mathrm{j}})$$

(2a)

$${\eta }_{\mathrm{ij},\mathrm{k}}={\epsilon }_{\mathrm{ij},\mathrm{k}}=\frac{1}{2}({\mathrm{u}}_{\mathrm{j},\mathrm{ik}}+{\mathrm{u}}_{\mathrm{i},\mathrm{jk}})$$

(2b)

$${\mathrm{E}}_{\mathrm{k}}=-{Ø}_{,\mathrm{k}}$$

(2c)

whereas ${\mathrm{u}}_{\mathrm{i}}$ and Ø represent the component of mechanical displacement and electrostatic potential.

In case of bulk materials, the constitutive equations can be written as:

$${\sigma }_{\mathrm{ij}}=\frac{\partial {\mathrm{U}}_{\mathrm{b}}}{\partial {\epsilon }_{\mathrm{ij}}}={\mathrm{C}}_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}}-{\mathrm{e}}_{\mathrm{kij}}{\mathrm{E}}_{\mathrm{k}}$$

(3a)

$${\tau }_{\mathrm{jkl}}=\frac{\partial {\mathrm{U}}_{\mathrm{b}}}{\partial {\eta }_{\mathrm{jkl}}}=-{\mu }_{\mathrm{ijkl}}{\mathrm{E}}_{\mathrm{i}}$$

(3b)

$${\mathrm{D}}_{\mathrm{i}}=-\frac{\partial {\mathrm{U}}_{\mathrm{b}}}{\partial {\mathrm{E}}_{\mathrm{i}}}={\chi }_{\mathrm{ij}}{\mathrm{E}}_{\mathrm{j}}+{\mathrm{e}}_{\mathrm{ijk}}{\epsilon }_{\mathrm{jk}}+{\mu }_{\mathrm{ijkl}}{\eta }_{\mathrm{jkl}}$$

(3c)

where ${\chi }_{\mathrm{ij}}$ , ${\mathrm{C}}_{\mathrm{ijkl}}$ , ${\mathrm{e}}_{\mathrm{ijk}}$ , and ${\mu }_{\mathrm{ijkl}}$ are the dielectric susceptibility, elastic, piezoelectric, and flexoelectric coefficient, respectively.

In case of piezoelectric continuum, the surface internal energy density function is deduced as ¹⁸ :

$${\mathrm{U}}_{\mathrm{s}}=\frac{1}{2}{\sigma }_{\alpha \beta }^{\mathrm{o}}{\epsilon }_{\alpha \beta }^{\mathrm{s}}+\frac{1}{2}{\sigma }_{\alpha \beta }^{\mathrm{s}}{\epsilon }_{\alpha \beta }^{\mathrm{s}}-\frac{1}{2}{\mathrm{D}}_{\gamma }^{\mathrm{s}}{\mathrm{E}}_{\gamma }^{\mathrm{s}}$$

(4)

where ${\sigma }_{\alpha \beta }^{\mathrm{o}}$ indicates the residual surface stress tensor of the second order, ${\epsilon }_{\alpha \beta }^{\mathrm{s}}$ and ${\sigma }_{\alpha \beta }^{\mathrm{s}}$ denote the surface strain and stress; ${\mathrm{D}}_{\gamma }^{\mathrm{s}}$ and ${\mathrm{E}}_{\gamma }^{\mathrm{s}}$ denote the surface electric flux density and field and it can be obtained as:

$${\epsilon }_{\alpha \beta }^{\mathrm{s}}=\frac{1}{2}({\mathrm{u}}_{\alpha ,\beta }^{\mathrm{s}}+{\mathrm{u}}_{\beta ,\alpha }^{\mathrm{s}})$$

(5a)

$${\mathrm{E}}_{\gamma }^{\mathrm{s}}=-{Ø}^{\mathrm{s}},\gamma$$

(5b)

whereas ${\mathrm{u}}^{\mathrm{s}}$ signifies the surface displacement and ${Ø}^{\mathrm{s}}$ signifies the surface electric potential. The constitutive expressions in case of surface effects (with superscript “s”) are approximately identical as that of constitutive relations for bulk material but some residual terms are also present.

Using equation (4), the linear constitutive relations are determined considering surface effect ¹⁸ :

$${\sigma }_{\alpha \beta }^{\mathrm{s}}=\frac{\partial {\mathrm{u}}_{\mathrm{s}}}{\partial {\epsilon }_{\alpha \beta }}={\tau }_{\alpha \beta }+{\mathrm{C}}_{\alpha \beta \gamma \mathrm{k}}^{\mathrm{s}}{\epsilon }_{\gamma k}^{\mathrm{s}}-{\mathrm{e}}_{\mathrm{k}\mathrm{\alpha }\mathrm{\beta }}^{\mathrm{s}}{\mathrm{E}}_{\mathrm{k}}^{\mathrm{s}}$$

(6a)

$${\mathrm{D}}_{\gamma }^{\mathrm{s}}=-\frac{\partial {\mathrm{u}}_{\mathrm{s}}}{\partial {\mathrm{E}}_{\gamma }^{\mathrm{s}}}={\chi }_{\gamma \mathrm{k}}^{\mathrm{s}}{\mathrm{E}}_{\mathrm{k}}^{\mathrm{s}}+{\mathrm{e}}_{\gamma \mathrm{k}\beta }^{\mathrm{s}}{\epsilon }_{\alpha \beta }^{\mathrm{s}}$$

(6b)

To predict the static response of NRs with a ratio of thickness to length $(h/L\ge 1/20)$ , thin nanobeam theory is utilized while thick nanobeam theory is utilized if this ratio decreases with shear deformation effects taken into consideration. Hence, the required ratio is $h/L<1/10$ for the thick hybrid rod. Without any mechanical load in the axial direction (Figure 3) the axial force becomes zero, hence, a displacement field can be obtained by using EB theory:

$$\mathrm{u}(\mathrm{x})=-\mathrm{z}\frac{\mathrm{dv}(\mathrm{x})}{\mathrm{dx}}=-\mathrm{z}{\mathrm{v}}^{\prime }(\mathrm{x})$$

(7)

Here, $\mathrm{u}(\mathrm{x})$ and $\mathrm{v}(\mathrm{x})$ denote the horizontal and vertical displacement of the NR.

Figure 3.

NRs subjected to UDL with different BCs: (a) FC, (b) FF, and (c) SS.

Therefore, the following nonzero strain and strain gradients are obtained by using equations (2) and (7).

$$\begin{array}{c}\multicolumn{1}{c}{{\epsilon }_{11}=-\mathrm{z}\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2}=-\mathrm{z}{\mathrm{v}}^{″};{\eta }_{111}=-\mathrm{z}\frac{{\mathrm{d}}^3\mathrm{v}}{\mathrm{d}{\mathrm{x}}^3}=-\mathrm{z}{\mathrm{v}}^{″\prime };}\\ \multicolumn{1}{c}{{\eta }_{311}=-\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2}=-{\mathrm{v}}^{″}}\end{array}$$

(8)

The internal energy density function for bulk material is again deduced by using equations (3) and (8) into equation (1) ⁶⁸ :

$$\begin{array}{c}\delta {\mathrm{U}}_{\mathrm{b}}=\underset{\mathrm{v}}{\int }({\sigma }_{11}\delta {\epsilon }_{11}+{\tau }_{311}\delta {\eta }_{311}+{\tau }_{111}\delta {\eta }_{111})\mathrm{dV}\hfill \\ =-\underset{0}{\overset{\mathrm{L}}{\int }}\mathrm{M}\mathrm{\delta }{\mathrm{v}}^{″}\mathrm{dx}-\underset{0}{\overset{\mathrm{L}}{\int }}{\mathrm{M}}_{\mathrm{h}}\delta {\mathrm{v}}^{″\prime }\mathrm{dx}-\underset{0}{\overset{\mathrm{L}}{\int }}\mathrm{P}\mathrm{\delta }{\mathrm{v}}^{″}\mathrm{dx}\hfill \end{array}$$

(9)

wherein,

$$\mathrm{M}=\underset{\mathrm{A}}{\int }{\sigma }_{11}\mathrm{zdA};\mathrm{P}=\underset{\mathrm{A}}{\int }{\tau }_{311}\mathrm{dA};\mathrm{and}{\mathrm{M}}_{\mathrm{h}}=\underset{\mathrm{A}}{\int }{\tau }_{111}\mathrm{zdA}$$

(10)

whereas $\underset{\mathrm{A}}{\int }\mathrm{dA}$ indicates the integration over the whole area “A” in case of bulk material.

Using equations (5)–(8) into equation (4), the internal energy density function is again re-expressed as:

$$\begin{array}{c}{\mathrm{U}}_{\mathrm{s}}=\underset{\mathrm{a}}{\int }\frac{{\mathrm{d}\mathrm{\sigma }}_{11}^{\mathrm{s}}}{\mathrm{dx}}\mathrm{z}\mathrm{\delta }{\mathrm{v}}^{\prime }\mathrm{da}-\underset{\mathrm{z}=\frac{\mathrm{H}}{2}}{\int }({\sigma }_{11}^{\mathrm{s}}\mathrm{k}{)}^{\mathrm{u}}\delta \mathrm{v}\mathrm{da}\hfill \\ \multicolumn{1}{c}{+\underset{\mathrm{z}=\frac{-\mathrm{H}}{2}}{\int }({\sigma }_{11}^{\mathrm{s}}\mathrm{k}{)}^{\mathrm{l}}\delta \mathrm{vda}}\\ =\underset{0}{\overset{\mathrm{L}}{\int }}\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{s}}}{\mathrm{dx}}\delta {\mathrm{v}}^{\prime }\mathrm{dx}-\underset{0}{\overset{\mathrm{L}}{\int }}{\mathrm{T}}_{\mathrm{z}}^{\mathrm{s}}\delta \mathrm{vdx}\hfill \end{array}$$

(11)

where $\underset{\mathrm{a}}{\int }\mathrm{da}$ indicates the integration over the small infinitesimal surface area “a.”

By using equation (6a), the axial surface stress ${\sigma }_{11}^{\mathrm{s}}$ can be expressed as:

$${\sigma }_{11}^{\mathrm{s}}={\tau }_0+{\mathrm{C}}_{11}^{\mathrm{s}}{\epsilon }_{\mathrm{x}}^{\mathrm{s}}-{\mathrm{e}}_{31}^{\mathrm{s}}{\mathrm{E}}_{\mathrm{z}}^{\mathrm{s}}.$$

(12a)

For surface effect, ${\sigma }_{11}^{\mathrm{s}}$ can be rewritten as:

$${\sigma }_{11}^{\mathrm{s}}={\tau }_0+{\mathrm{C}}_{11}^{\mathrm{s}}(-\mathrm{z}\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2})-{\mathrm{e}}_{31}^{\mathrm{s}}{\mathrm{E}}_{\mathrm{z}}^{\mathrm{s}},$$

(12b)

in which ${\tau }_0$ is the constant of residual surface stress.

wherein ${\mathrm{M}}_{\mathrm{s}}$ and ${\mathrm{T}}_{\mathrm{z}}^{\mathrm{s}}$ denote the bending moment (BM) and lateral loadings considering the surface effect and can be obtained as:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{M}}_{\mathrm{s}}=\underset{\mathrm{c}}{\int }{\sigma }_{11}^{\mathrm{s}}\mathrm{z}\mathrm{dC};{\mathrm{T}}_{\mathrm{z}}^{\mathrm{s}}=\underset{\mathrm{z}=\frac{\mathrm{H}}{2}}{\int }({\sigma }_{11}^{\mathrm{s}}\mathrm{k}{)}^{\mathrm{u}}\mathrm{da}}\\ \multicolumn{1}{c}{-\underset{\mathrm{z}=\frac{-\mathrm{H}}{2}}{\int }({\sigma }_{11}^{\mathrm{s}}\mathrm{k}{)}^{\mathrm{l}}\mathrm{da}}\end{array}$$

(13)

whereas $\mathrm{C}$ denotes the perimeter of NR and superscript “u” and “l” denote the upper and lower surface of NRs, respectively. Here, $\mathrm{k}$ denotes the curvature which is determined by $\mathrm{k}=\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2}$ when considering the EB theory.

If the NR is subjected to uniformly transverse load $\mathrm{q}\left(\mathrm{x}\right),$ end moment $\tilde{\mathrm{M}}$ and end force $\mathrm{Q}$ , the virtual work done induced due to the external forces is acquired as:

$$\delta \mathrm{W}=\underset{0}{\overset{\mathrm{L}}{\int }}\mathrm{q}\left(\mathrm{x}\right)\delta \mathrm{vdx}+\mathrm{Q}\mathrm{\delta }\mathrm{v}+\tilde{\mathrm{M}}\delta {\mathrm{v}}^{\prime }$$

(14)

The variation principle is expressed using equations (9), (11), and (14):

$$\delta \mathrm{W}-(\delta {\mathrm{U}}_{\mathrm{b}}+\delta {\mathrm{U}}_{\mathrm{s}})=0$$

(15a)

$$\begin{array}{c}\delta \mathrm{W}-(\delta {\mathrm{U}}_{\mathrm{b}}+\delta {\mathrm{U}}_{\mathrm{s}})=\underset{0}{\overset{\mathrm{L}}{\int }}\mathrm{q}\left(\mathrm{x}\right)\delta \mathrm{vdx}+\mathrm{Q}\delta \mathrm{v}+\tilde{\mathrm{M}}\delta {\mathrm{v}}^{\prime }-\hfill \\ (\underset{0}{\overset{L}{\int }}\mathrm{M}\delta v^{″}\mathrm{dx}+\underset{0}{\overset{L}{\int }}{\mathrm{M}}_{\mathrm{h}}\delta {\mathrm{v}}^{″\prime }\mathrm{dx}+\underset{0}{\overset{L}{\int }}\mathrm{P}\delta {\mathrm{v}}^{″}\mathrm{dx}\hfill \\ \multicolumn{1}{c}{+\underset{0}{\overset{\mathrm{L}}{\int }}\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{s}}}{\mathrm{dx}}\mathrm{\delta }{\mathrm{v}}^{\prime }\mathrm{dx}-\underset{0}{\overset{\mathrm{L}}{\int }}{\mathrm{T}}_{\mathrm{z}}^{\mathrm{s}}\mathrm{\delta }\mathrm{vdx})}\end{array}$$

(15b)

The equation (15b) is again re-formulated by using integration by parts:

$$\begin{array}{c}\delta \mathrm{W}-(\delta {\mathrm{U}}_{\mathrm{b}}+\delta {\mathrm{U}}_{\mathrm{s}})={\mathrm{M}}_{\mathrm{h}}\delta {\mathrm{v}}^{″}|\begin{array}{c}\mathrm{L}\hfill \\ \multicolumn{1}{c}{0}\end{array}+\mathrm{Q}\mathrm{\delta }\mathrm{v}+\tilde{\mathrm{M}}\delta {\mathrm{v}}^{\prime }\hfill \\ \multicolumn{1}{c}{-(\frac{\mathrm{dM}}{\mathrm{dx}}+\frac{\mathrm{dP}}{\mathrm{dx}}-\frac{{\mathrm{d}}^2{\mathrm{M}}_{\mathrm{h}}}{\mathrm{d}{\mathrm{x}}^2}+\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{s}}}{\mathrm{dx}})\delta \mathrm{v}|\begin{array}{c}\mathrm{L}\hfill \\ \multicolumn{1}{c}{0}\end{array}}\\ +\underset{0}{\overset{\mathrm{L}}{\int }}(\frac{{\mathrm{d}}^2\mathrm{M}}{\mathrm{d}{\mathrm{x}}^2}+\frac{{\mathrm{d}}^2\mathrm{P}}{\mathrm{d}{\mathrm{x}}^2}-\frac{{\mathrm{d}}^3{\mathrm{M}}_{\mathrm{h}}}{\mathrm{d}{\mathrm{x}}^3}+\frac{{\mathrm{d}}^2{\mathrm{M}}_{\mathrm{s}}}{\mathrm{d}{\mathrm{x}}^2}+\mathrm{q}\left(\mathrm{x}\right)+{\mathrm{T}}_{\mathrm{z}}^{\mathrm{s}})\hfill \\ \multicolumn{1}{c}{\delta \mathrm{vdx}+({\mathrm{M}}_{\mathrm{T}}\mathrm{P}-\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{h}}}{\mathrm{dx}})\delta {\mathrm{v}}^{\prime }|\begin{array}{c}\mathrm{L}\hfill \\ \multicolumn{1}{c}{0}\end{array}}\end{array}$$

(16)

Because of the uncertainty of $\delta \mathrm{v}$ in equation (16), one can get governing expression as:

$$\frac{{\mathrm{d}}^2\mathrm{M}}{\mathrm{d}{\mathrm{x}}^2}+\frac{{\mathrm{d}}^2\mathrm{P}}{\mathrm{d}{\mathrm{x}}^2}-\frac{{\mathrm{d}}^3{\mathrm{M}}_{\mathrm{h}}}{\mathrm{d}{\mathrm{x}}^3}+\frac{{\mathrm{d}}^2{\mathrm{M}}_{\mathrm{s}}}{\mathrm{d}{\mathrm{x}}^2}+\mathrm{q}\left(\mathrm{x}\right)+{\mathrm{T}}_{\mathrm{z}}^{\mathrm{s}}=0$$

(17)

The following BCs are applied at the ends of a GRHC NRs:

$${\mathrm{M}}_{\mathrm{h}}$$

$$(\mathrm{M}+\mathrm{F}-\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{h}}}{\mathrm{dx}})$$

$$\frac{\mathrm{d}}{\mathrm{dx}}(\mathrm{M}+\mathrm{F}-\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{h}}}{\mathrm{dx}}+{\mathrm{M}}_{\mathrm{s}})$$

For example, if cantilever NR subjected to end point load (P) then corresponding BCs are written as follows (x = 0 and L):

At x = 0,

$$\mathrm{v}\left(0\right)={\mathrm{v}}^{\prime }\left(0\right)=0$$

$${\mathrm{M}}_{\mathrm{h}}=0$$

At x = L,

$${\mathrm{M}}_{\mathrm{h}}=0$$

$$\mathrm{M}\left(\mathrm{L}\right)=\mathrm{M}+\mathrm{F}-\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{h}}}{\mathrm{dx}}=0,$$

$$\mathrm{Q}\left(\mathrm{L}\right)=\frac{\mathrm{dM}}{\mathrm{dx}}+\frac{\mathrm{dF}}{\mathrm{dx}}-\frac{{\mathrm{d}}^2{\mathrm{M}}_{\mathrm{h}}}{\mathrm{d}{\mathrm{x}}^2}+\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{s}}}{\mathrm{dx}}=-\mathrm{P}.$$

whereas M, ${\mathrm{M}}_{\mathrm{h}}$ , and F indicate the BM, higher-order BM, and higher-order axial couple.

If we consider the open-circuit case, then the electric field is obtained using the concept of electric flux density $({\mathrm{D}}_{\mathrm{y}}\to 0)$ which is zero on the surface.

$${\mathrm{E}}_3=-\frac{{\mathrm{e}}_{31}}{{\chi }_{33}}{\epsilon }_{11}-\frac{{\mu }_{31}}{{\chi }_{33}}{\eta }_{311}$$

(18)

Using equations (3) and (8) into equation (18), the BM and the higher-order axial couple can be reformulated as:

$$\mathrm{M}=-({\mathrm{C}}_{11}+{\mathrm{e}}_{31}\frac{{\mathrm{e}}_{31}}{{\chi }_{33}}){\mathrm{I}}_{\mathrm{y}}{\mathrm{v}}^{″}$$

(19a)

$$\mathrm{F}=-\left(\frac{{\mu }_{31}^2}{{\chi }_{33}}\right){\mathrm{v}}^{″}\mathrm{A}$$

(19b)

whereas $\mathrm{A}$ and ${\mathrm{I}}_{\mathrm{y}}$ signify the cross-sectional area and the moment of inertia of GRHC NR. Equation (13) is again re-formulated as:

$${\mathrm{M}}_{\mathrm{s}}=-({\mathrm{C}}_{11}^{\mathrm{s}}+{\mathrm{e}}_{31}^{\mathrm{s}}\frac{{\mathrm{e}}_{31}}{{\chi }_{33}}){\mathrm{I}}^{*}{\mathrm{v}}^{″}$$

(20)

wherein ${\mathrm{I}}^{*}$ denotes the perimeter moment of inertia.

For Ex. NR with a circular (diameter D) and a rectangular (height H, width B) cross-section, the perimeter moment of inertia can be obtained as ³⁰ :

$${\mathrm{I}}^{*}=\{\begin{array}{c}\frac{{\mathrm{H}}^3}{6}+\frac{{\mathrm{BH}}^2}{2}\\ \frac{{\pi \mathrm{D}}^3}{8}\end{array}$$

The curvatures give the equivalent magnitudes but opposite directions on the top and bottom surfaces of the NR. Hence, by ignoring the nonlinear effect caused due to term ( ${\mathrm{e}}_{31}^{\mathrm{s}}{\mathrm{E}}_{\mathrm{z}}^{\mathrm{s}}$ ) in equation (13), the lateral loading ( ${\mathrm{T}}_{\mathrm{z}}^{\mathrm{s}}$ ) can be again deduced as:

$${\mathrm{T}}_{\mathrm{z}}^{\mathrm{s}}=\underset{\mathrm{z}=\frac{\mathrm{H}}{2}}{\int }({\sigma }_{11}^{\mathrm{s}}\mathrm{k}{)}^{\mathrm{u}}\mathrm{da}-\underset{\mathrm{z}=\frac{-\mathrm{H}}{2}}{\int }{\left({\sigma }_{11}^{\mathrm{s}}\mathrm{k}\right)}^{\mathrm{l}}\mathrm{da}={\mathrm{S}}^{*}{\tau }_0\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2},$$

(21)

whereas ${\mathrm{S}}^{*}=\{\begin{array}{c}2B\hfill \\ \multicolumn{1}{c}{\pi D/2}\end{array}$

By considering the flexoelectric and surface effect in equation (17), the governing equation is formulated as ⁶⁹ :

$${\left(\mathrm{BR}\right)}^{\mathrm{eff}}\frac{{\mathrm{d}}^4\mathrm{v}}{\mathrm{d}{\mathrm{x}}^4}={\mathrm{S}}^{*}{\tau }_0\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2}-\mathrm{q}\left(\mathrm{x}\right)$$

(22)

wherein

$$\begin{array}{c}\multicolumn{1}{c}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}=({\mathrm{C}}_{11}+\frac{{\mathrm{e}}_{311}^2}{{\chi }_{33}})\mathrm{I}+\left(\frac{{\mu }_{31}^2}{{\chi }_{33}}\right)\mathrm{A}}\\ \multicolumn{1}{c}{+({\mathrm{C}}_{11}^{\mathrm{s}}+\frac{{\mathrm{e}}_{311}^{\mathrm{s}}{\mathrm{e}}_{311}}{{\chi }_{33}}){\mathrm{I}}^{*}}\end{array}$$

(23)

Static loading of the hybrid NRs

Equation (22) can be re-expressed to determine the nondimensional static bending deflection ( $\overline{\mathrm{v}}=\mathrm{v}/\mathrm{L}$ ) with respect to the length of NR $(\overline{x}=\mathrm{x}/\mathrm{L})$ as:

$$\frac{{\mathrm{d}}^4\overline{\mathrm{v}}}{\mathrm{d}{\mathrm{x}}^4}-\Gamma \frac{{\mathrm{d}}^2\overline{\mathrm{v}}}{\mathrm{d}{\mathrm{x}}^2}+\frac{\mathrm{q}\left(\mathrm{x}\right){\mathrm{L}}^3}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}}=0$$

(24)

where as $\Gamma =\frac{{\mathrm{S}}^{*}{\tau }_0{\mathrm{L}}^2}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}}$

To get the generalized solution, equation (24) is again simplified for applied point/UDL which can be determined as follows:

$$\overline{\mathrm{v}}={\mathrm{D}}_1{\mathrm{e}}^{\overline{x}\sqrt{\Gamma }}+{\mathrm{D}}_2{\mathrm{e}}^{-\overline{x}\sqrt{\Gamma }}+{\mathrm{D}}_3+{\mathrm{D}}_4\overline{x}+\frac{{\mathrm{q}}_0{\mathrm{L}}^3}{2\Gamma {\left(\mathrm{BR}\right)}^{\mathrm{eff}}}{\overline{x}}^2$$

(25)

whereas ${\mathrm{D}}_1-{\mathrm{D}}_4$ are the arbitrary constants determined from the respective necessary BCs mentioned in following Table 1 in equation (25). If the NR is subjected to point load, then we can assume UDL as zero $({\mathrm{q}}_0=0)$ and vice versa.

Table 1.

Loading and boundary conditions applied on different GRHC NRs to obtain arbitrary constants.

Loadings	Boundary conditions	Arbitrary constants ${\mathrm{D}}_1-{\mathrm{D}}_4$	Equations
GRHC NRs subjected to point load $({\mathbf{q}}_0=0)$
FC NR subjectedto point load P	$\overline{\mathrm{v}}\left(0\right)={\overline{\mathrm{v}}}^{\prime }\left(0\right)=\mathrm{M}\left(1\right)=0;$ $\mathrm{Q}(1)=-\mathrm{P}$	${\mathrm{D}}_1=\frac{\mathrm{P}{\mathrm{L}}^2}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}({\mathrm{e}}^{2\sqrt{\Gamma }}+1)};$ ${\mathrm{D}}_2=-\frac{\mathrm{P}{\mathrm{L}}^2{\mathrm{e}}^{2\sqrt{\Gamma }}}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}({\mathrm{e}}^{2\sqrt{\Gamma }}+1)}$ ${\mathrm{D}}_3=\frac{\mathrm{P}{\mathrm{L}}^2({\mathrm{e}}^{2\sqrt{\Gamma }}-1)}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}({\mathrm{e}}^{2\sqrt{\Gamma }}+1)};$ ${\mathrm{D}}_4=-\frac{\mathrm{P}{\mathrm{L}}^2}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\Gamma }$	(26)
FF NR is subjected tomidpoint load P	$\overline{\mathrm{v}}(0)={\overline{\mathrm{v}}}^{\prime }(0)={\overline{\mathrm{v}}}^{\prime }(1/2)=0;$ $\mathrm{Q}(1/2)=-\mathrm{P}/2$	${\mathrm{D}}_1=\frac{{\mathrm{PL}}^2{\mathrm{e}}^{\raisebox{1ex}{$-\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}\left({\mathrm{e}}^{\raisebox{1ex}{$-\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+1\right)};$ ${\mathrm{D}}_2=\frac{{\mathrm{PL}}^2}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}\left({\mathrm{e}}^{\raisebox{1ex}{$-\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+1\right)};$ ${\mathrm{D}}_3=-\frac{{\mathrm{PL}}^2\left(1-{\mathrm{e}}^{\raisebox{1ex}{$-\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}\left({\mathrm{e}}^{\raisebox{1ex}{$-\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+1\right)};$ ${\mathrm{D}}_4=-\frac{\mathrm{P}{\mathrm{L}}^2}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\Gamma }$	(27)
SS NR is subjected to midpoint load P	$\overline{\mathrm{v}}(0)={\overline{\mathrm{v}}}^{\prime }(1/2)=\mathrm{M}(0)=0;$ $\mathrm{Q}(1/2)=-\mathrm{P}/2$	${\mathrm{D}}_1=\frac{{\mathrm{PL}}^2}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}\left(e^{\raisebox{1ex}{$\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+e^{\raisebox{1ex}{$-\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)};$ ${\mathrm{D}}_2=\frac{{\mathrm{PL}}^2}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}\left(e^{\raisebox{1ex}{$\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+e^{\raisebox{1ex}{$-\sqrt{\Gamma }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)};$ ${\mathrm{D}}_3=0;{\mathrm{D}}_4=-\frac{\mathrm{P}{\mathrm{L}}^2}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\Gamma }$	(28)
GRHC NRs subjected to UDL $(\mathit{P}=0)$
FC NR is subjected to UDL	$\overline{\mathrm{v}}\left(0\right)={\overline{\mathrm{v}}}^{\prime }\left(0\right);$ $\mathrm{M}(1)=\mathrm{Q}(1)=0$	${\mathrm{D}}_1=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3(1-e^{-\sqrt{\Gamma }}\sqrt{\Gamma })}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^2({\mathrm{e}}^{\sqrt{\Gamma }}+{\mathrm{e}}^{-\sqrt{\Gamma }})};$ ${\mathrm{D}}_2=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3(1+{\mathrm{e}}^{\sqrt{\Gamma }}\sqrt{\Gamma })}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^2({\mathrm{e}}^{\sqrt{\Gamma }}+{\mathrm{e}}^{-\sqrt{\Gamma }})};$ ${\mathrm{D}}_3=\frac{{\mathrm{q}}_0{\mathrm{L}}^3(2+{\mathrm{e}}^{\sqrt{\Gamma }}\sqrt{\Gamma }-{\mathrm{e}}^{-\sqrt{\Gamma }}\sqrt{\Gamma })}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^2({\mathrm{e}}^{\sqrt{\Gamma }}+{\mathrm{e}}^{-\sqrt{\Gamma }})};$ ${\mathrm{D}}_4=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\Gamma }$	(29)
FF NR is subjected to UDL	$\overline{\mathrm{v}}\left(0\right)={\overline{\mathrm{v}}}^{\prime }\left(0\right)=0;$ $\overline{\mathrm{v}}(1)={\overline{\mathrm{v}}}^{\prime }(1)=0$	${\mathrm{D}}_1=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}({\mathrm{e}}^{\sqrt{\Gamma }}-1)};$ ${\mathrm{D}}_2=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3{\mathrm{e}}^{\sqrt{\Gamma }}}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}({\mathrm{e}}^{\sqrt{\Gamma }}-1)};$ ${\mathrm{D}}_3=\frac{{\mathrm{q}}_0{\mathrm{L}}^3({\mathrm{e}}^{\sqrt{\Gamma }}+1)}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^{\frac{3}{2}}({\mathrm{e}}^{\sqrt{\Gamma }}-1)};$ ${\mathrm{D}}_4=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\Gamma }$	(30)
SS NR is subjected to UDL	$\overline{\mathrm{v}}\left(0\right)=\overline{\mathrm{v}}\left(1\right)=0;$ $\mathrm{Q}(0)=\mathrm{Q}(1)=0$	${\mathrm{D}}_1=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3{\mathrm{e}}^{\overline{x}\sqrt{\Gamma }}}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^2(1+{\mathrm{e}}^{\sqrt{\Gamma }})};$ ${\mathrm{D}}_2=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3{\mathrm{e}}^{-\overline{x}\sqrt{\Gamma }}}{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\Gamma }^2(1+{\mathrm{e}}^{-\sqrt{\Gamma }})};$ ${\mathrm{D}}_3=\frac{{\mathrm{q}}_0{\mathrm{L}}^3{\overline{x}}^2}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\Gamma };{\mathrm{D}}_4=-\frac{{\mathrm{q}}_0{\mathrm{L}}^3\overline{x}}{2{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\Gamma }$	(31)

Finite element (FE) formulation considering the surface and flexoelectric effects

The proposed FE formulation considers the thin NRs subjected to point load and UDL under static loading condition, by considering no axial force and ignoring inertial terms equation (22) is re-expressed as:

$${\left(\mathrm{BR}\right)}^{\mathrm{eff}}\frac{{\mathrm{d}}^4\mathrm{v}}{\mathrm{d}{\mathrm{x}}^4}-\left({\tau }_0{\mathrm{s}}^{*}\right)\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2}\mathrm{q}\left(\mathrm{x}\right)=0$$

(32)

In FE formulation, Galerkin’s weighted residual method is applied to equation (32) which can be reformulated as:

$$\Psi =\underset{0}{\overset{L}{{\displaystyle \int }}}\left\{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\frac{{\mathrm{d}}^4\mathrm{v}}{{\mathrm{dx}}^4}-({\tau }_0{\mathrm{s}}^{*})\frac{{\mathrm{d}}^2\mathrm{v}}{{\mathrm{dx}}^2}+\mathrm{q}(\mathrm{x})\right\}\overline{\mathrm{v}}\mathrm{dx}=0$$

(33)

L denotes the NR length and $\overline{\mathrm{v}}$ denotes the weight function.

$$\Psi =\underset{0}{\overset{L}{{\displaystyle \int }}}\left\{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\frac{{\mathrm{d}}^2\mathrm{v}}{{\mathrm{dx}}^2}\frac{{\mathrm{d}}^2\overline{\mathrm{v}}}{{\mathrm{dx}}^2}+\left({\tau }_0{\mathrm{s}}^{*}\right)\frac{\mathrm{dv}}{\mathrm{dx}}\frac{\mathrm{d}\overline{\mathrm{v}}}{\mathrm{dx}}+\mathrm{q}\overline{\mathrm{v}}\right\}\mathrm{dx}+\left(\mathrm{M}\frac{\mathrm{d}\overline{\mathrm{v}}}{\mathrm{dx}}-\mathrm{Q}\overline{\mathrm{v}}\right)|\begin{array}{c}\mathrm{L}\\ 0\end{array}=0$$

(34)

where the BM and shear force is formulated as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{M}=-{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\left(\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2}\right);\mathrm{Q}=-{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\left(\frac{{\mathrm{d}}^3\mathrm{v}}{\mathrm{d}{\mathrm{x}}^3}\right)}\\ \multicolumn{1}{c}{+\left({\tau }_0{\mathrm{s}}^{*}\right)\frac{\mathrm{dv}}{\mathrm{dx}}}\end{array}$$

Equation (33) again expressed in weak form with application of integration by parts as:

$$\begin{array}{c}\multicolumn{1}{c}{\Psi =\underset{0}{\overset{\mathrm{L}}{\int }}\{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2}\frac{{\mathrm{d}}^2\overline{\mathrm{v}}}{\mathrm{d}{\mathrm{x}}^2}+\left({\tau }_0{\mathrm{s}}^{*}\right)\frac{\mathrm{dv}}{\mathrm{dx}}\frac{\mathrm{d}\overline{\mathrm{v}}}{\mathrm{dx}}+\mathrm{q}\overline{\mathrm{v}}\}\mathrm{dx}}\\ -\{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\frac{{\mathrm{d}}^2\mathrm{v}}{\mathrm{d}{\mathrm{x}}^2}\frac{\mathrm{d}\overline{\mathrm{v}}}{\mathrm{dx}}-{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\left(\frac{{\mathrm{d}}^3\mathrm{v}}{\mathrm{d}{\mathrm{x}}^3}\right)\overline{\mathrm{v}}+\left({\tau }_0{\mathrm{s}}^{*}\right)\frac{\mathrm{dv}}{\mathrm{dx}}\overline{\mathrm{v}}\}|\begin{array}{c}\mathrm{L}\hfill \\ \multicolumn{1}{c}{0}\end{array}=0.\hfill \end{array}$$

(35)

From the above weak form in equation (35), it is noted that the highest order of derivative of $\mathrm{v}$ is 3, that is, the required approximation function can be obtained by differentiating the equation thrice. Therefore, we can say that function is a third/cubic degree polynomial. Moreover, the highest order of derivative inside the integral is 2 in equation (35), and therefore, the overall approximation should be C ¹ continuous (1 is a non-negative integer). Here, Figure 4 represents the two noded NR elements with consideration of surface parameters having translational, $\mathrm{v}$ and rotational, Ø degrees of freedom at each node. In addition, the cubic function also satisfies the conditions of displacement and slope continuity at nodes.

Figure 4.

Two noded NR element.

Hence, the nodal displacement vector is interpreted as:

$${\mathrm{v}}^{\mathrm{e}}={\left[\begin{array}{cccc}{\mathrm{v}}_1\hfill & {Ø}_1\hfill & {\mathrm{v}}_2\hfill & {Ø}_2\hfill \end{array}\right]}^{\prime }$$

(36a)

Therefore, making use of shape function $\mathrm{N}\left(\mathrm{x}\right)$ , the vertical displacement is interpolated as:

$$\mathrm{v}=\sum N_i{v^e}_i,i=1,2,3,4.$$

(36b)

Here, ${\mathrm{N}}_1$ , ${\mathrm{N}}_2$ , ${\mathrm{N}}_3$ , and ${\mathrm{N}}_4$ indicate the shape functions for rod elements. These cubic shape functions are known as Hermitian cubic interpolation (or cubic spline) functions and the elements formulated using these are called Hermite elements. In case of NR element, ${\mathrm{N}}_1=1$ when evaluated at node 1 and ${\mathrm{N}}_1=0$ when calculated at node 2. Because ${\mathrm{N}}_2$ is related to ${Ø}_1$ , we have $\frac{\mathrm{d}{\mathrm{N}}_2}{\mathrm{dx}}=1$ when evaluated at node 1. Likewise, for node 2, ${\mathrm{N}}_3$ and ${\mathrm{N}}_4$ have equal outcomes. $\mathrm{Ni}\left(\mathrm{x}\right)$ is the shape function which can be obtained as follows:

$$\begin{array}{c}{\mathrm{N}}_1\left(\mathrm{x}\right)=1-\frac{3{\mathrm{x}}^2}{{\mathrm{L}}^2}+\frac{2{\mathrm{x}}^3}{{\mathrm{L}}^3};{\mathrm{N}}_2\left(\mathrm{x}\right)=\mathrm{x}-\frac{2{\mathrm{x}}^2}{\mathrm{L}}+\frac{{\mathrm{x}}^3}{{\mathrm{L}}^2}\hfill \\ {\mathrm{N}}_3\left(\mathrm{x}\right)=\frac{3{\mathrm{x}}^2}{{\mathrm{L}}^2}-\frac{2{\mathrm{x}}^3}{{\mathrm{L}}^3};{\mathrm{N}}_4\left(\mathrm{x}\right)=-\frac{{\mathrm{x}}^2}{\mathrm{L}}+\frac{{\mathrm{x}}^3}{{\mathrm{L}}^2}\hfill \end{array}$$

(37)

Making use of equation (37) into equation (35), we can get

$$\underset{0}{\overset{\mathrm{L}}{\int }}\{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\mathrm{N}{″}^{\mathrm{T}}\mathrm{N}″+{\tau }_0{\mathrm{s}}^{*}\mathrm{N}{\prime }^{\mathrm{T}}\mathrm{N}\prime \}{\mathrm{v}}^{\mathrm{e}}\mathrm{dx}=-\underset{0}{\overset{\mathrm{L}}{\int }}\mathrm{q}\left(\mathrm{x}\right)\mathrm{Ndx}$$

(38)

The element stiffness matrix is given by:

$${\mathrm{K}}_{\mathrm{e}}=\underset{0}{\overset{\mathrm{L}}{\int }}\left\{\{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}\mathrm{N}{″}^{\mathrm{T}}\mathrm{N}″+{\tau }_0{\mathrm{s}}^{*}\mathrm{N}{\prime }^{\mathrm{T}}\mathrm{N}\prime \}\mathrm{dx}\right\}$$

(39a)

$${\mathrm{K}}_{\mathrm{e}}=\underset{0}{\overset{\mathrm{L}}{\int }}\{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}{\left(\frac{{\mathrm{d}}^2\mathrm{N}}{\mathrm{d}{\mathrm{x}}^2}\right)}^{\mathrm{T}}\frac{{\mathrm{d}}^2\mathrm{N}}{\mathrm{d}{\mathrm{x}}^2}+\left({\tau }_0{\mathrm{s}}^{*}\right){\left(\frac{\mathrm{dN}}{\mathrm{dx}}\right)}^{\mathrm{T}}\frac{\mathrm{dN}}{\mathrm{dx}}\}\mathrm{dx}$$

(39b)

$${\mathrm{K}}_{\mathrm{e}}=\frac{{\left(\mathrm{BR}\right)}^{\mathrm{eff}}}{{\mathrm{L}}^3}\left[\begin{array}{cccc}12& 6\mathrm{L}& -12& 6\mathrm{L}\\ 6\mathrm{L}& 4{\mathrm{L}}^2& -6\mathrm{L}& 2{\mathrm{L}}^2\\ -12& -6\mathrm{L}& 12& -6\mathrm{L}\\ 6\mathrm{L}& 12& -6\mathrm{L}& 4{\mathrm{L}}^2\end{array}\right]+\frac{\left({\tau }_0{\mathrm{s}}^{*}\right)}{30\mathrm{L}}\left[\begin{array}{cccc}36& 3\mathrm{L}& -36& 3\mathrm{L}\\ 3\mathrm{L}& 4{\mathrm{L}}^2& -3\mathrm{L}& -{\mathrm{L}}^2\\ -36& -3\mathrm{L}& 36& -3\mathrm{L}\\ 3\mathrm{L}& -{\mathrm{L}}^2& -3\mathrm{L}& 4{\mathrm{L}}^2\end{array}\right]$$

(40)

It is noted that equation (40) is divided into two sections: the initial section is analogous to improve the stiffness matrix while another section is interrelated to stiffness matrix accounting for the surface residual stress. The negative or positive magnitude of residual surface stresses influences the value of effective bending rigidity.

Therefore, the generalized nodal force vector is obtained as:

$${\mathrm{f}}_{\mathrm{e}}=-\underset{0}{\overset{\mathrm{L}}{\int }}\mathrm{q}\left(\mathrm{x}\right){\mathrm{N}}^{\mathrm{T}}\mathrm{dx}$$

(41)

Here, we considered the generalized nodal force vector which contains the effect of the point load P and UDL.

By using the assembly of nodal force and element stiffness matrix, the equilibrium equation with global force vectors and element stiffness can be obtained as:

$$\mathrm{Kd}=\mathrm{f}$$

(42)

where f, d, and K signify the global force, displacement vector, and global stiffness matrix, respectively.

Results and discussion

The effective properties of GRHC were evaluated using a three-phase MOM model with consideration of graphene volume fraction which is considered 0.2 times volume fraction of PZT-5H which is taken from Shingare and Naskar ⁵⁷ and are enlisted in Table 2. These properties of GRHC are used to envisage the electrostatic response of hybrid nanorods (NR). In this, we derived the analytical solutions for GRHC NRs using EB theory considering surface and flexoelectric effect in conjunction with FE model based on Galerkin residual method to validate the results obtained from the analytical model.

Table 2.

Effective properties of GRHC.

Material	${\mathbf{C}}_{11}$ (GPa)	${\mathbf{e}}_{31}$ (C/m2)	${ϵ}_{33}\times {10}^{-9}$ (F/m)	Volume fraction
GRHC	112.43	-6.9337	3.264	${\mathrm{v}}_{\mathrm{g}}=0.2{\mathrm{v}}_{\mathrm{p}}$

The surface elastic coefficients for the current model are equal to the elastic coefficients of GRHC multiplied by its surface layer thickness, assumed as 1 nm.^29,31,70,71 Similarly, the coefficient of surface piezoelectricity can be calculated. We have taken the flexoelectric coefficient as ${\mu }_{31}\approx {10}^{-8}\mathrm{C}{\mathrm{m}}^{-1}$ for further calculations. It is significant to note in NEMS applications that the piezoelectric NRs are basic building blocks in the design of nanostructures while analyzing the mechanical performances with high potential. We considered nanorod having length, $(\mathrm{L}=20\mathrm{H})$ with a thickness $(\mathrm{H})$ , diameter $(\mathrm{D})$ , and width $(\mathrm{B}=0.5\mathrm{H})$ . The results are evaluated for NRs subjected to a point load (P) and UDL $({\mathrm{q}}_0)$ for a cantilever (FC), fixed-fixed (FF), and simply supported (SS) BCs using both the theoretical and numerical modeling. According to the practical point of view, these loading and BCs are extremely important in automotive, aircraft, mechanical, and spacecraft engineering structures. These structures are mainly subjected to penetration of external load (i.e. resemblance to point load) and UDL with all mentioned BCs. In addition to this, two types of the cross-section of NRs are considered such as rectangular and circular cross-sections. In case of FE analysis, a convergence study is also carried out, and it is observed that the obtained results converge for a reasonable number of elements and agree well with the analytical results for static bending deflections of NRs. The results of the convergence study for normalized deflection of FC nanorod subjected to endpoint load are summarized in Table 3. From this, it is found that the converged solutions for static bending deflections of cantilever NRs can be achieved with more than three finite elements while it is achieved with five elements for the remaining two BCs. Therefore, each FE simulation is performed using five finite elements.

Table 3.

Convergence study for normalized deflection of NRs with respect to the number of finite elements.

Number of finite elements	1	2	3	4	5
FE results	−7.9215	−7.9870	−8.0091	−8.0157	−8.0157

Effect of graphene nanofillers

Figure 5 shows the variation of normalized deflection over the length of hybrid FC GRHC NR when subjected to UDL considering the different combinations of graphene and PZT. The different combinations of graphene and PZT considered here are: ${\mathrm{v}}_{\mathrm{g}}=0.2{\mathrm{v}}_{\mathrm{p}},0.4{\mathrm{v}}_{\mathrm{p}},0.6{\mathrm{v}}_{\mathrm{p}},0.8{\mathrm{v}}_{\mathrm{p}}$ and ${\mathrm{v}}_{\mathrm{g}}={\mathrm{v}}_{\mathrm{p}}$ , whereas ${\mathrm{v}}_{\mathrm{g}}\mathrm{and}{\mathrm{v}}_{\mathrm{p}}$ denote the volume fraction of graphene and PZT, respectively. From this Figure 5, it can be noticed that the normalized deflection of NR decreases as the value of graphene content in hybrid composite NR increases. It is due to fact that the effective properties of GRHC were determined with consideration of strong covalent bond which offers interaction and in-plane stability of graphene as well as strong van der Wall forces formed between graphene and matrix. It is also reported that the large surface area of graphene aids in strong interaction with fiber and matrix which results in improved effective properties and obviously exhibits increased stiffness of GRHC nanorod. According to the experimental studies, it is clear that the appropriate graphene volume fraction in composite should be less than 10%. Therefore, in this study, the volume fraction of graphene is considered 0.2 times the volume fraction of PZT (i.e. ∼8% in composite). In the view of the experimental outlook, we have considered ${\mathrm{v}}_{\mathrm{g}}=0.2{\mathrm{v}}_{\mathrm{p}}$ for further calculation of static bending deflection which consists of graphene content less than 10% in hybrid GRHC NR.

Figure 5.

Normalized deflection of FC NR under UDL with different graphene percentages.

Effect of flexoelectricity and surface parameters on GRHC nanorods

In this Section, the flexoelectric and surface stress effects on the NRs with different BCs are discussed to investigate their elastic behavior (i.e. softening and stiffening). The deviation of normalized bending rigidity ${\left(\mathrm{BR}\right)}^{\mathrm{eff}}/{\left(\mathrm{BR}\right)}^0$ and ${\left(\mathrm{BR}\right)}^{\mathrm{s}}/{\left(\mathrm{BR}\right)}^0$ against rod thickness/diameter is demonstrated in Figure 6, whereas ${\left(\mathrm{BR}\right)}^{\mathrm{eff}},{\left(\mathrm{BR}\right)}^{\mathrm{s}}$ , and ${\left(\mathrm{BR}\right)}^0$ denote the bending rigidity of NRs accounting for the flexoelectricity, surface parameters, and without considering the effect of flexoelectric and surface parameters, respectively. According to equation (23), bending rigidity depends only on rod diameter/thickness and not on the in-plane dimensions of NRs. It is observed that the flexoelectricity effect on the normalized bending rigidity is size-dependent and it is more eminent for NRs with a smaller thickness/diameter. For example, the bending rigidity considering the effect of both flexoelectricity and surface stress and considering only flexoelectricity is 3.2 and 2.2 times the bending rigidity without consideration of surface parameters and flexoelectricity when the NR thickness is kept 10 nm while 1.89 and 1.05 times greater when the NR diameter is kept 10 nm. Such a difference is obvious and for envisaging the electrostatic response of nanostructures accurately, it cannot be ignored. As the thickness of NR increases, the flexoelectricity effect reduces promptly and as the normalized bending rigidity reaches unity, the effect of flexoelectricity again disappears.

Figure 6.

Normalized bending rigidity with respect to: (a) NR thickness and (b) NR diameter.

Figures 7 to 18 depict that the flexoelectricity and surface effect show a significant role in elastic response of NRs for different BCs and shows the clear difference during the time of comparison with the results of NR which considers pure flexoelectricity, combined surface and flexoelectricity effects, and without consideration of surface and flexoelectricity effects (classical or conventional NR). In case of SS and FF NRs, due to symmetry, only half of NRs deflections are shown. From this, it is observed that the deflection of NRs with consideration of both flexoelectricity and piezoelectricity is less when compared with results obtained considering only piezoelectricity. Besides the flexoelectricity effect, the surface stress effect is considered to examine the elastic response of deflection of NRs. In case of cantilever (FC) NR, this elastic behavior depends on the sign of surface stresses (positive ${\tau }_0>0$ and negative ${\tau }_0<0$ ). In contrast to FC NR, the deflections of SS and FF NRs are stiffened by the positive surface stresses. Simultaneously, it is softened by the negative surface stresses and found in good coherence with earlier conveyed results. ²⁵

Figure 7.

Normalized deflection over the length of SS rectangular NR.

Figure 8.

Normalized deflection over the length of SS circular NR.

Figure 9.

Normalized deflection over the length of FC rectangular NR.

Figure 10.

Normalized deflection over the length of FC circular NR.

Figure 11.

Normalized deflection over the length of FF rectangular NR.

Figure 12.

Normalized deflection over the length of FF circular NR.

Figure 13.

Normalized deflection over the length of SS rectangular NR.

Figure 14.

Normalized deflection over the length of SS circular NR.

Figure 15.

Normalized deflection over the length of FC rectangular NR.

Figure 16.

Normalized deflection over the length of FC circular NR.

Figure 17.

Normalized deflection over the length of FF rectangular NR.

Figure 18.

Normalized deflection over the length of FF circular NR.

It can be seen that SS and FF NRs show a stiffer elastic behavior while FC NR shows a softer elastic response for ${\tau }_0>0$ while the converse is correct for ${\tau }_0<0$ as compared to the corresponding NRs without considering the surface stress and flexoelectric effect (conventional NR). This may be due to fact that the different curvatures of FC, SS, and FF NRs under the point load and UDL. Hence, the deflection of FC NR exhibits the downward curvature (negative curvature) while the SS NR exhibits upward curvature (positive curvature) and FF NR exhibits both downward as well as upward curvature for the same load. The main phenomenon behind the softer elastic behavior of NR is that the value of imposed point/UDL is identical to the negative uniformly transverse load in the same direction because of residual surface stress in cantilevered NR; else, the enhancement in deflection of NR which is mechanically deformed opposed by this positive uniformly transverse load in case of SS and FF NRs. The effect of surface stress and flexoelectricity plays a significant role in the stiffening and softening response of NRs with different BCs under point and UDL loading conditions. The sign of residual surface stress becomes important in the stiffening and softening response of NRs whereas the flexoelectric effect often stiffens means reduces the deflection of the nanostructure. In case of FC NR, surface stress softens the elastic response though the flexoelectric effect helps to overcome this softness. In contrast to FC NR, the elastic behavior of SS and FF NRs stiffens by these surface stresses.

In FE calculations, these NRs are subjected to the same point load and UDL is defined in a downward direction. Also, the material properties and geometry used in the FE calculations are similar to those of the analytical model. If the load is applied downward then the FC NR bent with concave downward shape results in negative curvature. Therefore, the additional corresponding uniform transverse load will be applied by surface stresses effect and the opposite signs of mechanical load. While the SS and FF NRs are exposed to downward load, it bent with a concave upward shape and the additional uniform transverse loads improve the deflection which is mechanically deformed. As the sign of curvature of NRs is similar to the surface stress, the deflection is enriched, and the reverse is true when the sign of curvature of NRs is opposite to the surface stress. The FF NR exhibits stiffer behavior when compared to SS NR because FF NR exhibits both downward as well as upward curvature. This is attributed to the reason that the SS and FF NRs are stiffened by positive surface stress effect $({\tau }_0>0)$ while the FC NR softens by the same effect. Moreover, the results evaluated by the FE simulations are in great accordance with those of the analytical models. From the above results, in case of circular NRs, it is also obvious that the combined effect of flexoelectricity and surface shows the significant enhancement in the reduction of deflection of NRs.

It is evident from Figures 7 to 18 that due to the incorporation of graphene as nanofillers, surface stress, and flexoelectricity result in a reduction of deflection of GRHC NRs as that of deflection of conventional NRs that is, without surface and flexoelectric effects (classical NRs). The other most important reason behind the reduction in static bending deflection is nothing but enriched stiffness/bending rigidity of GRHC due to consideration of graphene nanofillers which exhibit excellent electro-thermo-mechanical properties. Despite the fact, the magnitude of maximum deflection of GRHC NRs under point load and UDL observed in the following order: FC > SS > FF while the deflection of NRs irrespective of edge support conditions and cross-section gives maximum deflection for point load as compared to UDL. For instance, the deflection of GRHC NRs with consideration of pure flexoelectricity and considering combined flexoelectricity and surface stress effects are enhanced significantly when compared to that of deflection of conventional NRs irrespective of all edge support conditions. But if we compare the results between pure flexoelectricity and combine effects, then the deflection of NRs is reduced significantly due to consideration of later case. From these Figures 5 to 18, it is noticed that the surface parameters and flexoelectricity have a remarkable influence for smaller thickness and diameter of NRs and hence, it must be accounted properly. Therefore, it is concluded that the results presented in this work by using the analytical and FE models are found in excellent coherence.

Such numerical consequences fundamentally open up the paths of potential exploitation and enrichment of the desired electromechanical responses in design engineering including the factors such as open- and short-circuit condition, strain/electric gradient, surface effects, electromechanical loading as well as inverse piezo- and flexo-electric effects. With the recent advances in nano-scale manufacturing and experimental capabilities, this article will offer the essential physical understandings in modeling the size-dependent electromechanical coupling in multifunctional materials, systems, and devices for applications in distributed sensors, actuators, active controllers, and energy harvesters.

Conclusions and perspective

This article explores the static bending deflection behavior of graphene-reinforced hybrid composite (GRHC) nanorods (NRs) with flexoelectric and surface effects. In this, we derived the closed-form solutions for GRHC NRs using EB theory and linear piezoelectricity theory by considering surface stress and the flexoelectric effect. Moreover, the theoretical FE model is derived using Galerkin residual method to validate the results obtained from the analytical model under the same loading and BCs. Based on this, the static bending deflection of GRHC NRs for different types of BCs is considered to investigate the role of flexoelectricity and surface effect. For instance, we considered cantilever (FC), fixed-fixed (FF), and simply supported (SS) nanorods subjected to point and UDL loadings. The reduction in static deflection of GRHC NRs is enhanced by considering flexoelectricity and surface effect over the deflection of GRHC conventional NRs. From these results, it can be concluded that the magnitude of maximum deflection of GRHC NRs under point load and UDL observed in the following order: FC > SS > FF while the deflection of NRs irrespective of edge support conditions and cross-section gives maximum deflection for point load as compared to UDL. For instance, the deflection of GRHC NRs with consideration of pure flexoelectricity and considering combined flexoelectricity and surface stress effects are enhanced significantly when compared to that of deflection of conventional NRs irrespective of all edge support conditions. This current study offers pathways for developing new proficient novel GRHC materials with enhanced control authority and offers a guideline for the design of nanodevices in NEMS applications and several other industries.

The flexoelectricity is found to be more dominant for thin structures and it cannot be ignored while modeling 1-D, 2-D, and 3-D composite nanostructures. The novelty of the current study considers the introduction of graphene nanofillers in evaluating the overall properties of hybrid piezoelectric composite using the representative volume element (RVE) technique in context to composite laminates to form the different structural elements (MEMS/NEMS). The overall conclusion of this work is that the developed analytical and numerical models herein may provide the theoretical base for investigating the electrical and mechanical response of energy harvesters in form of beam and wire/rod. In the future, the proposed computationally effective method for electromechanical analysis of composite structures can be exploited for numerous other materials as well as structural systems to study their global response.