Periodic solutions of a graphene based model in micro-electro-mechanical pullin device

Phase plane analysis of the nonlinear spring-mass equation arising in modeling vibrations of a lumped mass attached to a graphene sheet with a fixed end is presented. The nonlinear lumped-mass model takes into account the nonlinear behavior of the graphene by including the third-order elastic stiffness constant and the nonlinear electrostatic force. Standard pull-in voltages are computed. Graphic phase diagrams are used to demonstrate the conclusions. The nonlinear wave forms and the associated resonance frequencies are computed and presented graphically to demonstrate the effects of the nonlinear stiffness constant comparing with the corresponding linear model. The existence of periodic solutions of the model is proved analytically for physically admissible periodic solutions, and conditions for bifurcation points on a parameter associated with the third-order elastic stiffness constant are determined. c © 2017 University of West Bohemia. All rights reserved.


Introduction
The importance of graphene can be highlighted in the title "The Arms Race for Graphene is Officially On" of an article published in the Wall Stress Daily on May 21, 2014 [20].The potential applications of graphene are everywhere [16].Nanoscale engineering devices that use graphene as a material for basic components in resonators, switches, and valves are being developed in many industries [1,4,5,12,13,23,24].Understanding the mechanical responses of graphene structure elements subject to applied loads are clearly important, see e.g., [15], and above mentioned references.The use of graphene as an actuating cantilever beam subject to electric pull-in force can be a major improvement over the use of the traditional pull-in materials, since graphene is more durable, stronger, and corrosion resistant, see e.g., [1,4,12,13].Many researchers have been using the linear constitutive equation in their models for graphene devices, see [1,4,19,23].For example, a headphone equipped with a graphene membrane has become a real life application of graphene [16,23], however the mathematical model for the membrane reported in the work [23] considered graphene as a linear elastic material and a linear lumped parameter model based on Hooke's law was used to estimate the frequencies of vibrations.Nonetheless, the nonlinear mechanical behavior of graphene is well-known even for small strains and there are theoretical [2,11,15] and experimental [14] validations of the nonlinear mechanical behavior by using nonlinear constitutive equations.The widely accepted continuum mechanics based nonlinear constitutive equation for graphene based materials in one dimensional form is given by σ where ε is the axial strain, σ is the axial stress, E is the Young's modulus, and D = −E 2 /(4σ max ) the third-order elastic stiffness constant, and σ max the ultimate yield stress of the graphene.Equation ( 1) is valid for |σ| ≤ σ max .It is demonstrated in [2] and [3] that the material constant D depends on the direction considered on the graphene plane.
The nonlinear mechanical behavior of graphene presents some difficulties in their structural simulations and designs.Graphene based materials have been considered for applications in actuators by material scientists, see e.g., [10].The diagram in Fig. 1 below is used to display the corresponding physical model of this study which consists of a spring as a model for the graphene based material strip attached with a plate.This kind of model can be found in [10].A substrate under the plate is kept parallel to the plate with a distance-called the gap and denoted by g in the following.It is assumed that a voltage is applied to the system creating an electrostatic force to pull the plate towards the substrate.This is commonly considered as the simplest model for stability study of electrostatic pull-in devices.Many researchers do not take into account the third-order stiffness constant and use only the linear model in their works associated with this physical model, see, e.g., [16,23,24].A recent review paper by Zhang et.al., [22], provides the history and a detailed account of the work on stability study of electrostatic pull-in phenomenon based on lumped-mass and other models.Since we are unable to solve the mathematical equations analytically, numerical methods have been used for most of the results presented here.
The goal of this work is to extend the results in [9] on this model for free and mechanically forced vibrations to the analysis of forced vibrations by considering the effect of the third-order stiffness constant as well as the electrostatic force.Through our results, we demonstrate that the second order material constant D is an important factor in modeling the patterns of graphene in vibrations and stability of electrostatic pull-in devices just as reported in [9] for mechanically operated devices [22].
In Section 2, we present the lumped-mass nonlinear spring mathematical model for electrostatically forced axial vibrations, and then in Sections 3 to 6, we present analytically and numerically the pull-in voltages for two standard cases of applications, the phase diagrams at low and high energy levels, numerical wave solutions, amplitudes, and frequencies for specific parameters of the spring-mass model.The effects of the second order material constant D on the solutions are demonstrated graphically in comparing with the corresponding linear model.In Section 7, we establish a mathematical condition on a parameter depending upon D for existence of periodic solutions of the forced vibrations within certain energy levels, which can be used to determine critical length of the graphene strip in the model leading to instability of the system.

The nonlinear lumped-mass model
We assume that the cross-sectional area of the graphene sheet is A c , one end of the sheet is fixed, and the other end is attached with a flat plate with mass m.The plate is subjected to a force F causing axial displacement of the graphene based material strip, which is mechanically modeled as a nonlinear spring.The axial displacement from equilibrium of the plate is denoted by x, the engineering strain is approximately modeled as x/L, where L is the length of the graphene sheet.The restoring force of the spring is modeled by the restoring force (1).By the Newton's second law of motion, the vertical displacement variable x in the lumped-mass nonlinear spring model shown in Fig. 1 satisfies the following nonlinear wave equation where the constants E and D are the material constants appeared in (1).We wish to know how the solutions of nonlinear spring-mass model ( 2) responds to the electrostatic force of the form 2 and how important it is to include the third order stiffness constant D in the graphene's constitutive equation (1).In this case we have the following spring-mass equation: where ε 0 is the electric emissivity, V the applied voltage, and g the gap between the plate and the substrate, A c the cross-sectional area of graphene based material strip, A the area of the plate in the device illustrated by Fig. 1.For the derivation and physics and engineering aspect of models similar to (2) with the electrostatic force F elec , see e.g., [21] and [17].We remark that this model can be used for nonlinear materials similar to graphene with second order stress-strain responses.

The pull-in voltages
The pull-in of the system is defined as the result of a sudden displacement of plate when the electrostatic force overcomes the restoring force and the stiffness of the spring, causing the plate physically touch the substrate under it.The pull-in voltage is defined as the required voltage which causes the pull-in the spring-mass system.To find pull-in voltage for the model, we consider the corresponding static equilibrium equation of (3): which can be solved for the pull-in voltage: As it is pointed out in [7] that the pull-in occurs when the structure deflects to 2/3 of the initial gap g for a number of pull-in models.This statement has been experimentally and numerically validated and reported in [6,8,18].In this case, x = 2g/3 and the pull-in voltage of ( 3) is given by where α = − Dg EL = Eg 4σmaxL .It is also assumed that pull-in can occur for many traditional pull-in devices when the structure deflects to 1/3 of the initial gap.In this case, we also have pull-in for (3), see also, e.g., [5], for the assumptions leading to these practices.The pull-in voltages corresponding to various plate areas and the two standard gap sizes are listed in Table 1.Adopted from [13], we use E = 1 000 GPa, A c = 1 μm 2 , g = 2, L = 600 μm, and ε 0 = 8.85 × 10 −12 F/m for the permittivity of vacuum for numerical calculations of the pull-in voltages.The results provided in Table 1 appear to be consistent with the values found in [5] and [13], where D is not included in the models.

Phase diagrams
For plotting the phase diagram and the bifurcation analysis, we will work with the following dimensionless form of (3): For simplicity, we will still use x and t to denote the x and t in (4).Multiply both sides of (6) by ẋ and integrate the result relative to t, we then have the conservation of energy equation where y = ẋ and U 0 is the integration constant.First, we provide numerical solutions for some specific values of K, g, L and U 0 for plotting the phase diagram of (6) by using (7).By using the change of variables, we can convert the phase diagrams of (6) to obtain the phase diagrams of (3).For simplicity, we only show the phase diagrams of (3) for some specific sets of parameters.The numerical phase diagrams are provided at different energy levels defined by the values of U 0 to show the existence of periodic solutions and provide some further inside information about the existence of periodic solutions.We plot diagrams for the voltage of 30 V which is lower than the pull-in voltage of 37.27 V for the plate with area that equals to 160 000 μm 2 .Several phase diagrams are presented below to demonstrate the effect of the nonlinear term and the electrostatic force on the periodic solutions.Fig. 3a shows the phase diagrams for higher levels of the initial energy U .The impact of the nonlinear term α|x|x versus the effect of the electrostatic force term K/(g − x) 2 on the phase diagram for (3) is also illustrated by Fig. 3, in which (a) corresponds to the case when α|x|x dominates K/(g − x) 2 and (b) to the opposite case.Our calculations indicate that the nonlinear restoring force term has stronger effects than the electrostatic force on the periodic solutions.This seems to indicate that a full nonlinear analysis of the pull-in device is necessary for nonlinear materials with second order stress-strain responses.

Numerical wave forms
In this section, we present graphic presentations of the periodic wave solutions.We use α = 0.032 1, K = 0.038 23, x(0) = 0.3, x (0) = 1 and α = 0.0, K = 0.038 23, x(0) = 0.3, x (0) = 1, respectively, to obtain Fig. 4a for comparison of solutions of (3) with and without nonlinear stiffness effect.The nonlinear solution is slightly different than the linear solution for the same linear stiffness constant.Our numerical results also show how the nonlinear numerical solution changes when the magnitude of the nonlinear term is increased significantly (by increasing the values of α, see Fig. 4b.Hence, we can suspect that the period increases as the nonlinear term becomes bigger and frequency decreases as period goes up.It is clear that there is a bifurcation point, but it is not, for now, possible to find the general analytic expression for it.For particular values of parameters of the spring-mass system, numerical value of the bifurcation point may be computed.However, it is not practical to do this computation every time, so we will establish a theoretical result later in Section 7 to determine such a bifurcation point.Since we are unable to solve the mathematical equation (3) analytically, numerical methods are used.The numerical technique adopted to solve the nonlinear differential equation ( 3) is based on the Matlab function ode45, which uses the four stage Lobatto IIIA formula due to the good stability properties of the procedure.The second-order nonlinear differential equation ( 3) is transformed by using y = (x, ẋ) into a first order system of the form dỹ dt = f (t, y), which can be solved by following the s-stage Runge-Kutta method that is, solving a system of nonlinear equations where . ., N, are the mesh points on the time interval [0, T * ], h is the step size.The number T * > 0 is the chosen maximum time for the numerical simulation.Different algorithms, such as Lobatto IIIA, are obtained by choosing a specific sets of values for parameters α i , β i and γ i .See, e.g., the reference [18] for more details.

Periods and frequencies
Solving (7) for ẋ and integrate the separable first-order ODE, converting to the dimensional variables in (7), we then obtain the integral formula to determine period and frequency where x * denotes the smallest x-intercept of the energy equation ( 7) in the phase plane in (x, y) and T denotes the period of the corresponding periodic solution.We assume that U 0 = 1, V = 30 V, A = 160 000 μm 2 , and use the following sets of parameters: For the given three sets of values, we have the following values of the corresponding intercepts: By using these values, we then integrate numerically the integrals and have the following three numerical values of the periods for these three cases: 1) T 1 ∼ = 6.189 48, 2)T 2 ∼ = 6.367 34 and 3) T 3 ∼ = 6.854 31.It can be observed that as we increase the value of α from 0.003 21 to 0.128 4 (we increase the effect of the nonlinear term), the corresponding period increases from 6.189 48 s to 6.854 31 s.The frequencies can be found by the formula f = 1/T and we get f 1 ∼ = 0.161 56 s −1 , f 2 ∼ = 0.157 05 s −1 and f 3 ∼ = 0.145 89 s −1 .

Existence of periodic solutions and bifurcation analysis
In this section we discuss the existence of periodic solutions of (6) with initial conditions by studying the graphs of the energy equation (7) for various values of the initial condition and the parameter α.We assume that the initial conditions satisfy |y(0)| > √ 2K and x(0) = 0.Then, we have U 0 = 1 2 y(0) 2 − K > 0. The y-intercepts of the graph can be easily found to be ± 2(U 0 + K).We consider the function for finding a bifurcation point.Here, we demonstrate the existence of physically admissible periodic solutions of the forced vibrations within the gap g = 1 in the non-dimensional variable x and for certain initial energy levels U 0 through the following Lemma and Proposition: Lemma 1 There exists a physically admissible periodic solution if and only if f(x) = 0 for some x ∈ (0, 1).
Proof: We note that (7) defines a periodic solution of (6) if it defines a closed curve in the phase diagram.To be a physically admissible periodic solution it needs to be the closed curve around the equilibrium solution (0, 0) in the phase diagram.Since y = ± 2 f(x) and f (0) = K + U 0 > 0 we see that for the equation defines a closed curve around (0, 0) there must exist x 1 < 0 < x 2 < 1 such that f(x 1 ), f (x 2 ) ≤ 0. We note that for any x ∈ (0, 1) it holds that f (x) > f(−x) as K/(1 − x) > K/(1 − (−x)).Thus, the existence of x ∈ (0, 1) with f (x) = 0 immediately implies the existence of x < 0 satisfying f(x ) = 0, since f(x) is continuous and f(0) = K+U 0 .Hence, it is necessary and sufficient to look at if f (x) = 0 for some x ∈ (0, 1).
Proof: In view of the Lemma 1, for given K, U 0 > 0, we investigate the existence of a solution of f (x) = 0 in the interval (0, 1) depending on the parameter α ≥ 0. As f (0) > 0 and f is continuous on [0, 1), we see that it suffices to find x ∈ (0, 1) such that f (x) ≤ 0. Solving for α from f (x) ≤ 0, we get α ≤ H(x).Therefore, the existence of x ∈ (0, 1) such that f (x) ≤ 0 is equivalent to the existence of x ∈ (0, 1) such that α ≤ H(x).Clearly, such a non-negative α exists if and only if h(x) ≥ 0 for some x ∈ (0, 1), where h this is possible if and only if there exists a critical point of h(x) such that h(x 0 ) ≥ 0. Taking the derivative, we have h (x) = 2x−3x 2 +2U 0 , we see that h(x) has two critical points (1± √ 1 + 6U 0 )/3 of which only x 0 = (1+ √ 1 + 6U 0 )/3 is positive and it satisfies x 0 < 1 if U 0 < 1/2.Thus, we conclude that (7) has a closed curve, i.e. the equation ( 4) has a physically admissible periodic solution if and only if U 0 < 1/2 as well as h(x 0 ) ≥ 0. The latter is equivalent to H(x 0 ) ≥ 0 which finishes the proof.
We now assume that U 0 < 1/2 and H(x 0 ) > 0. From Proposition 1, we see that α ≥ H(x 0 ).Moreover, we see from Lemma 1 that if f (x) > 0 on (0, 1) then there are no admissible periodic solutions.Clearly, on (0, 1) and the right hand side function has the global minimum at x = 1/α on (0, ∞).So, if − 1 2α 2 + α 3α 3 + K + U 0 ≥ 0, we have no admissible periodic solutions.Although our numerical phase diagrams indicate the existence of periodic solutions with |x(t)| > g, t ≥ 0, they are considered extraneous and physically impossible since they would be below the line of the substrate in the model.Therefore, we omit the mathematical analysis of these cases.

Conclusion
A nonlinear spring-mass model for electrostatically forced axial vibration of a graphene in a simple pull-in device is presented following the work of [9].The existence of periodic solutions of free and forced vibrations is shown through phase plane analysis.Dependence of the periodic solutions to the model and bifurcation points on the parameter α = −Dg/(EL) is identified numerically.Pull-in voltages are calculated for different areas of the plate at two different gap sizes.The analysis of phase diagrams demonstrates that, for high levels of the initial energy, some unknown behaviors of the solutions are discovered which are distinct from those found in [9] which only considers mechanically forced vibrations.Numerical solutions of some periodic waves have been presented and the values of periods and frequencies for three different values of length of the graphene are shown.Comparisons of the solutions with and without the presence of the second-order stiffness constant D = −E 2 /(4σ max ) in the constitutive stress-strain equation for graphene are illustrated graphically and dramatic differences in these solutions are found.However, our results indicate that graphene behaves mechanically like all the nonlinear continuum materials with 3rd order elastic stiffness.Among these materials, no graphene-specific phenomena are found.

Fig. 4 .
Fig. 4. Comparison of waves of the linear model and the nonlinear model