<MNA/ODEs>

1. where <MNA/ODEs> G=conductance ce=constitutional eqs C=capacitance V=voltage Q=charge M=mass D=damping K=stiffness F=force q=displacement Electrical & Mechanical = Coupled MNA solution vector MNA matrix Excitation vector Multi-domain state vector Input coupling matrix System matrix 1st order ODE

2. Transient Analysis, (TA) SUGAR uses Matlab’s ODE solvers to integrate the 1st order system. Assuming the “mass” (G1) of the ode is constant, it’s stored into the solver. So we only need to return the RHS. Options may be chosen whereby stability and accuracy may be traded for speed. Mechanical example: Stored into ODE solver ODE integration

3. Gap Test Case

4. TA / Gap-Closing Actuator A) B) Transient response of a gap-closing actuator. A) shows a plot of displacement as a function of time. The voltage ramps from 5V at t=5usec to 12V at t=500usec, and then releases. As the voltage increases linearly during this time interval, the space between the gap decreases at a nonlinear rate due to electrostatic forces; likewise, the period of oscillation decreases. The amplitude of oscillations decrease exponentially due to the viscous layer of air between the device and substrate.

5. TA w/ Reduced Order Modeling Z. Bai ROM dim n ode dim N Reduced order modeling simulation (blue) superimposed onto our other transient solver (red) demonstrates its accuracy. The lower dimensionality of the reduced order system decreased simulation time by a factor 1/60. ROM via Krylov/Lanczos.

6. Sensitivity Analysis, (SA) L Schenato, W Wu Monte-Carlo Solves the above equation by drawing many samples from d probability distributions. It produces the most likely outcomes of performance variables. Ellipsoidal Calculus Finds the extreme bounds on performance parameters, i.e. the worst case scenario. EE245 project

7. Geometric Variation L Schenato, W Wu Worst case scenario via Ellipsoidal Calculus Performance possibilities via Monte Carlo [Knominal + DK] q = Fstatic

8. <Model Basics> Mechanical example 3D linear beam matrices Matrix coordinate transformation Model functions

9. EOM

10. Stiffness Matrix

11. Mass Matrix

12. The ith damping matrix is given by Di • = viscosity, D = fluid layer thickness, Di = damping matrix Couette Flow Damping Matrix

13. Coordinate Transformations

14. y x z oy y x z x x ox y y oz z z Rotation / Positioning node1 node2 Beam initiates along the X axis. The global XYZ and local xyz axes coincide. Think of the XY-plane is the substrate. Step 1: +Rotation about local y. Local x & z are repositioned in the global XY-plane. Step 3: +Rotation (twist) about local x. Local y & z are effected. Step 2: +Rotation about the new local z. Local y & x are effected. 3 positive rotations (right hand rule) are shown here. Important: rotations in SUGAR are performed in the order of y-z-x!

15. Model Functions function [output] = MF_beam(flag, R, params, x, t, nodes, varargin); switch(flag) case 'vars‘ output.dynamic = {1 {‘x’ ‘y’ ‘z’ ‘rx’ ‘ry’ ‘rz’}; 2 {‘x’ ‘y’ ‘z’ ‘rx’ ‘ry’ ‘rz’}}; case ‘M’ output = beam_matrices('M', params); case ‘D’ output = beam_matrices('D', params); case ‘K’ output = beam_matrices('K', params); case ‘F’ [F, dFdx] = compute_forces(flag, params, x, t); output = F; case ‘dFdx’ [F, dFdx] = compute_forces(flag, params, x, t); output = dFdx; ETC… E.g. Residual stress Thermal expansion Accelerating frames Electrostatic force Nonlinear effects Modelfunctions contain a model’s information. The output (e.g. stiffness, electrostatic force, stress) is determined by the flag.

16. <Accelerating Frames>

17. Understanding Inertial Forces Perpendicular to rotation axis only if node is moving; perpendicular to velocity If angular acceleration If substrate translates

18. Simple Gyro on a Rotating Substrate Z w 40us w(t)=1rad/sec E Y -X Ex(t) Ey(t) The gyro is first set in motion along the y-axis. The plots show the Ey and Exdisplacements of node E as functions of time. Midway through the simulation the left anchor is set to spin about the z-axis at 1 rad/sec. Though this angular velocity has almost no affect on Ey, Ex is significantly affected.

19. <Residual Stress>

20. Residual Stress Gauge 8.4um C Pan, W Hsu SUGAR simulation of residual stress. MEMS devices are often subject to residual stress effects which may affect device performance. Simulated deflection of this residual stress gauge is within 0.59% of measured data. Dy=8.35um

21. Residual Stress Gauge C Pan, W Hsu

22. ADXL-05 / BiCMOS Analog Devices Close-up view of the residual strain effects.

23. <Thermal Expansion>

24. T=600C T=150C Heatuator Simulation B Allen T=600C T=150C Dy=4.82u In the real device there’s a heat distribution along the hot and cold arms. Averaging the temperature along the beam produces the same linear thermal expansion. SUGAR is within 0.5% of the measured deflection done by B. Allen, Hilton Head 1998.

25. <Electrostatic Gap> This integral finds the equivalent nodal force and moments caused by the distributed surface forces, where H(x) is the Hermitian shape function. For this case, p(x) is the electrostatic pressure.

26. TA Pull-in (ramp) Abrupt V(t) = 11V to 15V, ramp. Nonlinear resonance, then pull-in. Smooth V(t) = 0V to 20V ramp. Pulled in before it got to 20V

27. <Digital to Analog Converter> R Yeh, KSJ Pister Using a subnet description for device building blocks, the nestlist description of this 4-bit MEMS DAC was reduced to just 14 lines text. SEM (left) and SUGAR visualization (right).

28. DAC Data (nonlinear beams) Output displacement vs digital input of a DAC. The nonlinear beam model matches measured data to within 5%.

29. <Mirror Scanner> KSJ Pister Two-degree-of-freedom optical scanner prototype. Mode 1 = 739Hz, mode 2 = 745Hz.

30. Projection Display Texas Instruments Electrostatic gap TI mirror for flat panel projection displays. Electrostatic gap electrodes (not visible) are underneath the plates. The mass is a little off due to the overlapping plates.

31. Mode Analysis of a Torsional Mirror BSAC c) a) d) b) The mode shapes and frequencies of modes 1, 3, 4 and 6 are shown in Figures a) through d). Respectively they are 15.5kHz, 31.1kHz, 41.7kHz, and 123kHz.

32. < Gyro Modes> Seshia, Howe Design by Inertial Sensors Group UCB

33. Mode Analysis of Accelerometer LLNL Design provided by: Jonathan Simon, LLNL a) b) c) a) schematic design. b) first mode shape corresponding to 27.73 Hz from SUGAR. c) third mode corresponding to 133.02 Hz from SUGAR, mMatches Simon’s analysis 132.46 Hz within 0.5%

34. Induced Current (Multimode resonator) R Brennen This demonstrates steady-state analysis applied to a multi-mode resonator. The left figure shows the Bode and phase plot of the current induced on the sensing comb as a function of the frequency of the voltage at the driving comb of the figure on the left. The measured modes are with in 5% of experimental frequency given by Brennen et. al.

35. y EI L x ds dy F0 d Y dx 1 2 L - D Y0 D <Nonlinear model> s = distance along the beam from node 1 to node 2 r(s) = radius of curvature at s y(s) = angle of the beam at s with respect to the horizontal M(s) = moment at s E = modulus of elasticity I = moment of inertia about the z-axis (out of plane) x and y are the in-plane horizontal and vertical axes with the origin at node 1. Fo = external force at node 2 L = beam length D = projected beam shortening

36. Nonlinear stiffness model To obtain nonlinear stiffness, we first assume that the curves can be approximated by a third order polynomial of the form were q stands for q, x, y, etc. Seeing that the solution has odd symmetry, we only need to keep constants B and D, which also eases iterative computations. Absorbing the material and geometric terms into B and D, respectively K1 and K2, we find that The coefficients of these polynomial curves are associated with the linear stiffnesses K1,i and the cubic nonlinearities K2,i. In order to maintain accuracy, it is applied in a continuous piecewise fashion by dividing the total physical range into, say, 4 intervals, were i(q)=[1,4]. (see plot)

37. Nonlinear Model Results vs Theory Circles are theory. The red, blue, & green curves are nonlinear model

38. <9-Node Plate Model> Excitation 1

39. 9-Node Plate Model, Excitation 2

40. Conclusion SUGAR is . . . Simple Text netlist description Matlab environment

41. Conclusion SUGAR is . . . Accurate Simulated results < 5% of actual device performance e.g. gap pull-in, Heatuator, residual stress, multimode resonator

42. Conclusion SUGAR is . . . Extensible User-definable modelfunctions Easy to add more Matlab functionality

43. Conclusion SUGAR is . . . FREE! www-bsac.eecs.berkeley.edu/~cfm cfm@bsac.eecs.berkeley.edu

44. some near-term goals CIF in/out, Squeeze film damping, Nonlinear stiffness, Curved beams, Plates, Buckling, Fringing fields, Trapezoidal beams, Initially loaded beams, Capacitive sensing, Friction, Collisions, Charge transfer, Noise analysis, Genetic Algorithms for design optimization, AutoMacroLevel Algorithms for large dynamic systems, Joints, Strain limiting, Heat transfer, Plasticity, Piezoelectric, Piezoresistive, Resistance(temperature), Thermoelectrics, Multiple time scale handling, Fluidic systems, Flight, Performance analysis, Assembly, Magnetics, Rigid, Bimetalic, Shape-memory, charge leakage, repulsion, Millennium web service, etc.