Effect of thermoelastic damping on nonlinear vibrations of a microbeam under electrostatic actuation using nonlinear normal modes

Abstract

A Galerkin-based nonlinear normal mode approach based on the concept of invariant manifolds is adopted to analyze the effect of thermoelastic damping on the nonlinear vibrations of a microbeam with mid-plane stretching and electrostatic actuation. The derivation of the mathematical model starts with two coupled equations: a general Euler–Bernoulli equation for microbeam vibrations, which includes thermal moment and thermal force terms due to the thermoelastic damping effect; and a thermal equation describing the generation of the thermal moment due to bending vibrations. Given the available information about the material properties and geometry of conventional microbeams, and a physical reasoning based on it about the time constant of the thermal equation, the two governing equations are integrated into a single nonlinear partial differential equation. Galerkin discretization is used to transform the governing equation into a system of coupled nonlinear ordinary differential equations that are uncoupled at linear order. The modal analysis is then followed by finding an invariant manifold leading to a nonlinear ordinary differential equation that governs the vibration behavior of the system. The results of the case study showed that the amount of thermoelastic damping is dependent on the excitation voltage, and the damping increases with increasing voltage. Furthermore, in the presence of thermoelastic damping, by applying the step voltage and starting the vibrations, the frequency increases over time until it reaches a certain limit. With a higher excitation voltage, the frequency increases at a higher rate, but to a lower final limit. In an undamped system, the frequency of oscillations is constant over time, but a higher excitation voltage causes lower frequency vibrations. The results of pull-in analysis showed that thermoelastic damping has no significant effect on static or dynamic pull-in voltage.

Keywords

thermoelastic damping nonlinear vibrations nonlinear normal modes microbeam electrostatic actuation

1. Introduction

Micro-electromechanical systems (MEMS) continue to make for an active area of research as a result of their reliability, low cost, and low power consumption. Therefore, it is still relevant to look for more efficient analysis methods aimed at accurate designs, especially when it comes to MEMS that exhibit nonlinear behavior. Many MEMS devices incorporate microbeams under various loading conditions. The electrostatic force is one such loading which is generated by a step DC voltage and is a source of nonlinearity in many MEMS devices including micro-switches. A MEMS switch can be described as a voltage-controlled electrostatic actuator that can also function as a mechanical relay (Gabriel M. Rebeiz et al., 2013), with applications ranging from switching between the transmit and receive paths in wireless communications and radar systems to signal routing in phones and phased-array antennas (Tilmans et al., 2003).

Thermoelastic damping (TED) is classified as an intrinsic dissipation mechanism in MEMS (Younis, 2011; Tilmans et al., 1992), as it occurs inside the material of the structure. Generally, intrinsic losses become important only when extrinsic ones are minimized or eliminated (Younis, 2011). It has been shown experimentally that TED can be a dominant source of energy loss in MEMS at room temperature (Roszhart, 1990). In a vibrating structure, vibration-induced elastic deformations create strain gradients. As a result, regions in compression become hotter and those in tension cooler. The resulting irreversible heat flow within the structure increases entropy. Therefore, part of the vibrational energy is converted to the internal energy of the structure, and no longer contributes to the mechanical vibrations.

Long before MEMS devices became practical, Zener (Zener, 1937) pioneered TED research by deriving an analytical approximation to the quality factor of slender thermoelastic beams of rectangular cross-section. Lifshitz and Roukes (Lifshitz and Roukes, 2000) offered a refined solution to a similar problem in the context of MEMS by deriving an exact expression for TED in beams under flexure. The change in frequency due to TED is important in frequency-agile applications, such as micro-resonators (Guo and Rogerson, 2003), which convert the measured quantity to a frequency-shift (Tilmans et al., 1992). Guo and Rogerson (Guo and Rogerson, 2003) studied the effect of thermoelastic coupling on a micro-resonator and showed that the frequency-shift ratio due to TED is much larger than that due to air-damping. Nayfeh and Younis (Nayfeh and Younis, 2004) questioned the application of the classical TED theory to micro-plates. They decoupled the thermal equation from the vibration equation of a thermoelastically damped microplate and used a perturbation method to obtain an analytical expression for the quality factor. De and Aluru (De and Aluru, 2006) also found the classical theory inadequate in the presence of nonlinear electrostatic forces in MEMS, and refined it by taking into account the higher-order harmonics of the excitation frequency. Wong et al. (Wong et al., 2006) utilized the classical thermoelastic model to the in-plane vibration of uniform rings and explored the relationship between ring geometry, scale, and quality factor. MEMS resonators

The fundamental merits of modal analysis as a powerful technique for studying linear systems motivated its extension to nonlinear modal analysis. More specifically, attempts were made to generalize the concept of eigenfunctions in linear vibratory systems, known as mode shapes, to nonlinear normal modes (NNMs). NNMs, like their linear counterparts, are intrinsic phenomena of nonlinear systems and are independent of external excitation or initial conditions.

However, the concept of superposition does not hold in the usual sense for NNMs: The general response of a nonlinear system cannot be constructed by the simple superposition of its NNMs. This translates into more involved computations and less tractable equations as far as exact analytical solutions are concerned. However, once the equation of motion in a single NNM is obtained, one can achieve the same accuracy as that obtained using several linear modes, but with significantly reduced computational cost (Boivin et al., 1995). Furthermore, as inherent features of a vibrating system, NNMs can provide insight into its dynamics and lend themselves easily to parametric studies (Kerschen, 2014) and therefore serve as an efficient design tool. Rosenberg (Rosenberg, 1962) was the first to make an attempt at a formal definition of NNMs. Although his work was later improved upon by Rand (Rand, 1971), Manevich and Mikhlin (Manevich and Mikhlin, 1972), and Vakakis (Vakakis, 1991), all those efforts were limited to conservatives systems. The next theoretical leap was made by Shaw and Pierre (Shaw and Pierre, 1991, 1993), who proposed a new definition of NNMs that can be applied to non-conservative, gyroscopic, and piece-wise linear systems. Their method utilized concepts in the theory of invariant manifolds and contained all previous methods as special cases (Shaw and Pierre, 1991). They defined a motion in a normal mode as one where all displacements and velocities are functionally related to a single arbitrarily selected displacement and velocity pair as the reference. They also developed the method to be applicable to continuous systems, using two approaches. The first approach tackled the original partial differential equations (PDEs) directly (Shaw and Pierre, 1994; Hsieh et al., 1994) and was based on the selection of a single suitable material point on the structure as the reference point, whose dynamics determined the entire displacement and velocity fields. The second approach (Shaw, 1994) discretized the governing PDEs into a set of ordinary differential equations (ODEs) using a Galerkin projection and then applied the invariant manifold technique developed in Shaw and Pierre (1993). The latter approach involved selection of a single linear mode (modal coordinate), rather than a material point, as the reference. Boivin et al. (Boivin et al., 1995) compared the two approaches and reported the latter to be more reliable and more practical. Nayfeh (Nayfeh, 1995) compared another approach for constructing NNMs based on the method of multiple scales with the two approaches based on invariant manifolds and reported the multiple scales approach to be the simplest and least computationally demanding. However, Shaw (Shaw, 1994) indicated a discrepancy between the results obtained from the perturbation and invariant manifold approaches and argued that since the invariant manifold approach is based on the fundamental principle of dynamic invariance, it yields the correct results, and perturbation strategies should be used with caution.

The pull-in voltage has been the subject of extensive research as an important design parameter in the development of many electrostatic devices (Magrab, 2012). The static pull-in voltage is the minimum value for which the electrostatic force pulling the two electrodes together meets the following two conditions: first, it is balanced statically (i.e., in the absence of inertia forces) by the restoring force of the deformable electrode; and second, it is large enough for the electromechanical stiffness of the system to vanish. If the latter condition is met dynamically (i.e., during motion), then the corresponding voltage is referred to as the dynamic pull-in voltage. It has been shown both theoretically (Krylov and Maimon, 2004) and experimentally (Sattler et al., 2002) that the voltage at which dynamic pull-in occurs is lower than the corresponding static/quasi-static value for the same device under a given step voltage.

The effects of TED on the forced vibrations of MEMS resonators have already been the subject of numerous investigations (Rezazadeh et al., 2012; Fathalilou and Rezazadeh, 2016; Bostani and Karami-Mohammadi, 2017; Sun et al., 2006). The effects of TED on the transient vibrations of microbeams have also been the subject of some studies. Haddadzadeh Hendou and Karmi-Mohammadi (Haddadzadeh Hendou and Karami-Mohammadi, 2014) applied an NNM method to investigate the effect of TED on the transient vibrations of a cantilever microbeam, with large deflections as the source of nonlinearity. They applied the NNM method directly to the governing PDEs.

The present study examines the effects of TED on the vibrational and dynamic behavior of a clamped-clamped microbeam under the excitation of a DC electrostatic voltage. The nonlinearity is caused by electrostatic excitation and mid-plane stretching. A Galerkin-based NNM approach based on the concept of invariant manifolds is adopted.

The rest of the text is organized as follows. Section 2 includes the definition of the problem and its characteristics and assumptions. Section 3 begins with a general Euler–Bernoulli equation for microbeam vibrations that includes the effects of thermal moment and thermal force due to TED. Then, a thermal equation that describes the relationship between thermal moment and bending vibrations is added and the appropriate final form is obtained. Based on a physical argument about the time constant of the thermal equation, these two equations are integrated and a single governing equation is obtained as a nonlinear PDE, which is subsequently written in a non-dimensional form. In Section 4, an NNM analysis based on the concept of invariant manifolds is presented. First, discretization is carried out using the Galerkin method to arrive at a system of coupled nonlinear equations that are linearly uncoupled. The process of modal analysis is then followed by finding the invariant manifold and leads to a nonlinear ODE governing the vibration behavior of the system. Finally, the equation is solved to arrive at the complete solution. A case study is presented in Section 5, where the solution results are reported and the effects of TED on the dynamic and vibrational behavior of the microbeam are discussed. The paper closes on some concluding remarks in Section 6.

2. Problem statement

As shown in Figure 1, a clamped-clamped Euler–Bernoulli microbeam of length $l$ , width $b$ , and thickness $h$ with TED and moderately large deflection is considered. The beam deflection $w (x, t)$ is assumed to be large enough to give rise to vibration-induced tensile force, and hence mid-plane stretching due to the axially fixed ends; but small enough to justify the linear curvature assumption. The beam is also subjected to an excitation force along its length generated by an electrostatic voltage. The problem involved three sources of nonlinearity, namely, the nonlinear dependence of the electrostatic force on beam deflection, the geometric nonlinearity due to mid-plane stretching, and the nonlinear TED effect. As a result of the mutual dependence of the excitation force and deflection, the system is self-excited upon the application of a step voltage.

A 2-D temperature profile $θ (x, z, t)$ is used to model the thermoelastic behavior of the beam, with negligible thermal gradients along the beam length compared to those along the thickness direction. Temperature variations $θ (x, z, t) = T - T_{0}$ at any material point are assumed to be insignificant compared to the reference (ambient) temperature $T_{0}$ , at which the beam is unstrained and unstressed; in other words, $T - T_{0} ≪ T_{0}$ . Furthermore, with the beam vibrating in vacuum, zero heat flux is considered at the top and bottom surfaces of the beam (i.e., adiabatic thermal boundary conditions). It is also assumed that any thermal disturbance propagated at infinite speed allowing for the “second sound” effect to be neglected.

3. Governing equations

The formulation starts by considering the general Euler–Bernoulli equation for a microbeam with mid-plane stretching and electrostatic excitation, which also includes terms related to the TED effect through the thermal moment (Nayfeh and Pai, 2004):

\begin{array}{l} ρ A \frac{\partial^{2} w}{\partial t^{2}} + \frac{\partial^{2}}{\partial x^{2}} (E I \frac{\partial^{2} w}{\partial x^{2}} + M^{T}) + \frac{\partial^{2} w}{\partial x^{2}} \int_{0}^{l} \frac{1}{l} [N^{T} - \frac{\partial M^{T}}{\partial x} \frac{\partial w}{\partial x} - \frac{E A}{2} {(\frac{\partial w}{\partial x})}^{2}] d x \\ = \frac{ε_{0} b V^{2}}{2 {(g - w)}^{2}} \end{array}

(1)

where

M^{T} = E α b \int_{- \frac{h}{2}}^{\frac{h}{2}} z θ d z

(2)

and

N^{T} = E α b \int_{- \frac{h}{2}}^{\frac{h}{2}} θ d z

(3)

are the thermal moment and thermal axial force, respectively;

w = w (x, t)

represents the deflection of the beam and

θ = θ (x, z, t)

denotes the temperature variations;

ρ

is the mass density,

E

is Young’s modulus, and

α

is the coefficient of thermal expansion;

A = b h

and

I = b h^{3} / 12

are the area and the area moment of inertia of the cross-section of the beam;

ε_{0}

is the permittivity of free space,

V

represents the DC step voltage,

g

denotes the gap between the two electrodes, and the term on the right-hand side is the electrostatic force.

The thermal equation with thermoelastic coupling is given by (Lifshitz and Roukes, 2000)

c_{1} \frac{\partial θ}{\partial t} = χ \frac{\partial^{2} θ}{\partial z^{2}} + c_{2} z \frac{\partial}{\partial t} (\frac{\partial^{2} w}{\partial x^{2}}),

(4)

where

c_{1} = 1 + 2 Δ_{E} \frac{1 + ν}{1 - 2 ν}, c_{2} = \frac{Δ_{E}}{α};

(5)

ν

is Poisson’s ratio and

Δ_{E} = \frac{E α^{2} T_{0}}{ρ c_{m}}

(6)

is the relaxation strength (

Δ_{E} ≪ 1

) representing the relative difference between the adiabatic and isothermal values of Young’s modulus (Wong et al., 2006).

c_{m}

in equation (6) is the specific heat capacity of the material. The thermal diffusivity of the beam denoted by

χ

in equation (4) is given by

χ = \frac{k}{ρ c_{m}},

(7)

where

k

is the thermal conductivity. The thermoelastic equation (4) is described as being one-way coupled, meaning that the stress field that leads to temperature variations is itself unaffected by the changes in temperature (Prabhakar and Vengallatore, 2008). The adiabatic boundary conditions at

z = \pm \frac{h}{2}

can be written as

{\frac{\partial θ}{\partial z} |}_{z = \pm \frac{h}{2}} = 0 .

(8)

Multiplying equation (4) by

z

, integrating from

z = - \frac{h}{2}

z = \frac{h}{2}

, and using equation (2), give

c_{1} \frac{\partial M^{T}}{\partial t} = χ \int_{- \frac{h}{2}}^{\frac{h}{2}} z \frac{\partial^{2} θ}{\partial z^{2}} d z + c_{2} \frac{I}{b} \frac{\partial}{\partial t} (\frac{\partial^{2} w}{\partial x^{2}}) .

(9)

At any material point within the beam, the temperature varies depending on elongation or compression at that point caused by bending vibrations. This mechanism of temperature change suggests a sinusoidal distribution as a reasonable choice for the temperature profile along the thickness of the beam (Sun et al., 2006). A sinusoidal profile that also satisfies the adiabatic boundary conditions (equation (8)) is given by

θ (x, z, t) = Θ (x, t) \sin (\frac{π}{h} z) .

(10)

The following relationship is then observed between $θ (x, z, t)$ and its second derivative with respect to $z$ :

\frac{\partial^{2} θ}{θ z^{2}} = - {(\frac{π}{h})}^{2} Θ (x, t) \sin (π z) = - {(\frac{π}{h})}^{2} θ .

(11)

One can now substitute $- {(\frac{π}{h})}^{2} θ$ for $\frac{\partial^{2} θ}{θ z^{2}}$ from equation (11) into equation (9). The result is

c_{1} \frac{\partial M^{T}}{\partial t} = - {(\frac{π}{h})}^{2} χ \int_{- \frac{h}{2}}^{\frac{h}{2}} z θ d z + c_{2} \frac{I}{b} \frac{\partial}{\partial t} (\frac{\partial^{2} w}{\partial x^{2}}) .

(12)

Multiplying equation (12) by

E α b

, substituting

M^{T}

from equation (2), and rearranging the result yield

\frac{E α b h^{2} c_{1}}{π^{2} χ} \frac{\partial M^{T}}{\partial t} + M^{T} = \frac{E α h^{2} c_{2} I}{π^{2} χ} \frac{\partial}{\partial t} (\frac{\partial^{2} w}{\partial x^{2}}) .

(13)

The coupled equations (1) and (13) in two unknowns $w$ and $M^{T}$ describe the vibrational behavior of the microbeam with electrostatic excitation and TED. $N^{T}$ can be eliminated from equation (1) by plugging the assumed temperature profile given by equation (10) into equation (3).

While the left-hand side of equation (13) shows the inherent dynamics of $M^{T}$ with the time constant

τ = \frac{E α b h^{2} c_{1}}{π^{2} χ},

(14)

the right-hand side can be considered as an excitation which is a function of

w

By examining the practical values of material properties and geometry of conventional silicon microbeams, it is concluded that the time constant $τ$ in equation (13) is very small and can be ignored. Therefore, equation (13) can be written as

M^{T} = γ \frac{\partial}{\partial t} (\frac{\partial^{2} w}{\partial x^{2}}), γ = \frac{E α h^{2} c_{2} I}{π^{2} χ} = \frac{E^{2} h^{2} α^{2} T_{0} I}{π^{2} k} .

(15)

By substituting $M^{T}$ from equation (15) in equation (1), a single equation governing the microbeam behavior is obtained:

\begin{array}{l} ρ A \frac{\partial^{2} w}{\partial t^{2}} + E I \frac{\partial^{4} w}{\partial x^{4}} + γ \frac{\partial}{\partial t} (\frac{\partial^{4} w}{\partial x^{4}}) + \frac{\partial^{2} w}{\partial x^{2}} \int_{0}^{l} [- \frac{γ}{l} (\frac{\partial^{4} w}{\partial t \partial x^{3}}) \frac{\partial w}{\partial x} - \frac{E A}{2 l} {(\frac{\partial w}{\partial x})}^{2}] d x \\ = \frac{ε_{0} b V^{2}}{2 {(g - w)}^{2}} . \end{array}

(16)

Next, by introducing the dimensionless parameters

w = \frac{w^{*}}{g}, x = \frac{x^{*}}{l}, t = \frac{t^{*}}{l^{2}} \sqrt{\frac{E I}{ρ A}},

(17)

one can rewrite the governing equation (16) in the following non-dimensional form

\frac{\partial^{2} w}{\partial t^{2}} + \frac{\partial^{4} w}{\partial x^{4}} + b_{1} \frac{\partial}{\partial t} (\frac{\partial^{4} w}{\partial x^{4}}) + \frac{\partial^{2} w}{\partial x^{2}} \int_{0}^{1} [b_{2} (\frac{\partial^{4} w}{\partial t \partial x^{3}}) \frac{\partial w}{\partial x} + b_{3} {(\frac{\partial w}{\partial x})}^{2}] d x = \frac{b_{4}}{{(1 - w)}^{2}},

(18)

where the asterisk superscripts are dropped to simplify the notation, and the constant parameters $b_{i}$ are given by

\begin{array}{l} b_{1} = \sqrt{\frac{E I}{ρ A}} \frac{E h^{2} α^{2} T_{0}}{π^{2} l^{2} k}, b_{2} = - \frac{g^{2}}{l^{4}} \sqrt{\frac{E I}{ρ A}} \frac{E h^{2} α^{2} T_{0}}{π^{2} k} = - \frac{g^{2}}{l^{2}} b_{1}, \\ b_{3} = - \frac{A g^{2}}{2 I}, b_{4} = \frac{6 ε_{0} l^{4}}{E h^{3} g^{3}} V^{2} . \end{array}

(19)

The non-dimensional electrostatic force on the right-hand side of equation (16) can be approximated by a third-degree polynomial (Xie et al., 2003)

F_{e} = \frac{b_{4}}{{(1 - w)}^{2}} = a_{0} V^{2} + a_{1} V^{2} w + a_{2} V^{2} w^{2} + a_{3} V^{2} w^{3} .

(20)

where coefficients

a_{0}

a_{3}

are constant coefficients calculated by curve fitting.

Substituting equation (20) into equation (16), one obtains

\begin{array}{l} \frac{\partial^{2} w}{\partial t^{2}} + \frac{\partial^{4} w}{\partial x^{4}} + b_{1} \frac{\partial}{\partial t} (\frac{\partial^{4} w}{\partial x^{4}}) + \frac{\partial^{2} w}{\partial^{2} x} \int_{0}^{1} [b_{2} (\frac{\partial^{4} w}{\partial t \partial x^{3}}) \frac{\partial w}{\partial x} + b_{3} {(\frac{\partial w}{\partial x})}^{2}] d x \\ - a_{0} V^{2} - a_{1} V^{2} w - a_{2} V^{2} w^{2} - a_{3} V^{2} w^{3} = 0 . \end{array}

(21)

Equation (21) is the final form of the equation of motion that is used in the NNM analysis in the following section. The last two terms on the left-hand side of equation (21) can be considered as two nonlinear springs with negative stiffness. The other nonlinear terms, that is, the first and second terms within the integral, are due to the thermal moment and the nonlinear stretching effect, respectively.

4. Nonlinear normal mode analysis

4.1. Galerkin discretization

Before applying the Galerkin procedure, it is convenient to use operator notation in subsequent computations by rewriting equation (21) in the form

\frac{\partial^{2} w}{\partial t^{2}} + D [\frac{\partial w}{\partial t}] + L [w (x, t)] + N [w (x, t)] = 0,

(22)

where

D [.]

and

L [.]

represent the linear damping and linear spatial operators, respectively, defined by

L [.] = (\frac{\partial^{4}}{\partial x^{4}} + c_{L}) [.], c_{L} = - a_{1} V^{2}

(23)

and

D [.] = b_{1} \frac{\partial^{4}}{\partial x^{4}} [.]

(24)

The spatial and temporal operator $N [w (x, t)]$ in equation (22) holds all the nonlinear terms as well as the non-homogenous term in equation (21):

N [w (x, t)] = - a_{0} V^{2} + \frac{\partial^{2} w}{\partial x^{2}} \int_{0}^{1} [b_{2} (\frac{\partial^{4} w}{\partial t \partial x^{3}}) \frac{\partial w}{\partial x} + b_{3} {(\frac{\partial w}{\partial x})}^{2}] d x - a_{2} V^{2} w^{2} - a_{3} V^{2} w^{3}

(25)

The linear undamped system corresponding to the nonlinear system represented by equation (21) can be written as

\frac{\partial^{2} w}{\partial t^{2}} + L [w (x, t)] = 0,

(26)

whose eigenfunctions and eigenvalues are given by (Hagedorn and DasGupta, 2007)

\begin{array}{l} W_{j} (x) = [\cosh (β_{j} x) - \cos (β_{j} x)] + \frac{\cos (β_{j}) - \cosh (β_{j})}{\sinh (β_{j}) - \sin (β_{j})} [\sinh (β_{j} x) - \sin (β_{j} x)], \\ j = 1,2,3, \dots, \end{array}

(27)

and

ω_{j}^{2} = β_{j}^{4} - a_{1} V^{2},

(28)

respectively, where

β_{j} = j π + \frac{π}{2}

(29)

are the roots of the approximate characteristic equation corresponding to the linear system equation (26). Since

L [.]

is a self-adjoint operator, the corresponding eigenfunctions are orthogonal (Nayfeh and Nayfeh, 1994), that is,

〈 W_{i}, W_{j} 〉 = δ_{i j},

(30)

where

〈, 〉

represents the inner product define in the usual manner, and

δ_{i j}

is the Kronecker delta. Using equation (30) and the eigenvalue problem derived from equation (26), one can show that

〈 W_{i}, L [W_{j}] 〉 = ω_{j}^{2} δ_{i j} .

(31)

Equations (23) and (24) can be combined to show that the following relationship holds between the damping and linear operators:

D [\cdot] = b_{1} L [\cdot] + c [\cdot], c = - c_{L} b_{1} .

(32)

This is a case of proportional damping, which offers the advantage of decoupling the discretized equations with respect to the damping term when the solution of the damped system is expanded in terms of eigenfunctions of the undamped system (Hagedorn and DasGupta, 2007). Furthermore, one can use equations (26), (30), and (31) to show that

〈 W_{i}, D [W_{j}] 〉 = 2 ζ_{j} ω_{j}^{2} δ_{i j},

(33)

where

ζ_{j} = \frac{b_{1} (ω_{j}^{2} - c_{L})}{2 ω_{j}} .

(34)

The Galerkin discretization procedure starts by writing the solution $w (x, t)$ to equation (21) as an infinite series in linear mode shapes $W_{j}$ :

w (x, t) = \sum_{j = 1}^{\infty} q_{j} (t) W_{j} (x)

(35)

where

q_{j} (t)

is the contribution of the

j

th linear mode to the general response. Plugging equation (35) into equation (22) results in

\sum_{j = 1}^{\infty} {\ddot{q}}_{j} W_{j} + D [\sum_{j = 1}^{\infty} {\dot{q}}_{j} W_{j}] + L [\sum_{j = 1}^{\infty} q_{j} W_{j}] + N [\sum_{j = 1}^{\infty} q_{j} W_{j}] = 0 .

(36)

The inner product of equation (36) with $W_{i}$ simplifies to

\sum_{j = 1}^{\infty} {\ddot{q}}_{j} \underset{δ_{i j}}{\underset{⏟}{〈 W_{i}, W_{j} 〉}} + \sum_{j = 1}^{\infty} {\dot{q}}_{j} \underset{2 ζ_{j} ω_{j} δ_{i j}}{\underset{⏟}{〈 W_{i}, D [W_{j}] 〉}} + \sum_{j = 1}^{\infty} q_{j} \underset{ω_{j}^{2} δ_{i j}}{\underset{⏟}{〈 W_{i}, L [W_{j}] 〉}} + \underset{G_{i} (q_{1}, q_{2}, . . .; {\dot{q}}_{1}, {\dot{q}}_{2}, . . .)}{\underset{⏟}{〈 W_{i}, N [\sum_{j = 1}^{\infty} q_{j} W_{j}] 〉}} = 0 .

(37)

Using the properties of the corresponding linear system given by equations (30), (31), and (33), one obtains

{\ddot{q}}_{i} + 2 ζ_{i} ω_{i} {\dot{q}}_{i} + ω_{i}^{2} q_{i} + G_{i} (q, \dot{q}) = 0, i = 1,2,3, \dots,

(38)

where

G_{i} (q, \dot{q}) = 〈 W_{i}, N [\sum_{j = 1}^{\infty} W_{j} q_{j}] 〉

(39)

and

q = q_{j},

(40)

\dot{q} = {\dot{q}}_{j}, j = 1,2,3, \dots .

(41)

Equation (38) are the discretized equations of motion, which are uncoupled at linear order (Shaw, 1994) with the coupling term $G_{i} (q, \dot{q})$ holding all the nonlinear terms as well as the non-homogenous term.

4.2. Invariant manifold derivation

As the first step in finding the invariant manifold, equation (38) is transformed into the state space form as

{\dot{q}}_{i} = p_{i},

(42)

{\dot{p}}_{i} = - 2 ζ_{i} ω_{i} p_{i} - ω_{i}^{2} q_{i} - G_{i} (q, \dot{q}), i = 1,2,3, \dots .

(43)

Based on the invariant manifold approach presented in Shaw (1994), the next step is to select one coordinate pair $(q_{k}, p_{k})$ as a master coordinate, and describe all other pairs $(q_{i}, p_{i})$ as functions of $(q_{k}, p_{k})$ :

q_{i} (t) = Q_{i k} (q_{k}, p_{k}),

(44)

p_{i} (t) = P_{i k} (q_{k}, p_{k}), i = 1,2,3, \dots, i \neq k .

(45)

By plugging equations (44) and (45) into equations (42) and (43), the following sets of PDEs are obtained:

\frac{\partial Q_{i k}}{\partial q_{k}} p_{k} + \frac{\partial Q_{i k}}{\partial p_{k}} (- 2 ζ_{k} ω_{k} p_{k} - ω_{k}^{2} q_{k} - G_{k} (Q_{k}, P_{k})) = P_{i k},

(46)

\begin{array}{l} \frac{\partial P_{i k}}{\partial q_{k}} p_{k} + \frac{\partial P_{i k}}{\partial p_{k}} (- 2 ζ_{k} ω_{k} p_{k} - ω_{k}^{2} q_{k} - G_{k} (Q_{k}, P_{k})) \\ = - 2 ζ_{i} ω_{i} P_{i k} - ω_{i}^{2} Q_{i k} - G_{i} (Q_{k}, P_{k}), i = 1,2,3, \dots, i \neq k, \end{array}

(47)

where

Q_{k} = Q_{i k} (q_{k}, p_{k}),

(48)

P_{k} = P_{i k} (q_{k}, p_{k}), i = 1,2,3, \dots, i \neq k .

(49)

The reason for writing equations (48) and (49) for $i \neq k$ is that they are trivially satisfied for $i = k$ . The solution to these equations govern the geometry of the invariant manifold that constraints the dynamics of the system when vibrating in the $k$ th NNM. Since these equations are generally as difficult to solve analytically as the original equation (38), an approximate solution near the origin is sought in terms of a series expansion. For this purpose, the nonlinearities represented by $G_{i} (q, \dot{q})$ can, to cubic order in the $q_{j}$ ’s and ${\dot{q}}_{j}$ ’s, be expressed as

G_{i} (q, \dot{q}) = γ_{1 i k} + γ_{2 i k} q_{k}^{2} + γ_{3 i k} q_{k}^{3} + γ_{4 i k} q_{k}^{3} + γ_{5 i k} p_{k} q_{k}^{2},

(50)

where

\begin{array}{l} γ_{1 i k} = - a_{0} V^{2} \int_{0}^{1} W_{i} d x, \\ γ_{2 i k} = - a_{2} V^{2} \int_{0}^{1} W_{i} W_{k}^{2} d x, \\ γ_{3 i k} = - a_{3} V^{2} \int_{0}^{1} W_{i} W_{k}^{3} d x, \\ γ_{4 i k} = b_{3} [\int_{0}^{1} (W_{k}^{″}) (W_{i}) d x] \int_{0}^{1} {(W_{k}^{'})}^{2} d x, \\ γ_{5 i k} = b_{2} [\int_{0}^{1} (W_{k}^{″}) (W_{i}) d x] \int_{0}^{1} W_{k}^{‴} W_{k}^{'} d x . \end{array}

(51)

The solution proceeds by approximating the functions $Q_{i k}$ and $P_{i k}$ as nonlinear polynomials. The nonlinearities in equation (50) justify the expansion of the $Q_{i k}$ and $P_{i k}$ series as follows:

Q_{i k} = Κ_{1 i k} + Κ_{2 i k} q_{k}^{2} + Κ_{3 i k} p_{k}^{2} + Κ_{4 i k} q_{k}^{3} + Κ_{5 i k} q_{k}^{2} p_{k} + Κ_{6 i k} q_{k} p_{k}^{2} + Κ_{7 i k} p_{k}^{3} + \dots

(52)

P_{i k} = Π_{1 i k} q_{k} p_{k} + Π_{2 i k} q_{k}^{3} + Π_{3 i k} q_{k}^{2} p_{k} + Π_{4 i k} q_{k} p_{k}^{2} + Π_{5 i k} p_{k}^{3} + \dots

(53)

i = 1,2,3, \dots, i \neq k

where only the non-zero terms are retained. Substituting equations (52) and (53) into equations (46) and (47) and collecting like power terms lead to

Κ_{1 i k} = - \frac{γ_{1 i k}}{ω_{j}^{2}},

(54)

Κ_{2 i k} = \frac{γ_{2 i k} (2 ω_{k}^{2} - ω_{i}^{2})}{ω_{i}^{2} (ω_{i}^{2} - 4 ω_{k}^{2})},

(55)

Κ_{3 i k} = \frac{2 γ_{2 i k}}{ω_{i}^{2} (ω_{i}^{2} - 4 ω_{k}^{2})},

(56)

Κ_{4 i k} = \frac{(7 ω_{k}^{2} - ω_{i}^{2}) (γ_{3 i k} + γ_{4 i k})}{(9 ω_{k}^{2} - ω_{i}^{2}) (ω_{k}^{2} - ω_{i}^{2})},

(57)

Κ_{5 i k} = \frac{(3 ω_{k}^{2} - ω_{i}^{2}) γ_{5 i k}}{(9 ω_{k}^{2} - ω_{i}^{2}) (ω_{k}^{2} - ω_{i}^{2})},

(58)

Κ_{6 i k} = \frac{6 (γ_{3 i k} + γ_{4 i k})}{(9 ω_{k}^{2} - ω_{i}^{2}) (ω_{k}^{2} - ω_{i}^{2})},

(59)

Κ_{7 i k} = \frac{2 γ_{5 i k}}{(9 ω_{k}^{2} - ω_{i}^{2}) (ω_{k}^{2} - ω_{i}^{2})},

(60)

Π_{1 i k} = - \frac{2 γ_{2 i k}}{ω_{i}^{2} - 4 ω_{k}^{2}},

(61)

Π_{2 i k} = \frac{γ_{5 i k} ω_{k}^{2} (ω_{i}^{2} - 3 ω_{k}^{2})}{(9 ω_{k}^{2} - ω_{i}^{2}) (ω_{k}^{2} - ω_{i}^{2})},

(62)

Π_{3 i k} = \frac{3 (γ_{3 i k} + γ_{4 i k}) (3 ω_{k}^{2} - ω_{i}^{2})}{(9 ω_{k}^{2} - ω_{i}^{2}) (ω_{k}^{2} - ω_{i}^{2})},

(63)

Π_{4 i k} = - \frac{2 γ_{5 i k} ω_{i}^{2}}{(9 ω_{k}^{2} - ω_{i}^{2}) (ω_{k}^{2} - ω_{i}^{2})},

(64)

Π_{5 i k} = \frac{6 (γ_{3 i k} + γ_{4 i k})}{(9 ω_{k}^{2} - ω_{i}^{2}) (ω_{k}^{2} - ω_{i}^{2})} .

(65)

4.3. Vibration response

From equation (35) and then equations (44) and (52), the vibration response in terms of the $k$ th NNM is given by

\begin{array}{l} w_{k} (x, t) = W_{k} (x) q_{k} (t) \\ + \sum_{i \neq k} W_{i} (x) (Κ_{1 i k} + Κ_{2 i k} q_{k}^{2} + Κ_{3 i k} p_{k}^{2} + Κ_{4 i k} q_{k}^{3} + Κ_{5 i k} q_{k}^{2} p_{k} + Κ_{6 i k} q_{k} p_{k}^{2} + Κ_{7 i k} p_{k}^{3} + \dots), \end{array}

(66)

where the governing equation of oscillator $q_{k} (t)$ is obtained by inserting equation (50) into equation (38). The result is

{\ddot{q}}_{k} + 2 ζ_{k} ω_{k} {\dot{q}}_{k} + w_{k}^{2} q_{k} + γ_{1 k k} + γ_{2 k k} q_{k}^{2} + (γ_{3 k k} + γ_{4 k k}) q_{k}^{3} + γ_{5 k k} q_{k}^{2} {\dot{q}}_{k} = 0 .

(67)

The velocity field is obtained by taking the first time derivative of equation (66):

\begin{array}{l} {\dot{w}}_{k} (x, t) = W_{k} (x) p_{k} (t) \\ + \sum_{i \neq k} W_{i} (x) (Π_{1 i k} q_{k} p_{k} + Π_{2 i k} q_{k}^{3} + Π_{3 i k} q_{k}^{2} p_{k} + Π_{4 i k} q_{k} p_{k}^{2} + Π_{5 i k} p_{k}^{3} + \dots) \end{array}

(68)

The response $w_{k} (t)$ of the system vibrating in the $k$ th NNM is calculated by solving equation (67) to obtain $q_{k} (t)$ and substituting it into equation (66).

Furthermore, by setting all the time derivatives in equation (67) equal to zero, one arrives at the following equilibrium equation for the modal coordinate $q_{k}$ :

w_{k}^{2} q_{k} + γ_{1 k k} + γ_{2 k k} q_{k}^{2} + (γ_{3 k k} + γ_{4 k k}) q_{k}^{3} = 0

(69)

Equation (69) is used for static pull-in analysis.

5. Case study results and discussion

The NNM method presented in Section 4 was implemented to perform a case study on a silicone microbeam with typical specifications as listed in Table 1. The first three modes were considered in Galerkin discretization and then NNM analysis was done by choosing

q_{1}

(

k = 1

) as the master coordinate in equations (44) and (45).

Table 1.

Parameter values.

Symbol	Value	Unit
$l$	$610 \times 10^{- 6}$	$[m]$
$b$	$40 \times 10^{- 6}$	$[m]$
$h$	$2.2 \times 10^{- 6}$	$[m]$
$g$	$2.3 \times 10^{- 6}$	$[m]$
$E$	$1.49 \times 10^{11}$	$[Pa]$
$ρ$	$2330$	$[k g / m^{3}]$
$ε_{0}$	$8.854 \times 10^{- 12}$	$[F / m]$
$α$	$2.5 \times 10^{- 6}$	$[1 / K]$
$k$	$1.7$	$[W / (m K)]$
$c_{m}$	$711$	$[J / (k g K)]$
$T_{0}$	$298$	$[K]$
$ν$	$0.2$	$[-]$
$V$	$2$	$[V]$

The coefficients

a_{0}

a_{3}

in the approximate electrostatic force described by equation (20) were calculated by running a least-squares routine. The values are given in Table 2.

Table 2.

Coefficients of the polynomial approximation to the electrostatic force.

Constant	Value
$a_{0}$	0.371804626618547
$a_{1}$	1.11254646236951
$a_{2}$	−1.59148160027687
$a_{3}$	7.83645091887688

5.1. Static pull-in

A static pull-in analysis was performed by solving equation (69) for the master coordinate $q_{1}$ . The absence of a time derivative in the equation shows that damping has no effect on the static pull-in voltage. Figure 2 shows the equilibrium state as a function of applied voltage. The discontinuity of the curve indicates a static pull-in voltage of about 14.5 V, which is consistent with the result reported in Xie et al. (2003). The microbeam becomes unstable for voltages exceeding this value.

5.2. Undamped microbeam, vibrations, and stability

Figure 3 shows the initial portion of the vibration time history of the microbeam at midpoint without TED for different voltages; the phase portraits are shown in Figure 4. Both figures indicate a dynamic pull-in voltage of about 13.457 V. It is observed that for voltages below this value, a higher voltage results in a larger amplitude but lower frequency. As the voltage exceeds the dynamic pull-in voltage, the behavior of the system changes from periodic to non-periodic and the response becomes unstable.

In order to analyze the effect of voltage on frequency for the undamped microbeam, the frequency of vibration was numerically determined from the time response for different voltages. The results shown in Figure 5 indicate both the dynamic pull-in and the nonlinear dependence of frequency on the applied voltage.

5.3. Microbeam with TED, vibrations, and stability

The phase portraits over a period of about 0.1 s for the microbeam with TED are shown in Figures 6 –8 separately for different voltages. For any voltage below the dynamic pull-in voltage, the phase trajectory shows a decrease in amplitude and follows a spiral path inward. The phase trajectory for a voltage value slightly above the dynamic pull-in voltage is also shown in Figure 8. This curve diverges out immediately from a periodic motion and becomes unstable.

Figure 9 shows the effect of voltage on frequency for the microbeam with TED. It is observed that unlike the case of the undamped microbeam, the frequency changes over time due to the TED effect. The results are therefore presented as separate curves for different voltages. It can be seen that the initial frequency values for each voltage are identical to the corresponding values in Figure 5. As the TED takes effect and diminishes the amplitude, the frequency of vibration increases. This increase is observed to be greater at higher voltages and larger amplitudes, because larger-amplitude vibrations invoke a stronger TED effect.

Another point of interest was the effect of voltage on the linear TED ratio $ζ_{m}$ . This was studied by using equations (34) and (28) to plot $ζ_{m}$ as a function of applied voltage. The resulting curve shown in Figure 10 indicates about 32% increase in $ζ_{m}$ when the applied voltage was increased from 2 V to 14 V. Furthermore, the effect of the nonlinear damping term, that is, the last term on the left-hand side of equation (67), was found to be insignificant for all voltages compared to the effect of linear damping represented by the second term in the same equation.

6. Conclusion

A Galerkin-based NNM approach based on the concept of invariant manifolds was implemented to analyze the vibrations of an Euler–Bernoulli microbeam under a DC step voltage. The mathematical model initially consisted of two coupled equations, namely, the thermal and vibration equations. A simplification argument about the time constant of the thermal equation was used to combine the two equations into one. A case study was carried out on a silicon microbeam with typical specifications. The first three modes were considered in Galerkin discretization and the analysis was developed by choosing $q_{1}$ ( $k = 1$ ) as the master coordinate. The results showed that the amount of TED is dependent on the excitation voltage and the damping increases with increasing voltage. Furthermore, in the presence of TED, after the initiation of vibrations, the frequency increases over time until it reaches a certain limit. With a higher excitation voltage, the frequency increases at a higher rate, but to a lower final limit. In an undamped system, the frequency of oscillations is constant over time, but a higher excitation voltage causes lower frequency vibrations. The results of pull-in analysis showed that TED has no significant effect on static or dynamic pull-in voltage.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Ardeshir Karami Mohammadi

References

Boivin

Pierre

Shaw

(1995) Non-linear normal modes, invariance, and modal dynamics approximations of non-linear systems. Nonlinear Dynamics 8: 315–346.

Bostani

Karami-Mohammadi

(2017) Thermoelastic damping in microbeam resonators based on modified strain gradient elasticity and generalized thermoelasticity theories. Acta Mechanica 229: 173–192.

Aluru

(2006) Theory of thermoelastic damping in electrostatically actuated microstructures. Physical Review B 74: 144305.

Fathalilou

Rezazadeh

(2016) Effects of the length scale parameter on the thermoelastic damping of a microbeam considering the couple stress theory. International Journal of Applied Mechanics 08(6): 1650083.

Guo

Rogerson

(2003) Thermoelastic coupling effect on a micro-machined beam resonator. Mechanics Research Communications 30: 513–518.

Hagedorn

DasGupta

(2007) Vibrations and Waves in Continuous Mechanical Systems. John Wiley & Sons.

Haddadzadeh Hendou

Karami-Mohammadi

(2014) Transient analysis of nonlinear Euler–Bernoulli micro-beam with thermoelastic damping, via nonlinear normal modes. Journal of Sound and Vibration 333: 6224–6236.

Shang-Rou

Shaw

Pierre

(1994) Normal modes for large amplitude vibration of a cantilever beam. International Journal of Solids and Structures 31(14): 1981–2014.

Kerschen

(2014) Modal Analysis of Nonlinear Mechanical Systems. Springer.

10.

Krylov

Maimon

(2004) Pull-in dynamics of an elastic beam actuated by continuously distributed electrostatic force. Journal of Vibration and Acoustics 126(3): 332–342.

11.

Lifshitz

Roukes

(2000) Thermoelastic damping in micro- and nanomechanical systems. Physical Review B 61(8): 5600–5609.

12.

Magrab

(2012) Vibrations of Elastic Systems - with Applications to MEMS and NEMS. Springer.

13.

Manevich

Mikhlin

(1972) On periodic solutions close to rectilinear normal vibration modes. Journal of Applied Mathematics and Mechanics 36: 988–994.

14.

Nayfeh

(1995) On direct methods for constructing nonlinear normal modes of continuous systems. Journal of Vibration and Control 1: 389–430.

15.

Nayfeh

(1994) On nonlinear modes of continuous systems. Journal of Vibration and Acoustics 116: 129–136.

16.

Nayfeh

Pai

(2004) Linear and Nonlinear Structural Mechanics. John Wiley & Sons, Inc.

17.

Nayfeh

Younis

(2004) Modeling and simulations of thermoelastic damping in microplates. Journal of Micromechanics and Microengineering 14: 1711–1717.

18.

Prabhakar

Vengallatore

(2008) Theory of thermoelastic damping in micromechanical resonators with two-dimensional heat conduction. Journal of Microelectromechanical Systems 17(2): 494–502.

19.

Rand

(1971) A higher-order approximation for non-linear normal modes in two-degree-of-freedom systems. International Journal of Non-Linear Mechanics 6: 545–547.

20.

Rebeiz

Patel

Han

, et al. (2013) The search for a reliable MEMS switch? IEEE Microwave Magazine 14(1): 57–67.

21.

Rezazadeh

Vahdat

Tayefeh-rezaei

, et al. (2012) Thermoelastic damping in a micro-beam resonator using modified couple stress theory. Acta Mechanica 223: 1137–1152.

22.

Rosenberg

(1962) The normal modes of nonlinear n-Degree-of-Freedom systems. Journal of Applied Mechanics 29(1): 7–14.

23.

Roszhart

(1990) The effect of thermoelastic internal friction on the Q of micromachined silicon resonators.IEEE 4th Technical Digest on Solid-State Sensor and Actuator Workshop. Hilton Head, SC: IEEE.

24.

Sattler

Plötz

Fattinger

, et al. (2002) Modeling of an electrostatic torsional actuator: demonstrated with an RF MEMS switch. Sensors and Actuators A: Physical 97–98: 337–346.

25.

Shaw

(1994) An invariant manifold approach to nonlinear normal modes of oscillation. Journal of Nonlinear Science 4: 419–448.

26.

Shaw

Pierre

(1991) Non-linear normal modes and invariant manifolds. Journal of Sound and Vibration 150(1): 170–173.

27.

Shaw

Pierre

(1993) Normal modes for non-linear vibratory systems. Journal of Sound and Vibration 164(1): 85–124.

28.

Shaw

Pierre

(1994) Normal modes of vibration for non-linear continuous systems. Journal of Sound and Vibration 169(3): 319–347.

29.

Sun

Fang

Soh

(2006) Thermoelastic damping in micro-beam resonators. International Journal of Solids and Structures 43: 3213–3229.

30.

Tilmans

HAC

Elwenspoek

Fluitman

JHJ

(1992) Micro resonant force gauges. Sensors and Actuators A: Physical 30: 35–53.

31.

Tilmans

HAC

Raedt

Beyne

(2003) MEMS for wireless communications: from RF-MEMS components to RF-MEMS-SiP. Journal of Micromechanics and Microengineering 13: S139–S163.

32.

Vakakis

(1991) Analysis and Identification of Linear and Nonlinear Normal Modes in Vibrating Systems. California Institute of Technology.

33.

Wong

Fox

CHJ

McWilliam

(2006) Thermoelastic damping of the in-plane vibration of thin silicon rings. Journal of Sound and Vibration 293: 266–285.

34.

Xie

Lee

Lim

(2003) Nonlinear dynamic analysis of MEMS switches by nonlinear modal analysis. Nonlinear Dynamics 31: 243–256.

35.

Younis

(2011) MEMS Linear and Nonlinear Statics and Dynamics. Springer Science+Business Media.

36.

Zener

(1937) Internal friction in solids I. Theory of internal friction in reeds. Physical Review 52: 230–235.