## 3 Solution Procedure
### Natural Frequency Analysis  
Firstly, the free vibration analysis is performed to obtain the natural frequency of the structure. Setting the base excitation \(w_{b}\), damping coefficient \(c_{v}\) and voltage potential \(V\) to be zero, neglecting the nonlinear terms, and assuming harmonic excitation, the lateral displacement is assumed in the form at the prebuckling state:  
\[w=W(t)\sin\Big{(}\frac{\pi xx}{L}\Big{)}\text{, prebuckling} \tag{39}\]By substituting Equation (39) into Equation (29), applying the [PERSON] method, one obtains:  
\[\begin{array}{c}-L^{2}\pi^{2}P\Big{[}L^{2}+(\epsilon_{0}a)^{2}\pi^{2}\Big{]}We^{ i\omega_{n}t}+\pi^{4}bB_{11}\big{(}L^{2}+l^{2}\pi^{2}\big{)}We^{i\omega_{n}t} \\ +L^{4}\Big{\{}I_{0}\Big{[}L^{2}+(\epsilon_{0}a)^{2}\pi^{2}\Big{]}+\frac{4m_{t} L}{\pi}\Big{\}}(j\omega_{n})^{2}We^{i\omega_{n}t}=0\end{array} \tag{40}\]  
Equation (40) can be written as:  
\[\ddot{W}+\omega_{pre}^{2}W=0 \tag{41}\]  
where \(\omega_{pre}\) is so-called linear natural frequency at prebuckling state and:  
\[\omega_{\rm pre}=\sqrt{\frac{\pi^{5}bB_{11}(\pi^{2}l^{2}+L^{2})-\pi^{3}L^{2}P \Big{[}\pi^{2}(\epsilon_{0}a)^{2}+L^{2}\Big{]}}{L^{4}\Big{\{}\pi l_{0}\Big{[} \pi^{2}(\epsilon_{0}a)^{2}+L^{2}\Big{]}+4 Lm_{s}\Big{\}}}},\,\,\text{prebuckling} \tag{42}\]  
It is seen in Equation (42) that the natural frequency is related to the compressive load \(P\).  
When the beam is buckled, the lateral displacement is in the from:  
\[w=W_{s}\sin\Big{(}\frac{\pi x}{L}\Big{)}+W(t)e^{i\omega_{n}t}\sin\Big{(}\frac{ \pi x}{L}\Big{)},\,\,\text{postbuckling} \tag{43}\]  
where \(W_{s}\sin(\frac{\pi x}{L})\) is the static vertical displacement, and the expression of \(W_{s}\) is provided in Equation (37) By substituting Equation (43) into Equation (29), neglecting the nonlinear terms, and applying the [PERSON] method, one obtains:  
\[\begin{array}{c}L^{2}\pi^{2}P\Big{[}L^{2}+(\epsilon_{0}a)^{2}\pi^{2}\Big{]} We^{i\omega_{n}t}-\pi^{4}bB_{11}\big{(}L^{2}+l^{2}\pi^{2}\big{)}We^{i\omega_{n}t} \\ +\frac{l^{4}}{2}\Big{\{}I_{0}\Big{[}L^{2}+(\epsilon_{0}a)^{2}\pi^{2}\Big{]}+ \frac{4m_{t}L}{\pi}\Big{\}}(j\omega_{n})^{2}We^{i\omega_{n}t}=0\end{array} \tag{44}\]  
Equation (44) can be written as:  
\[\ddot{W}+\omega_{post}^{2}W=0 \tag{45}\]  
where \(\omega_{post}\) is the natural frequency at postbuckling state and:  
\[\omega_{\rm post}=\sqrt{2}\sqrt{\frac{\pi^{3}L^{2}P\Big{(}\pi^{2}(\epsilon_{0 }a)^{2}+L^{2}\Big{)}-\pi^{5}bB_{11}(\pi^{2}l^{2}+L^{2})}{L^{4}\Big{[}\pi l_{0} \Big{(}\pi^{2}(\epsilon_{0}a)^{2}+L^{2}\Big{)}+4 Lm_{s}\Big{]}}},\,\,\text{ postbuckling} \tag{46}\]  
In Equation (46) the porosity related coefficient \(B_{11}\) is shown in Equation (36), and the other porosity related coefficient \(I_{0}\) is calculated as:  
\[I_{0}=b\int_{-h_{t}/2}^{h_{t}/2}\rho(z)dz=\left\{\begin{array}{l}bh_{t}\rho_ {\rm max}\sqrt{1-\chi_{0}\alpha}\,\,\,\text{for distribution 1}\\ bht_{t}\rho_{\rm max}\Big{[}1-\frac{2h_{c}}{\pi h_{t}}\chi_{m}\sin\Big{(}\frac {\pi h_{t}}{2h_{c}}\Big{)}\Big{]}\,\,\,\text{for distribution 2}\\ bht_{t}\rho_{\rm max}\Big{[}1-\frac{2\chi_{2}h_{t}}{\pi h_{t}}\chi_{m}\sin \Big{(}\frac{\pi h_{t}}{4h_{c}}\Big{)}\Big{]}\,\,\,\text{for distribution 3}\end{array}\right. \tag{47}\]