Optimal design of a bistable switch

Michael P. Brenner, Jeffrey H. Lang, Jian Li, Jin Qiu, and Alexander H. Slocum
Division of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138; and Departments of Electrical Engineering and Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139

Edited by George C. Papanicolaou, Stanford University, Stanford, CA, and approved June 6, 2003 (received for review March 14, 2003)

Determining optimally designed structures is important for diverse elds of science and engineering. Here we describe a procedure for calculating the optimal design of a switch and apply the method to a bistable microelectromechanical system relay switch. The approach focuses on characterizing the unstable transition state connecting the two stable equilibria to control the force displacements. Small modications in component shape lead to a substantial improvement in device operation. Fabrication of the optimized devices conrms the predictions.

eveloping the ability to predict and quantify optimal designs for devices of increasing complexity is an important direction for mathematical and computational research. As technologies are scaled down, device optimization within the allowable design space can make or break its functionality. Many of the fabrication technologies associated with small-scale device manufacture (e.g., deep reactive ion etching) allow significant freedom in designing the shapes of components, at minimal additional fabrication cost; thus, there is potentially significant room for improving device performance by correctly choosing component shapes. This article presents a methodology for optimizing the component shapes of a switch; the method is applied to optimizing a bistable microelectromechanical system (MEMS) switch made from elastic beams. The difficulty of the optimization is that the device involves the continuum of states of the switch, connecting the two stable equilibria. Traditional methods in structural optimization find shapes to optimize the response of a single state of the system (1, 2), including designing the strongest column (3, 4), the tallest column (5), or compliant structures (610). The methods of structural optimization have also been applied to designing minimum drag structures, including the fastest falling shape at very low velocities (11, 12) and the shape of airplane wings with a desired pressure distribution (13, 14). The central idea of the method for optimizing a switch is that the force displacement characteristics of the switch can be manipulated by focusing on the unstable transition state connecting the two stable equilibrium states. Although the force displacement characteristics of the switch are highly nonlinear, the properties of this transition state depend only on the linear modes of the structures comprising the switch. Therefore the optimization reduces to a linear problem, approachable using traditional methods: The central calculation of the optimization reduces to a variation of Kellers strongest column problem (3). The organization of this article is as follows: in the next section, we describe the bistable MEMS relay switch to be optimized, together with a mathematical model for the device. Then, we discuss numerical calculations of the force displacement characteristics for the switch, together with a mathematical analysis that captures its essential features. Then, we describe the optimization and present optimal designs for several different applications of the switch. Finally, we discuss fabrication and testing. The Bistable Switch The bistable switch that we will optimize (15, 16) operates with two initially curved beams that are clamped together in the center and at the ends. In this initial position there is no stress in the beams. When the structure is pushed in the center it

deflects and then snaps to a second stable position, as shown in Fig. 1. One application of the device (fabricated on a silicon wafer) attaches a metallic crossbar to the bottom of the switch, so that when the beam is pushed from one state to another, the switch can be closed. The advantage of this switch over other alternatives (1720) is that it uses no latches, hinges, or residual stress to achieve its bistability. In its original implementation (15), the beam had constant thickness h, with a backbone shape w 0(x) (d 2)(1 cos(2 x L)) (Table 1). Here d is the transverse distance of the center of the curved backbone relative to the ends, and L is the length of the beam. The thickness of the beams along their length is h, so that the moment of inertia of the beam is I bh 3 12, where b is the breadth of the beam (given by the thickness of the silicon wafer). These design choices fix the force-displacement characteristics of the switch, which follow from the equation for the deflection w of the beam under an applied force f. Working in dimensionless units, with the beam deflection w scaled in units of d, the horizontal scale in units of L, and the applied force in units of Eh 3b 12d L 3, where E is the Youngs modulus, we have d2 dx 2
2

d2 w w0 dx 2 6Q 2
1

d 2w dx 2 dw 0 dx
2

[1]

dy
0

dw dx

[2]

where Q d h. The function (x) is the (dimensionless) moment of inertia, and the dimensionless tension in the beam is 2 TL 2 (Eh 3b 12), with T the dimensional tension. The second equation reflects that the tension is determined by the relative extension of the beam. The boundary conditions for a clamped beam are w(0) w (0) w(1) w (1) 0. Fig. 2A shows the force as a function of displacement w(1 2) w 0(1 2) of the center of the switch, for a force applied at the center of the beam: f f (x 1 2), assuming Q 2 10. Fig. 2B shows the corresponding tension in the beam during its deformation. The noteworthy features of these curves are: (i) bistability, there are two stable equilibria separated by an unstable transition state; (ii) the maximum force f push( 757, in units of EId L 3) required to activate the switch is twice that necessary 370); (iii) the for returning to its original configuration f pop( stroke required for actuation is determined by the distance over which a positive force must be applied to flip the switch, given by the location of the unstable transition state * 1.33; 2 and (iv) the tension in the beam at the transition state is * 2 8 2. 16 2 and in final state f Note that the magnitude of the forces f pop and f push are almost 3 orders of magnitude greater than that expected from the dimensional estimate EId L 3. In optimizing the switch, we want to design the beam to both decrease the absolute magnitude of
This paper was submitted directly (Track II) to the PNAS ofce.
To

whom correspondence should be addressed. E-mail: brenner@[Link].

[Link] cgi doi 10.1073 pnas.1531507100

PNAS

August 19, 2003

vol. 100

no. 17

96639667

APPLIED MATHEMATICS

Fig. 1. Schematic for the operation of the relay switch. Cantilevered starting zippers (Top) initiate the zipping by closing the gap and then allowing the beam to zip closed when the middle electrode is charged (Middle), causing the switch to close. Turning on the top electrode (Bottom) causes the top starting zippers to engage, and the switch opens. Fig. 2. (A) Force displacement characteristics for the switch, assuming the force is applied in the center of the switch. The ratio of the maximum force to the minimum force is 2. (B) Tension 2 (2 )2 in the beam as a function of the displacement in the center. Through the middle part of the cycle the tension is constant with 2 16 2.

f pop and f push and manipulate their relative magnitudes. In particular, since in the envisioned application, the force f pop is the relay contact force, while f push is the actuating force, we want as much contact force as possible for a given actuation force. Analysis of the Switch Optimizing the force displacement characteristics of this switch requires first understanding which features of the beam shape determine the force displacement characteristics. What determines the magnitude of the forces f push and f pop? What determines the ratio of the forces? What determines the location of the transition state and the required stroke? Can the location of the stroke and the magnitude of the applied forces be tailored for a specific application? The central idea of the optimal design analysis is that the magnitudes of the forces and the stroke are almost completely determined by the linear properties of the three equilibrium states; this follows from bifurcation analysis, which provides a correspondence between the eigenfunctions of a linearized problem and the critical points of an associated nonlinear problem (21). The linear properties can be easily manipulated by changing (x), which allows optimization of the force displacement response. Changing the overall shape of the backbone of the beam w 0(x) has a much smaller effect on the essential properties of the switch.
Equilibrium States. To begin, we characterize the three equilibrium states of the switch. The first unstressed stable equilibrium is w w 0. The second stable equilibrium follows from writing 2 2 w Aw 0; Eq. 2 then implies A 1 O( 4Q 4 0 (6Q ) 0 3). The third equilibrium is the unstable transition state, which controls the path for passing from the unstressed stable equilibrium to the stressed stable equilibrium. The beam tension in Table 1. List of symbols
w(x) w0(x) (x)
2 n

the transition state must correspond to a free eigenvalue. The lowest energy deformation mode that is consistent with the boundary conditions has eigenvalue 2 and eigenfunction w 2. w 0 Bw 2. The solution in the transition state is therefore w t Eq. 1 then implies that 1
2 0 2 1

[3]

and B is determined by inserting w t in Eq. 2 and setting 2. Eq. 3 is an exact result for any beam cross-sectional shape (x). 1 3. A uniform beam has 0 2 , 2 4 , implying that Thus the displacement of the center of the beam at the transition state is * 1 4 3, exactly as observed in numerical simulations (Fig. 2).
Forces. We now turn to characterize the force displacement

relations, close to each of the equilibrium points; the information we obtain will be used subsequently for optimization. Fig. 2 demonstrates that a beam with a uniform cross section has a linear force-displacement response around the transition state. This linearity is a direct result of the constant tension in the beam (Fig. 2 Lower) over this wide range. A constant tension 2 implies that the force displacement characteristics follow from w w0
2 2w

f 2

1 2 .

[4]

wn *

Displacement of beam Backbone shape of beam Dimensionless moment of inertia Dimensionless tension nth free eigenvalue nth free eigenmode Displacement of beam center Location of transition state Force ratio

The solution satisfying the boundary conditions is w w0 B w 2 (f 2 3)w f, where is given by Eq. 3, B is determined from 2 sin( 2 x). Therefore, near the transition Eq. 2 and w f (x) 2x state for the uniform beam the force displacement relationship is given by f 4
2 2

[5]

Assuming that this linear regime extrapolates to the stable 16 3(4 ) 2 EId L 3 equilibria, Eq. 5 predicts f push 8 3(4 ) 2EId L 3 421EId L 3; 842EId L 3 and f pop these predictions agree with the simulation of Fig. 2, with an error of 5%. The results show that the prefactors in f pop and
Brenner et al.

9664

[Link] cgi doi 10.1073 pnas.1531507100

f push are so high because of 2 (4 )2 157. Decreasing the 2 absolute magnitude of the forces thus requires decreasing 2. Moreover, we see that the force ratio R f push f pop is controlled by the linear modes of the beam, despite the fact that the force displacement characteristics of the beam are highly nonlinear. For the uniform beam the force ratio is R f push fpop * * 2
2 2 0 0 2 2

shape, we need to derive a relation calculating the gradient of the eigenvalue ratio 2 0 with respect to changes in the beam shape (x). late the gradient of a function C({ i}) of the eigenvalues { i} with respect to changes in (x). Since C i , we i 1 (dC d i ) must compute the change in the eigenvalues caused by changing . If 3 , then the ith eigenfunction (eigenvalue) changes from w i 3 w i w i( i 3 i i ). Hence wi
2 i

Computing Gradients. First, we describe how to efficiently calcu-

2.

[6]

Nonuniform Cross Section. For nonuniform cross sections, the tension is not constant near the transition state: If we assume that w 2 at the transition state, then there is a 2 and w solvability condition for a solution to Eq. 4 to exist near the transition state: There is a solution when f 0 only if w 2(1 2) 0. Although this condition is satisfied for a uniform beam, it is not satisfied for general (x). Hence, for general (x), the tension near the transition state depends on the applied force f. Taking , the solvability 2 condition implies that fw 2(1 2) 1w 2. There is thus a 0 2 nonlinear relation between the force and the displacement. However, the change in the tension fw 2(1 2) is small when the force is small; numerical experience shows that w 2(1 2) is always 10 3 for the optimal nonuniform shapes. The tension is thus nearly constant on the entire branch through the transition state. Therefore, we can still approximate the force displacement relation through the transition state as f c( * ), where * 1 * is the location of the transition state and c is a constant depending on (x). Hence, Eq. 3 for * works for general (x), so the force ratio is given by Eq. 6. Furthermore, 2 we expect that c 2 , so that the magnitude of the slope of the force-displacement curve through the transition point is controlled by 2. We note that the magnitude of the forces through the stable equilibria are controlled by 0. Thus, these forces can be decreased by finding beam shapes which decrease 0.

wi

i wi

wi .

[7]

A solution only exists if the right side of this equation is orthogonal to w w i. Hence we need T iw i This implies that wi 2
i

wi

w idx

0.

. w i2dx

[8]

This formula explicitly gives the gradient of C with respect to ; the numerical cost of computing this gradient is identical to that for computing the eigenvalues and eigenfunctions. Note that this calculation assumes that the eigenvalues are simple; this assumption is validated a posteriori by the calculations that follow. Taking C({ i})
Optimal Designs. Now we use Eq. 8 for optimizing the force ratio.
2
2 0

Optimization We are now in a position to optimize the force displacement characteristics of the switch. The results of the previous section demonstrate that the absolute and relative magnitudes of f push and f pop are determined by the free eigenmodes of the structure: the force ratio is controlled by the ratio of 2 0, through Eq. 6; the magnitude of the forces is controlled by 0 (for the force displacement characteristics through each of the stable equilibria), and 2 for the force displacement characteristics through the transition state. In this section we will choose (x) to achieve an optimal design. It should be understood that the behavior of 2 and 0 is not independent, when considered as functionals of the beam shape (x). However, we will see in what follows that the freedom of changing (x) is sufficient to produce greatly improved designs. For the optimizations, we take the backbone shape w 0 to correspond to the lowest eigenmode of the structure. The original application of the switch as a component of a relay circuit had the requirements that: (a) f pop 20 mN, so that the relay contact would reliably close; (b) minimize both f push and f pop, for ease of both actuation and resetting the switch to its initial configuration; (c) the strain in the device cannot exceed the yield strain; (d) the switch must maintain bistability; and (e) fabrication constraints require that the thickness of the beam must be larger than a minimum feature size, e.g., 10 m. These conditions translate into finding (x) so that the force ratio R is near unity, minimizing the magnitude of the forces, while simultaneously maintaining the strain and bistability constraints. Formula 6 implies that the force ratio tends to unity as 2 0 3 . We therefore seek a beam shape (x) that maximizes 2 0 while maintaining the other constraints. To the optimal beam
Brenner et al.

0,

Eq. 8 implies that 1

w2 2 w2 2
0, 2,

2 0

w0 2 w0 2 [9]

w 0, w 2 ,

1 2 where wi 2 0 dx wi . To find a beam shape with a prespecified 2 0, we therefore use the following procedure: Starting with a candidate (x), we numerically calculate the eigenvalues and eigenfunctions [ i, w i]; then Eq. 9 implies that if we choose [ 0, 2, w 0, w 2] (with a small number), the new beam shape will have a larger 2 0. Iterating this procedure and enforcing any additional constraints produces a beam shape with maximal 2 0. Fig. 3 shows the evolution of (x) during the optimization. In these calculations, was chosen to ensure that 2 and 0 do not change by 10% from iteration to iteration. As the optimization proceeds, the minimum thickness of the beam decreases. Fig. 3 Lower shows the evolution of the predicted force ratio (Eq. 6) as a function of the minimum . Each of the structures with different min solves an optimal design problem: with decreasing min, the maximum strain in the structure during its operation increases; thus we expect that the optimal structure will be determined by the strain constraint. We will investigate how this maximum strain depends on min in what follows. Fig. 4 shows the actual thickness profiles of the beam for all three of the optimized shapes, corresponding to min 0.5, 0.25, and 0.15, respectively. The modulations in the thickness are extremely modest (note the 103-fold difference in the scales of

PNAS

August 19, 2003

vol. 100

no. 17

9665

APPLIED MATHEMATICS

Fig. 3. (Upper) Evolution of the moment of inertia of the beam during the optimization. The initial (x) is uniform; with increasing 2 0 the beam becomes modulated. Note that the beam thickness is exaggerated relative to its length by 103. (Lower) Evolution of the R ( 2 0)2 (( 2 0)2 2) as a function of min the minimum moment of inertia of the beam.

Fig. 5. (Top) Force displacement characteristics for the structure with min 0.5. R 1.43. (Middle) Force displacement characteristics for the structure with min 0.25. R 1.18. (Bottom) Force displacement characteristics for the structure with min 0.125. R 1.05.

the vertical and horizontal axis) but the change in the predicted force ratio is substantial. The eigenvalue ratio for these three shapes is 2 2 6, 9, and 90, respectively, corresponding to 2 0 predicted force ratios of R 1.5, 1.29, and 1.02, respectively. Fig. 5 shows the calculated force displacement characteristics of the optimized beam for the three optimized structures, with minimum thickness min 0.5, 0.25, and 0.125, respectively. The experimental force ratios in all three of these structures tend to be smaller than those predicted. For min 0.5, 0.25, and 0.125, the force ratios are R 1.43, 1.18, and 1.05, respectively. Additionally, for the min 0.125 design, the maximum force is smaller than that for the uniform beam: f push 718, compared with f push 800 for the uniform beam. Thus the optimized designs perform better than the specifications they were designed for! In all cases, the location of the transition state is exactly in accordance with the present analysis: for min 0.5, 0.25, and 0.125, respectively, the transition state moved from 4 3 3

8.65 7.65 1.13. Moreover, the force displacement characteristics are indeed linear through the transition state as anticipated above. The slope of this linear region tends to be slightly higher than that anticipated from the uniform beam calculation. The improvement in the force ratio occurs because the force displacement relation around the two stable equilibria is weaker than that of the uniform beam; a smaller force is required to maintain a given displacement near these stable equilibria. This causes the range of displacements around the transition state where the force is linear to decrease. As min decreases, both of these effects are more pronounced; the net result is to create both smaller overall forces and smaller force ratios than anticipated by the optimization procedure. This occurs because although our optimization procedure is designed to increase

Fig. 4. Optimized thickness prole of the beam for min 0.5, 0.25, and 0.15, with R of 1.5, 1.29, and 1.02, respectively. Note that the beam thickness is exaggerated relative to its length by 103.

Fig. 6. Experiments for the optimized beam with min 0.5, 0.25, and 0.125, respectively. The fabricated beam has length 15 mm, minimum thickness 20 m, and d 250 m. The optimization analysis designed these structures to have R 1.43, 1.18, and 1.05, respectively, compared with the measured values R 1.64, 1.37, and 1.28.

9666

[Link] cgi doi 10.1073 pnas.1531507100

Brenner et al.

2 0, it simultaneously decreases 0; this controls the strength of the force-displacement response around the stable equilibria. Finally, we address the maximum strain during the operation of the switch. The strain is given by zw , where z is the distance from the centerline of the beam. We are interested in computing the maximum strain max max(hw ), where h h(x) is the beam thickness. If we assume that the breadth b of the 1/3 , and the maximum strain structure is constant, then h max(hw ). In dimensional units, max h 0d L 2, where max h 0 is the thickness of the beam before optimization, and depends on min. We computed the maximum strain in each of the three optimized designs and found that 11.8, 13.8, 20.0, and 34 for min 1, 0.5, 0.25, and 0.125, respectively. As anticipated, the maximum strain increases with decreasing min. The beam shape that should be used for a given application depends largely on the strain constraint. Since the force ratio for 0.125 1 nearly achieves the optimum force ratio for min minimal strain, we present this as our optimal design.

profile of the beam in the transparency mask differs slightly from the predicted profile. Nonetheless, the experiments demonstrate the major point, that the predicted structures are systematically better than those originally fabricated, even without taking great care to control fabrication error. By designing beams with even smaller min we could presumably achieve force ratio even closer to unity even with fabrication error. As a final remark, it is worth noting that the calculations presented herein provide a prediction for the change in the force ratio given a certain level of fabrication error. Eq. 6 implies that the change in the force ratio R 2R x x, where x 2 0 . Eq. 9 states x . If we assume that the fabrication error is a zero mean white noise with variance , we therefore have 2 . For the optimal design with R R 2 0.125, this min gives R R 19.2 . A fabrication error of 5 m thus implies 0.3R 0.31. This the error in the force ratio should be R is consistent with the observed error (comparing the measured force ratio to the theoretical prediction). Summary In this article we have described a method for optimizing a bistable switch. The central idea is to use the unstable transition state connecting the two stable equilibria to tune the force displacement relation of the switch. By focusing on the properties of the transition state, the seemingly nonlinear problem of optimizing the nonlinear force-displacement characteristics of a switch is recast as a linear problem, whose solution can be addressed with standard methods from optimal structural design. The method produced a designed structure with the optimal force ratio, without any substantial penalty (e.g., the maximum strain that the device withstands, etc.) The specific problem we have described herein raises many general questions worthy of further study: Is it always possible to design the nonlinear force displacement properties of a switch by manipulating the properties of the transition state? Can switch design always be reduced to the analysis of a linear problem, or are some switches inherently nonlinear?
This research was supported by the National Science Foundation Division of Mathematical Sciences; National Science Foundation Division of Design, Manufacturing, and Industrial Innovation Award 9900792; and a research grant from ABB Research Ltd. (Zurich, Switzerland).
11. Pironneau, O. (1984) Optimal Shape Design for Elliptic Systems (Springer, New York). 12. Pironneau, O. (1974) J. Fluid Mech. 64, 97110. 13. Jameson, A. (2001) Prog. Aero. Sci. 37, 197243. 14. Glowinski, R. & Pironneau, O. (1975) J. Fluid Mech. 72, 385389. 15. Qiu, J., Lang, J. H. & Slocum, A. H. (2001) MEMS Digest, 353356. 16. Jensen, B. D., Parkinson, M. B., Kurabayashi, K., Howell, L. L. & Baker, M. S. (2001) Proc. Am. Soc. Mech. Eng. Intl. Mech. Eng. Congress, 17. 17. Hoffman, M., Kopka, P. & Voges, E. (1999) IEEE J. Selected Top. Quantum Electronics 5, 4651. 18. Zhou, S., Sun, X. & Carr, W. N. (1999) J. Micromech. Microeng. 9, 4550. 19. Jensen, B., Howell, L. & Salmon, S. G. (1999) J. Mech. Des. 121, 416423. 20. Saif, M. T. A. (2000) J. MEMS 9, 157170. 21. Guckenheimer, J. & Holmes, P. (1983) Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields (Springer, New York).

Experimental Tests The optimal structures that minimize the force ratio were fabricated by using deep reactive ion etching of a silicon wafer (15), and the force-displacement characteristics were measured by using a novel instrument for measuring stiffness of MEMS devices.** The results are summarized in Fig. 6 for the beam with 0.5, 0.25, 0.125 respectively. min The measured force ratios for min 0.5, 0.25, and 0.125 are, respectively, R 1.64, 1.37, and 1.28, compared with the predicted force ratios R 1.43, 1.18, and 1.05. The measured force ratios are higher than the predictions, though their relative values are in line with the predictions. The reason that the measurements are systematically higher than the predictions is that there are errors in the fabrication process that cause differences between the fabricated shape and the designed shape. Scanning electron microscopy shows that the fabricated structures tended to be overetched in the middle; in addition, the
There is a chance that designs with smaller min will require smaller overall forces to actuate; since increasing the force ratio further does have the advantage of further decreasing 0. **Li, J., Qiu, J., Lang, J. & Slocum, A., Proceedings of the 10th International Conference on Precision Engineering, July 18 20, 2001, Yokohama, Japan.

1. Banichuk, N. V. (1983) Problems and Methods of Optimal Structural Design (Plenum, New York). 2. Bendsoe, M. P. (1988) Comp. Methods Appl. Mech. Eng. 71, 197224. 3. Keller, J. B. (1960) Arch. Rat. Mech. Anal. 5, 275285. 4. Cox, S. J. (1992) Math. Intell. 14, 1624. 5. Keller, J. B. & Niordson, F. I. (1966) J. Math. Mech. 16, 433446. 6. Sigmund, O. (2000) Philos. Trans R. Soc. London A 358, 211227. 7. Frecker, M. I., Ananthasuresh, G., Nishiwaki, S., Kikuchi, N. & Kota, S. (1997) Trans. Am. Soc. Mech. Eng. 119, 238. 8. Frecker, M., Ananthasuresh, G. K., Nishiwaki, S., Kikuchi, N. & Kota, S. (1997) J. Mech. Des. Trans. Am. Soc. Mech. Eng. 119, 238245. 9. Kota, S., Hetrick, J., Li, Z. & Saggere, L. (1999) IEEE Am. Soc. Mech. Eng. Trans. Mech. 4, 396408. 10. Kohn, R. V. & Strang, G. (1986) Commun. Pure Appl. Math. 39, 113.

Brenner et al.

PNAS

August 19, 2003

vol. 100

no. 17

9667

APPLIED MATHEMATICS