## 2 Mathematical Formulation
Considering the NSGT with flexoelectric effect, the constitutive equations for the nanobeam are written as [18]:  
\[\Big{[}1-(e_{0}a)^{2}\
abla^{2}\Big{]}t_{xx}=C_{11}(1-l^{2}\
abla^{2}) \varepsilon_{xx}-e_{311}E_{z} \tag{7}\]  
\[\Big{[}1-(e_{0}a)^{2}\
abla^{2}\Big{]}\sigma_{xxz}=-\mu_{3113}E_{z} \tag{8}\]  
Figure 1:0  
where \(e_{0}a\) and \(l\) are the nonlocal parameter and the strain gradient constant, respectively. \(e_{311}\), \(\mu_{3113}\) and \(a_{33}\) are the piezoelectric coefficient, the electric field-strain gradient coupling coefficient and the dielectric constant, respectively. \(D_{z}\) is the electric displacement. The variation in the electric entropy can be calculated as [18]:  
Figure 1:1  
By substituting Equation (6) into Equation (10), and performing the calculation, one obtains:  
\[ \begin{array}{l}\delta\iiint_{V}Hd\Omega=\iiint_{V}\big{(}t_{xx}\delta \varepsilon_{xx}-D_{z}\delta E_{z}+\sigma_{xxz}\delta\varepsilon_{xx,z}\big{)} d\Omega\\ =\int_{0}^{L}\Big{(}N_{xx}\delta\frac{\partial u}{\partial x}+N_{xx}\frac{ \partial w}{\partial x}\delta\frac{\partial w}{\partial x}-M_{xx}\delta \frac{\partial^{2}w}{\partial x^{2}}-G_{xx}\delta\frac{\partial^{2}w}{ \partial x^{2}}+\frac{B_{xx}}{h_{p}}\delta V\Big{)}dx\\ =\int_{0}^{L}\Big{\{}-N_{xx,x}\delta u+\frac{B_{xx}}{h_{p}}\delta V-\left[ \frac{\partial}{\partial x}\left(N_{xx}\frac{\partial w}{\partial x}\right) +M_{xx,xx}+G_{xx,xx}\right]\delta w\Big{\}}dx\end{array}  \tag{11}\]  
where  
\[ \begin{array}{l}N_{xx}=b\int_{-h_{l}/2}^{ht/2}t_{xx}dz,\,M_{xx}=b\int_{-h_{l} /2}^{ht/2}zt_{xx}dz,\\ G_{xx}=b\int_{-h_{l}/2}^{ht_{l}/2}\sigma_{xxz}dz,\,B_{xx}=b\int_{-h_{l}/2}^{ht _{l}/2}D_{z}dz\end{array}  \tag{12}\]  
By combining Equation (12) with Equations (7)-(9), one has:  
\[ \begin{array}{l}\Big{[}1-(e_{0}a)^{2}\
abla^{2}\Big{]}N_{xx}=bA_{11}\big{(}1 -l^{2}\
abla^{2}\big{)}\Big{[}\frac{\partial u}{\partial x}+\frac{1}{2}\Big{(} \frac{\partial w}{\partial x}\Big{)}^{2}\Big{]}+2 be_{311}V\\ \Big{[}1-(e_{0}a)^{2}\
abla^{2}\Big{]}M_{xx}=-bB_{11}\big{(}1-l^{2}\
abla^{2} \big{)}\frac{\partial^{2}w}{\partial x^{2}}\\ \Big{[}1-(e_{0}a)^{2}\
abla^{2}\Big{]}G_{xx}=2b\mu_{3113}V\\ \Big{[}1-(e_{0}a)^{2}\
abla^{2}\Big{]}B_{xx}=-2 ba_{33}V+2 bh_{p}\bigg{\{}e_{311 }\bigg{[}\frac{\partial u}{\partial x}+\frac{1}{2}\Big{(}\frac{\partial w}{ \partial x}\Big{)}^{2}\bigg{]}-\mu_{3113}\frac{\partial^{2}w}{\partial x^{2 }}\bigg{\}}\end{array}  \tag{13}\]  
where \(A_{11}=\int_{-h_{l}/2}^{h_{l}/2}C_{11}dz\), \(B_{11}=\int_{-h_{l}/2}^{h_{l}/2}z^{2}C_{11}dz\). Neglecting the rotational inertia of the beam, the total kinetic energy of the beam is written as:  
\[ \begin{array}{l}T=\frac{1}{2}\iiint_{V}\rho(z)\Big{[}\dot{u}^{2}+\big{(}\dot{w} +\dot{w}_{b}\big{)}^{2}\Big{]}d\Omega+\frac{1}{2}m_{s}\big{[}\dot{w}(L/2)+\dot{w }_{b}\big{]}^{2}\\ =\frac{1}{2}\int_{0}^{L}l_{0}\Big{[}\dot{u}^{2}+\dot{w}^{2}+\dot{w}_{b}^{2}+2 \dot{w}\dot{w}_{b}\Big{]}dx+\frac{1}{2}m_{s}\big{[}\dot{w}(L/2)+\dot{w}_{b}\big{]} ^{2}\end{array}  \tag{14}\]  
where \(l_{0}=b\int_{-h_{0}/2}^{h_{0}/2}\rho(z)dz\). The variation in the total work performed by the compressive axial load \(P\), damping force and the electric charge output \(Q(t)\) is expressed as:  
\[\delta W=P\int_{0}^{L}\delta\Bigg{[}\frac{\partial u}{\partial x}+\frac{1}{2} \bigg{(}\frac{\partial w}{\partial x}\bigg{)}^{2}\Bigg{]}dx-\int_{0}^{L}c_{v }\dot{w}\delta wdx+2Q(t)\delta V(t) \tag{15}\]  
where \(c_{v}\) is the viscous damping coefficient. The Hamilton principle is utilized as:  
\[\delta\int_{t_{1}}^{t_{2}}\Big{(}T-\iiint_{V}Hd\Omega+W\Big{)}dt=0 \tag{16}\]  
By substituting Equations (11)-(15) into Equation (16), and performing the calculation, one obtains:  
\[