Thermoelastic vibrations of a Timoshenko microbeam based on the modified couple stress theory

The dependence of the quality factor of nonlinear microbeam resonators under thermoelastic damping for Timoshenko beams with regard to geometric nonlinearity has been studied. The constructed mathematical model is based on the modified couple stress theory which implies prediction of size-dependent effects in microbeam resonators. The Hamilton principle has yielded coupled nonlinear thermoelastic PDEs governing dynamics of the Timoshenko microbeams for both plane stresses and plane deformations. Nonlinear thermoelastic vibrations are studied analytically and numerically and quality factors of the resonators versus geometric and material microbeam properties are estimated. Results are presented for gold microbeams for different ambient temperatures and different beam thicknesses, and they are compared with results yielded by the classical theory of elasticity in linear/nonlinear cases.


Introduction
Thermoelastic damping/dissipation (TED) occurs in all elastic materials subjected to cycling deformations, and particularly in the case when the period of an exciting cycle is close to the material relaxation [1][2][3][4]. For instance, in the case of vibrations of an elastic beam major part of mechanical work is transformed into elastic energy, but part of the work is converted into thermal energy. During bending one side of the structure is extended and cooled, while the other side is compressed and heated. Temperature non-homogeneity described so far means the transfer of heat from the more heated part to the more cooled part of the structure. According to the second law of thermodynamics, the latter heat transfer causes the increase in entropy and finally the loss of mechanical energy in the form of heat. There are various mechanisms of energy loss, such as thermoelastic damping [1,[5][6][7], energy loss in supports [8], internal energy loss [9] and air damping [10].
In almost all macroscopic systems with basic damping sources, thermoelastic coupling plays a secondary role because it is an order of magnitude lower than the damping caused by external and internal viscous effects [11]. However, in the case of MEMS/NEMS devices and resonators, it plays a crucial role. It should be emphasized that micromechanical resonators with high quality factors are employed in a wide range of MEMS applications, such as sensors [12][13][14][15][16], electric filters [17,18], accelerometers [19], giros [7], atomic force microscopes [20,21], MEMS type switches, and various other mechanisms [22][23][24]. From a technical point of view, the evaluation of thermoelastic damping (TED) of those structural elements at the design stage plays an important role in maintaining the required engineering characteristics. TED belongs to fundamental sources of eigen damping in the microelectromechanical systems (MEMS) [25] and nanoelectromechanical systems (NEMS) [3,26] working in the vacuum regime. Evoy et al. [27] and Duwel et al. [7] have shown experimentally that TED is the dominant source of damping in MEMS and NEMS devices.
Zener [1,28,29] belongs to the first who defined the occurrence of TED and showed its crucial role in dissipation exhibited by bending resonators. Landau and Lifshitz [30] estimated analytically the TED coefficients of elastic vibrations. The second widely employed analytical model has been developed by Lifshitz and Roukes [2]. In that model, a change in the resonance frequency associated with TED was observed, and thus, the classical Zener model was improved. The mentioned authors have also derived an analytical expression for quality factor Q in microbeam resonators and studied the influence of numerous geometric parameters on their quality factors. Nayfeh and Younis [31] constructed the analytical expression for Q of a microplate of arbitrary forms for different boundary conditions and studied the occurrence of TED. Sun et al. [32][33][34], Sharma et al. [35] and Rezazadeh et al. [36] investigated thermoelastic damping in capacity microbeam resonators using the hyperbolic heat trans-fer model. Besides, Vahdat and Rezazadeh [25] investigated the effect of axial and residual stresses on thermoelastic damping in a capacity microbeam resonator.
Berry [37] obtained experimentally Q factor as a function of frequency. Duwel et al. [7] presented experimental data illustrating the importance of TED in the resonance MEMS sensors. They also investigated the influence of various materials on TED behaviour. Zhang and Turner [37] developed a theory for TED for micro-and sub-micro longitudinally vibrating beams. Vengallatore [39] studied TED in laminated composite micromechanical beam resonators. Prabhakar and Vengallatore [40] constructed an exact theory useful for TED computation in asymmetric two-layer micromechanical beam resonators. Yi [41] investigated the effect of geometric parameters on TED in the MEMS resonators. Analytical methods to study TED in nanobeams are presented in reference [42].
On the other hand, most of the models described in the literature use Euler-Bernoulli's beam theory which does not apply to a small ratio of beam length to its thickness. Massalas and Kalpakidis [43] investigated coupled thermoelastic vibrations of simply supported beams. They presented analytical solutions for the coupled problem of thermoelasticity of beams made from homogeneous and isotropic materials applicable for the Euler-Bernoulli and Timoshenko beam models. Shieh [44] studied TED of a simply supported Timoshenko beam with a circular cross section. He used 2D heat transfer equations and obtained damping coefficients. In reference [45], the authors proposed an analytical solution for the resonator quality factor under thermoelastic damping in beams based on the Timoshenko theory using the Lifshitz and Roukes method [2] for the Euler-Bernoulli beams. Further investigations could be carried out as part of generalized beam theories [46,47].
There are known attempts to solve the problems of thermoelastic damping numerically. Prevost and Tao [48] employed FEM (Finite Element Method) to solve a coupled thermoelastic problem. The authors used the second-order heat transfer PDE instead of the Fourier equation. Prabhakar and Vengallatore [49] matched 1D Euler-Bernoulli beam theory and 2D heat transfer PDE and obtained the quality factor using the introduced theory of series. Silver and Peterson [50] investigated elastic thermodynamic damping for the beam using FEM and also employed the method of perturbation. Sun et al. [51] obtained TED of microbeam resonators using the method of finite Fourier transform and modal analysis. Lepage [52] applied FEM to estimate the quality factor of TED of the Euler-Bernoulli beams. He employed the cubic approximation of temperature variation along beam thickness. Guo et al. [53] used 2D FEM to study MEMS resonators. De and Aluru [54] employed a combined method of finite and boundary elements to validate the results of their TED-related modified theory.
In MEMS structures made of metals and polymers in which the smallest length of elements was in the microrange, the size-dependent phenomena were observed experimentally in many systems [55][56][57][58][59]. Effects of the dimensional dependence can be divided into two categories. The first one is implied by the boundary effect coupled with a molecular layer located on the specimen surfaces. Brezny and Green [60] and Onck et al. [61] observed that the effective Young modulus and resistance against compression of some materials essentially decreased while decreasing the size of the specimen. The second one known as the Kosser effect, or micropolar effect [62], occurs due to violation of the classical mechanical statements. Owing to the lack of eigen size of length, the classical mechanical theories cannot give proper interpretation and prediction of the observed size-dependent phenomena.
More recently theories of continuum of a higher order for prediction of the mentioned size-dependent relations were developed [63]. In the 1960s, some of researchers proposed a couple stress theory of elasticity which can be classified as the non-classical theory [64]. The theory made it possible to interpret the sizedependent effects using two higher-order material constants in state equations. A simplification introduced by Yang et al. [65] to the scale complex relations of the couple stress theory of elasticity resulted in a modified coupled stress theory which is suitable because of effects using only one scale parameter of the material length. Recently, many researchers investigated nonclassical theorems of continuum in order to formulate and study the size-dependent mechanical behaviour of beams and plates.
Tsiatas [58] proposed a novel model of Kirchhoff plates based on a modified couple stress theory (MCST) which can take into account the size-dependent effects. Park et al. [66] proposed the Euler-Bernoulli beam model based on MCST and they investigated the influence of the length material parameter on static mechanical microbeam properties.
Kong et al. [67,68] solved the dynamical problem of Euler-Bernoulli beams based on the MCST and the gradient theory of elasticity, and they demonstrated that eigenfrequencies of the microbeams obtained with the help of MCST are higher than those predicted by the classical theory of Euler-Bernoulli beams. Using MCST, Kahrobaiyan et al. [69] and Asghari et al. [70] studied analytically the static and dynamic behaviour of functionally graded microbeams and material of atomic force microscope microcantilevers.
Ma et al. [71] used Hamilton principle and relations of the non-local Eringer's theory to develop non-local theories of the Euler-Bernoulli, Timoshenko, Reddy and Levinson beams.
Besides, Wang et al. [63] proposed a microscale Timoshenko beam model based on the theory of gradient deformation. Using the Hamilton principle, they derived the governing equations, initial and boundary conditions. The model took into account Poisson's effects, including three length material parameters, and consequently, it exhibited size effects.
In recent years, a series of works have been devoted to the study of MCST of beam and plates also functionally graded taking into account shear deformations [36,[72][73][74][75][76][77]. Reddy [73] and Reddy and Kim [74] derived equations of beams and plates based on MCST.
However, in general, there are not many works reported in the existing literature devoted to the study of size-dependent effect in the problems of coupled thermoelasticity of microbeams. Guo and Rogerson [78] investigated the influence of the size on thermoelastic coupling in clamped elastic prismatic beams. Sun et al. [51] developed a theory for the coupled thermoelastic microbeam problem using the hyperbolic model of heat transfer. They demonstrated how the size effect generated by the thermoelastic coupling disappears when thickness of the microbeam exceeds its critical value depending on the material properties and boundary conditions. Rezazadeh et al. [79] derived analytical solutions to estimate the quality factor of TED obtained on the basis of MCST for the plane stress and plane deformation states of a linear Euler-Bernoulli beam. They found that when the beam thickness is close to the material length parameter, the results obtained by MCST are different compared to the results obtained using the classical beam theory.
In order to explain size dependence for microbeams, Taati et al. [80] combined the MCST and hyperbolic model of heat transfer which differs from the classical Fourier model. They showed that if one uses the hyperbolic heat transfer model, then MCST shows more deviations from results compared to the classical theory.
In reference [81], the size-dependent model was proposed using the time of thermal relaxation for temporal analysis of thermoelasticity of Timoshenko microbeams. Besides, the size-dependent effects of a simply supported beam were studied under constant impulse force where two beam ends were exposed to constant ambient temperature.
Since eigenfrequencies of microbeam vibrations change when accounting for MCST, it may be expected that TED forecast by MCST will also differ from the results yielded by the classical theory of beams. Therefore, the influence of the size length material parameter on quality factor Q of microbeam resonators still requires a deeper analysis (in particular, when the characteristic size is compared with the internal size parameter of material length).
The microbeam resonators often work in nonlinear regimes with large amplitude while harvesting energy [82]. A key role in the studies of beam resonators plays energy dissipation in micro-and nanoscales [1,3,4,83,84]. In particular, TED belongs to one of the principal mechanisms of energy dissipation for large amount of the micro/nanomaterial resonators [2][3][4]85,86]. However, as it has already been mentioned, almost all published results deal with linearized regimes of vibrations of the microbeam resonators taking into account only vibrations with small amplitude.
Nonlinearity in micro-resonators may occur due to various reasons including large deflections (geometric nonlinearity or material nonlinearity). Geometric nonlinearity plays a key role in the microbeam resonator [87][88][89]. Though it was observed in high-frequency MEMS/NEMS resonators [82,[90][91][92], to our knowledge there are no reference studies on large deflection impact on thermoelastic dissipation. There are a few papers devoted to the analysis and quantitative estimation of the influence of thermoelastic damping in MEMS. Majority of the mentioned works deal with small amplitude linear vibrations [31,33,93]. For instance, in the work [86] though large static deflection generated by electrostatic forces under the action of constant stress is studied, thermoelastic investigation is based on the linearized vibration of small amplitude in the neighbourhood of static equilibrium with large deflection. An analogous approach can be found in ref-erence [54] where large static deflection generated by electrostatic input affects the thermoelastic damping of linear small amplitude vibrations around static equilibrium. A series of analogous examples are reported in paper [42], where TEDs of resonators made from carbon tubes and vibrating around equilibrium were studied. Despite the thermoelastic dissipation analysis within linear theory, there are several studies devoted to nonlinear dynamics under large deflections of beam MEMS resonators [20,21,94,95]. Méndez et al. [95] studied large deflection nonlinearity at the resonant frequency localization and investigated the velocity of damping of cantilever microbeams but omitted the effect of large amplitude vibration on TED. In reference [11], the investigated thermoelastic damping with an account of nonlinear effects of the circle plate modelled with the use of the von Kármán hypothesis was analysed. Employing the Kantorovich method and the methods of perturbation, a formula describing the resonator quality factor was obtained.
The authors of reference [96] studied nonlinear effects under large deflections and under thermoelastic dissipation of the microbeam resonators for adiabatic/isothermic beam surfaces. The models of thermoelasticity and thermoelastic dissipation were formulated for the case of large amplitude vibrations for the Euler-Bernoulli beam without the size-dependent effects. Mendez et al. [95] carried out numerical modelling for MEMS cantilever, and a difference between linear and nonlinear approximations was illustrated and discussed. It was shown that for large displacements, the beam becomes stiffer (frequency increase) and its vibrations are damped faster. Those observations made it possible to improve fabrication of high-frequency MEMS resonators. In the work [97], the importance of nonlinear MEMS behaviour with a strong coupling between mechanical and electric fields was demonstrated. The variational procedure we used without the help of entropy inequality has been successfully validated by recently published papers that allow us to get clues to understand the nonlinear/chaotic behaviour of structural members including size-dependent effects. Namely, in paper [98], regular and chaotic vibrations of flexible curvilinear beams with (and without) the action of temperature and electric fields are studied. There are numerous results regarding the transition scenarios from regular to chaotic vibrations including the occurrence of hyper-hyper chaos and deep chaos. In reference [99], a mathematical model of two-layer beams coupled by boundary conditions in a stationary temperature field with regard to geometric nonlinearity is studied. Phase synchronization of beam vibrations versus both character and intensity of the applied temperature field is investigated. A mathematical model describing the nonlinear dynamics of a plate-beam system is proposed in [100]. The model takes account of coupling between temperature and deformation fields as well as external mechanical, temperature, and noise excitations. Both the proof of existence and uniqueness of the solution to the problem and the proof of convergence of the Faedo-Galerkin method used to solve the problem are given. In the work by Krysko et al. [101], a mathematical model of flexible shallow nano-shells is considered based on the modified couple stress theory in a higher approximation. The Fourier power spectra, various wavelets-type spectra, phase and modal portraits, as well as signs of the LEs (Lyapunov exponents) are investigated.
The overview of the state of the art of the MEMS/ NEMS beam vibrations carried out so far yields the following conclusions: (i) thermoelastic damping of microbeam vibrations with large amplitudes plays a crucial role in a proper fabrication of the resonators; (ii) the influence of geometric nonlinearity on the values of frequency, damping and on quality factors of the resonators is often omitted in the available literature; (iii) both thermoelastic damping and large displacements may have an important effect not only on the high-frequency resonators but also on the MEMS/NEMS devices [40,82,91,92]; (iv) in general, the thermoelastic PDEs governing the dynamics of microbeams were obtained based on the classical theory of continuum. There are no sufficient results following from the nonlinear model of coupled thermoelasticity of the Timoshenko microbeams based on the MCST.
Our work is aimed at filling the gaps that exist in the research of MEMS/NEMS resonators. In this work, we employ the MCST to study quality factor of microbeam resonators with an account of thermoelastic damping and geometric nonlinearity. The associated mathematical model is constructed which, in contrary to the classical theory, consists of the internal scale length material parameter allowing for a reliable prediction of the sizedependent effect in microbeam resonators. The Hamil-ton principle yields coupled nonlinear thermoelastic equations of motion of the Timoshenko microbeams in the case of plane stress/deflection state. Based on the obtained equations, both analytical and numerical investigations of the nonlinear thermoelastic vibrations of the microbeams are carried out and then, the values of quality factors of the resonators versus geometric and material microbeam properties are estimated. The results are presented for gold microbeams at different ambient temperatures and beam thicknesses which are compared with the results obtained based on the classical theory of elasticity in both linear and nonlinear cases.
It should be mentioned that some authors of this paper have earlier investigated the influence of temperature field on the statics and dynamics of nonlinear Euler-Bernoulli nanobeams [102] fabricated based on the optimization method [103].
The paper is organized in the following way. Formulation of the problem with reference to the stress and deformation fields, plane deformation conditions and the governing beam dynamics equations within a coupled problem of thermoelasticity with regard to the Timoshenko beams is presented in Sect. 2. Section 3 refers to the quality factor of resonators with thermoelastic damping. Numerical and analytical results for a rectangular cross section beam resonator are given in Sect. 4. The obtained results as well as novel features of the paper are briefly described in Sect. 5.

Formulation of the problem
We introduce a rectangular system of coordinates to study a microbeam as shown in Fig. 1. The origin of coordinates is fixed on the left microbeam end. The microbeam has length L, thickness h, width b, and its transversal cross section parallel to the Y O Z plane is denoted by A (see Fig. 1b).

Stress and deformation fields
Heat deformations occur in an elastic body due to a linear heat extension coefficient. Therefore, the total deformation composed of mechanical and heat components is as follows [79] (1) By θ = T − T 0 , we denote the temperature change with respect to the initial temperature T 0 , and we assume that the shear deformations generated by the temperature change can be neglected. Then, in the case of a linear, isotropic and homogeneous body, heat (T ) and mechanical (M) deformations are governed by the following formula [100] where α temperature coefficient of a linear extension, σ i j stress tensor, δ i j Kronecker's symbol, E Young modulus, and ν is the Poisson coefficient. Owing to (1) and (2), the total deformation field is as follows

(i) Plane stress state condition
In the case when beam thickness (in direction of the O Z axis) and its width (in direction of the OY axis) is sufficiently small in comparison to the beam length (in direction of the O X axis), then assumption of the plane stress state implies that all components of the stress tensor in directions of the axes OY and O Z are equal zero (σ 12 = σ 22 = σ 32 = 0, σ 23 = σ 33 = 0). However, in the case of the Timoshenko beams we have σ 13 = 0 [100]. In the latter case, components of the deformations tensor can be recast to the following simple form A stress tensor under the plane stress state assumption has only two components different from zero, which is expressed by the components of the deformation tensor in the following way where k s is correction coefficient introduced due to the assumption of non-homogeneous shear deformation with respect to the transverse beam cross section and it depends on the forms of the beam cross sections, G is the shear modulus.
Owing to the known Cauchy-Green deformations the kinematic relations of the theory of Timoshenko beams can be presented in the following form where u (x, t) , w (x, t) and ψ (x, t) are the axial displacement of the beam middle line, transversal beam deflection and the angle of rotation of the transversal load with respect to the vertical direction, respectively. However, in this work, we consider the geometrically nonlinear Kármán model, i.e. beam with small deformations and rotations, but with relatively large displacements w that cause the occurrence of geometric nonlinearity. Parameters ψ, ∂u/∂ x, ∂ψ/∂x under the mentioned conditions are very small. Employing the Timoshenko kinematic hypothesis and taking into account (4)- (7), the nonzero deformation and stress components can be defined via the displacements field in the following way is sufficiently small in comparison with its length (in direction of the O X axis), but beam width (in direction of the OY axis) is sufficiently large, then based on the assumption of the plane deformation state one may conclude that components of the stress tensor in direction of the axis O Z and deformation components in direction of the axis OY are equal zero (σ 23 = σ 33 = 0, ε 12 = ε 22 = ε 32 = 0). The remaining nonzero components of the deformation and stress tensors can be expressed by the field of displacements (7) and temperature in the following way

Beam equation of motion based on the modified couple stress theory
In the modified couple stress theory [65], the accumulated energy of deformation U in a linear elastic body occupying space Ω is defined as follows where V is the volume of elastic matter; m i j are the components of a deviator part of the symmetric tensor momentum of higher order and χ i j are the components of symmetric curvature tensor [56], which are defined in the following way where l is the material length parameter characterizing the influence of a higher-order moment. Substituting (7) into (11) yields the following nonzero components of ϕ i , χ i j , m i j : In the case of plane stress state after substitution of (8) and (12) into (10), the full deformation energy of the bended beam takes the following form After computation of integrals with regard to z and a few simple transformations one gets where 3 12 , Ebαθdz, It should be emphasized that in order to get a plane deformable state, it is sufficiently to substitute E bỹ E = E/ 1 − ν 2 and α byα = α (1 + ν) in formula (15).
Kinetic beam energy vibration K takes the following form External work associated with the forces, moments and stress on the boundary plane is governed by the formula (17) whereN ,V ,M σ andM m denote the axial resulting forms of normal stress σ 11 , transversal resulting force of shear stress, resulting moment of normal stresses σ 11 and the resulting moment about the OY axis due to action m 12 in the transversal cross sections, respectively.
The Hamilton principle yields the following formula which after changing variables u, w and ψ, integration by parts and comparison of the terms standing by δu, δw and δψ to zero, yields the following resulting equations of motion and the following boundary conditions 2.3 Coupled problem of thermoelasticity for Timoshenko beams Assuming small temperature variation, a coupled heat transfer equation for a thermoelastic body can be recast to the following form [104] λθ ,ii = ρc v ∂θ ∂t where λ heat transfer coefficient, c v heat capacity of the unit material volume.
Recall that the heat transfer equation is coupled with beam deformation by heat production through compression/extension of the beam. It is assumed that the volume deformation (trace of the deformation tensor) is the only source of heat production and hence there is a lack of curvature tensor in Eq. (21). In the case of plane stress state (8), the trace of deformation tensor ε ii can be presented in the following form Substituting (22) into (21) and assuming that heat gradients along the O Z axis are essentially larger than the gradient along the remaining O X and OY axes, Eq. (21) takes the following form Let us introduce the coefficients of heat diffusion D = λ/ρc v and relaxation R = Eα 2 T 0 /ρc v . Then, Eq. (23) can be presented as follows In the case of plane deformation the trace of the deformation tensor is yielded by relations (9) and takes the following form Substituting (25) into (21) and taking into account heat gradients only along the O Z axis, the following equation is derived which can be rewritten in terms of heat diffusion and relaxation: As boundary conditions for PDEs (24) and (27), we take the thermal isolation condition with regard to the upper and lower microbeam surfaces, i.e. ∂θ/∂z = 0 for z = ±h/2.

Quality factor of the resonator with thermoelastic damping
We consider the plane stress state. In what follows we neglect the quality = 2Eα 2 T 0 (1 + ν)/(1 − 2ν) ρc v with comparison to one (1) [2], and we recast Eq. (24) to the following form It follows from (28) that temperature θ induced in the deformed beam and expressed in terms of displacements u, ψ, w can be divided into two parts: bending deformation and extension deformations, i.e. we have where θ 1 (x, t) and θ 2 (x, t) are defined by the following two equations Let the transversal beam cross section be symmetric with respect to its neutral axis z = 0. Therefore, since non-homogeneous terms in Eqs. (30), (31) (secondary terms in these equations) are symmetric and anti-symmetric with regard to z = 0, respectively, functions θ 1 (x, z, t) = θ 1 (x, t) f 1 (z) should be even and function θ 2 (x, z, t) = θ 2 (x, t) f 2 (z) should be odd with respect to z.
Before integration along the area of the cross section A Eq. (30), (31), we multiply Eq. (31) by z, and we get where four constants G i are as follows If in the first equation of (19) we neglect longitudinal force G (x, t) and longitudinal inertial force ρ A∂ 2 u/∂t 2 , then it can be double integrated with respect to x, and we get In the case of the beam with fixed ends (u (0, t) = u (L , t) = 0), we get C 2 (t) = 0, and hence Therefore, longitudinal force N does not depend on x and remains constant along the beam In what follows, we consider only coupled thermoelastic vibrations, and hence, we take F (x, t) = C (x, t) = 0 in (19) andN =V =M σ =M m = 0 in (20). In these conditions and taking into account (37), the following two equations of motion are obtained and the boundary conditions follow Let us express N and M T by temperature θ 1 (x, t) and θ 2 (x, t) as follows As a result, we obtain the closed system of Eqs. (32), (33) and (38)-(40) for the nonlinear analysis of TED of the Timoshenko microbeams.
Our goal is to derive a non-dimensional counterpart form of the governing PDEs. For this purpose, we introduce the following non-dimensional parameters Formulas (40) yield and hence Eqs. (32) and (33) take the following nondimensional form (bars over non-dimensional quantities are omitted) ∂θ 1 (x, t) ∂t whereas Eqs. (38) and (39) take the following counterpart non-dimensional form We consider simple harmonic vibrations in the form where ω stands for circular frequency, and w 0 (x) , ψ 0 (x) are approximated by a linear vibration mode. Observe that a large beam deflection may cause the frequency to depend on the amplitude of vibrations. Therefore, owing to relations (43), (44) and (42a), (42b), other remaining functions have the following forms Substituting (47) into (43) yields whereas integration along the beam length gives 1 τ We have which is substituted into (42a). Therefore where Function S ω can be presented in the following form We proceed in the same way with Eqs. (42b) and (44) to get Hence is substituted into (42b) to yield where Again, we find real and imaginary parts of D ω in the following way A thermoelastic dissipation can be quantified by a ratio of the loss of mechanical energy per a cycle of vibration to the general conserved energy of deformation, and it can be estimated using a required amount of energy for each cycle to keep a periodic harmonic motion [2,85,86]. The required amount of energy for an infinitely small beam element dx located in point x over period T = 2π/ω is defined by the following formula Substituting (51), (56) for N and M into (59), one obtains [2,85,86,96] Energy loss extension Energy loss bending and it is clear that imaginary parts of the complex axial (S ω ) and bending (D ω ) stiffness define the mechanical energy loss. General generated energy of deformation due to extension bending, shear deformation and higher-order moments is given by the following nondimensional formula Therefore, employing Eqs. (60)-(62), general formula describing thermoelastic dissipation or inverse quality factor Q is as follows

Numerical and analytical results for a rectangular cross section beam resonator
We study the beam resonator with rectangular cross section shown in Fig. 1. Functions f 1 (z) , f 2 (z) that occur in (29) can be yielded by the adiabatic condition applied to the upper and lower microbeam surfaces, i.e. we have d f /dz = 0 for z = ±1/2 in the nondimensional coordinates. Note that the adiabatic condition yields G 2 = 0, which is a consequence of the general theorem on divergence that holds for arbitrary forms of the beam cross section. As we have already mentioned, function f 1 (z) representing the temperature field generated by mechanical extension should be an even function of z. If one develops f 1 (z) into a series with regard to z up to the second order, it is visible that a coefficient standing by the second order should be equal to zero to fit the adiabatic state of the beam surface. Therefore, since deformation along extension does not depend on z, it implies f 1 (z) = const. It follows from the definition of f 1 (z) that the constant can be chosen arbitrarily, and hence it is allowed to take f 1 (z) = 1. Therefore, we obtain the values of coefficients in Eq. (43) In the case of function f 2 (z) describing the change in the temperature field generated by the bending deformation, which is odd, for the rectangular cross section one gets Temperature profile f 2 (z) along beam thickness is described by the following formula [20] and In further computations, we take into account (65). In order to study the influence of the length scalar parameter on TED, we consider beams made from gold, which are extensively employed in MEMS. Material and geometric properties of the resonator are given in Table 1. The values of the length material parameter l for different beam thicknesses can be found in references [105,106].

Analytical solution for the linear size-dependent beam
We take into account undamped modes of beam vibrations as the eigenfunctions (see [107]). It is assumed that the effect of large deflection in the vibrational regime plays the secondary role, and therefore is not meaningful and the vibrational beam regime can be approximated in a linear way modelled by a beam clamped from both sides [20,108,109]. Substituting (47) into (45) for θ 1 (x) = θ 2 (x) = 0 and N = 0 (linear case), the following system of governing PDEs is obtained The latter system of PDEs can be recast to the following counterpart form The following notation is introduced In the case of classical (F = 0) beam clamped on its both sides, equations (70) follow from (67)-(69) Solutions to (70) in the following forms are searched Formula (66) yields The characteristic equation has the following roots and α and β can be easily found (see (71) and (72)) or equivalently Boundary conditions (46) for the clamped beam (cc) follow w| x=0,1 = ψ| x=0,1 = 0.
Hence, the following eigenfunctions are derived where δ = M 2 sin β − M 1 shα (chα − cos β) , H α , and the frequency equation is as follows Quantity a in (79) is equal to the maximum value of the eigenfunction in square brackets of the first formula of (79), whereas W 0 is the value of non-dimensional amplitude of beam deflection. Solving Eq. (81) and defining a, functions w 0 (x) , ψ 0 (x) in (79) are substituted to (63) in order to get quality factor Q of the resonator. As a case study, we consider a beam made from gold (see Table 1).
Graphs of the fundamental vibration mode w 0 (x) , ψ 0 (x) for W 0 = a = 1.589 for the case of beam clamped on both sides (c-c) for L/ h = 100 are shown in Figs. 2 and 3.
The algebraic equation for frequency determination (81) is solved with MathCad. Figure 4a presents the dependence of linear non-dimensional frequencies of the fundamental vibration mode for the classical linear Timoshenko beam made from gold (F = 0) for boundary conditions (c-c), and for different values of its length L as well as for thickness h = 2 μm versus the beam length: 1-analytical solution of Eq. (81); 2-numerical solution. Reduction of the problem of infinite dimension is reduced to the Cauchy problem based on the method of finite differences of the secondorder accuracy. The Cauchy problem is solved using the sixth-order Runge-Kutta method. Figure 4b shows the dependence of quality factor of the beam resonator estimated by formulas (63) for W 0 = 10 −6 (linear case) versus beam length L.
In the case of simply supported classical beam (s-s), the boundary conditions take the following form The frequency equation is obtained based on (71), (72) and takes the following form sin β = 0.
In the case of the fundamental mode β = π . It follows from (80) that the eigenfunctions take the following form Formula (66) for F = 0 yields which allows us to define the frequency of the fundamental mode of vibrations ω 2 F=0 for (s-s) Timoshenko beam.
The dependence of linear non-dimensional frequencies of the fundamental mode of the classical linear Timoshenko beam (F = 0) for the (s-s) boundary condition in dependence of its length is shown in Fig. 5a where analytical/numerical solution is marked by 1/2. Figure 5b illustrates the dependence of beam quality factor versus beam length L obtained from formulas (63) for the linear case.
In order to obtain linear frequencies in the modified couple stress theory (F = 0), we employ Eq. (66) and notations (69) to get the following sixth-order characteristic equation The equation can be recast to the following form [110] where and s 2 is the solution to the following cubic equation The analysis of coefficients F = 1 4 l 2 G A E I in (41) leads to a conclusion that in the case of gold it has the order of 10 −1 for l = 10 −6 m, whereas the remaining coefficients P and H are of several thousand times larger in the interval L/ h = 25 − 100. If we neglect the first term in Eq. (85), since a 2 1 ≈ 10 −2 , hence they are much smaller than coefficients a 1 a 3 , a 4 a 2 , a 2 4 standing by other terms. So, one gets The analysis shows that q 1,2 are real, q 3,4 are complex conjugated, and q 5,6 are real. Therefore We assume the following form of the unknown functions Substituting (89) into (66), the following coefficients are defined We consider as earlier two types of the boundary conditions: clamped-clamped beam (c-c) with boundary conditions w| x=0,1 = w x=0,1 = ψ| 0,1 = 0 and the simple-simple supported beam (s-s) with boundary conditions w| x=0,1 = w x=0,1 = ψ 0,1 = 0.
In the case of (c-c) boundary conditions, the associated frequency equation takes the following form where The obtained Eq. (91) is strongly transcendent and to find its solution is not an easy task. Therefore, in order to estimate the frequencies for the modified couple stress theory (F = 0), we employ an approximate method. Namely, using modes (79) as the vibration modes for F = 0, one may also approximately compute fundamental frequency F = 0 via the Rayleigh method. The eigenfrequency is defined by comparison of the formulas governing kinetic and potential energies of the vibrating beam. As follows from (62), under the lack of temperature field Re [S ω ] = Re [D ω ] = 1 and the values of kinetic and potential energies of the linear beams carrying out vibrations on the fundamental frequencies ω F =0 are defined by the following formulas where functions w(x), ψ (x) are defined by (89) and (c-c) beam boundary conditions. The exact value of frequency ω F =0 , owing to the Rayleigh relation, is as follows If we substitute w 0 (x), ψ 0 (x) from (79) for F = 0, instead of w(x), ψ (x), we obtain In other words, the eigenfrequencies for F = 0 are computed directly from the classical case for F = 0. In reference [111] a similar approach has been employed to find approximate value ω 2 F =0 using the Bernoulli hypothesis ∂w 0 /∂ x = −ψ 0 . If we neglect Poisson's effect, (96) yields However, in this case, the shear deformations disappear, which does not hold for the Timoshenko beam. Relation (96) is free from this statement, and allows us to get more exact values for the Timoshenko beam using functions (79). Figure 6a illustrates the dependence of linear nondimensional eigenfrequencies of the fundamental vibrational mode for the non-classical linear Timoshenko beam (F = 0.260417) for the boundary condition (cc) versus its length L: 1-analytical solution obtained from (96); 2-analytical solution yielded by (97); 3numerical method. Figure 6b shows the dependence of the quality factor of the beam resonator on the beam length estimated via formulas (63) (linear case). Figure 6a shows that the values of ω 2 F =0 obtained from Eqs. (96) and (97) practically coincide only for both values of L, i.e. if the Timoshenko model is close to the Bernoulli-Euler model. Therefore, we employ formula (96) to detect the size-dependent frequencies In the case of the simply supported beam the analysis is simpler. If F = 0, then only the form of the coefficient M 2 = Fπ 4 + H π 2 − ω 2 / −Fπ 3 + H π is changed. Frequency equation for F = 0 takes the following form After defining the frequencies and substituting the eigenfunctions w = sin π x, ψ = M 2 cos π x into relations (63), the quality factor Q is defined in the considered linear case. Figure 7a presents the dependence of the linear non-dimensional eigenfrequencies of the fundamental mode for the non-classical linear (s-s) Timoshenko beam versus its length L: 1-analytical solution based on (98); 2-numerical solution. Figure 7b shows the dependence of the beam quality resonator obtained by formulas (63) in the linear case and for the frequencies estimated from equation (98). The influence of vibration amplitude on the resonance frequencies can be estimated from equations (45) with regard to TED. For this purpose we employ the Rayleigh ratio which is analogous to that [see (93) If functions w(x), ψ (x) are the exact eigenfunctions for the first vibration mode for the given boundary conditions, the Rayleigh ratio yields an exact value of nonlinear frequency ω N L , which follows where δ = M 2 sin β − M 1 shα (chα − cos β) , for F = 0, then one may estimate that Here, as before, W 0 /a plays the role of vibration amplitude. Therefore, nonlinear eigenfrequencies can be approximately obtained from linear frequencies using linear modes of vibrations (102). The quality factor of the beam resonator is estimated based on formulas (63) for the nonlinear case and the frequencies obtained from (104).
Dependence of nonlinear frequencies ω 2 N L versus amplitude w obtained from formula (104) is shown in Fig. 8a. Figure 8b presents quality factor Q N L of the beam resonator under the lack of the size beam dependence (F = 0) versus w for three values of beam length L of the clamped beam.
Analogous results for F = 0.260417 are given for different w in Fig. 9.
Substituting (105) into (104) for a = 1, one may find the nonlinear frequency ω 2 N L for different values of the vibration amplitude W 0 . The quality factor of the beam resonator is computed by formula (63) based on the obtained frequencies ω 2 N L and the linear vibration mode (105). Figure 10a presents the dependence of nonlinear frequencies ω 2 N L obtained via formula (104) vs. the vibration amplitude w, and in Fig. 10b the quality factor Q N L of the beam resonator with the lack of the size dependence (F = 0) is given for three values of the beam length. Analogous results for F = 0.260417 are presented in Fig. 11.
In order to check the reliability of computation of eigenfrequencies, the obtained analytical results have been compared with numerical solutions. Figure 12 presents the results of comparison for the beam (c-c) for F = 0 (a) and for F = 0.260417 (b), i.e. taking into account the size-dependent behaviour obtained numerically (2) and analytically (1).  Figure 13 presents the results of comparison of the eigenfrequencies computed for the (s-s) beam for F = 0 (without the size-dependent behaviour) and for F = 0.260417 (b) with the size-dependent behaviour.

Concluding remarks
In this paper, we have mainly studied the influence of von Kármán nonlinearity on the values of frequency, thermoelastic damping and quality factors on MEMS/NEMS Timoshenko beam resonators based on the modified couple stress theory. Both numerical and analytical studies of the derived governing equations show the significance of geometric nonlinearity and thermoelastic nonlinear damping on the size effects of other dynamic features of the studied micro-and nanobeams.
The most important conclusions of our study are summarized in the following three points.
1. Results reported in Figs. 4, 5, 6, 7 with sizedependent behaviour (F = 0) and in the case of linear vibrations exhibit increase in the eigenfrequencies and increase in the quality factor of the beam resonator. 2. In both tested boundary conditions (s-s) and (cc) occurrence of the nonlinearity implies an essential increase in the frequency of the fundamental vibration mode as well of the quality factor of the resonator (Figs. 8, 9, 10, 11). 3. The analysis of the obtained results shows that in the case of the thermoelastic damping, it is necessary to take into account the nonlinear behaviour as well as the size-dependent effects of the microbeams. Both mentioned features are crucial for the quality factor of the beam resonator.