Accurate Simulation of Parametrically Excited Micromirrors via Direct Computation of the Electrostatic Stiffness.

Electrostatically actuated torsional micromirrors are key elements in Micro-Opto-Electro- Mechanical-Systems. When forced by means of in-plane comb-fingers, the dynamics of the main torsional response is known to be strongly non-linear and governed by parametric resonance. Here, in order to also trace unstable branches of the mirror response, we implement a simplified continuation method with arc-length control and propose an innovative technique based on Finite Elements and the concepts of material derivative in order to compute the electrostatic stiffness; i.e., the derivative of the torque with respect to the torsional angle, as required by the continuation approach.

1. Introduction

In the family of Micro-Opto-Electro-Mechanical-Systems (MOEMS), electrostatically-actuated torsional micromirrors are currently widely diffused in laser micro-scanners (e.g., [1]), optical shutters, micro-spectrometers, micro lenses, and pico-projectors. Indeed, electrostatic actuators can be easily fabricated using micromachining techniques, and driving voltages are compatible with standard integrated circuit (IC) technology.

Electrostatic micromirrors have been analysed in several recent papers [2,3,4,5,6,7], where it has been shown that they are governed by parametric resonance—a known phenomenon in the field of non-linear dynamics affecting, for example, the stability of ships, the forced motion of bridges, and Faraday surface wave patterns on water [8,9,10,11,12,13]. The importance of parametric amplification at the microscale has also been widely recognized for different applications (e.g., [14,15,16,17,18,19,20,21,22]).

The 2D layout and a scanning electron microscopy (SEM) image of the silicon-on-insulator (SOI) micromirror is depicted in Figure 1. Similar components are utilized in the Intel Real Sense technology [23]. The central circular reflecting surface is attached to the substrate via two coaxial beams acting as torsional springs. Four sets of 33 fingers each are anchored to the trapezoidal regions directly attached to the mirror. These plates, interdigitated with their stator counterparts, form a comb drive structure, providing the electrostatic actuation mechanism. During operation, sensing of the opening angle is performed via the same comb drive electrodes.

Figure 1

Resonant micro-mirror: SEM image of the device and layout. Courtesy of STMicroelctronics.

Following the discussion published in [7], the mirror is treated as a rigid body hinged in its centre and connected to the substrate via the two elastically deformable springs. The torsional response is hence governed by the simple 1D model: where ng m is the inertia around the torsional axis, N m s is a damping coefficient, N m/rad is the torsional stiffness of the springs, and is the electrostatic torque due to a unit voltage bias. The voltage V is given in volts, and time in Equation (1) is measured in microseconds. Since pF/m, C is measured in microns and is a purely geometrical feature. The torsional eigenfrequency of the mirror is Hz. The dissipation term, in particular, is a delicate issue. Due to the presence of comb-finger actuation and of a large gap between the mirror plate and the substrate, the main contributions to dissipation are from shear flow in the comb fingers and the transport of mass induced by the large rotation of the mirror. Preliminary results have been obtained on a similar structure in [24]; however, to increase the accuracy, in the present work the constant B has been utilized and calibrated starting from the known maximum opening angle of the mirror in operative conditions.

It is worth stressing that, due to symmetry, the electrostatic torque around the torsional axis vanishes for ; i.e., in the rest configuration, . However, this trivial solution becomes unstable for some combinations of the input voltage and frequency, and triggers the mirror rotation. Indeed, in [7], it is shown that Equation (1) reduces—in the case of small rotation, negligible dissipation, and sinusoidal excitation at frequency f—to the well-known Mathieu’s equation [8,16]: where is a non-dimensional time, , and is proportianal to the square of the maximum voltage. Equation (2) admits a nontrivial response only in specific regions of the plane called instability tongues, which emanate from axis. In particular, the first and largest instability tongue occurs when the system is excited by a periodic forcing function at frequency , which is totally different from classical phenomena of harmonic resonance.

While the functioning of similar devices has already been extensively investigated in the literature, we focus here on a specific topic pertaining to the numerical simulation of stable and unstable branches of their response. A complete characterization of the mirror requires refined numerical techniques like the continuation approach that is briefly described in Section 2. This tool rests on an accurate evaluation of the electrostatic forcing terms and of their derivative with respect to the rotation angle. The latter term is typically computed naively with finite difference, often inducing convergence issues in the procedure. On the contrary, in Section 3 we propose a new direct procedure based on classical finite elements. Numerical results are finally discussed in Section 4.

2. Numerical Solution via a “Continuation” Approach

The numerical simulation of Equation (1) can be performed with several techniques. The brute force approach—which closely reproduces the actual operation of the mirror—performs a sweep over the frequencies of interest, and for each frequency simulates a sufficient number of cycles by direct integration in time to reach a steady state; the amplitude is then recorded and the next frequency is addressed, using the final amplitude and phase of the previous analysis as initial conditions. This is a very robust technique; however, it only permits the simulation of the stable branches of the amplitude vs. frequency response.

On the contrary, the continuation approach [12,25,26] with arc length control is more versatile, and is adopted herein. The model in Equation (1) is rewritten as a first-order non-autonomous differential system of equations in terms of the fictitious time , with :

In Equation (3), the prime denotes differentiation w.r.t. , and the voltage V is a generic function of , with . We limit our attention to periodic forcing functions with and to periodic solutions, so that our problem can be rewritten in condensed form as: where is a constant matrix, collects the two unknown functions, and contains the forcing function. It is worth stressing that, according to our assumptions: (i) In the interval T of analysis, and might contain and cycles, respectively, with in general. This is essential for the simulation of different instability tongues; (ii) The system could be generalized and transformed into an autonomous one by adding a nonlinear oscillator generating the desired periodic forcing, as done for instance in the software AUTO [25]. This option is not implemented herein.

Let us suppose now that is a known solution of the system. The simplest choice consists of taking T as continuation parameter: we fix and solve Equation (4) for through an iterative Newton–Raphson procedure using a suitable initial guess. However, this classical parameter continuation fails in the presence of an unstable branch as is clear from Figure 2. Indeed, by imposing an increment (and hence ) at the peak, the solution would jump to the stable solution, completely missing the unstable dashed branches. For this reason, it is customary to introduce an arc-length control in which is part of the unknowns, the abscissa s along the solution branch is taken as the continuation parameter, and a new constraint is added. A typical choice is:

Figure 2

Sinusoidal excitation in the range of 5 kHz, = 55 V; experimental upward and downward sweep (discrete symbols) and numerical continuation (continuous and dashed line).

An alternative is the Keller’s pseudo arc-length method [25] in which the increment is sought such that its projection along a specific direction (typically the tangent to the manifold) equals a fixed arc-length .

In our simple implementation, we enforce Equation (4) in a weak manner, fixing a suitable space of vector test functions :

The segment [0-1] is partitioned into N equal elements, and the unknown function is discretized over each element with quadratic Lagrangian shape functions. On the contrary, the space of test functions is selected as the space of piecewise Legendre orthogonal polynomials . If Equation (5) is integrated numerically over each element with a two-point Gauss–Legendre quadrature rule, this is equivalent to the method of orthogonal collocation, in which Equation (4) is collocated at the two zeros of in each element.

Equation (5) is solved iteratively by means of a Newton–Raphson procedure. In the generic iterate, given an estimate for (and hence for ), a small correction is sought such that is the solution of the linearized system:

The implementation of Equation (7) requires the computation of the derivative of the electrostatic torque with respect to the rotation angle . A commonly adopted strategy is to numerically compute using some dedicated software and then differentiate it using standard finite differences. However, this approach often lacks the required accuracy for the iterative procedure to converge. In the following section, we detail an alternative direct formulation, based on the notion of material derivative.

3. Derivatives of Electrostatic Force and Torque via Direct Finite Element Computation

In this work, a finite element method (FEM) has been selected to compute the electrostatic force and torque directly by post-processing the potential field of the classical electrostatic problem.

The rigid fingers of the mirror are collected in two groups: the “stator” and the ”shuttle” —the former being fixed and the latter being movable. The stator and the shuttle are immersed in the infinite space truncated at , where homogeneous Neumann conditions are enforced. Without any loss of generality, we will assume that on and on , and we will express the solution as . The unknown function is governed by the following set of equations:

If denotes the space of sufficiently continuous functions respecting Dirichlet boundary conditions, the weak form of Equation (8) is:

Find such that:

3.1. Direct Computation of Force and Torque

A straightforward way of computing forces and couples by a simple post-processing of flows from the notion of material derivative applied to the total potential energy [27,28,29].

Let us assume that a change of configuration is generated by the variation of a given parameter p (e.g., the rotation angle of the shuttle), and is associated to the fictitious velocity . Following classical texts of optimization [27], the material derivative (also called total or particle derivative) of a scalar f, of a gradient and volume element due to are, respectively:

The restriction to of any admissible must be a rigid body motion , while must vanish on . The following expression for the virtual velocity will be adopted:

Now: where , . Let denote the material derivative of the total potential energy W. If and are the force and the torque exerted on : where is the electrostatic energy for a unit voltage bias. Applying Equation (10):

However, the first term vanishes, since: due to Equation (9) written with . Indeed, is an admissible test function that vanishes on both and . At equilibrium, , and hence: yielding the final general equation: since . In particular, the force is obtained setting and the torque setting where . It is worth stressing that Equations (18) and (19) for the force and the torque do not involve the material derivative of the potential. However, this is required in order to differentiate and .

3.2. Material Derivative of the Potential

Since the material derivative of the weak form Equation (9) vanishes identically: we have:

It is worth stressing that the discretization of Equation (21) leads to the same matrix as the primal problem, which can be efficiently factorized only once and solved for different rhs vectors.

Let and denote two vectors such that:

The components of and can be computed from Equation (21), fixing six independent choices for and .

3.3. Differentiation of Force and Torque

The procedure can be iterated by applying material derivatives to Equations (18) and (19). At this level, is required. It turns out that , with:

Similarly, , with:

4. Numerical Results

While Equations (23) and (26) are general and could be applied to any rigid-body movement of the structure, for the case of interest herein, only is required. The virtual velocity is a rotation around the torsional axis of unit vector , hence , , and the material derivative of the potential is . Moreover, . Eventually, setting :

One of the four sets of comb fingers has been simulated with a custom FEM code implementing quadratic tetrahedra, coupled with the mesh generator GMSH [30]. The surface mesh of the fingers and a clipping of the volume mesh are collected in Figure 3. The configuration considered corresponds to a rotation of rad.

Figure 3

Surface mesh of one set of comb fingers (a) and clipping of the volume mesh (b).

The potential distribution is plotted in Figure 4. The large-scale analyses require specific hardware, but eventually compare well with the approach based on integral equations utilized in [7].

Figure 4

Potential distribution in the comb. The shuttle is set to (red) and stator is set to (blue).

The computed electrostatic torque and its derivative are shown in Figure 5 for one set of comb fingers as a function of . Assuming a sinusoidal excitation at frequency f: , it is worth noting that the driving term will contain the harmonic. Thus, if f is close to the natural frequency of the mirror ( kHz), this will be excited with kHz activating the first instability tongue. Figure 2 presents the numerical results compared with experimental data (taken from [7]), assuming a sinusoidal excitation around 5 kHz with V.

Figure 5

Electrostatic torque (a) and derivative of electrostatic torque (b) for one set of comb-fingers with respect the torsional angle .

The trivial solution always exists for any f, but becomes unstable (red dashed line) between approximately 5380 Hz and 5420 Hz. The numerical simulation predicts the existence of two non-trivial branches: one is stable (to the right of the peak) and one is unstable (to the left of the peak). This is reflected in experimental data obtained controlling f, both in an upward and in a downward sweep. In the former case, f is increased starting from a low frequency, where is the only solution and is stable. However, when Hz the mirror becomes unstable and the system jumps to the stable non-trivial branch which is next followed up to the intersection with the trivial solution at Hz. Similar remarks hold when a downward sweep is performed, with the difference that the stable branch is followed until the peak is reached and next suddenly drops to zero, “switching off’ the torsional movement.

The mirror response strongly depends on the voltage V. According to the theoretical analysis of the Mathieu equation in the presence of dissipation, a threshold exists such that below this value only the trivial response exists. This is confirmed by the numerical simulations: for V, the maximum aperture at the peak is 0.7 deg, and no solution is found when is further decreased.

5. Discussion and Conclusions

We have developed a new approach based on finite elements allowing the computation of the electrostatic force and torque exerted on solid bodies, but also the electrostatic stiffness; i.e., their derivative with respect to parameters governing the movement of the shuttle. Although the applicability of these results is quite general, our investigation has been motivated by the need to simulate phenomena of parametric resonance in electrostatically actuated micro-mirrors. The highly non-linear response of similar devices typically presents both stable and unstable branches which can be simulated using “continuation approaches”. In these techniques, starting from a known solution corresponding to a given value of a parameter, the non-linear governing equations are linearized following a Newton–Raphson procedure in order to compute the system response to a variation of the parameter. The electrostatic stiffness appears in the first-order expansion of the equation, and must be carefully evaluated in order not to spoil the quadratic convergence of the numerical scheme.

This degradation of convergence (or even lack thereof) typically arises when the torque is computed directly (or even only the capacitance), while the stiffness is obtained indirectly by applying finite difference formulas to the torque curve. This procedure does not guarantee the required accuracy. The aim of the proposed approach is to also directly compute the stiffness terms, with a precision comparable to the torque computation.

While forces and torques are based on a simple post-processing of the scalar potential of the primal problem, the computation of the electrostatic stiffness requires solving for the (material) derivative of the potential. However, this latter problem leads to the same system matrix as the primal one, thus inducing a minimal increase of computational cost.

Alternative approaches are indeed available, like those based on integral equations. However, they require dedicated codes of great complexity, while the present proposal could be easily configured as an add-on to any commercial finite element code. Clearly, a possible limitation is represented by the need of a 3D volume mesh, which leads to large-scale problems. The availability of low-cost multicore processors with large amounts of dedicated memory is, however, rapidly pushing the limit, making the solution of realistic structures on standard hardware possible.