- Email: [email protected]
An iterative method for nonlinear dynamical system of an electrostatically actuated micro-cantilever
Physics Letters A 374 (2010) 3455–3459
Contents lists available at ScienceDirect
Physics Letters A www.elsevier.com/locate/pla
An iterative method for nonlinear dynamical system of an electrostatically actuated micro-cantilever Y.M. Chen a,∗ , G. Meng a , J.K. Liu b a b
State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China Department of Mechanics, Sun Yat-sen University, Guangzhou 510275, China
a r t i c l e
i n f o
Article history: Received 27 April 2010 Received in revised form 12 June 2010 Accepted 30 June 2010 Available online 8 July 2010 Communicated by A.R. Bishop
a b s t r a c t This study proposes an iterative method, in which all iterations are linear, for solving the nonlinear dynamical system of an electrostatically actuated micro-cantilever in MEMS. The presented method is further improved, so that the CPU time needed is significantly reduced. This method provides solutions in excellent agreement with numerical ones. © 2010 Elsevier B.V. All rights reserved.
Keywords: MEMS Micro-cantilever Nonlinear Iterative method
1. Introduction Most micro-electro-mechanical (MEMS) inherently contains nonlinearities, such as intrinsic and exterior nonlinearities arising from coupling of different domains [1,2]. Also, there exist mechanical nonlinearities, i.e., large deformations, surface contact, creep phenomena, time-dependent masses and nonlinear damping effects, etc. [3,4]. Effective nonlinear dynamic analysis becomes an increasingly important task in MEMS research and manufacturing. The nonlinear dynamical behaviors of micro-cantilever based instrument in MEMS under various loading conditions have stimulated the curiosities and interests of many researchers [5–9]. For instance, the nonlinear vibrations of an electrostatically actuated micro-cantilever based device in MEMS were investigated through a simplified mass-spring-damping model subjected to nonlinear electrostatic force [10–13]. Complex nonlinear terms resulted from electrostatic force and from squeeze film damping make it difficult to address the system directly. Zhang et al. [14] suggested an approximate approach via expanding the nonlinearities into Taylor series and retaining only the first one or two terms to simplify the considered system. The harmonic balance method was then employed to analyze the simplified system, and the responsefrequency relationship was investigated extensively. However, the attained results have not been validated by neither exact nor numerical solutions. It is worthwhile, therefore, to propose some new approaches to solve the aforementioned system.
*
Corresponding author. Tel.: +86 021 62818562; fax: +86 021 54747990. E-mail address: [email protected] (Y.M. Chen).
0375-9601/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2010.06.068
Predicting solutions of nonlinear dynamical systems via analytical or semi-analytical methods has been of fundamental interest for many years. The KBM methods [15] and the method of multiple scales (MMS) [16] are among the most widely-used ones. Both the two approaches are based on series expansions of solutions in a small parameter, i.e., the perturbation quantity. They are not only capable of providing periodic and/or limit cycle solutions, but also can analyze instantaneous vibration responses. When applying the KBM method and the MMS to obtain periodic/limit cycle solutions, the procedure of eliminating the so-called secular terms has to be implemented. Moreover, this procedure must be repeated in seeking every approximation. Often, equations deduced to eliminate secular terms are nonlinear, which can result into additional difficulty in solving them. This study aims at seeking highly accurate solutions of the nonlinear dynamical system of an electrostatically actuated microcantilever in MEMS. The above mentioned idea of expanding the nonlinearities into Taylor series [14] is adopted, in order to avoid the difficulty in handling complex nonlinearities in the frequency domain. Then, an iterative method is proposed. Numerical examples are taken to validate this approach. Periodic solutions are obtained very precisely. Also, this study analyzes the effects of the nonlinear stiffness, the squeeze film damping and the applied voltage on the amplitudes of vibration. 2. Dynamical systems The electrostatically actuated micro-cantilever in MEMS as shown in Fig. 1 is 4.5 um × 80 um × 200 um in dimensions [10].
3456
Y.M. Chen et al. / Physics Letters A 374 (2010) 3455–3459
Fig. 1. A simplified dynamical model of the micro-cantilever in MEMS.
The governing equation of motion for the micro-structure can be described as [10,14]
m y¨ + c y˙ + ky = F E (t )
(1)
where the superscript denotes the differentiation with respect to time t, y the vertical displacement of the micro-cantilever relative to the origin of the fixed plate, m the mass, k and c the effective spring stiffness and damping coefficient of the simplified system, respectively. According to the parallel plate theory, the fringe effects at the edges of the plates are ignored [17], hence F E , the electrostatic force between the capacitor plates (the fixed plate and the movable plate) generated by applying a voltage V (t ), can be expressed by
FE =
(2)
2 (d − y )2
where ε0 is the absolute dielectric constant of vacuum, ε0 = 8.5 × 10−12 N/m, A the overlapping area between the two plates, and d is the gap between them. Other parameters are given as m = 3.5 × 10−11 kg, k = 0.17 N/m, c = 1.78 × 10−6 kg/s and A = 1.6 × 10−9 m2 [14], the values of d varying. Introducing the following dimensionless variables
ω0 = k/m,
τ = ω0t
n = 1, 2, . . . .
In [14], the authors discussed the resonance responses of the original system using both the harmonic balance method and numerical integration technique to study the cases of n = 1 and n = 2. Expanding both nonlinear terms in the left and the right-hand side of (4), we construct the following iterative algorithm
x¨ n+1 + (ζ + ζ2 )˙xn+1 + xn+1
i =0
x¨ + ζ x˙ + x = T V 2 (τ )/(1 − x)2 (3) √ where ζ = c / km, T = ε0 A /(2kd3 ). In Eq. (3) and the following paper, the superscript (.) denotes the differentiation with respect
τ.
When the applied voltage, V (τ ), includes alternating current (ac) voltage, V 0 cos(ωτ ), and polarization voltage, V P , V (τ ) = V P + V 0 cos(ωτ ) can be considered for analyzing the nonlinearities of the system. In addition, MEMS devices are often characterized by structures that are a few microns in size, separated by micronsized gaps. At these sizes, air damping dominates over other dissipation. Squeeze film damping may be used to represent the air damping experienced by the moving plates [18]. When considering a cubic nonlinear spring stiffness term k1 x3 and the effect of squeeze film damping, system (3) is given as
ζ2 T [ V P + V 0 cos(ωτ )]2 x¨ + ζ + x˙ + k1 x3 + x = 3 (1 − x) (1 − x)2
(5)
n +1 2 = T V P + V 0 cos(ωτ ) (i + 1)xni
then one can rewrite Eqs. (1) and (2) as
to
n 2 ( i + 1) x i ,
x¨ + ζ x˙ + x = T V P + V 0 cos(ωτ )
i =0
ε0 A V 2 (t )
x = y /d,
Fig. 2. Comparisons between numerical solutions and results provided by the presented algorithm, where heavy dots denote the former and solid lines the latter. Parameter values are given as V P = 0.5, V 0 = 3, d = 1.2 μm, k3 = c 2 = 0 and ω = 1.
(4)
√
where ζ2 = c 2 /(d3 km ) is a non-dimensional scale, c 2 is the squeeze film damping (Pa μm4 s). 3. An iterative method In order to avoid the difficulty in directly dealing with the nonlinear terms in Eq. (3), i.e., 1/(1 − x)2 , Zhang et al. [14] expanded it as Maclaurin series, so that the harmonic balance method can be applied. Without considering the nonlinear spring or the squeeze film damping, retaining a given number of terms under the assumption that |x| < 1 results into the following modified systems
− x˙ n
n +1
αi xni − βn xn3 , n = 0, 1, 2, . . .
(6)
i =0
with a starting solution governed by
2
x¨ 0 + (ζ + ζ2 )˙x0 + x0 = T V P + V 0 cos(ωτ )
(7)
where
αi = βi =
i = 0, ζn (i + 1)(i + 2)/2, i = 1, 2, . . . ,
0,
0, k1 ,
i = 0, 1 , 2 , i = 3, 4, . . . .
Since Eqs. (6) and (7) are both linear, it is rather easy to practice the presented iterative algorithm. A conjecture is drawn without proof, that as long as series {xi , i = 0, 1, 2, . . .} converges, it must converge to one solution of Eq. (4) if |x| < 1 satisfies. Using the iterative algorithm proposed above, theoretically, analytical solutions can be obtained because all iterative equations are linear. However, they can become rapidly complicated as iteration proceeding. Fortunately, it is rather easy to obtain semi-analytical solutions. The attained results attained by the presented algorithm are validated through comparison with numerical solutions. Note that, all the numerical solutions in this study are obtained via directly integrating Eq. (4) using the fourth-order Runge–Kutta method, with the initial values as x(0) = 0, x˙ (0) = 0. Fig. 2 shows a phase plane for system (3), namely without nonlinear spring or squeeze film damping considered. The solutions provided by the
Y.M. Chen et al. / Physics Letters A 374 (2010) 3455–3459
Fig. 3. Comparisons between numerical solutions and results provided by the presented algorithm, where heavy dots denote the former and solid lines the latter. Parameter values are given as V P = 1, V 0 = 5, d = 1.5 μm, k3 = 0.5, c 2 = 0 and ω = 1.
3457
Fig. 5. Amplitudes of harmonics in iteration equation. Parameter values are given as V P = 1, V 0 = 2, d = 1 μm, k3 = 0.5, c 2 = 2 and ω = 1.
to excellent precision. Those imply feasibility and potential in implementing the method to analyze nonlinear dynamical behaviors. 4. An improved approach for saving computational effort Due to the high-order power like xnn included in the nth iteration equation, the number of harmonics in the iteration approximations increases acceleratedly. Denote the number of harmonics in xn as Γn , then we can estimate the counterpart for xn+1 as Γn+1 = (n + 1)Γn . The starting solution x0 provided by Eq. (7) has the first two harmonics (Γ0 = 2), hence Γ1 = 4, Γ2 = 12, Γ3 = 48, Γ4 = 240 and Γn = 2(n + 1)!, etc. In most cases, vibrations, no matter how complex they are, are dominated by several lower harmonics. Therefore, the harmonics accumulate so rapidly that it seems to be unworthy to improve the accuracy slightly with such a vast increasing of computational cost. The iteration solutions can be obtained with a given number (N) of harmonics while neglecting all higher-order harmonics than the Nth. Denote the (n + 1)th-order approximation as Fig. 4. Comparisons between numerical solutions and results provided by the presented algorithm, where heavy dots denote the former and solid lines the latter. Parameter values are given as V P = 0.5, V 0 = 2, d = 1 μm, k3 = 0.5, c 2 = 2 and ω = 0.5.
presented method converge to the numerical one as n increasing, so that the 10th-order approximation is nearly the same as the numerical one. Also observed is the rapid convergence of the approximations to the numerical solution. After less than 10 steps are performed, iteration solutions can be provided very accurately. In fact, solutions can be obtained to any desired accuracy without additional difficulty. When the nonlinear spring and/or squeeze film damping are taken into account, the proposed algorithm can still provide very accurate solutions, as shown in Figs. 3 and 4. Both the figures show clearly rapid convergence rates of the approximations as n increasing. Actually, the iteration solutions are rather accurate even after 5 steps. In addition, the vibrations denoted by these phase planes are nonharmonic, namely the solutions cannot be expressed by the first harmonic alone with a reasonable precision. That is probably because high harmonics contribute to the responses to an innegligible extent. The phase plane shown in Fig. 3 even has a loop. Even though, the presented iterative method can still track them
xn+1 =
M
cni +1 cos(ωτ ) + sni +1 sin(ωτ )
i =0
where M denotes the number of harmonics. Once M > N, retain the first to the Nth harmonics and rewrite xn+1 as
xn+1 =
N
cni +1 cos(ωτ ) + sni +1 sin(ωτ ) .
(8)
i =0
Interestingly, as long as a solution provided by algorithms (6) and (8) together is convergent, it must converge to one harmonic balance solution. To test this, letting xn = xn+1 and substituting Eq. (8) into (4), we have
S n = x¨ n+1 + (ζ + ζ2 )˙xn+1 + xn+1 n +1 2 − T V P + V 0 cos(ωτ ) (i + 1)xni +1 i =0
− x˙ n
n +1 i =0
ai xni +1 − βn xn3+1
3458
Y.M. Chen et al. / Physics Letters A 374 (2010) 3455–3459
Fig. 6. Ratios of CPU time and of amplitude, for the algorithm without neglecting harmonics and that with 5 harmonics retained. Parameter values are given as V P = 1, V 0 = 2, d0 = 1 μm, k3 = 0.5, c 2 = 2 and ω = 1.
Fig. 7. The amplitudes of solutions, where heavy dots denote numerical solutions and solid lines the 20th-order approximations provided by the improved algorithm when retaining 10 harmonics. Parameter values are given as V P = 1, d0 = 2, and k3 = c 2 = 0.
N (n+1)
=
C in cos(i ωτ ) + S ni sin(i ωτ ) .
i =0
Fig. 5 shows the convergence (vs n) of harmonic amplitudes of in the iteration equations. It is noticeable that the summation of the absolute amplitudes of the first 5 harmonics converges to 0 by a logarithmic rate, whereas the amplitudes of the 6th and of the 7th harmonics approach constants, respectively, though very small they are. It is implied that the attained solutions comply with the harmonic balancing principle. In other words, the attained iteration solutions converge to one of harmonic balance solutions. Denote the amplitudes as an of iteration solutions (xn ) without neglecting harmonics, and its counterpart as a¯ n for algorithm with N harmonics retained. Also, assuming the respective CPU time needed in solving xn as T n and T¯ n , respectively. The ratios vs n of time and of amplitude are plotted in Fig. 6. As one can see, as n increases, the ratio of amplitudes approaches 1 very quickly. That means the two types of solutions almost have the same pre-
Fig. 8. The amplitudes of solutions, where heavy dots denote numerical solutions and solid lines the 20th-order approximations provided by the improved algorithm when retaining 10 harmonics. Parameter values are given as V P = 1, V 0 = 3, d0 = 1.5, and c 2 = 0.
Fig. 9. The amplitudes of solutions, where heavy dots denote numerical solutions and solid lines the 20th-order approximations provided by the improved algorithm when retaining 10 harmonics. Parameter values are given as V P = 0.5, V 0 = 2, d0 = 1 and k3 = 0.5.
cision. As for the CPU time, however, there is a large discrepancy between the one for (6) alone and that for (6) and (8) together. The time spent on solving x7 without neglecting higher harmonics is more than 105 times of that with 5 low harmonics retained. Furthermore, the ratio can increase even more rapidly as n increasing. Accordingly, it is very worthwhile using the proposed approach for saving computational effort. The solutions are obtained for different values of ω and plotted in Figs. 7–9. Also plotted are the corresponding numerical solutions. All amplitude–frequency curves show the fine agreement between the solutions attained by the improved algorithm and those by numerical method. As is shown in these figures, the amplitudes of vibration decrease when ω increases from 0. Noticeably, they decrease much more rapidly when ω < 1 than when ω > 1. The amplitudes are positively related to the applied ac voltage V 0 , namely the larger V 0 is the bigger the amplitude becomes.
Y.M. Chen et al. / Physics Letters A 374 (2010) 3455–3459
3459
That is because the increasing of V 0 results in the growing of the external forces. On the other hand, the amplitudes are negatively related to both the nonlinear stiffness coefficient and to the squeeze film damping coefficient, respectively. That means positive damping and stronger nonlinear restoring forces can suppress the vibrations of the considered micro-structure.
(iii) When the KBM method, the MMS and the Mickens iteration method are used to obtain periodic solutions, an additional procedure must be repeatedly adopted to eliminate secular terms if higher approximations are needed. In the presented approach, this procedure is avoided as no secular terms will arise in all iteration steps.
5. Conclusions and remarks
Acknowledgements
By expanding the nonlinear terms into series step by step, we have proposed a new iterative algorithm for solving the nonlinear dynamical system of an electrostatically actuated micro-cantilever. Numerical examples validate the proposed approach, with excellent approximations for periodic solutions are obtained. This approach was further improved by restricting the number of harmonics in every iteration solutions as a given number N. It turns out that the improved approach can save a large amount of computational effort, at the same time can keep a high accuracy. Also, the frequency-response curves are obtained and addressed detailly. The feasibility and efficiency, demonstrated by numerical examples, imply that we could expect the presented approach can be applicable in more nonlinear vibration problems, especially those arising in MEMS. The iteration algorithm has been introduced to solve nonlinear dynamical systems by Mickens [19] since 1980s. Compared the presented approach with the KBM method and the MMS, as well as the Mickens iteration method [19], several remarks are listed below for a better understanding.
This work is supported by the National Natural Science Foundation of China (10772202, 10102023, 10972241), Doctoral Program Foundation of Ministry of Education of China (20090171110042), China Postdoctoral Science Foundation (20090460628) and Shanghai Postdoctoral Science Foundation (10R21413600).
(i) The presented approach and the Mickens iteration method are only effective for attaining periodic solutions. Both the KBM method and the MMS can be used to seek not only periodic solutions but also instantaneous ones. (ii) The KBM method and the MMS require the perturbation quantity, which should be a small parameter. The iteration methods can be applicable in nonlinear dynamical systems without any small/large parameter.
References [1] S.D. Senturia, Sens. Actuators A 67 (1998) 1. [2] S.E. Lyshevski, in: Proceeding of the 40th IEEE Conference on Decision and Control Orlando, FL, USA, 1997, pp. 4681–4686. [3] S.K. De, N.R. Aluru, J. Microelectromech. Syst. 15 (2006) 355. [4] M. Ashhab, M.V. Salapaka, M. Dahleh, I. Mezic, Automatica 35 (1999) 1663. [5] A. Passiana, G. Muralidharana, A. Mehtaa, H. Simpsonb, T.L. Ferrella, T. Thundat, Ultramicroscopy 97 (2003) 391. [6] Y. Ahn, H. Guckel, J.D. Zook, J. Micromech. Microeng. 11 (2001) 70. [7] E.K. Chan, R.W. Dutton, J. Microelectromech. Syst. 9 (2000) 321. [8] S.N. Mahmoodi, N. Jalili, Sens. Actuators A 150 (2009) 131. [9] Y.M. Fu, J. Zhang, Acta Mech. Sin. 25 (2009) 211. [10] S. Liu, A. Davidson, Q. Liu, in: Proceedings of the IEEE, Transducers’03, the 12th International Conference on Solid State Sensors Actuators and Microsystems, Boston, 8–12 June 2003, pp. 1092–1095. [11] S. Liu, A. Davidson, Q. Liu, J. Micromech. Microeng. 14 (2004) 1064. [12] W.M. Zhang, G. Meng, D. Chen, Sensors 7 (2007) 760. [13] G. Meng, W.M. Zhang, H. Huang, H.G. Li, D. Chen, Chaos Solitons Fractals 40 (2009) 538. [14] W.M. Zhang, G. Meng, Sens. Actuators A 119 (2005) 291. [15] N.N. Bogoliubov, Y.A. Mitropolsky, Asymptotic Methods in the Theory of Nonlinear Oscillations, Gordon and Breach Inc., New York, 1961. [16] A.H. Nayfeh, D.T. Mook, John Wiley & Sons, New York, 1979. [17] R.H. Price, J.E. Wood, S.C. Jacobsen, Sens. Actuators A 20 (1989) 107. [18] T. Veijola, H. Kuisma, J. Lahdenpera, Sens. Actuators A 66 (1998) 83. [19] R.E. Mickens, J. Sound Vib. 116 (1987) 185.
Contents lists available at ScienceDirect
Physics Letters A www.elsevier.com/locate/pla
An iterative method for nonlinear dynamical system of an electrostatically actuated micro-cantilever Y.M. Chen a,∗ , G. Meng a , J.K. Liu b a b
State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China Department of Mechanics, Sun Yat-sen University, Guangzhou 510275, China
a r t i c l e
i n f o
Article history: Received 27 April 2010 Received in revised form 12 June 2010 Accepted 30 June 2010 Available online 8 July 2010 Communicated by A.R. Bishop
a b s t r a c t This study proposes an iterative method, in which all iterations are linear, for solving the nonlinear dynamical system of an electrostatically actuated micro-cantilever in MEMS. The presented method is further improved, so that the CPU time needed is significantly reduced. This method provides solutions in excellent agreement with numerical ones. © 2010 Elsevier B.V. All rights reserved.
Keywords: MEMS Micro-cantilever Nonlinear Iterative method
1. Introduction Most micro-electro-mechanical (MEMS) inherently contains nonlinearities, such as intrinsic and exterior nonlinearities arising from coupling of different domains [1,2]. Also, there exist mechanical nonlinearities, i.e., large deformations, surface contact, creep phenomena, time-dependent masses and nonlinear damping effects, etc. [3,4]. Effective nonlinear dynamic analysis becomes an increasingly important task in MEMS research and manufacturing. The nonlinear dynamical behaviors of micro-cantilever based instrument in MEMS under various loading conditions have stimulated the curiosities and interests of many researchers [5–9]. For instance, the nonlinear vibrations of an electrostatically actuated micro-cantilever based device in MEMS were investigated through a simplified mass-spring-damping model subjected to nonlinear electrostatic force [10–13]. Complex nonlinear terms resulted from electrostatic force and from squeeze film damping make it difficult to address the system directly. Zhang et al. [14] suggested an approximate approach via expanding the nonlinearities into Taylor series and retaining only the first one or two terms to simplify the considered system. The harmonic balance method was then employed to analyze the simplified system, and the responsefrequency relationship was investigated extensively. However, the attained results have not been validated by neither exact nor numerical solutions. It is worthwhile, therefore, to propose some new approaches to solve the aforementioned system.
*
Corresponding author. Tel.: +86 021 62818562; fax: +86 021 54747990. E-mail address: [email protected] (Y.M. Chen).
0375-9601/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2010.06.068
Predicting solutions of nonlinear dynamical systems via analytical or semi-analytical methods has been of fundamental interest for many years. The KBM methods [15] and the method of multiple scales (MMS) [16] are among the most widely-used ones. Both the two approaches are based on series expansions of solutions in a small parameter, i.e., the perturbation quantity. They are not only capable of providing periodic and/or limit cycle solutions, but also can analyze instantaneous vibration responses. When applying the KBM method and the MMS to obtain periodic/limit cycle solutions, the procedure of eliminating the so-called secular terms has to be implemented. Moreover, this procedure must be repeated in seeking every approximation. Often, equations deduced to eliminate secular terms are nonlinear, which can result into additional difficulty in solving them. This study aims at seeking highly accurate solutions of the nonlinear dynamical system of an electrostatically actuated microcantilever in MEMS. The above mentioned idea of expanding the nonlinearities into Taylor series [14] is adopted, in order to avoid the difficulty in handling complex nonlinearities in the frequency domain. Then, an iterative method is proposed. Numerical examples are taken to validate this approach. Periodic solutions are obtained very precisely. Also, this study analyzes the effects of the nonlinear stiffness, the squeeze film damping and the applied voltage on the amplitudes of vibration. 2. Dynamical systems The electrostatically actuated micro-cantilever in MEMS as shown in Fig. 1 is 4.5 um × 80 um × 200 um in dimensions [10].
3456
Y.M. Chen et al. / Physics Letters A 374 (2010) 3455–3459
Fig. 1. A simplified dynamical model of the micro-cantilever in MEMS.
The governing equation of motion for the micro-structure can be described as [10,14]
m y¨ + c y˙ + ky = F E (t )
(1)
where the superscript denotes the differentiation with respect to time t, y the vertical displacement of the micro-cantilever relative to the origin of the fixed plate, m the mass, k and c the effective spring stiffness and damping coefficient of the simplified system, respectively. According to the parallel plate theory, the fringe effects at the edges of the plates are ignored [17], hence F E , the electrostatic force between the capacitor plates (the fixed plate and the movable plate) generated by applying a voltage V (t ), can be expressed by
FE =
(2)
2 (d − y )2
where ε0 is the absolute dielectric constant of vacuum, ε0 = 8.5 × 10−12 N/m, A the overlapping area between the two plates, and d is the gap between them. Other parameters are given as m = 3.5 × 10−11 kg, k = 0.17 N/m, c = 1.78 × 10−6 kg/s and A = 1.6 × 10−9 m2 [14], the values of d varying. Introducing the following dimensionless variables
ω0 = k/m,
τ = ω0t
n = 1, 2, . . . .
In [14], the authors discussed the resonance responses of the original system using both the harmonic balance method and numerical integration technique to study the cases of n = 1 and n = 2. Expanding both nonlinear terms in the left and the right-hand side of (4), we construct the following iterative algorithm
x¨ n+1 + (ζ + ζ2 )˙xn+1 + xn+1
i =0
x¨ + ζ x˙ + x = T V 2 (τ )/(1 − x)2 (3) √ where ζ = c / km, T = ε0 A /(2kd3 ). In Eq. (3) and the following paper, the superscript (.) denotes the differentiation with respect
τ.
When the applied voltage, V (τ ), includes alternating current (ac) voltage, V 0 cos(ωτ ), and polarization voltage, V P , V (τ ) = V P + V 0 cos(ωτ ) can be considered for analyzing the nonlinearities of the system. In addition, MEMS devices are often characterized by structures that are a few microns in size, separated by micronsized gaps. At these sizes, air damping dominates over other dissipation. Squeeze film damping may be used to represent the air damping experienced by the moving plates [18]. When considering a cubic nonlinear spring stiffness term k1 x3 and the effect of squeeze film damping, system (3) is given as
ζ2 T [ V P + V 0 cos(ωτ )]2 x¨ + ζ + x˙ + k1 x3 + x = 3 (1 − x) (1 − x)2
(5)
n +1 2 = T V P + V 0 cos(ωτ ) (i + 1)xni
then one can rewrite Eqs. (1) and (2) as
to
n 2 ( i + 1) x i ,
x¨ + ζ x˙ + x = T V P + V 0 cos(ωτ )
i =0
ε0 A V 2 (t )
x = y /d,
Fig. 2. Comparisons between numerical solutions and results provided by the presented algorithm, where heavy dots denote the former and solid lines the latter. Parameter values are given as V P = 0.5, V 0 = 3, d = 1.2 μm, k3 = c 2 = 0 and ω = 1.
(4)
√
where ζ2 = c 2 /(d3 km ) is a non-dimensional scale, c 2 is the squeeze film damping (Pa μm4 s). 3. An iterative method In order to avoid the difficulty in directly dealing with the nonlinear terms in Eq. (3), i.e., 1/(1 − x)2 , Zhang et al. [14] expanded it as Maclaurin series, so that the harmonic balance method can be applied. Without considering the nonlinear spring or the squeeze film damping, retaining a given number of terms under the assumption that |x| < 1 results into the following modified systems
− x˙ n
n +1
αi xni − βn xn3 , n = 0, 1, 2, . . .
(6)
i =0
with a starting solution governed by
2
x¨ 0 + (ζ + ζ2 )˙x0 + x0 = T V P + V 0 cos(ωτ )
(7)
where
αi = βi =
i = 0, ζn (i + 1)(i + 2)/2, i = 1, 2, . . . ,
0,
0, k1 ,
i = 0, 1 , 2 , i = 3, 4, . . . .
Since Eqs. (6) and (7) are both linear, it is rather easy to practice the presented iterative algorithm. A conjecture is drawn without proof, that as long as series {xi , i = 0, 1, 2, . . .} converges, it must converge to one solution of Eq. (4) if |x| < 1 satisfies. Using the iterative algorithm proposed above, theoretically, analytical solutions can be obtained because all iterative equations are linear. However, they can become rapidly complicated as iteration proceeding. Fortunately, it is rather easy to obtain semi-analytical solutions. The attained results attained by the presented algorithm are validated through comparison with numerical solutions. Note that, all the numerical solutions in this study are obtained via directly integrating Eq. (4) using the fourth-order Runge–Kutta method, with the initial values as x(0) = 0, x˙ (0) = 0. Fig. 2 shows a phase plane for system (3), namely without nonlinear spring or squeeze film damping considered. The solutions provided by the
Y.M. Chen et al. / Physics Letters A 374 (2010) 3455–3459
Fig. 3. Comparisons between numerical solutions and results provided by the presented algorithm, where heavy dots denote the former and solid lines the latter. Parameter values are given as V P = 1, V 0 = 5, d = 1.5 μm, k3 = 0.5, c 2 = 0 and ω = 1.
3457
Fig. 5. Amplitudes of harmonics in iteration equation. Parameter values are given as V P = 1, V 0 = 2, d = 1 μm, k3 = 0.5, c 2 = 2 and ω = 1.
to excellent precision. Those imply feasibility and potential in implementing the method to analyze nonlinear dynamical behaviors. 4. An improved approach for saving computational effort Due to the high-order power like xnn included in the nth iteration equation, the number of harmonics in the iteration approximations increases acceleratedly. Denote the number of harmonics in xn as Γn , then we can estimate the counterpart for xn+1 as Γn+1 = (n + 1)Γn . The starting solution x0 provided by Eq. (7) has the first two harmonics (Γ0 = 2), hence Γ1 = 4, Γ2 = 12, Γ3 = 48, Γ4 = 240 and Γn = 2(n + 1)!, etc. In most cases, vibrations, no matter how complex they are, are dominated by several lower harmonics. Therefore, the harmonics accumulate so rapidly that it seems to be unworthy to improve the accuracy slightly with such a vast increasing of computational cost. The iteration solutions can be obtained with a given number (N) of harmonics while neglecting all higher-order harmonics than the Nth. Denote the (n + 1)th-order approximation as Fig. 4. Comparisons between numerical solutions and results provided by the presented algorithm, where heavy dots denote the former and solid lines the latter. Parameter values are given as V P = 0.5, V 0 = 2, d = 1 μm, k3 = 0.5, c 2 = 2 and ω = 0.5.
presented method converge to the numerical one as n increasing, so that the 10th-order approximation is nearly the same as the numerical one. Also observed is the rapid convergence of the approximations to the numerical solution. After less than 10 steps are performed, iteration solutions can be provided very accurately. In fact, solutions can be obtained to any desired accuracy without additional difficulty. When the nonlinear spring and/or squeeze film damping are taken into account, the proposed algorithm can still provide very accurate solutions, as shown in Figs. 3 and 4. Both the figures show clearly rapid convergence rates of the approximations as n increasing. Actually, the iteration solutions are rather accurate even after 5 steps. In addition, the vibrations denoted by these phase planes are nonharmonic, namely the solutions cannot be expressed by the first harmonic alone with a reasonable precision. That is probably because high harmonics contribute to the responses to an innegligible extent. The phase plane shown in Fig. 3 even has a loop. Even though, the presented iterative method can still track them
xn+1 =
M
cni +1 cos(ωτ ) + sni +1 sin(ωτ )
i =0
where M denotes the number of harmonics. Once M > N, retain the first to the Nth harmonics and rewrite xn+1 as
xn+1 =
N
cni +1 cos(ωτ ) + sni +1 sin(ωτ ) .
(8)
i =0
Interestingly, as long as a solution provided by algorithms (6) and (8) together is convergent, it must converge to one harmonic balance solution. To test this, letting xn = xn+1 and substituting Eq. (8) into (4), we have
S n = x¨ n+1 + (ζ + ζ2 )˙xn+1 + xn+1 n +1 2 − T V P + V 0 cos(ωτ ) (i + 1)xni +1 i =0
− x˙ n
n +1 i =0
ai xni +1 − βn xn3+1
3458
Y.M. Chen et al. / Physics Letters A 374 (2010) 3455–3459
Fig. 6. Ratios of CPU time and of amplitude, for the algorithm without neglecting harmonics and that with 5 harmonics retained. Parameter values are given as V P = 1, V 0 = 2, d0 = 1 μm, k3 = 0.5, c 2 = 2 and ω = 1.
Fig. 7. The amplitudes of solutions, where heavy dots denote numerical solutions and solid lines the 20th-order approximations provided by the improved algorithm when retaining 10 harmonics. Parameter values are given as V P = 1, d0 = 2, and k3 = c 2 = 0.
N (n+1)
=
C in cos(i ωτ ) + S ni sin(i ωτ ) .
i =0
Fig. 5 shows the convergence (vs n) of harmonic amplitudes of in the iteration equations. It is noticeable that the summation of the absolute amplitudes of the first 5 harmonics converges to 0 by a logarithmic rate, whereas the amplitudes of the 6th and of the 7th harmonics approach constants, respectively, though very small they are. It is implied that the attained solutions comply with the harmonic balancing principle. In other words, the attained iteration solutions converge to one of harmonic balance solutions. Denote the amplitudes as an of iteration solutions (xn ) without neglecting harmonics, and its counterpart as a¯ n for algorithm with N harmonics retained. Also, assuming the respective CPU time needed in solving xn as T n and T¯ n , respectively. The ratios vs n of time and of amplitude are plotted in Fig. 6. As one can see, as n increases, the ratio of amplitudes approaches 1 very quickly. That means the two types of solutions almost have the same pre-
Fig. 8. The amplitudes of solutions, where heavy dots denote numerical solutions and solid lines the 20th-order approximations provided by the improved algorithm when retaining 10 harmonics. Parameter values are given as V P = 1, V 0 = 3, d0 = 1.5, and c 2 = 0.
Fig. 9. The amplitudes of solutions, where heavy dots denote numerical solutions and solid lines the 20th-order approximations provided by the improved algorithm when retaining 10 harmonics. Parameter values are given as V P = 0.5, V 0 = 2, d0 = 1 and k3 = 0.5.
cision. As for the CPU time, however, there is a large discrepancy between the one for (6) alone and that for (6) and (8) together. The time spent on solving x7 without neglecting higher harmonics is more than 105 times of that with 5 low harmonics retained. Furthermore, the ratio can increase even more rapidly as n increasing. Accordingly, it is very worthwhile using the proposed approach for saving computational effort. The solutions are obtained for different values of ω and plotted in Figs. 7–9. Also plotted are the corresponding numerical solutions. All amplitude–frequency curves show the fine agreement between the solutions attained by the improved algorithm and those by numerical method. As is shown in these figures, the amplitudes of vibration decrease when ω increases from 0. Noticeably, they decrease much more rapidly when ω < 1 than when ω > 1. The amplitudes are positively related to the applied ac voltage V 0 , namely the larger V 0 is the bigger the amplitude becomes.
Y.M. Chen et al. / Physics Letters A 374 (2010) 3455–3459
3459
That is because the increasing of V 0 results in the growing of the external forces. On the other hand, the amplitudes are negatively related to both the nonlinear stiffness coefficient and to the squeeze film damping coefficient, respectively. That means positive damping and stronger nonlinear restoring forces can suppress the vibrations of the considered micro-structure.
(iii) When the KBM method, the MMS and the Mickens iteration method are used to obtain periodic solutions, an additional procedure must be repeatedly adopted to eliminate secular terms if higher approximations are needed. In the presented approach, this procedure is avoided as no secular terms will arise in all iteration steps.
5. Conclusions and remarks
Acknowledgements
By expanding the nonlinear terms into series step by step, we have proposed a new iterative algorithm for solving the nonlinear dynamical system of an electrostatically actuated micro-cantilever. Numerical examples validate the proposed approach, with excellent approximations for periodic solutions are obtained. This approach was further improved by restricting the number of harmonics in every iteration solutions as a given number N. It turns out that the improved approach can save a large amount of computational effort, at the same time can keep a high accuracy. Also, the frequency-response curves are obtained and addressed detailly. The feasibility and efficiency, demonstrated by numerical examples, imply that we could expect the presented approach can be applicable in more nonlinear vibration problems, especially those arising in MEMS. The iteration algorithm has been introduced to solve nonlinear dynamical systems by Mickens [19] since 1980s. Compared the presented approach with the KBM method and the MMS, as well as the Mickens iteration method [19], several remarks are listed below for a better understanding.
This work is supported by the National Natural Science Foundation of China (10772202, 10102023, 10972241), Doctoral Program Foundation of Ministry of Education of China (20090171110042), China Postdoctoral Science Foundation (20090460628) and Shanghai Postdoctoral Science Foundation (10R21413600).
(i) The presented approach and the Mickens iteration method are only effective for attaining periodic solutions. Both the KBM method and the MMS can be used to seek not only periodic solutions but also instantaneous ones. (ii) The KBM method and the MMS require the perturbation quantity, which should be a small parameter. The iteration methods can be applicable in nonlinear dynamical systems without any small/large parameter.
References [1] S.D. Senturia, Sens. Actuators A 67 (1998) 1. [2] S.E. Lyshevski, in: Proceeding of the 40th IEEE Conference on Decision and Control Orlando, FL, USA, 1997, pp. 4681–4686. [3] S.K. De, N.R. Aluru, J. Microelectromech. Syst. 15 (2006) 355. [4] M. Ashhab, M.V. Salapaka, M. Dahleh, I. Mezic, Automatica 35 (1999) 1663. [5] A. Passiana, G. Muralidharana, A. Mehtaa, H. Simpsonb, T.L. Ferrella, T. Thundat, Ultramicroscopy 97 (2003) 391. [6] Y. Ahn, H. Guckel, J.D. Zook, J. Micromech. Microeng. 11 (2001) 70. [7] E.K. Chan, R.W. Dutton, J. Microelectromech. Syst. 9 (2000) 321. [8] S.N. Mahmoodi, N. Jalili, Sens. Actuators A 150 (2009) 131. [9] Y.M. Fu, J. Zhang, Acta Mech. Sin. 25 (2009) 211. [10] S. Liu, A. Davidson, Q. Liu, in: Proceedings of the IEEE, Transducers’03, the 12th International Conference on Solid State Sensors Actuators and Microsystems, Boston, 8–12 June 2003, pp. 1092–1095. [11] S. Liu, A. Davidson, Q. Liu, J. Micromech. Microeng. 14 (2004) 1064. [12] W.M. Zhang, G. Meng, D. Chen, Sensors 7 (2007) 760. [13] G. Meng, W.M. Zhang, H. Huang, H.G. Li, D. Chen, Chaos Solitons Fractals 40 (2009) 538. [14] W.M. Zhang, G. Meng, Sens. Actuators A 119 (2005) 291. [15] N.N. Bogoliubov, Y.A. Mitropolsky, Asymptotic Methods in the Theory of Nonlinear Oscillations, Gordon and Breach Inc., New York, 1961. [16] A.H. Nayfeh, D.T. Mook, John Wiley & Sons, New York, 1979. [17] R.H. Price, J.E. Wood, S.C. Jacobsen, Sens. Actuators A 20 (1989) 107. [18] T. Veijola, H. Kuisma, J. Lahdenpera, Sens. Actuators A 66 (1998) 83. [19] R.E. Mickens, J. Sound Vib. 116 (1987) 185.