Nonlinear dynamic stability of piezoelectric thermoelastic electromechanical resonators

This research work deals with analyzing instability and nonlinear behaviors of piezoelectric thermal nano-bridges. An adjustable thermo-elastic model with the ability to control stability conditions is developed to examine the system behavior at different temperatures. To increase the performance range and improve system characteristics, a piezovoltage is applied and a spring is connected to the sliding end of the deformable beam as design parameters. The partial differential equations (PDEs) are derived using the extended Hamilton’s principle and Galerkin decomposition is implemented to discretize the nonlinear equations, which are solved via a computational method called the step-by-step linearization method (SSLM). To improve the accuracy of the solution, the number of mode shapes and the size of voltage increments are analyzed and sufficient values are employed in the solution. The validity of the formulation and solution method is verified with experimental, analytical, and numerical data for several cases. Finally, the vibration and eigenvalue problem of the actuated nano-manipulator subjected to electrostatic and Casimir attractions are investigated. It is concluded that the fringing-fields correction changes the system frequency, static equilibrium, and pull-in characteristics significantly. The results are expected to be instrumental in the analysis, design, and operation of numerous adjustable advanced nano-systems.

Recently, micro/nano-electromechanical systems (M/NEMS) have attracted numerous attentions due to their ultra-small dimensions, favorable mechanical and electrical properties, and low power consumptions. These tiny structures represent the most promising contenders in numerous micro/nano-based applications in several fields. As a result, modelling and analysis of pull-in characteristics in the micro-scale have found significant industrial applications. However, as the pull-in instability restricts the operational range of tiny devices, obtaining system characteristics becomes important in numerous micro/nano-structures [1][2][3][4][5][6][7][8][9][10][11] .
In general, classical theories are not able to model the real behaviors of electromechanical systems in the nano-scale. That is because several effects, especially intermolecular interactions play significant roles in the submicron-scale. To predict the accurate response of different NEMS structures, several theoretical models have been developed and modified. For example, electric fringing-fields corrections (FFC) 12,13 , dispersion forces (Casimir and vdW) 14,15 , size-dependent theories (couple stress (CST), strain gradient (SGT), and so on) 16,17 , and surface layer elasticity models 18,19 have been taken into account so far. Each of these models and theories has special applications based on the miniature system properties and circumstances. It has been demonstrated that without taking the fringing corrections and dispersion effects into account, the structure will behave harder 14,20,21 . As a result, the predicted NEMS responses may not be accurate enough. Moreover, the impacts of size-dependency and surface layer energy can be a hardening or softening factor at different conditions in submicron structures [21][22][23] . The role of all these effects on the behavior and stability of nano-systems will be considered and discussed by a set of parametric studies in this paper.
In order to be able to overcome the problems of conducting conscientious experimental examinations, researchers have been competing to improve the modelling tools capable of predicting the behavior of different types of NEMS. One of the most important issues is the electrostatic force, which acts an essential role in the pull-in instability of micro/nano-systems. On the other hand, it has been demonstrated that the electric FFC is an effective nonlinearity in investigating actuated small-scale systems 24,25 . Recently, the impacts of FFC on the response of several fully-clamped nano-beams have extensively been examined 20 . Ouakad 20 found that FFC

theoretical Model
The schematic of a piezoelectric nano-beam fabricated from a movable rectangular arm and a stationary plate as the substrate is represented in Fig. 1. The deformable electrode deflects toward the substrate due to the attractions of Casimir and electric forces. Furthermore, the piezobeam is actuated by the direct current V P as the polarization voltage, which is applied along the thickness of the nano-beam and leads to a load in the longitudinal direction 33 . Such a configuration as an adjustable system is useful in several nano-bridges, tunable filters, thermal gates, cooling devices, and variable switches 31,35,36 . In this case, the boundary conditions (BCs) of the structure are such that the beam does not tolerate any traction along its neutral axis. There is not any rotation or vertical displacement at the right end.
Strain, potential, and kinetic energy. Based on the thermo-elasticity mechanics, the strain energy of a thermo-sensitive beam U t is where A, E, ΔΤ, and α t are the rectangular cross-sectional area, bulk Young's modulus of the beam or modulus of elasticity, temperature variation, and thermal expansion coefficient, respectively. Furthermore, w is the beam displacement along z coordinate axis. It should be mentioned that all the symbols used in this work are defined in the Nomenclature. The nonlinear curvature effect should be considered when the deformable beam undergoes distortion due to the external force. Therefore, we obtain the strain as xx where s and ξ xx are the actual length arched beam during deflection and beam axial strain at its neutral axis, respectively. It should be mentioned that the energy of the linear spring can be given as where U k , K s , and u are the spring potential energy, beam displacement along x coordinate axis, and linked spring stiffness, respectively. The nonlinear curvature ζ is stated as 37 where θ is the angle of arched element. By taking the nonlinear curvature of the piezobeam into account, both of axial strain ε xx and stress σ xx can be expressed as 19,38 ; all other 0 , where e 31 is the transversal piezoelectric coefficient (−6.4 C.m −2 ).
To model the electrostatic response of the piezobeam, the energy of the piezobeam U p is stated as 38 where h, L, and V P are the beam thickness, beam length, and the piezovoltage between top and bottom surfaces, respectively. Moreover, the moment of the nano-beam by considering bulk and surface layer is expressed as where b, E s , τ 0 , and ε xx s are the beam width, surface Young's modulus around the beam, surface strain around the beam (beam surface layer strain), and residual stress of the beam surface area, respectively. Therefore, the energy of the nano-beam is derived as where U m and E eff I eff are the mechanical strain energy and effective Young's modulus and second moment of area of beam (bulk and surface layer), respectively. The energy due to the surface layer tension U s is obtained as Due to the nonclassical modified couple stress theory (MCST), the energy of a beam in the submicron-scale is given by 17 where U l , ℓ, and ν are the strain energy due to material size, material nano-scale factor or size-dependent parameter according to MCST, and Poisson's ratio, respectively. According to another nonclassical SGT 17 , the energy of a nano-beam is 39 www.nature.com/scientificreports www.nature.com/scientificreports/ Furthermore, the kinetic energy of the beam is where ρ is the beam density. external work. In general, the work carried out by possible attractions can be given by ext L ext 0 where W ext and F ext are the external work and distributed external force (per unit length), respectively.
Considering two parallel plates (moveable electrode and substrate) without considering FFC, the electric attraction can be estimated as 12,25  where G, V, and ε 0 are the initial separation or gap, external applied DC voltage or electric potential difference, and electrical permittivity constant of free space vacuum (8.854 × 10 −12 F.m −1 ), respectively. Although, Parallel-Plates (PP) model for electric force is the most common configuration in numerous applications; NEMS designers are turning the focus to the electric FFC. In order to complement this effect, different corrections can be suggested to the pure electric expression. Taking into account the simplest as well as the most extensively used FFC in the literature, which is known as Palmer's (PM) model, the electrostatic attraction can be given by 12 In addition, Mejis-Fokkema (MF) model is another famous one, which can adjust the electrostatic force as follows 13 : where the effects of both beam width and thickness have been reflected. On the other hand, Casimir regime is important in nano-devices, which can be disregarded beyond this scale. Casimir force is given by 40 The studied nano-manipulator is under both Casimir F cas and electric F els forces by considering FFC, so the performed work by the mentioned forces is given by ( ) By replacing the potential, kinetic, and strain energy plus the performed work by other (corrected electrostatic as well as Casimir) forces considering the related supplementary relations into Eq. (20), using the variational approach, the following PDE for the piezoelectric size-dependent nano-beam can be obtained  The boundary conditions (BC) of the considered system are The mentioned BC are not sufficient for the strain gradient model. Accordingly, they are written as For the sake of the scalability of the problem, the nonlinear equations of motion are considered in the dimensionless form as (1 )

Solution Methodology
In general, analytical techniques are not able to solve the nonlinear differential equations of excited thermal piezo NEMS. For that reason, the governing equation of motion is better to be linearized and solved by numerical methods. In this research, a computational method called SSLM is expended to estimate threshold characteristics. In the following, the term ϖ will be considered as a linear combination of a number of modes or a set of basic (2020) 10:2982 | https://doi.org/10.1038/s41598-020-59836-0 www.nature.com/scientificreports www.nature.com/scientificreports/ functions, which are defined independently from each other. In this condition, the approximate solution of the studied structure will be constructed as where Α i is the amplitude factor of the beam transverse displacement and ϕ is the non-normalized system mode shapes vector includes a set of harmonic functions ς i , e.g. as following for MCST case 42 where terms B i signifies the displacement scaling and the values δ i are the characteristic equation roots. For SGT also different solutions have been presented in several papers 39,43 .
Having the Eq. (27) and replacing the Eq. (26) into (25), multiplying both sides by ϕ(χ), ultimately the governing equation is converted as where beam stiffness expressions K i (i = 0, 1, 2, 3, 4) and the beam inertia expression m are obtained as  Having the Eq. (28), the total stiffness of the nano-beam is expressed as voltage. In this case, the response of nano-devices is achievable by considering Eq. (33) as At the unstable conditions, the system stiffness will become singular or det(K) equals zero when the instability happens 44 . This method prepares a useful and practical way of computing the instability condition of the structures.
The nonlinear motion equations are solved numerically. In this research work, SSLM is employed 45 , and instability characteristics will be determined. Because of the nonlinear treatment of considered nano-beam, linearizing initial states may be a source of errors. To minimize the possible errors, SSLM is considered to increase the external electrostatic actuation slightly. Due to sufficient small increments of the external actuation, the beam deflection will also increase slightly. As a result, we will be able to model the nonlinear behavior of the thermal nano-switch accurately. Moreover, by accounting Taylor's expansions, accuracy necessities are obtained. Considering Eq. (33), the equation of equilibrium at the ith step is derived as (2020) 10:2982 | https://doi.org/10.1038/s41598-020-59836-0 www.nature.com/scientificreports www.nature.com/scientificreports/ As an example, the term n i is achieved from Eq. (30) as  2 where "det" refers to the determinant operator and the terms K and m are presented in Eqs. (31) and (29), respectively. Accordingly, by considering the dynamic governing equation (Eq. (28)) as well as Eq. (36), the i-th step will be derived as (1 ) (2(1 ) ) (37)

Results and Discussion
The taken step length of the solution procedure can meaningfully influence the accuracy of obtaining responses. Firstly, we will try to determine an appropriate length to increase the external electrostatic actuation in SSLM. For this purpose, responses of the prepared system in ref. 35 considering different voltage lengths are listed in Table 1.
It can be concluded that the obtained responses will be converged by taking an adequately small step length into account.
In order to make validation, we have compared the threshold voltages with the experimental 46 , analytical 26 , and numerical 27 results reported in the literature for beams with different length ( Table 2). The width and thickness of the deformable beam are 50 and 1 μm, respectively, and the initial gap is 3 μm. Furthermore, the modulus of elasticity and Poison's ratio are 169 GPa and 0.6, respectively. As it can be seen, the results obtained from our model show a good correlation with experiments validating its high accuracy.
As another validation, a fixed-sliding micro-beam in the attendance of electrostatic is modeled and pull-in characteristics are obtained and compared with a micro-structure with similar BCs (Table 3). Koochi et al. 35 examined a fixed-sliding micro-beam using the Modified Adomain decomposition and numerical solution as well as a lumped model, which assumes the uniform external forces along the beam. The available theoretical results, which are provided in Table 4, illustrate that the responses of the performed analysis are in an appropriate agreement with available data in the open literature.
The dimensions and constants employed to model the nano-system are mentioned in Table 4. More details and explanations will be discussed in the subsequent section. www.nature.com/scientificreports www.nature.com/scientificreports/ Static analysis. The relationships between the mid-beam deflections versus the voltage will be displayed in this section. Based on the figures, mid-beam deflections increase with an increase in the value of the applied voltage. In general, the actuated nano-beam is not able to endure relatively considerable voltages. Note that all reported results presented in this section (Static analysis) are according to MCST unless specified otherwise. Figure 2 shows the static defection variation with respect to different expressions of the fringing correction when the DC voltage increases from zero. Two well-known models for FFC, which are PM and MF models in addition to the electrostatic PP model for the electric actuation, are assumed to correct the distributed electric force. It is recognized that the impacts of the FFC on the system stability are remarkable. The figures display that considering the FFC makes the electrode softer and causes the decline of instability voltages. Furthermore, it can be seen that the effects of FFC are more significant by taking MF than PM model. The obtained results completely agree with the results of clamped-clamped nano-beams for different types of electric FFC reported by Ouakad 20 . He found that the responses of MF model are generally in better agreement with the finite element method (FEM) and experimental data. Figure 3 displays the deflection of the electrode, including the thermal variations and FFC (K S = 1 N/m). In this case, we consider the actuated system with (F els,PP and F els,MF ) and without (F els ) fringing effects at the temperatures of 263 K and 283 K. It can be concluded that the effects of FFC on the stability are remarkable at different temperatures or when the system operates in thermally fluctuating conditions. The effect of the linear spring on the structural behavior should be explained with respect to the sign of the thermal expansion coefficient (Table 4). It is also seen that the nano-system becomes more stable (i.e., higher pull-in voltage is obtained) with increasing the temperature. This is due to the negative coefficient of the thermal expansion (i.e., the beam starts shrinking while heated) that causes increasing the tensile strength by applying a load in the x-direction. Moreover, the curves illustrate that an increase in the ambient temperature makes the system harder and increases the contributions of electric force corrections. It means that the effect of FFC becomes greater at higher temperatures and differences between instability voltages will increase.
In Fig. 4, the impact of Casimir force with (considering MF model) and without taking FFC (assuming PP model) on the system behavior are presented. Generally, the instability condition or pull-in phenomenon happens with a delay in nano-systems without considering Casimir force. In addition, when there is no electric force, the    www.nature.com/scientificreports www.nature.com/scientificreports/ deformable electrode endures a primary displacement due to Casimir effect. On the other hand, the influence of this force on the critical applied voltage becomes more obvious in the Parallel-Plates model. Accordingly, differences between the instability voltages decrease slightly with considering electric FFC. Moreover, by using the relation π = c h cL h G E /20 cas p 2 4 3 5 , it can be understood that the impact of Casimir force enhances with an increase in the nano-beam length, unlike the thickness. This also specifies that the impact of Casimir on the system behavior is more considerable for slender electrodes.
Comparing the obtained graphs confirms that FFC is vital to be considered in evaluating the behavior of nano-switches. It can be computed in a relative error of almost 24-25% in predicting the instability parameters of the MEMS when neglecting this effect. Therefore, we have found reasonable agreement between the responses of this dimensionless nano-system and the numerically calculated values for arch MEMS in Ref. 20 . Figure 5 illustrates the variations of the mid-beam at different amounts of spring stiffness vs. the external voltage. Here, as the spring stiffness tends to zero, there exists no thermal force in this configuration. Therefore, the thermal effect is meaningless without considering the spring in such beams with the fixed-sliding BCs (since a fixed-sliding beam like a cantilever one is capable of moving in the longitudinal axis, the temperature variation results in changing its length and it will be free of stress). On the other hand, as the spring stiffness tends to ∞, the responses become close to those for a clamped-clamped block. The influence of spring stiffness is more considerable at lower temperatures due to the difference between the pull-in voltages of blue curves, which should be explained regarding the thermal expansion coefficient. Finally, as a practical point, it is notable that the influence of temperature variations increases while the stiffness enhances, so the difference between the pull-in voltages becomes greater in devices with doubly clamped BCs. Figure 6 shows the relationship between the curves and the excitation including/excluding the effects of FFC using different size theories. To investigate the impact of higher-order material size parameters, the achieved results have been compared using MCST and SGT. First, note that in this research, the higher-order bending rigidity of material is assumed 21 nm. As a result, the modified couple stress length-scale parameter is 15 nm 17,47 . Assuming three equal strain gradient length-scale parameters 39 , these parameters will be equal to 9 nm 43,48 . The  www.nature.com/scientificreports www.nature.com/scientificreports/ results demonstrate that both MCST and SGT make the beam to be harder and cause an increase in the threshold voltage.
In Fig. 7, it can be found that the piezovoltage affects the system response significantly by changing the structural behavior. Based on the application, the piezovoltage can compress or stretch the beam by considering the polarity. The results show that the positive actuation increases the instability voltage by increasing the beam   www.nature.com/scientificreports www.nature.com/scientificreports/ stiffness where the negative polarity produces a compression force in the x-direction. It is comprehended that piezoexcitations shift the equilibrium manifold and consequently the instability voltage. The results agree well with the reported ones in clamped-clamped piezobased bridges 32,33 .
The relationships between the mid-beam vs. the electrostatic excitation by considering different FFC models with both positive and negative piezoexcitation are shown in Fig. 8. It can be observed that the impact of piezoelectricity changes by taking different models to consider FFC. Generally, as the piezoelectric voltage increases, the difference between the results of FFC models becomes more considerable. Consequently, differences between the instability voltages in FFC models are increased by considering the positive piezoelectric excitation. On the other hand, differences between the results of piezovoltage will increase as the negative polarity considered.
It is interesting to note that by tuning the value of piezoexcitation, the bridge behavior is changeable, so the pull-in voltage is controllable. As a result, we can increase the stability range or decrease the critical voltage by taking positive or negative piezoelectric actuation. Therefore, improving the performance of several nano-switches are possible by considering piezomaterials and applying the piezovoltage as an important design parameter in numerous applications.

Dynamic analysis.
Here, the vibrations of the rectangular electromechanical nano-system will be analyzed.
It should be noted that for more accurate investigations of each parameter, these assumptions are taken into account unless specified otherwise (T = 273 K, λ = 0, η = 0, ι = 0, ι i = 0, V P = 0, φ = 0, and c cas = 0). In these cases, the system frequency decreases with an increase in the external voltage at first. After that, the system response tends to 0 promptly, when the DC voltage goes to a maximum amount (near the critical value). Figure 9 shows the relationships between the mid-beam and the voltage parameter by considering all PP, PM, and MF models. It can be concluded, taking FFC affects pull-in characteristics more than system frequencies. In addition, the impacts of the electrostatic force corrections on the nano-system frequency become more obvious with increasing the voltage gradually. The responses for system frequencies are qualitatively agreeable by accounting numerical results of Ouakad 20 , including the electric fringing-fields effect (all models). Finally, note that without applying the DC voltage, there is no difference between the response frequencies with/without  www.nature.com/scientificreports www.nature.com/scientificreports/ accounting FFC (Fig. 9). It means that the start points of curves are the same in this figure unlike other figures of dynamic behavior.
The impact of the residual tension in the surface layer on the nano-system frequency is shown in Fig. 10. It is comprehended that system frequencies and instability voltages change with variations of the tension. Here, response frequencies and critical voltages increase by considering positive surface stress. On the other hand, the values for the residual surface stress nano-materials may be negative. In this case, the results are exactly the opposite of the positive one, so considering the negative stress makes the fixed-sliding beams soften. Moreover, using the relation λ , it is recognized the influence of residual surface tension reduces with decreasing the nano-beam length and increasing the thickness. Therefore, this molecular surface effect will be more dominant in the slender nano-structures. It is necessary to mention that the effect of FFC becomes more significant by considering positive term; however, the difference between the positive and negative stress decreases, including FFC. Figure 11 displays the effect of different nonclassical size theories (MCST and SGT) on the dynamic behavior and response frequencies of the nano-beam. The curves illustrate that by considering the length-scale parameters, we have a stiffer nano-system. Therefore, both system frequencies and dynamic voltage will be increased by taking nonclassical theories into account. Moreover, differences between the results of SGT and classical theory is more considerable. Ultimately, from relations of material size parameters, it is found that the influence of size parameters becomes weak by increasing the beam thickness.
The obtained results of this research allow a better understanding of the instability of the nano-beam based tiny devices and can guide scientists to improve their general performance, appropriately. Moreover, the outcomes of this study explain the causes of differences obtained in the literature when comparing the responses of the nano-system without considering modified models including nonclassical theories, small-scale effects, and external force corrections to experimental data. There is still plenty of exploration to be discovered regarding the reliability of electro-thermo-mechanical nano-systems. Forthcoming investigations need to be established, including modelling the impacts of out-of-plane electrodes arrangement on the behavior of nano-devices with different BCs by considering both size and surface effects, which should be verified experimentally.  www.nature.com/scientificreports www.nature.com/scientificreports/ conclusion In order to investigate the pull-in phenomenon of adjustable intelligent systems, a nano-beam with a sliding end connected to spring was taken into account and a thermomechanical model was presented. In the model, the impacts of the surface layer, molecular forces, size-dependency, piezoelectricity, and temperature variations were included with consideration of the nonlinear curvature. Two famous models for the electric fringing corrections, which are PM and MF models, were assumed to correct the distributed electrostatic force. The PDEs were derived by the extended Hamilton's principle. Galerkin decomposition was implemented to discretize the nonlinear equations. To improve the accuracy of responses by spending affordable computational calculations, the number of mode shapes and the size of voltage increments were optimized. Afterward, the structural behavior of the nano-system under Casimir and corrected electric attractions were analyzed numerically. Finally, the system vibrations with regard to the environment temperature, surface layer, and material size theories (MCST and SGT) were studied.
The obtained results reveal that the classical theory predicts lower values for instability voltages and system frequencies. It was also shown that considering the electric FFC is essential for accurate modelling of the nano-switches behavior. Moreover, the effect of MF model is more remarkable than the other one. In addition, this effect causes a significant change at higher temperatures as well as in cases of positive residual surface stress and positive piezoelectric voltage. As a practical point, the effect of temperature variation increases while the stiffness of the connected spring enhances, so differences between the critical voltage at different temperatures becomes greater in clamped-clamped beams. Finally, it should be mentioned that the piezoexcitation could be employed as a design parameter to improve the controllability and increase the performance range of adjustable smart devices.