On Bending of Piezoelectrically Layered Perforated Nanobeams Embedded in an Elastic Foundation with Flexoelectricity

Abstract

Analysis of the electromechanical-size-dependent bending of piezoelectric composite structural components with flexoelectricity has been considered by many researchers because of the developments of nanotechnology and the applicability of piezoelectric composite nanobeam structures in Micro/Nano-Electro-Mechanical Systems (MEMS/NEMS). Therefore, the work investigates the size-dependent electromechanical bending of piezoelectrically layered perforated nanobeams resting on elastic foundations including the flexoelectric effect. Within the framework of the modified nonlocal strain gradient elasticity theory, both the microstructure and nonlocality effects are captured. The governing equilibrium equations including piezoelectric and flexoelectric effects are derived using Hamilton’s principle. Closed forms for the non-classical electromechanical bending profiles are derived. The accuracy of the proposed methodology is verified by comparing the obtained results with the available corresponding results in the literature within a 0.3% maximum deviation. Parametric studies are conducted to explore effects of perforation parameters, elastic foundation parameters, geometric dimensions, nonclassical parameters, flexoelectric parameters, as well as the piezoelectric parameters on the bending behavior of piezoelectrically layered perforated nanobeams. The obtained results demonstrate that incorporation of the nondimensional elastic foundation parameters, K_p = 2 and K_w = 20, results in a reduction in the relative percentage reduction in the maximum nondimensional mechanical transverse deflection due to increasing the perforation filling ratio from 0.2 to 1 from 199.86% to 91.83% for a point load and 89.39% for a uniformly distributed load. On the other hand, with K_p = 5 and K_w = 50, the relative percentage difference of the electromechanical bending deflection due to increasing the piezoelectric coefficient, e₃₁₁, reaches about 8.7% for a point load and 8.5% for a uniformly distributed load at a beam aspect ratio of 50. Thus, the electromechanical as well as mechanical behaviors could be improved by controlling these parameters. The proposed methodology and the obtained results are supportive in many industrial and engineering applications, i.e., MEMS/NEMS.

1. Introduction

Nowadays, from both the fundamental and application sides, a new mechanical-electrocoupling effect (i.e., direct flexoelectric effect) between the strain gradient and electric polarization, or vice versa, has attracted a fair amount of attention [1,2]. The flexoelectric effect of materials usually intensifies at the microscale even under small deformation owing to the intrinsic scaling of the strain gradient. Due to such a feature, flexoelectricity has attracted much attention in the past decades in a variety of materials such as crystalline ceramics and ferroelectric thin films to develop their electromechanical model [3].

Baroudi et al. [4] developed an analytical model to study static and dynamic responses of piezoelectric flexoelectric strain gradient nanobeams with bidirectional polarization and a self-field effect. Sedighi [5] represented the influence of the nonlinear amplitude of vibration on the dynamic pull-in instability and fundamental frequency of actuated microbeams. Kheibari and Beni [6] studied the free vibration of piezoelectric nanotubes by using Love’s cylindrical thin-shell model. Kaghazian et al. [7] studied the free vibration of a piezoelectric nonlocal Euler–Bernoulli nanobeam. Moory-Shirbani et al. [8] investigated experimentally and numerically a piezoelectrically actuated multilayered imperfect microbeam subjected to an applied electric potential. Chen et al. [9] exploited modified strain gradient theory in analyzing the dynamic responses of a flexoelectricity functionally graded (FG) porous piezoelectric sandwich Euler–Bernoulli nanobeam. Gupta et al. [10] presented the surface and flexoelectric effects on the response of the piezoelectric and dielectric properties of a Boron-Nitride-reinforced nanocomposite. Van Minh and Van Ke [11] developed comprehensive studies on mechanical responses of piezoelectric nanoplate foundations with the flexoelectric effect. Qian and Wang [12] examined the electromechanical coupling bandgap of a locally resonant piezoelectric/elastic phononic crystal double-layer nonlocal nanobeam. Zhao et al. [13] exploited the differential quadrature method to present the effect of porosity and flexoelectricity on the static bending and free vibration of axially FG piezoelectric nanobeams. Malikan and Eremeyev [14] studied the impact of flexoelectricity on a piezoelectric nanobeam involving internal viscoelasticity. Malikan and Eremeyev [15] evaluated the nonlinear bending of a piezo-flexomagnetic strain gradient nanobeam based on an analytical-numerical solution. Monaco et al. [16] developed a trigonometric solution for the bending response of magneto-electro-elastic nanoplates in a hygrothermal environment. Naderi et al. [17] investigated the size-dependent piezoelectric load on the vibration, buckling, and energy harvesting response of local and nonlocal piezoelectric nanobeams. Hao-nan et al. [18] examined the vibration response of rotating FG piezoelectric nonlocal nanobeams by using the differential quadrature method. Ren et al. [19] developed an integral nonlocal gradient piezoelectric model to investigate the bending of an FG piezoelectric nanobeam. Li and Jin [20] explored the electromechanical vibration and induced stresses of a piezoelectric nanobeam with a symmetrical FGM core under low-velocity impact. Ansari et al. [21] developed a numerical study for free vibration analysis of piezoelectric Bernoulli–Euler nonlocal nanobeams considering flexoelectric effects. Alam and Mishra [22] exploited a boundary layer solution to study the post-critical thermo-electromechanical stability of nonlocal strain gradient FG piezoelectric cylindrical shells. Boyina and Piska [23] studied the wave propagation in viscoelastic Timoshenko nanobeams under surface and magnetic field effects.

In the frame of nonlocal strain gradient theory, Lu et al. [24] obtained analytical solutions for natural frequencies of a size-dependent sinusoidal shear deformation nanobeam. She et at. [25] predicted the wave propagation of FG porous nanobeams based on Reddy’s higher-order shear deformation theory. Hashemian et al. [26] studied bending and buckling responses of Timoshenko and higher-order nanobeams. Karami et al. [27] investigated the vibration response of a 2D tapered porous FG Timoshenko beam including temperature and porosity influences on the material properties. Ducottet and El Baroudi [28] derived a mathematical model and frequency equation to study the radial vibration of nanospheres with nonlocal strain gradient theory. Hai et al. [29] studied the vibrational behavior of a sandwich honeycomb piezoelectric actuators microplate resting on the Pasternak elastic foundation. Vahidi-Moghaddam et al. [30] used the Galerkin approach and perturbation technique to examine the nonlinear forced vibrations of nonlocal strain gradient microbeams. Moradi-Dastjerdi and Behdinan [31] evaluated the stress waves in a lightweight substrate plate actuated with piezoelectric layers under sinusoidal pressures. Eghbali et al. [32] investigated vibrations of a cracked nonlocal strain gradient piezoelectric Euler–Bernoulli nanobeam to the effects of Euler–Bernoulli beam theory. Yang et al. [33] investigated the vibration response of unified nanobeams based on two-phase local/nonlocal strain and stress gradient theory, as well as surface elasticity theory. Gui and Wu [34] examined the buckling stability of thermo-magneto-electro-elastic cylindrical nanoshells surrounded by an elastic medium. Zhang et al. [35] exploited Mathieu–Hill equations with the aid of the Navier solution to study the dynamic stability of magneto-electro-thermo-elastic cylindrical nanoshells resting on the Winkler–Pasternak elastic foundations.

Based on the listed review and the authors’ knowledge, investigation of the electromechanical size-dependent bending of piezoelectrically layered perforated nanobeams embedded in a two-parameter elastic foundation with flexoelectricity based on the nonlocal strain gradient theory has not been handled elsewhere. The present work fills this gap. The rest of the paper has been organized as follows: Section 2 presents the mathematical model including equivalent geometrical parameters, displacement relations, basic constitutive equations, electromechanical relations, and modified nonlocal strain gradient theory with flexoelectricity, the energy formulation, and the governing electromechanical equilibrium equations. The analytical solution methodology and verification are proved in Section 3. Verification of the developed procedure and main effects of perforation parameters, elastic foundation parameters, geometric dimensions, nonclassical parameters, flexoelectric parameters, and piezoelectric parameters on the bending behavior of piezoelectrically layered perforated nanobeams is shown and discussed in Section 4. Conclusions and main points are highlighted in Section 5. Future work recommendations are offered in Section 6.

2. Mathematical Model

A piezoelectric layered nanobeam with a perforated core resting on an elastic foundation is presented in Figure 1. The beam is composed of two similar piezoelectric face sheets with the same height h_p and material characteristics with elasticity modulus E_p. On the other hand, the perforated core is made of a homogenous linear elastic material with Young’s modulus E_c and height h_c. All layers have the same beam length L and width W_b. The overall beam thickness is h, and h = h_c + 2h_p. The polarization direction of both piezoelectric face sheets is assumed upward. The two face sheets are exposed to two opposite-direction electric fields. The x-axis is directed through the nanobeam length direction, while the z-axis is through the thickness direction.

2.1. Equivalent Geometrical Parameters

As shown in Figure 1, the number of holes throughout the beam cross-section is referred to as N. Thus, the filling ratio of the beam can be expressed as [36,37,38,39]

$$\alpha =\frac{t_s}{l_s}, 0\le \alpha \le 1, \alpha =\left\{\begin{array}{c}1 \mathrm{Fully\; filled} \left(\mathrm{solid\; beam}\right)\\ 0 \mathrm{Completely\; perforated} \left(\mathrm{Artifitial\; case}\right)\end{array}\right.$$

(1)

where t_s and l_s are the length between two holes and the distance between centers of two holes, respectively. The equivalent bending and shear stiffness ratio are given by [40,41,42]

$${\left(EI\right)}_{eq}={\left(EI\right)}_s\left\{\frac{\alpha \left(N+1\right)\left(N^2+2N+{\alpha }^2\right)}{\left(1-{\alpha }^2+{\alpha }^3\right)N^3+3\alpha N^2+\left(3+2\alpha -3{\alpha }^2+{\alpha }^3\right){\alpha }^{2}N+{\alpha }^3}\right\}$$

(2)

$${\left(GA\right)}_{eq}={\left(EA\right)}_s\left(\frac{\left(1+N\right){\alpha }^3}{2N}\right)$$

(3)

where (EI)_eq and (EI)_s are, respectively, the equivalent bending stiffness of perforated and solid beams, and (GA)_eq and (GA)_s are the equivalent shear stiffness of perforated and solid beams, respectively.

2.2. Displacement and Basic Mechanical Relations

Considering the geometry of the composite beam shown in Figure 1, the displacement field according to Euler Bernoulli beam theory (EBBT) is given by [43]

$$\left\{\begin{array}{c}{u}_x\left(x,z,t\right)\\ u_z\left(x,z,t\right)\end{array}\right\}=\left\{\begin{array}{c}-z\frac{\partial w\left(x,t\right)}{\partial x}\\ w\left(x,t\right)\end{array}\right\},\hspace{1em}{u}_y\left(x,z,t\right)=0$$

(4)

where u_x(x, y, z), u_y(x, y, z), and u_z(x, y, z) are the displacement functions in x, y, and z directions, respectively. $w\left(x,t\right)$ refers to the transverse displacement.

The only nonzero strain component based on EBBT is expressed as [44,45]

$${\epsilon }_{xx}\left(x,z,t\right)=\frac{\partial u_x\left(x,z,t\right)}{\partial x}=-z\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}$$

(5)

Considering an isotropic linear elastic material constitutive behavior for the perforated core, the constitutive relation based on the EBBT can be expressed as [46]

$$\left\{\begin{array}{c}{\sigma }_{xx}\left(x,z,t\right)\\ {\sigma }_{yy\left(x,z,t\right)}\\ {\sigma }_{zz}\left(x,z,t\right)\end{array}\right\}=\left\{\begin{array}{c}\widehat{E}{\epsilon }_{xx}\left(x,z,t\right)\\ \lambda {\epsilon }_{xx}\left(x,z,t\right)\\ \lambda {\epsilon }_{xx}\left(x,z,t\right)\end{array}\right\}=\left\{\begin{array}{c}\frac{\left(1-\nu \right)E}{\left(1+\nu \right)\left(1-2\nu \right)}\left(-z\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)\\ \lambda \left(-z\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)\\ \left(\frac{\nu }{1-\nu }\right){\sigma }_{xx}\left(x,z,t\right)\end{array}\right\}$$

(6a)

with $\widehat{E}$ is the equivalent modulus of elasticity, $\widehat{E}=2\mu +\lambda$; E is Young’s modulus; v is Poisson’s ratio; σ_xx (x, y, z), σ_yy (x, y, z), and σ_zz (x, y, z) are the components of the Cauchy normal stress tensor. λ and μ are Lame’s constants in classical elasticity, which are related to the elasticity modulus (E) and Poisson’s ratio (v) as

$$\mu =\frac{E}{2\left(1+\nu \right)}, \lambda =\frac{\nu E}{\left(1+\nu \right)\left(1-2\nu \right)}$$

(6b)

2.3. Electric Field and the Electromechanical Relations

Incorporating the piezoelectric effect, the electric field is given by [47,48]

$$E_z\left(x,z,t\right)=-\frac{\partial {\varphi }_z\left(z,t\right)}{\partial z}$$

(7)

where ${\varphi }_z\left(z,t\right)$ refers to the electric potential field. The electric enthalpy energy density function is expressed as [49,50]

$$H=-\frac{1}{2}{a}_{kl}{E}_k{E}_l+\frac{1}{2}{c}_{ijkl}{\epsilon }_{ij}{\epsilon }_{kl}-e_{ijk} E_k {\epsilon }_{ij}-u_{ijkl} E_i {\epsilon }_{jk,l}$$

(8)

where $a_{kl}, c_{ijkl}, e_{ijk}, \mathrm{and} u_{ijkl}$ are, respectively, the second-order permittivity tensor, the fourth-order elasticity tensor, the piezoelectric coefficient tensor, and the electric field strain gradient coupling coefficient, which denotes the higher-order electromechanical coupling induced by strain gradients. In terms of the electric enthalpy energy density, the electric displacement is expressed as [48,51]

$$D_i=-\frac{\partial H}{\partial E_i}=a_{ij}{E}_j+e_{ijk} {\epsilon }_{jk}+u_{ijkl} {\epsilon }_{jk,l}$$

(9)

Based on the Gaussian theorem, the following condition is verified:

$$\frac{\partial D_z}{\partial z}=0\Rightarrow a_{33}\frac{\partial E_z\left(x,z,t\right)}{\partial z}+e_{311} \frac{\partial {\epsilon }_{xx}\left(x,z,t\right)}{\partial z}+u_{3111}\frac{ {\partial }^2{\epsilon }_{xx}\left(x,z,t\right)}{\partial x\partial z}+u_{3113}\frac{ {\partial }^2{\epsilon }_{xx}\left(x,z,t\right)}{\partial z^2}=0$$

(10)

Substituting Equation (7) into Equation (10) yields

$$a_{33}\frac{\partial E_z\left(x,z,t\right)}{\partial z}-e_{311} \left(\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)-u_{3111}\left(\frac{{\partial }^{3}w\left(x,t\right)}{\partial x^3}\right)=0$$

(11a)

The partial derivative of the electric field is expressed as

$$\frac{\partial E_z\left(x,z,t\right)}{\partial z}=\frac{1}{a_{33}}\left[e_{311} \left(\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)+u_{3111}\left(\frac{{\partial }^{3}w\left(x,t\right)}{\partial x^3}\right)\right]$$

(11b)

The electric field can be obtained by integrating Equation (11) with respect to z as

$$E_z\left(x,z,t\right)-E_{z0}=\frac{z}{a_{33}}\left[e_{311} \left(\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)+u_{3111}\left(\frac{{\partial }^{3}w\left(x,t\right)}{\partial x^3}\right)\right]$$

(12)

Equation (12) can be rewritten as

$$E_z\left(x,z,t\right)=E_{z0}+\frac{z}{a_{33}}\left[e_{311} \left(\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)+u_{3111}\left(\frac{{\partial }^{3}w\left(x,t\right)}{\partial x^3}\right)\right]=E_{z0}-\left[\left(\frac{e_{311}}{a_{33}}\right) {\epsilon }_{xx}\left(x,z,t\right)+\left(\frac{u_{3111}}{a_{33}}\right)\left(\frac{\partial {\epsilon }_{xx}\left(x,z,t\right)}{\partial x}\right)\right]$$

(13)

in which $E_{z0}$ denotes the initial electric field in the z-direction.

2.4. The Modified Nonlocal Strain Gradient Theory with Flexoelectricity Effect

To incorporate the nonclassical nonlocal as well as the strain gradient effects, the modified nonlocal strain gradient theory is adopted. The “nonlocal strain gradient theory” is a consistent model to describe the small-scale mechanical response of nanostructures. Based on the nonlocal strain gradient theory (NSGT), the stress field ${\sigma }_{ij}^t$ accounts for not only the nonlocal elastic stress field but also the strain gradient stress field [52]. Introducing the flexoelectricity effect, the following nonclassical constitutive relations for the components of the stress, the higher gradient strains, and the nonlocal electric potential based on the EBBT are given as:

$$\left(1-{\left(e_{0}a\right)}^2\frac{{\partial }^2}{\partial x^2}\right)\left\{\begin{array}{c}{\sigma }_{xx}^t\\ {\sigma }_{113}\\ {\sigma }_{111}\\ D_z\end{array}\right\}=\left\{\begin{array}{c}\left\{E\left(1-l^2\frac{{\partial }^2}{\partial x^2}\right)+\frac{e_{311}^2}{a_{33}}\right\}{\epsilon }_{xx}+\frac{e_{311}}{a_{33}}{\mu }_{3111}\frac{\partial {\epsilon }_{xx}}{\partial x}-e_{311}{E}_{z0}\\ \frac{{\mu }_{3113}}{a_{33}}\left(e_{311}{\epsilon }_{xx}+{\mu }_{3111}\frac{\partial {\epsilon }_{xx}}{\partial x}\right)-{\mu }_{3113}{E}_{z0}\\ \frac{{\mu }_{3111}}{a_{33}}\left(e_{311}{\epsilon }_{xx}+{\mu }_{3111}\frac{\partial {\epsilon }_{xx}}{\partial x}\right)-{\mu }_{3111}{E}_{z0}\\ a_{33}{E}_{z0}+{\mu }_{3113}\frac{\partial {\epsilon }_{xx}}{\partial z}\end{array}\right\}$$

(14)

However, while studying higher-order metaphysic effects such as flexoelectric or flexomagnetic properties, there is a problem. The right-side term in NSGT is the strain gradient, obviously. On the other hand, the term “σ₁₁₁ ε_xx,x” is also a strain gradient located in the energy relation. The problem is that these two terms act the same and both are the strain gradient. This occurs for flexoelectric and flexomagnetic effects when using NSGT. This manner is called the “double effect” for the strain gradient [53].

2.5. Energy Formulation and the Governing Equation

Based on the EBBT, the strain energy (${\pi }_s$) is expressed as [52]

$${\pi }_s={{\displaystyle \int }}_{\mathsf{\Omega }}Hd\mathsf{\Omega }=\frac{1}{2}{{\displaystyle \int }}_{\mathsf{\Omega }}\left({\sigma }_{xx}^t{\epsilon }_{xx}-D_z{E}_z+{\sigma }_{111}{\epsilon }_{xx,x}+{\sigma }_{113}{\epsilon }_{xx,z}\right)d\mathsf{\Omega }=\frac{1}{2}{{\displaystyle \int }}_0^L\left(\overline{M}{\epsilon }_{xx0}+{\overline{M}}^{\prime }{\epsilon }_{xx0,x}\right)dx$$

(15)

with

$$\begin{array}{ccc}[{\epsilon }_{xx0}& \overline{M}& {\overline{M}}^{\prime }]\end{array}=\begin{array}{ccc}[-\frac{{\partial }^{2}w}{\partial x^2}& M+\frac{e_{311}}{a_{33}}{M}_D+N_{113 }& \frac{{\mu }_{3111}}{a_{33}}{M}_D+N_{111}]\end{array}$$

(16)

Then, the strain energy can be evaluated by

$${\pi }_s=\frac{1}{2}{{\displaystyle \int }}_0^L\left[\left(M+\frac{e_{311}}{a_{33}}{M}_D+N_{113}\right)\left(-\frac{{\partial }^{2}w}{\partial x^2}\right)+\left(\frac{{\mu }_{3111}}{a_{33}}{M}_D+N_{111}\right)\left(-\frac{{\partial }^{3}w}{\partial x^3}\right)\right]d\mathrm{x}$$

(17)

Recalling the modified nonlocal strain gradient theory, the following constitutive relation can be written.

$$\left(1-{\left(e_{0}a\right)}^2\frac{{\partial }^2}{\partial x^2}\right)\left\{\begin{array}{c}M\\ N_{113}\\ M_{111}\\ M_D\end{array}\right\}=\left\{\begin{array}{c}\left(1-l^2\frac{{\partial }^2}{\partial x^2}\right)\left\{{\left(E_c{I}_c\right)}_{eq}-E_p{I}_p\right\}{\epsilon }_{xx0}+\frac{e_{311}^2}{a_{33}}{I}_p{\epsilon }_{xx}+\frac{e_{311}{\mu }_{3111}{I}_p}{a_{33}}\frac{\partial {\epsilon }_{xx0}}{\partial x}-e_{311}{I}_1{E}_{z0}\\ 0\\ w_b\left\{\frac{{\mu }_{3111}}{a_{33}}{I}_p\left(e_{311}{\epsilon }_{xx0}+{\mu }_{3111}\frac{\partial {\epsilon }_{xx0}}{\partial x}\right)-{\mu }_{3111}{I}_1{E}_{z0}\right\}\\ w_b{a}_{33}{I}_1{E}_{z0}\end{array}\right\}$$

(18)

where

$$\left[\begin{array}{ccccc}{I}_c& I_p& I_1& A_c& A_p\end{array}\right]=\left[\begin{array}{ccccc}\frac{h_c^3}{12}& \left(\frac{h^3}{12}-I_c\right)& \left(\frac{h^2}{4}-\frac{h_c^2}{4}\right)& h_c{w}_b& h_p{w}_b\end{array}\right]$$

(19)

The work done by the external load and the elastic foundation forces is expressed as

$$W_{ex}=-w_b{{\displaystyle \int }}_0^L\left\{q\left(x,t\right)D_s\left(x-x_p\right)+\left[k_{w}w\left(x,t\right)-k_p\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right]\right\}w\left(x,t\right)dx$$

(20)

where D_s(.) is the Dirac delta function; x is the abscissa, measured from the left end of the beam; q(x,t) is the applied external shear force; kw and kp are, respectively, Winkler and shear parameters of the elastic foundation.

Applying Hamilton’s principle, the electromechanical equilibrium equation can be expressed as [54]

$$-{{\displaystyle \int }}_{t_1}^{t_2}\left[\delta {\pi }_s+\delta W_{ex}\right]dt=0$$

(21)

where $\mathsf{\Omega }$ refers to the volume of the composite beam, and $\delta {\pi }_s$ is the variation in the strain energy, which is evaluated as

$$\delta {\pi }_s=\delta {{\displaystyle \int }}_{\mathsf{\Omega }}Hd\mathsf{\Omega }={{\displaystyle \int }}_0^L\left(\overline{M}\delta {\epsilon }_{xx0}+{\overline{M}}^{\prime }\delta {\epsilon }_{xx0,x}\right)d\mathrm{x}={{\displaystyle \int }}_0^L\left[\left(M+\frac{e_{311}}{a_{33}}{M}_D+N_{113}\right)\delta \left(-\frac{{\partial }^{2}w}{\partial x^2}\right)+\left(\frac{{\mu }_{3111}}{a_{33}}{M}_D+N_{111}\right)\delta \left(-\frac{{\partial }^{3}w}{\partial x^3}\right)\right]d\mathrm{x}$$

(22)

The variation in the external work done, $\delta W_{ex}$, is given by

$$\delta W_{ex}=w_b{{\displaystyle \int }}_0^L\left\{q\left(x,t\right)D_s\left(x-x_p\right)+\left[k_{w}w\left(x,t\right)-k_p\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right]\right\}\delta w\left(x,t\right)dx$$

(23)

Substituting Equations (22) and (23) into Equation (21) and evaluating the integrals, the governing electromechanical nonclassical equilibrium equation is derived as

$$\frac{{\partial }^2\overline{M}}{\partial x^2}-\frac{{\partial }^3{\overline{M}}^{\prime }}{\partial x^3}=\left\{q\left(x,t\right)+\left[k_{w}w\left(x,t\right)-k_p\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right]\right\}$$

(24)

For convenience, the governing electromechanical equilibrium equation is expressed in terms of normalized variables. To do this, the following nondimensional quantities are proposed:

$$\left\{\begin{array}{c}W\\ X\\ \overline{q}\end{array}\right\}=\left\{\begin{array}{c}\frac{w}{h}\\ \frac{x}{L}\\ \frac{q{L}^3}{\left({\left[E_c{I}_c\right]}_{eq}+E_p{I}_p\right)}\end{array}\right\} \mathrm{and} \left\{\begin{array}{c}dx\\ \frac{dw}{dX}\\ \frac{d^{4}w}{d{x}^4}\\ \frac{d^{6}w}{d{x}^6}\end{array}\right\}=\left\{\begin{array}{c}LdX\\ \left(\frac{h}{L}\right)\frac{dW}{dX}\\ \left(\frac{h}{L^4}\right)\frac{d^{4}W}{d{X}^4}\\ \left(\frac{h}{L^6}\right)\frac{d^{6}W}{d{X}^6}\end{array}\right\}$$

(25)

Substituting Equations (16) and (25) into Equation (24), the following electromechanical equilibrium equation is obtained:

$$-\left[\left\{{\left(E_c{I}_c\right)}_{eq}+E_p{I}_p\right\}+I_p\frac{e_{311}^2}{a_{33}}\right]\left(\frac{h}{L^4}\right)\left[\frac{d^{4}W\left(X\right)}{d{X}^4}\right]+\left(l^2\left[{\left(E_c{I}_c\right)}_{eq}-E_p{I}_p\right]+I_p\frac{{\mu }_{3111}^2}{a_{33}}\right)\left(\frac{h}{L^6}\right)\frac{{\partial }^{6}W\left(X\right)}{\partial X^6}=\left(1-\frac{{\left(e_{0}a\right)}^2}{L^2}\frac{{\partial }^2}{\partial X^2}\right)\left[q\left(X\right)+\left[h{k}_{w}W\left(X\right)-\frac{h{k}_p}{L^2}\frac{d^{2}W\left(X\right)}{d{X}^2}\right]\right]$$

(26)

3. Analytical Solution Methodology

Within this section, an analytical solutions procedure for the nonclassical electrochemical equilibrium equation is developed. To study and analyze the nonclassical bending behavior of piezoelectrically layered perforated nanobeams embedded in an elastic foundation incorporating the flexoelectric effect based on the nonlocal strain gradient theory, the solution of the equilibrium equation, Equation (26), should be provided. Based on Fourier expansion, closed forms for bending profiles for simply supported piezoelectrically layered perforated nanobeams are derived for different external loading patterns. The boundary conditions associated with simply supported beams are expressed as follows:

$${\left.W\left(X\right)\right|}_{at X=0}={\left.W\left(X\right)\right|}_{at X=1}={\left.\frac{d^{2}W\left(X\right)}{d{X}^2}\right|}_{at X=0}={\left.\frac{d^{2}W\left(X\right)}{d{X}^2}\right|}_{at X=1}=0$$

(27)

The normalized transverse bending deflection and its derivatives with respect to X can be expressed by the Fourier series as

$$\left\{\begin{array}{c}W\left(X\right)\\ \frac{d^{2}W\left(X\right)}{d{X}^2}\\ \frac{d^{4}W\left(X\right)}{d{X}^4}\\ \frac{d^{6}W\left(X\right)}{d{X}^6}\end{array}\right\}=\left\{\begin{array}{c}{\displaystyle \sum }_{n=1}^{\infty }{W}_n\sin \left(n\pi X\right)\\ -{\displaystyle \sum }_{n=1}^{\infty }{\left(n\pi \right)}^2{W}_n\sin \left(n\pi X\right)\\ {\displaystyle \sum }_{n=1}^{\infty }{\left(n\pi \right)}^4{W}_n\sin \left(n\pi X\right)\\ -{\displaystyle \sum }_{n=1}^{\infty }{\left(n\pi \right)}^6{W}_n\sin \left(n\pi X\right)\end{array}\right\}$$

(28)

On the other hand, the applied external shear force is expanded using the Fourier series as

$$\left\{\begin{array}{c}q\left(X\right)\\ \frac{d^{2}q\left(X\right)}{d{X}^2}\end{array}\right\}=\left\{\begin{array}{c}{\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{Q}_n\sin \left(n\pi X\right), \\ -{\displaystyle {\displaystyle \sum }_{n=1}^{\infty }}{\left(n\pi \right)}^2{Q}_n\sin \left(n\pi X\right) \end{array}\right\}$$

(29)

where Q_n refers to the Fourier expansion coefficients, which depends on the applied loading pattern. For a uniformly distributed shear force [30,40]:

$$Q_n=\left\{\begin{array}{c}\frac{2q_0}{\left(n\pi \right)}\left(1-\cos n\pi \right) \mathrm{Uniformly\; distributed\; load\; with\; intensity} q_0\\ \frac{2P}{L}\sin \left(n\pi {\overline{x}}_p\right) \mathrm{Concentrated\; load\; of\; intensity\; PD}\left(x-x_p\right) \mathrm{applied\; at\; a\; distance} {\overline{x}}_p \end{array}\right.$$

(30)

Substituting Equations (27)–(30) into Equation (26) yields

$$\left[{{\displaystyle \sum }}_{n=1}^{\infty }\left\{\left[\left\{{\left(E_c{I}_c\right)}_{eq}+E_p{I}_p\right\}+I_p\frac{e_{311}^2}{a_{33}}\right]\left(\frac{h}{L^4}\right){\left(n\pi \right)}^4+\left(l^2\left[{\left(E_c{I}_c\right)}_{eq}-E_p{I}_p\right]+I_p\frac{{\mu }_{3111}^2}{a_{33}}\right)\left(\frac{h}{L^6}\right){\left(n\pi \right)}^6+\left\{\left(k_{w}h+\frac{h{k}_p}{L^2}{\left(n\pi \right)}^2\right)\left(1+\frac{{\left(e_{0}a\right)}^2}{L^2}{\left(n\pi \right)}^2\right)\right\}\right\}{W}_n\sin \left(n\pi X\right)\right]=-{{\displaystyle \sum }}_{n=1}^{\infty }\left(1+\frac{{\left(e_{0}a\right)}^2}{L^2}\right)Q_n\sin \left(n\pi X\right)$$

(31)

Simplifying Equation (31), the normalized nonclassical electromechanical bending deflection profile,$W^{NELECM}\left(X\right)$, is derived as

$$W^{NELECM}\left(X\right)={\displaystyle \sum }_{n=1}^{\infty }\frac{-\left(1+\frac{{\left(e_{0}a\right)}^2}{L^2}{\left(n\pi \right)}^2\right)Q_n\sin \left(n\pi X\right)}{\left\{\left(\frac{h}{L^4}\right){\left(n\pi \right)}^4\left(\left[\left\{{\left(E_c{I}_c\right)}_{eq}+E_p{I}_p\right\}+I_p\frac{e_{311}^2}{a_{33}}\right]+\left(l^2\left[{\left(E_c{I}_c\right)}_{eq}+E_p{I}_p\right]+I_p\frac{{\mu }_{3111}^2}{a_{33}}\right){\left(\frac{n\pi }{L}\right)}^2\right)+\left(1+\frac{{\left(e_{0}a\right)}^2}{L^2}{\left(n\pi \right)}^2\right)\left(k_{w}h+\frac{h{k}_p}{L^2}{\left(n\pi \right)}^2\right)\right\}}$$

(32)

Neglecting the piezoelectric and flexoelectricity effects, the mechanical-size-dependent transverse bending deflection profile, $W^{NMEC}\left(X\right)$, can be expressed as

$$W^{NMEC}\left(X\right)={\displaystyle \sum }_{n=1}^{\infty }\frac{-\left(1+\frac{{\left(e_{0}a\right)}^2}{L^2}{\left(n\pi \right)}^2\right)Q_n\sin \left(n\pi X\right)}{\left\{\left(\frac{h}{L^4}\right){\left(n\pi \right)}^4\left(\left[\left\{{\left(E_c{I}_c\right)}_{eq}\right\}\right]+\left(l^2\left[{\left(E_c{I}_c\right)}_{eq}\right]\right){\left(\frac{n\pi }{L}\right)}^2\right)+\left(1+\frac{{\left(e_{0}a\right)}^2}{L^2}{\left(n\pi \right)}^2\right)\left(k_{w}h+\frac{h{k}_p}{L^2}{\left(n\pi \right)}^2\right)\right\}}$$

(33)

Neglecting the nonlocal effect, e₀a = 0, as well as the piezoelectric and flexoelectricity effects, the nonclassical mechanical-size-dependent bending based on the microstructure effect can be expressed as

$${\left[W\left(X\right)\right]}_{e_{0}a=0}={{\displaystyle \sum }}_{n=1}^{\infty }\frac{-Q_n\sin \left(n\pi X\right)}{\left\{\left(\frac{h}{L^4}\right){\left(n\pi \right)}^4\left(\left[\left\{{\left(E_c{I}_c\right)}_{eq}\right\}\right]+\left(l^2\left[{\left(E_c{I}_c\right)}_{eq}\right]\right){\left(\frac{n\pi }{L}\right)}^2\right)+\left(k_{w}h+\frac{h{k}_p}{L^2}{\left(n\pi \right)}^2\right)\right\}}$$

(34)

Neglecting the nonclassical effects, the classical electromechanical transverse bending profile $W^{CELECM}\left(X\right)$ can be expressed as

$$W^{CELECM}\left(X\right)={\displaystyle \sum }_{n=1}^{\infty }\frac{-Q_n\sin \left(n\pi X\right)}{\left\{\left(\frac{h}{L^4}\right){\left(n\pi \right)}^4\left(\left[\left\{{\left(E_c{I}_c\right)}_{eq}+E_p{I}_p\right\}+I_p\frac{e_{311}^2}{a_{33}}\right]+\left(I_p\frac{{\mu }_{3111}^2}{a_{33}}\right){\left(\frac{n\pi }{L}\right)}^2\right)+\left(k_{w}h+\frac{h{k}_p}{L^2}{\left(n\pi \right)}^2\right)\right\}}$$

(35)

Neglecting the nonclassical as well as the piezoelectric and flexoelectricity effects, the classical mechanical transverse bending profile, $W^{CMECES}\left(X\right)$, with an elastic foundation can be expressed as

$$W^{CMECES}\left(X\right)={\displaystyle \sum }_{n=1}^{\infty }\frac{-Q_n\sin \left(n\pi X\right)}{\left\{\left(\frac{h}{L^4}\right){\left(n\pi \right)}^4{\left(E_c{I}_c\right)}_{eq}+\left(k_{w}h+\frac{h{k}_p}{L^2}{\left(n\pi \right)}^2\right)\right\}}$$

(36)

Neglecting the nonclassical as well as the piezoelectric and flexoelectricity and elastic foundation effects, the classical mechanical transverse bending profile, $W^{CMECRS}\left(X\right)$, without an elastic foundation can be expressed as

$$W^{CMECRS}\left(X\right)={\displaystyle \sum }_{n=1}^{\infty }\frac{-Q_n\sin \left(n\pi X\right)}{\left(\frac{h}{L^4}\right){\left(n\pi \right)}^4{\left(E_c{I}_c\right)}_{eq}}$$

(37)

4. Numerical Results and Discussion

Within this section, the accuracy of the developed procedure is checked and verified. Furthermore, the applicability of the analytical procedure is demonstrated. Numerical experiments are performed to explore the design variables effects on the electromechanical as well as mechanical bending behaviors of piezoelectrically layered perforated nanobeams under different loading profiles.

4.1. Verification of the Developed Methodology

To check the accuracy of the proposed methodology to accurately detect the bending response of nanobeams under different loading conditions, consider a simply supported nanobeam with the following material and geometrical characteristics [55]: the beam is assumed to be made of an elastic material with a modulus of elasticity E = 210 GPa and Poisson’s ratio ν = 0.24. The geometrical variable of the beam is as follows: beam length L = 50 nm; beam cross-sectional area b × H = 5 × 5 nm. To consider a solid beam, the perforation filling ratio is taken as α = 1. The beam is subjected to a central concentrated load of intensity, P = 50 nN, while the distributed load intensity is taken as Q = P/L N/m. The classical and nonclassical bending behaviors of this problem previously discussed by the same problem have been analyzed by Şimşek [55] and Hashemian et al. [26] for solid EBNBT. The developed approach is applied to detect the bending profile for EBNBT with a perforation filling ratio of α = 1. Comparison of the nondimensional deflection, W/H, for EBNB under a central point load and uniformly distributed load are shown in Figure 2 and Figure 3, respectively. An excellent agreement is noticed with the corresponding results reported in Şimşek [55] and Hashemian et al. [26].

We seek a deeper verification of the proposed analytical procedure to check the accuracy of the analytical procedure to efficiently investigate the static behavior of the beam structure embedded in an elastic foundation. For this comparison, the following nondimensional maximum central deflection parameter is defined: ${\Delta }_{max}=w_{max}\left(\frac{E_{eq}I}{q{L}^4}\right)$, where E_eq is the equivalent elasticity modulus, which can be expressed as $E_{eq}=\lambda +2G$, where λ is Lame’s constant, $\lambda =\frac{\nu E}{\left(1+\nu \right)\left(1-2\nu \right)}$, and G is the rigidity modulus, $G=\frac{E}{2\left(1+\nu \right)}$. q is the intensity of the distributed load and L is the beam length. The nondimensional elastic foundation parameters are defined as $K_w=\frac{k_w{L}^4}{E_{eq}I} \mathrm{and} K_p=\frac{k_p{L}^2}{E_{eq}I}.$ Comparisons of the nondimensional maximum central deflection parameter, ${\overline{W}}_{max}$, for a simply supported beam (SS) under a uniformly distributed load for different values of the nondimensional elastic foundation parameters, K_w and K_p, for beam aspect ratio, L/H = 120, are depicted in

Table 1. Comparison of the estimated maximum classical nondimensional central transverse deflection ${\Delta }_{max}$ for simply supported (SS) beam embedded in two-parameter elastic foundation under uniformly distributed load of intensity q at beam aspect ratio L/H = 120 at different values of the elastic foundation parameters and Poisson’s ratio ν = 0.3.

Foundation Parameters		${\mathsf{\Delta }}_{\mathit{m}\mathit{a}\mathit{x}}={\mathit{w}}_{\mathit{m}\mathit{a}\mathit{x}}\left(\frac{\mathit{E}\mathit{I}}{\mathit{q}{\mathit{L}}^4}\right)\times {10}^2$
Kw	${\mathit{K}}_{\mathit{p}}$	Present	Ref. [56]	Analytical Ref. [56]	Ref. [57]	% Error
0	0	1.302083	1.302290	1.302290	1.3033	0.093378347
	10	0.644771	0.644827	0.644827	0.6457	0.143874864
	25	0.366091	0.366111	0.366111	0.3671	0.274856987
10	0	1.180396	1.180567	1.180567	1.1814	0.084983917
	10	0.613275	0.613325	0.613326	0.6141	0.134342941
	25	0.355649	0.355668	0.355668	0.3566	0.266685362
102	0	0.64002	0.640074	0.640074	0.6403	0.043729502
	10	0.425557	0.425582	0.425582	0.4261	0.127434874
	25	0.282834	0.282846	0.282846	0.2836	0.270098731

. Excellent agreement with the results reported by [56,57] is observed.

4.2. Numerical Experiments and Discussion

To demonstrate the effectiveness of the developed procedure to efficiently explore the classical and nonclassical electromechanical as well as the mechanical bending behavior of a piezoelectrically layered perforated nanobeam embedded in an elastic medium under uniformly or concentrated loads, consider a simply supported beam with geometrical and material characteristics shown in

Table 2. Geometrical and material characteristics of the piezoelectric composite nanobeam.

Parameters	Thickness, h (nm)	Length, L (nm)	Width, wb(nm)	Young’s Modulus E (GPa)	Mass Density ρ (Kg/m3)	Poisson’s Ratio, ν	e311 (C/m2)	μ3111 (C/m)	a33 N/(m2.K)	l
Elastic core	3	100	5	0.130	1380	0.24	----	----	------	--
Piezoelectric Layer	1	100	5	132	7500	0.27	−4.1	5 × 10−8	7.124 × 10−9	--

. Both the uniformly distributed load q with intensity 1 N/m and central point load P = q × L are considered and discussed. Numerical experiments are performed by applying the developed procedure to explore the influences of the perforation configuration parameters, the elastic foundation parameters, the strain gradient parameter, the nonlocal parameter, the piezoelectric coefficient, and the electric-field–strain-gradient coupling coefficient on the classical and size-dependent electromechanical as well as the mechanical bending behaviors.

To investigate the electromechanical as well as the mechanical bending behaviors, the following nondimensional variables are defined: for the electromechanical behavior, ${\Delta }_{maxq}=w_{max}\left(\frac{{\left(E_{e}I\right)}_{ceq}+{\left(EI\right)}_p}{q{L}^4}\right)\times {10}^2, {\Delta }_{maxp}=w_{max}\left(\frac{{\left(E_{e}I\right)}_{ceq}+{\left(EI\right)}_p}{p{L}^3}\right)\times {10}^2, K_w=\frac{k_w{L}^4}{{\left(E_{e}I\right)}_{ceq}+{\left(EI\right)}_p}, \mathrm{and} K_p=\frac{k_p{L}^2}{{\left(E_{e}I\right)}_{ceq}+{\left(EI\right)}_p}$, and for the mechanical bending behavior, h_p = 0, ${\Delta }_{maxq}=w_{max}\left(\frac{{\left(EI\right)}_{ceq}}{q{L}^4}\right)\times {10}^2, {\Delta }_{maxp}=w_{max}\left(\frac{{\left(E_{e}I\right)}_{ceq}}{p{L}^3}\right)\times {10}^2, K_w=\frac{k_w{L}^4}{{\left(E_{e}I\right)}_{ceq}}, \mathrm{and} K_p=\frac{k_p{L}^2}{{\left(E_{e}I\right)}_{ceq}}.$

The dependency of the nonclassical and classical electromechanical and mechanical bending deflection, Δ_max, on the nondimensional elastic foundation parameter K_w at different values of perforation filling ratio is illustrated in Figure 4. It is observed that the maximum nondimensional transverse deflection nonlinearly decreases with the increase in the nondimensional foundation parameter K_w due to decreasing the overall system flexibility for classical and nonclassical electromechanical and mechanical behaviors. Incorporating the piezoelectric effect results in a greater stiffening effect, leading to a greater decrease in the maximum nondimensional transverse deflection compared with the corresponding mechanical behavior. Considering the applied loading pattern, larger values of the maximum nondimensional bending deflection are detected under the action of the applied point load compared with the corresponding uniformly distributed load. On the other hand, regarding the effect of filling ratio, it may be seen that the mechanical behavior is significantly influenced by the filling ratio while it has no significant effect on the electromechanical bending behavior. Increasing the perforation filling ratio results in smaller values of the maximum nondimensional transverse deflection. Additionally, it may be noticed that incorporating the nonlocality effect leads to a softening effect and, thus, larger values of the nondimensional bending deflection compared with the corresponding classical behavior. On the other hand, introduction of the nonclassical strain gradient effect results in a stiffening effect, so smaller values of the nonclassical bending deflection are detected compared with the corresponding classical behavior.

Variations in the maximum nondimensional transverse deflection Δ_max with the nondimensional elastic foundation parameter K_p for nonclassical and classical electromechanical and mechanical behavior under uniformly distributed and point loads at different values of perforation filling ratio are depicted in Figure 5. It may be observed that almost the same qualitative response is observed; the maximum nondimensional transverse deflection nonlinearly decreases with the increase in the nondimensional elastic foundation parameter, K_p.

For deep investigation of the influence of the elastic foundation parameter K_w on the maximum nondimensional transverse deflection, the relative percentage difference of the maximum transverse deflection due to increasing the nondimensional elastic foundation parameter Kw is given by ${\left(\%D{\Delta }_{max}\right)}_{K_w}=100\times abs\left[\frac{{\left({\Delta }_{max}\right)}_{K_{w=0}}-{\left({\Delta }_{max}\right)}_{K_{w=200}}}{{\left({\Delta }_{max}\right)}_{K_w=0}}\right].$ The relative percentage difference ${\left(\%D{\Delta }_{max}\right)}_{K_w}$ for the electromechanical and mechanical nonclassical and classical behaviors at different values of α and nonclassical parameters for point and uniformly distributed loading is depicted in

Table 3. The relative percentage difference in the transverse deflection, ${\left(\%D{\Delta }_{max}\right)}_{K_w}$, due to increasing the nondimensional elastic foundation parameter, K_w, for electromechanical and mechanical nonclassical and classical behaviors for different values of filling ratio and nonclassical parameters under point and uniformly distributed loads at perforation filling ratio, α = 0.25 and 0.5, at N = 4 and beam aspect ratio L/h = 20, e₃₁₁ = −4.1 C/m², μ₃₁₁₁ = 5 × 10⁻⁸ C/m, a₃₃ = 7.124 × 10⁻⁹ C/m² K, K_w = 0 to 200 and K_p = 0.

%DΔmax	$\frac{\mathit{e}\mathit{a}}{\mathit{h}}=4 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=0$				$\frac{\mathit{e}\mathit{a}}{\mathit{h}}=0 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=4$
	Point Load		Distributed Load		Point Load		Distributed Load
	α = 0.25	α = 0.5	α = 0.25	α = 0.5	α = 0.25	α = 0.5	α = 0.25	α = 0.5
NCL_Elect	69.60735	69.60589	74.53625	74.53472	59.03702	59.03509	59.34576	59.34389
CL_Elect	65.98284	65.98116	67.14135	67.13964	65.98284	65.98116	67.14135	67.13964
NCL_Mec	82.49071	74.88125	88.2882	80.53019	78.0187	66.70369	78.4219	67.05033
CL_Mech	82.31208	72.84932	83.74434	74.14721	82.31208	72.84932	83.74434	74.14721

. It is observed that increasing the filling ratio significantly affects this relative percentage difference for the mechanical behavior, while it has an insignificant effect on the electromechanical behavior. Additionally, the uniformly distributed loading produces larger values of the relative percentage difference compared with the corresponding cases of the concentrated point loads for both electromechanical and mechanical behaviors.

The relative percentage difference of the maximum nondimensional transverse deflection due to increasing the nondimensional elastic foundation parameter K_p is given by ${\left(\%D{\Delta }_{max}\right)}_{K_p}=100\times abs\left[\frac{{\left({\Delta }_{max}\right)}_{K_{p=0}}-{\left({\Delta }_{max}\right)}_{K_{p=20}}}{{\left({\Delta }_{max}\right)}_{K_p=0}}\right].$ 

Table 4. The relative percentage difference in the transverse deflection, ${\left(\%D{\Delta }_{max}\right)}_{K_p}$ _x, due to increasing the nondimensional elastic foundation parameter Kp for electromechanical and mechanical nonclassical and classical behaviors for different values of filling ratio and nonclassical parameters under point and uniformly distributed loads at perforation filling ratio, α = 0.25 and 0.5, at N = 4 and beam aspect ratio L/h = 20, e₃₁₁ = −4.1 C/m², μ₃₁₁₁ = 5 × 10⁻⁸ C/m, a₃₃ = 7.124 × 10⁻⁹ C/m² K, K_p = 0 to 20 and K_w = 0.

%DΔmax	$\frac{\mathit{e}\mathit{a}}{\mathit{h}}=4 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=0$				$\frac{\mathit{e}\mathit{a}}{\mathit{h}}=0 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=4$
	Point Load		Distributed Load		Point Load		Distributed Load
	α = 0.25	α = 0.5	α = 0.25	α = 0.5	α = 0.25	α = 0.5	α = 0.25	α = 0.5
NCL_Elect	71.8572	71.8557	73.8245	73.8230	58.7368	58.7349	59.0230	59.0211
CL_Elect	65.8932	65.8915	66.7878	66.7861	65.8932	65.8915	66.7878	66.7861
NCL_Mec	86.1856	77.9658	87.5840	79.8012	77.8320	66.4352	78.1861	66.7520
CL_Mech	82.5244	72.8639	83.4433	73.8111	82.5244	72.8639	83.4433	73.8111

illustrates the relative percentage difference ${\left(\%D{\Delta }_{max}\right)}_{K_p}$ for the electromechanical and mechanical nonclassical and classical behaviors at different values of α and nonclassical parameters for point and uniformly distributed loading conditions. It may be noticed that the relative percentage difference ${\left(\%D{\Delta }_{max}\right)}_{K_p}$ is more sensitive to the nondimensional foundation parameter K_p compared to K_w. Larger values of this relative percentage difference are obtained compared to those obtained due to K_w.

The dependency of the nondimensional bending deflection on the perforation filling ratio for electromechanical and mechanical nonclassical and classical behaviors under uniformly distributed and point loads is shown in Figure 6. It may be noticed that, due to decreasing the overall beam flexibility, the maximum nondimensional mechanical bending deflection decreases with the increase in the perforation filling ratio, while there is no significant effect of the filling ratio on the electromechanical bending behavior.

For quantitative investigation of the significant effect of the perforation filling ratio α on the maximum nondimensional transverse deflection, the relative percentage difference of the maximum nondimensional transverse deflection due to increasing α at a constant value of the number of hole rows N is given by ${\left(\%D{\Delta }_{max}\right)}_{\alpha }=100\times abs\left[\frac{{\left({\Delta }_{max}\right)}_{\mathsf{\alpha }=0.2}-{\left({\Delta }_{max}\right)}_{\mathsf{\alpha }=1}}{{\left({\Delta }_{max}\right)}_{\mathsf{\alpha }=1}}\right].$ 

Table 5. The relative percentage difference in the transverse deflection ${\left(\%D{\Delta }_{max}\right)}_{\alpha }$ due to increasing the perforation filling ratio α for electromechanical and mechanical nonclassical and classical behaviors for different values of the nondimensional elastic foundation and nonclassical parameters under point and uniformly distributed loads at perforation filling ratio, α = 0.2 and 0.5, at N = 4 and beam aspect ratio L/h = 20, e₃₁₁ = −4.1 C/m², μ₃₁₁₁ = 5 × 10⁻⁸ C/m, a₃₃ = 7.124 × 10⁻⁹ C/m² K, and α = 0.2 to 1.

%DΔmax	$\frac{\mathit{e}\mathit{a}}{\mathit{h}}=4 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=0$				$\frac{\mathit{e}\mathit{a}}{\mathit{h}}=0 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=4$
	Point Load		Distributed Load		Point Load		Distributed Load
	Kw = Kp = 0	Kw = 20, Kp = 2	Kw = Kp = 0	Kw = 20, Kp = 2	Kw = Kp = 0	Kw = 20, Kp = 2	Kw = Kp = 0	Kw = 20, Kp = 2
NCL_Elect	0.0162	0.0107	0.0163	0.0103	0.0164	0.0127	0.0163	0.0127
CL_Elect	0.0163	0.0117	0.0162	0.0116	0.0163	0.0117	0.0162	0.0116
NCL_Mec	199.8589	81.4321	199.8589	72.6340	199.8589	106.9395	199.8589	106.3326
CL_Mech	199.8588	91.8268	199.8590	89.3870	199.8588	91.8268	199.8590	89.3870

shows the relative percentage difference ${\left(\%D{\Delta }_{max}\right)}_{\mathsf{\alpha }}$ for the electromechanical and mechanical nonclassical and classical behaviors at different values of the elastic foundation and nonclassical parameters for point and uniformly distributed loading conditions. It is seen that the filling ratio significantly affects the mechanical transverse deflection, while it has an insignificant influence on the electromechanical behavior. This influence decreases with the incorporation of the elastic foundation parameters.

The variation in the maximum nondimensional transverse bending deflection with the number of holes throughout the beam cross-section N for electromechanical and mechanical bending behaviors is depicted in Figure 7.

Keeping the value of perforation filling ratio constant, investigating the quantitative effect of the number of hole rows N on the maximum nondimensional transverse deflection, the relative percentage difference of the maximum nondimensional transverse deflection due to increasing N is defined as ${\left(\%D{\Delta }_{max}\right)}_N=100\times abs\left[\frac{{\left({\Delta }_{max}\right)}_{\mathrm{N}=10}-{\left({\Delta }_{max}\right)}_{\mathrm{N}=1}}{{\left({\Delta }_{max}\right)}_{\mathrm{N}=1}}\right].$ The relative percentage difference ${\left(\%D{\Delta }_{max}\right)}_N$ for electromechanical and mechanical nonclassical and classical behaviors at different values of the nondimensional elastic foundation and nonclassical parameters for point and uniformly distributed loading conditions is presented in

Table 6. The relative percentage difference in the transverse deflection ${\left(\%D{\Delta }_{max}\right)}_N$ due to increasing the number of hole rows N for electromechanical and mechanical nonclassical and classical behaviors for different values of the nondimensional elastic foundation and nonclassical parameters under point and uniformly distributed loads at perforation filling ratio, α = 0.5, and beam aspect ratio L/h = 20, e₃₁₁ = −4.1 C/m², μ₃₁₁₁ = 5 × 10⁻⁸ C/m, a₃₃ = 7.124 × 10⁻⁹ C/m² K, and N = 1 to 10.

%DΔmax	$\frac{\mathit{e}\mathit{a}}{\mathit{h}}=4 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=0$				$\frac{\mathit{e}\mathit{a}}{\mathit{h}}=0 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=4$
	Point Load		Distributed Load		Point Load		Distributed Load
	Kw = Kp = 0	Kw = 20, Kp = 2	Kw = Kp = 0	Kw = 20, Kp = 2	Kw = Kp = 0	Kw = 20, Kp = 2	Kw = Kp = 0	Kw = 20, Kp = 2
NCL_Elect	0.0081	0.0053	0.0081	0.0051	0.0081	0.0063	0.0082	0.0063
CL_Elect	0.0081	0.0058	0.0080	0.0057	0.0081	0.0058	0.0080	0.0057
NCL_Mec	34.1643	22.5597	34.1643	21.3813	34.1643	26.3062	34.1643	26.24783
CL_Mech	34.1643	24.2761	34.1643	24.0080	34.1643	24.2761	34.1643	24.0080

. It may be noted that significant values of the relative percentage difference are investigated for the mechanical behaviors for both types of loading conditions while insignificant values are detected for the electromechanical bending behavior.

It may be seen that the electromechanical bending behavior is independent of the number of holes while the mechanical bending behavior is significantly affected; larger values of the maximum nondimensional bending deflection are obtained with the increase in the number of holes due to increasing the system flexibility.

The dependency of the nonclassical electromechanical and mechanical bending behavior on the nondimensional strain gradient parameter is illustrated in Figure 8. It is seen that the maximum nondimensional transverse deflection is nonlinearly dependent on the nondimensional strain gradient parameter. Increasing the nondimensional strain gradient parameter results in decreasing the maximum nondimensional bending deflection due to increasing the overall system stiffness. Additionally, decreasing the beam thickness leads to a greater decrease in the strain gradient effect; a greater decrease in the nonclassical electromechanical and mechanical bending behavior is observed at small values of beam thickness. It is also observed that the classical bending behavior is detected at (e₀a/h) = (l/h).

Variations in the maximum nondimensional bending deflection with the nondimensional nonlocal parameter for electromechanical and mechanical behavior are shown in Figure 9. It may be noticed that the nonclassical electromechanical and mechanical bending behavior is also nonlinearly dependent on the nondimensional nonlocal parameter. Increasing the nondimensional nonlocal parameter produces larger values of the maximum nondimensional bending deflection due to the associated softening effect. In addition, a more significant effect is noticed at smaller beam thickness.

To investigate the influence of the strain gradient and the nonlocal parameters on the nonclassical electromechanical as well as the mechanical bending behaviors, the relative percentage difference of the maximum nondimensional transverse deflection due to increasing the nondimensional strain gradient parameter (l/h) is defined as ${\left(\%D{\Delta }_{max}\right)}_{l/h}=100\times abs\left[\frac{{\left({\Delta }_{max}\right)}_{l/h=4}-{\left({\Delta }_{max}\right)}_{l/h=1}}{{\left({\Delta }_{max}\right)}_{l/h=1}}\right]$, while it is defined due to increasing the nondimensional nonlocal parameter (ea/h) as ${\left(\%D{\Delta }_{max}\right)}_{ea/h}=100\times abs\left[\frac{{\left({\Delta }_{max}\right)}_{ea/h=4}-{\left({\Delta }_{max}\right)}_{ea/h=1}}{{\left({\Delta }_{max}\right)}_{ea/h=1}}\right]$. The relative percentage difference, ${\left(\%D{\Delta }_{max}\right)}_{l/h}$ and ${\left(\%D{\Delta }_{max}\right)}_{ea/h}$, for nonclassical electromechanical and mechanical behaviors at different values of the beam aspect ratio for point and uniformly distributed loading conditions is presented in

Table 7. The relative percentage difference in the transverse deflection due to increasing the nondimensional strain gradient parameter (l/h), ${\left(\%D{\Delta }_{max}\right)}_{l/h}$, and the nondimensional nonlocal parameter (ea/h), ${\left(\%D{\Delta }_{max}\right)}_{ea/h}$, for electromechanical and mechanical nonclassical behaviors under point and uniformly distributed loads at beam aspect ratio L/h = 10, 20, e₃₁₁ = −4.1 C/m², μ₃₁₁₁ = 5 × 10⁻⁸ C/m, a₃₃ = 7.124 × 10⁻⁹ C/m² K, K_w = 50, K_p = 5, α = 0.5, and N = 4.

%DΔmax	${\left(\mathbf{\%}\mathit{D}{\mathsf{\Delta }}_{\mathit{m}\mathit{a}\mathit{x}}\right)}_{\left(\frac{\mathit{l}}{\mathit{h}}\right)}, \frac{\mathit{e}\mathit{a}}{\mathit{h}}=1 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=1 \mathbf{to} 4$				${\left(\mathbf{\%}\mathit{D}{\mathsf{\Delta }}_{\mathit{m}\mathit{a}\mathit{x}}\right)}_{\left(\frac{\mathit{e}\mathit{a}}{\mathit{h}}\right)}, \frac{\mathit{e}\mathit{a}}{\mathit{h}}=1 \mathbf{to} 4 \mathbf{and} \frac{\mathit{l}}{\mathit{h}}=1$
	Point Load		Distributed Load		Point Load		Distributed Load
	L/h = 10	L/h = 20	L/h = 10	L/h = 20	L/h = 10	L/h = 20	L/h = 10	L/h = 20
NCL_Elect	41.13188	16.64304	39.5017	14.6687	53.9718	22.5435	37.1580	13.8975
NCL_Mec	37.550	14.85231	35.4802	12.4886	46.1735	20.4469	28.285	10.7989

. It may be noted that larger values of the nonclassical electromechanical and mechanical behaviors are detected at smaller values of the beam aspect ratio.

Variations in the maximum normalized transverse deflection Δ_max with the piezoelectric coefficient e₃₁₁ for nonclassical and classical electromechanical and mechanical behaviors under uniformly distributed and point loads at different values of the nonclassical and elastic foundation parameters are depicted in Figure 10. It may be noticed that increasing the absolute value of the piezoelectric coefficient results in the stiffening effect. This leads to decreasing the maximum nondimensional electromechanical bending deflection. In addition, it is noticed that larger values of the nondimensional bending deflection are detected at smaller beam aspect ratios. Additionally, a more significant difference between the nonclassical and classical bending behavior is observed due to the point load compared with the corresponding cases of the uniformly distributed load.

The dependency of the maximum normalized transverse deflection Δ_max on the electric-field–strain-gradient coupling coefficient μ₃₁₁₁ for nonclassical and classical electromechanical and mechanical behaviors under uniformly distributed and point loads at different values of the nonclassical and elastic foundation parameters is illustrated in Figure 11. It may be noticed that the electromechanical bending deflection is insignificantly affected by the electric-field–strain-gradient coupling coefficient; almost constant electromechanical bending behavior is detected.

The quantitative influences of the piezoelectric coefficient e₃₁₁ as well as the electric-field–strain-gradient coupling coefficient μ₃₁₁₁ on the electromechanical bending behavior could be investigated by defining the relative percentage difference of the maximum nondimensional transverse deflection due to increasing the piezoelectric coefficient e₃₁₁ as ${\left(\%D{\Delta }_{max}\right)}_{e311}=100\times abs\left[\frac{{\left({\Delta }_{max}\right)}_{e311=15}-{\left({\Delta }_{max}\right)}_{e311=0}}{{\left({\Delta }_{max}\right)}_{e311=0}}\right]$, while it is defined due to increasing the electric-field–strain-gradient coupling coefficient μ₃₁₁₁ as ${\left(\%D{\Delta }_{max}\right)}_{\mu 3111}=100\times abs\left[\frac{{\left({\Delta }_{max}\right)}_{\mu 3111=10}-{\left({\Delta }_{max}\right)}_{\mu 3111=0}}{{\left({\Delta }_{max}\right)}_{\mu 3111=0}}\right]$. The relative percentage difference, ${\left(\%D{\Delta }_{max}\right)}_{e311}$ and ${\left(\%D{\Delta }_{max}\right)}_{\mu 3111}$, for nonclassical and classical electromechanical bending behaviors at different values of the beam aspect ratio for point and uniformly distributed loading conditions is shown in

Table 8. The relative percentage difference in the transverse deflection due to increasing the piezoelectric coefficient e₃₁₁ and the electric-field–strain-gradient coupling coefficient (μ₃₁₁₁) for electromechanical and mechanical nonclassical and classical behaviors for different values of filling ratio under point and uniformly distributed loads at beam aspect ratio L/h = 20, 50, abs(e₃₁₁) = −(0 to 15) C/m², abs (μ₃₁₁₁) = (0 to 10) × 10⁻⁸ C/m, a₃₃ = 7.124 × 10⁻⁹ C/m² K, K_w = 50, K_p = 5, α = 0.5, N = 4, $\frac{ea}{h}=4 \mathrm{and} \frac{l}{h}=1$.

%DΔmax	${\left(\mathbf{\%}\mathit{D}{\mathsf{\Delta }}_{\mathit{m}\mathit{a}\mathit{x}}\right)}_{\left(\mathit{e}311\right)},$				${\left(\mathbf{\%}\mathit{D}{\mathsf{\Delta }}_{\mathit{m}\mathit{a}\mathit{x}}\right)}_{\left(\mathit{\mu }3111\right)}$
	Point Load		Distributed Load		Point Load		Distributed Load
	L/h = 20	L/h = 50	L/h = 20	L/h = 50	L/h = 20	L/h = 50	L/h = 20	L/h = 50
NCL_Elect	7.3857	8.5263	7.0698	8.2201	0.9274	0.1629	0.2627	0.0547
CL_Elect	8.6976	8.7156	8.4846	8.4896	0.6451	0.1120	0.3654	0.0585

. It is seen that larger values of the relative percentage differences are detected due to increasing the piezoelectric coefficient, e₃₁₁, especially at larger values of beam aspect ratio, while insignificant values are observed due to increasing the electric-field–strain-gradient coupling coefficient, μ₃₁₁₁.

5. Conclusions

A nonclassical size-dependent model is developed to explore the electromechanical bending response of a perforated nanobeam layered with two piezoelectric face sheets embedded in an elastic foundation including the flexoelectric effect. A squared perforation pattern with regular order is considered for the perforated core. The equivalent geometrical variables of the perforated core are expressed in explicit forms. In the context of 3D continuum mechanics, all kinematics and kinetics relations are expressed based on the Euler Bernoulli beam theory (EBBT). The Hamilton Principle is adopted to derive the electromechanical equilibrium equation, which is analytically solved for the bending behavior of a simply supported beam structure under uniform and concentrated loads. Numerical experiments are performed to investigate influences of the perforation parameters, elastic foundation parameters, nonclassical parameters, and the piezoelectric and flexoelectric parameters on the bending performance of piezoelectrically layered perforated nanobeams. The following remarks are derived:

✓

The electromechanical bending is affected by the piezoelectric coefficient, e_311. It decreases with the increase in the absolute value of the piezoelectric coefficient. Due to increasing the absolute value of the piezoelectric coefficient from 0 to 15, the relative percentage difference of the nondimensional maximum transverse deflection reaches about 8.7% for a point load and about 8.5% for a uniformly distributed load at a beam slenderness ratio of 50.

✓

On the contrary, the electromechanical bending behavior is slightly affected by the electric-field–strain-gradient coupling coefficient, μ3111; almost constant electromechanical bending behavior is obtained with increasing μ3111. With the increase in the absolute value of the electric-field–strain-gradient coupling coefficient, μ3111, from 0 to 10, the relative percentage difference of the nondimensional maximum transverse deflection reaches about 0.9274% for a point load and about 0.3654% for a uniformly distributed load.

✓

The elastic foundation parameters have a significant effect on the electromechanical and mechanical bending behavior. The bending deflection is nonlinearly dependent on the elastic foundation parameters, K_w and K_p. The elastic foundation introduces a stiffening effect and, thus, smaller values of the bending deflection are produced with increasing these parameters. Increasing the nondimensional foundation parameter, K_w, from 0 to 200 produces a relative percentage difference in the maximum transverse deflection of 74.54% for the electromechanical behavior and 88.29% for the mechanical behavior at a perforation filling ratio of 0.25. Furthermore, increasing the nondimensional foundation parameter K_p from 0 to 20 results in a relative percentage difference of the maximum transverse deflection of 73.8254% for the electromechanical behavior and 87.584% for the mechanical behavior.

✓

The perforation configuration parameters greatly affect the mechanical bending behavior. Increasing the perforation filling ratio at a constant number of holes results in a decrease in the mechanical bending deflection. The relative percentage difference in the maximum mechanical transverse deflection reaches about 199.86% due to increasing the filling ratio from 0.2 to 1 without the elastic foundation, while this percentage difference is decreased to 106.94% for the nonclassical mechanical behavior and 91.83% for the classical mechanical behavior under the point load with elastic foundation parameters K_w = 20 and K_p = 2. Furthermore, increasing the number of holes increases the mechanical bending deflection. The relative percentage difference of the maximum transverse deflection reaches 34.1643% without the elastic foundation, while this difference is decreased to 24.28% with elastic foundation constants K_w = 20 and K_p = 2.

✓

Utilizing the modified nonlocal strain gradient provides system flexibility to control stiffening and softening effects. Incorporating the nonlocal effect produces a softening effect, thus increasing the system flexibility. On the other hand, introduction of the strain gradient effect results in a stiffening effect, which leads to smaller values of the bending deflection. Due to increasing the nondimensional strain gradient parameter (l/h) from 1 to 4, the relative percentage difference of the maximum transverse deflection reaches 41.132% for the electromechanical behavior and 37.6% for the mechanical behavior at a beam aspect ratio of 10. Further increasing the nondimensional nonlocality parameter (ea/h) from 1 to 4 produces a relative percentage difference of 53.97% for the electromechanical behavior and 46.17% for the mechanical behavior.

✓

The applied loading pattern significantly affects the bending performance; larger values of the transverse bending deflection are detected with the point load compared with that obtained by the corresponding uniformly distributed loading patterns.

6. Future Work Recommendations

✓

The proposed methodology could be extended to include more nonclassical phenomena such as surface stress effects.

✓

To increase the flexibility of the developed procedure and improve the mechanical performance, a viscoelastic effect could be included in the bulk material as well as the elastic foundation.