Homotopy Analysis Method for Nonlinear Dynamical System of an Electrostatically Actuated Microcantilever

The homotopy analysis method HAM is employed to propose an approach for solving the nonlinear dynamical system of an electrostatically actuated micro-cantilever in MEMS. There are two relative merits of the presented HAM compared with some usual procedures of the HAM. First, a new auxiliary linear operator is constructed. This operator makes it unnecessary to eliminate any secular terms. Furthermore, all the deformation equations are purely linear. Numerical examples show the excellent agreement of the attained solutions with numerical ones. The respective effects of applied voltage, cubic nonlinear stiffness, gap distance, and squeeze film damping on vibration responses are analyzed detailedly.


Introduction
Over the last couple of decades, Liao 1-4 described a nonlinear analytical technique which does not require small parameters and thus can be applied to solve nonlinear problems without small or large parameters.This technique is based on homotopy, which is an important part of topology, called the homotopy analysis method HAM .Its main idea is to construct a class of homotopy in a rather general form by introducing an auxiliary parameter.This parameter can provide us with a convenient way to control the convergence of approximation series and adjust convergence rate and region when necessary.A systematical description of this method was presented in 5 .Liao 5 also studied the convergence properties of the HAM and proved that as long as an HAM series is convergent, it must converge to one solution of the considered problem.

Dynamical System
The electrostatically actuated microcantilever in MEMS as shown in Figure 1  where the superscript denotes the differentiation with respect to time t, y the vertical displacement of the microcantilever 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, if the fringe effects at the edges of the plates are ignored 21 , the electrostatic force F E generated by applying a voltage V t between the capacitor plates the fixed plate and the movable plate can be expressed by 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 μm is the gap between them.The other parameter values are given as m 3.5 × 10 −11 kg, and k 0.17 N/m, c 1.78 × 10 −6 kg/s and A 1.6 × 10 −9 m 2 20 , while the value of d varies.
Introducing the following dimensionless variables then one can rewrite 2.1 and 2.2 as where ζ c/ √ km, T ε 0 A/ 2kd 3 , and the superscript • refers to the differentiation with respect to t 1 .
When the applied voltage, V t 1 , includes alternating current ac voltage, V 0 cos ωt 1 , and polarization voltage, V P , V t 1 V P V 0 cos ωt 1 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 micron-sized 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 22 .When considering the effect of squeeze film damping and a cubic nonlinear spring stiffness k 3 x 3 , system 2.4 is given as where km is a nondimensional parameter, and c 2 is the squeeze film damping Pa•μm 4 •s .

Mathematical Formulation
In this section, the HAM is applied to solve 2.5 .First of all, introducing a new time scale τ ωt 1 , and multiplying 2.5 by 1 − x 3 , one can rewrite 2.5 as where the superscript denotes the differentiation with respect to τ.Compared with 2.5 , it is very easy to expand 3.1 once x is expressed as an HAM series.

Basic Idea of Homotopy Analysis Method
In the frame of the HAM, the first procedure is to choose an auxiliary linear operator and a nonlinear operator, respectively.A new equation, usually called as the zeroth-order deformation equation, is constructed as where L denotes the chosen auxiliary operator, N the nonlinear one, p ∈ 0, 1 the embedding parameter and h a nonzero auxiliary parameter.The HAM is based on the continuous variation of u τ, p .When the embedding parameter p increases from 0 to 1, u τ, p varies from an initial guess to one exact solution of the considered problem N u τ, p 0.

3.3
According to 3.1 , it can be given as

3.4
When p 0, 3.2 becomes a linear one, whose exact solution is evident.When p 1, it is the same as 3.1 provided that x τ u τ, 1 .Expanding u τ, p into a Taylor series u τ, p ∞ i 0 u i τ p i .

3.5
Note that u i τ is dependent upon h.In fact, the auxiliary parameter h is introduced to make it convenient to control and adjust the convergence of the HAM series i.e., 3.5 .As long as series 3.5 is convergent at p 1, one can obtain the mth-order HAM solution as x τ m i 0 u i τ .

3.6
Substituting 3.4 and 3.5 into 3.2 and equating the coefficients of p i to 0, one obtains where 3.9 when n 0, ϑ n 1, otherwise ϑ n 0; z i corresponds to the expansion of 1 − u 3 , that is, 1, 2, . ... Equations 3.7 and 3.8 are always linear for every positive integer n.Importantly, the nonlinear equation has so far been successfully transformed into a series of linear ones.These linear equations can be solved step by step.That is core idea and also the most outstanding advantage of the HAM.

A Widely Used Linear Operator
One of key procedures of the HAM is to choose an appropriate linear operator.According to Liao 5 , it should be chosen on the basis of the so-called rule of solution expression.Periodic/limit cycle solutions can be expressed as Fourier series in τ Note that the HAM has been widely used to obtain periodic and/or limit cycle solutions of nonlinear oscillators 23-29 .Particularly, this technique has found applications in microsystems 28 and microbeams subjected applied voltages 29 , respectively.To the best of out knowledge, most researchers 23-29 adopted such a following linear operator as If 3.11 is employed in this study, 3.7 becomes one homogenous equation in u 0 .This equation has a series of solutions that has to be determined using the first-order deformation equation, that is, 3.8 with n 0. To this end, an initial guess for u 0 τ is usually given as u 0 c 0,1 cos τ s 0,1 sin τ, 3.12 where c 0,1 and s 0,1 are prior unknown constants.Due to the rule of solution expression, u n can be given as truncated Fourier series.Substituting u n into R n , one has where ϕ n refers to the highest harmonic.Considering 3.7 and 3.13 , cos τ and sin τ can result in the so-called secular terms τ sin τ and τ cos τ, respectively.These terms do not comply with the rule of solution expression.One usual procedure to eliminate the secular terms is equating both C n,1 and S n,1 to 0, that is, When n ≥ 1, these equations are all linear.When n 0, however, 3.14 is nonlinear.Actually, it is essentially the first-order harmonic balance equations 25, 26 governing c 0,1 and s 0,1 .
Well-known, it is a major shortcoming of the harmonic balance method that nonlinear algebraic equations have to be solved.

A New Linear Operator
Considering that 3.1 contains powers like x 6 , there will be c 6 0,1 and s 6 0,1 in 3.14 with n 0. These nonlinearities can make it rather difficult to determine the initial solution.It is worthy to find another linear auxiliary operator to avoid this obstacle.Consider the following linear operator

3.15
Substitution of 3.15 into 3.7 and 3.8 results into L u 0 τ − T V P V 0 cos τ 2 0, 3.16 L u n 1 − σ n u n hR n u 0 , u 1 , . . ., u n , n 0, 1, 2, . .., 3.17 where σ 0 0, otherwise σ n 1 when n ≥ 1. Equations 3.16 and 3.17 provide us with a purely linear algorithm, without an additional task of eliminating the secular terms as long as ζ ζ 2 / 0. Since 3.16 is an inhomogeneous equation whose exact solution can be determined analytically, there is also no need to choose an initial guess for u 0 .

A Means for Saving Computational Effort
The number of harmonics in the HAM approximations increases acceleratedly, because of the existence of the powers like x 6 in 3.1 .Denote the number of harmonics in u n as Γ n .The starting solution u 0 provided by 3.16 contains the first two harmonics, that is, Γ 0 2. In the equation governing u 1 , the highest powers is u 6 0 , hence one knows that Γ 1 12.As n increases further, Γ 2 22 and Γ 3 32 can be obtained.In most cases, periodic or limit cycle oscillations are dominated by several lower harmonics.High harmonics contribute to the improvement of the precision quite limitedly, compared with the lower ones.Therefore, it seems to be unworthy to improve the accuracy slightly with such a vast increasing of computational cost.
Therefore, the high-order HAM approximations are obtained with a given number N of harmonics, while neglecting all harmonics higher than the Nth one.Denote the nthorder HAM approximation as s n,i sin iτ , n 0, 1, 2, . . ., 3.18 where M increases with n increasing.Once M > N, retain the first to the Nth harmonics in u n , and rewrite it as c n,i cos iτ s n.i sin iτ , n 0, 1, 2, . . . .

3.19
When seeking higher-order HAM approximations, 3.19 but not 3.18 is employed.Through this means, an HAM solution with N harmonics can be obtained.Interestingly, as long as this solution converges, it must converge to one harmonic balance solution, as revealed in 30 .

Results and Discussion
Theoretically, analytical solutions of u i can be obtained by solving 3.16 and 3.17 step by step, since all equations are linear.However, their expressions can rapidly become complicated as n increasing.It is rather easy to obtain semi-analytical solutions when all parameters are given.The HAM solutions are compared with numerical results.All the numerical solutions are obtained via integrating 2.5 using the RK method, with the initial values as x 0 0 and ẋ 0 0. Let T N,n denote the CPU time needed to obtain the nth-order HAM approximation u n with N harmonics, and T n is the counterpart when no harmonic is neglected.In addition, let A x be the maximum value amplitude of the nth-order HAM solution i.e., n i 0 u i with N harmonics retained, and B x the counterpart.From Figure 2, one can see the convergence of A x and B x to the numerical one as n increasing.It implies that the accuracy will not be lost by retaining only the first 10 harmonics N 10 .Also shown in Figure 2 is the ratio between T N,n and T n .It rises rapidly when n increases.Therefore, a lot of computational effort can be saved without losing substantial accuracy, through neglecting some high harmonics.Note that, all the following HAM solutions are obtained by employing 3.19 rather than 3.18 .
A preliminary procedure to obtain an HAM solution consists in the choosing of the auxiliary parameter h.Liao 23 suggested an h-curve approach for choosing a proper value of h, which can guarantee both the convergence of HAM series and a relatively rapid convergence rate.Figure 3 shows the relative errors versus varying h.These errors refer to the percentage-wise differences between the amplitudes A x of the 20th-order HAM solutions with 10 harmonics and the numerical results.According to the figure, an appropriate range for h can be chosen from −2 to −0.5 or so.When the absolute value of h is too large, for instance more than 2, it can probably make the HAM series divergent.On the other hand, as long as h approaches 0 enough, the HAM series converge by a relatively slow convergence rate.Figure 4 shows a phase plane of system 2.4 .The solutions provided by the HAM approach the numerical one when m increases.The 20th-order HAM solution is nearly the same as the numerical one.In fact, solutions can be obtained to any desired accuracy without additional difficulty.
When a nonlinear spring and/or squeeze film damping are taken into account, the proposed algorithm can still provide very accurate solutions, as shown in Figures 5 and 6   Both the figures show clearly the convergence of the HAM solutions to the numerical ones.It is worthy of emphasizing that, the vibrations described by these phase planes are nonharmonic, namely the solutions cannot be approximated by the first harmonic alone.That is probably because high harmonics contribute to the responses to an nonnegligible extent.
The phase plane shown in Figure 6 even has a loop.Even though, the presented HAM can still track them with excellent precision.These illustrative examples demonstrate the feasibility and efficiency of the presented approach in analyzing nonlinear dynamical behaviors.
Frequency-response curve is of substantial importance for dynamic analysis, when external excitations are included.Figures 7, 8 of periodic solutions versus varying ω, for different values of V P , V 0 and d, respectively.The 20th-order HAM solutions are in excellent agreement with the numerical ones.The amplitudes decrease with ω increasing for all the three cases.According to 2.5 , V P and V 0 correspond to the magnitude of the external excitation, respectively.Also, d is connected with the external excitation through T ε 0 A/ 2kd 3 , which implies that T increases when d decreases.That is why the amplitude increases when V P V 0 increases or d decreases, as shown in these figures.respectively, the larger k 3 or c 2 is, the smaller the amplitudes become.According to 2.5 , k 3 x 3 is the nonlinear restoring force, hence k 3 x 2 can be considered as an equivalent linear stiffness coefficient.Its value keeps nonnegative though varying.The increasing of k 3 implies an increasing of the stiffness, so that the magnitude of response decreases.As concerning c 2 , it denotes the coefficient of a positive nonlinear damping.Generally speaking, positive damping can suppress the vibration to some extent.It is worthy to point out, the influence of squeeze film damping on vibration amplitudes reduces rapidly when ω increases.When ω is relative large ω > 1.5 , the decreasing of c 2 makes a very little difference of the amplitudes.

Conclusions and Remarks
Based on the homotopy analysis method, we have proposed an approach to solve the nonlinear dynamical system of an electrostatically actuated microcantilever.This approach was further improved by restricting the number of harmonics as a given number N. It turns out that the improved approach can save a large amount of computational effort without reducing accuracy.Numerical examples validate the proposed approach.Approximations for periodic solutions are obtained very precisely.Using the presented approach, the frequencyresponse curves are obtained and discussed in details.The feasibility and efficiency, illustrated by the numerical examples, imply that we could expect the presented approach be applicable in more nonlinear vibration problems, especially those arising in MEMS.
As mentioned above, we have presented an auxiliary operator 3.15 that is different from another usually-adopted one 3.11 .It is advantageous to use 3.15 rather than 3.11 for the considered problem.First, there is no need to choose an initial guess for the zerothorder HAM approximation if one employs linear operator 3.15 .Second, if 3.11 is adopted, a series of algebraic equations have to be solved from time to time to eliminate secular terms when seeking every order HAM approximation.In many occasions, the first equation is nonlinear.In our procedures, there is no requirement of doing so.Last but not least, the proposed procedures are purely linear.

Figure 1 :
Figure 1: A simplified dynamical model of the microcantilever in MEMS.

Figure 2 :
Figure 2: a the HAM solutions; b the ratio between T N,n and T n , N 10, and h −1.Parameter values are given as V P 1, V 0 2, d 0 1, ω 0.5, k 3 1, and c 2 1.
Figure4shows a phase plane of system 2.4 .The solutions provided by the HAM approach the numerical one when m increases.The 20th-order HAM solution is nearly the same as the numerical one.In fact, solutions can be obtained to any desired accuracy without additional difficulty.When a nonlinear spring and/or squeeze film damping are taken into account, the proposed algorithm can still provide very accurate solutions, as shown in Figures5 and 6.

Figure 4 :
Figure 4: Comparison between numerical solutions and HAM results h −1, where heavy dots denote the former and solid lines the latter.Parameter values are given as V P 1, V 0 3, d 1.2, k 3 c 2 0, and ω 1.

Figure 5 :
Figure 5: Comparison between numerical solutions and HAM solutions h −1, 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, k 3 0.5, c 2 0, and ω 1.

Figure 6 :Figure 7 :
Figure 6: Comparison between numerical solutions and HAM solutions h −1, 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, k 3 0.5, c 2 2, and ω 0.5.