Parametric Design of a MEMS Accelerometer
ME 128 – Project 2
Professor Liwei Lin Spring 2006
Scott Moura SID 15905638
April 5, 2006
ME128 – ComputerAided Mechanical Design Spring 2006
Name: Project:
Scott Moura #2 Design of a MEMS Accelerometer
Introduction: Theory: Code Verification: Summary of Results: Results: FEM Results: Conclusions: Computer Source Code
10 20 10 10 10 10 20 10
____________ ____________ ____________ ____________ ____________ ____________ ____________ ____________
Total:
100
____________
April 5, 2006
Table of Contents
Table of Contents..................................................................................................1 Introduction ...........................................................................................................2 Nomenclature .......................................................................................................3 Variables ...........................................................................................................3 Subscripts .........................................................................................................3 Background...........................................................................................................4 Theory...................................................................................................................5 Relevant Equations ...........................................................................................6 Provided Derivations .........................................................................................6 Constraints & Engineering Goals ......................................................................8 Genetic Algorithm Optimization Method................................................................9 Code Verification.................................................................................................11 Summary of Results............................................................................................12 Theoretical Results .............................................................................................13 FEM Results .......................................................................................................16 Case 1 – Sensor Motion Resolution ................................................................17 Case 2 – Acceleration Survival........................................................................19 Case 3 – Maximum Displacement...................................................................21 Discussion & Conclusions...................................................................................23 Theoretical Analysis vs. FEA...........................................................................23 Error due to Dimension Magnitude..................................................................24 Anchor/Proof Mass Clearance.........................................................................24 Die Area and Acceleration Resolution.............................................................25 ANSYS Complications.....................................................................................25 Summary Conclusions ....................................................................................27 Source Code .......................................................................................................28 gaopt.m ...........................................................................................................28 enginprop.m ....................................................................................................28 LOCALden.m...................................................................................................28
1
Introduction
The objective of this investigation is to find the optimum design for a typical MEMS (Micro ElectroMechanical Systems) accelerometer, which satisfies a set of given constraints. Due to the complex nature of the problem, a genetic
algorithm (GA) is developed for optimization. The GA attempts to minimize the die area while satisfying all other engineering goals. Four major dimensions (L1, L2, L3, ym) are determined from this optimization. The optimal design from the theoretically derived genetic algorithm is compared to hand derived calculations and finite element analysis in order to ascertain its accuracy and verify the results.
The genetic algorithm, developed in MATLAB, utilizes concepts from evolution in order to ascertain the best performing design. Although computationally intense, this method is very useful for complex problems in which system of differential equations are difficult to write. All the details for implementing this method are provided in the optimization method section.
The software used for finite element analysis (FEA) is SolidWorks with the builtin COSMOSWorks finite element method (FEM) analysis package. A three
dimensional model of the best design is created and analyzed under three different loading conditions. These cases represent the limiting conditions of operation, and are therefore of interest for evaluating resolution and survival. Unlike the theoretical calculations, this method considers the mass of the beams, whereas the MATLAB computations assume the beam masses are negligible compared to the proof mass size.
2
Nomenclature
Variables
A b E F h I k L m M x λ ∧ ∏ θ Φ ρ σ Die Area, [m2] Beam Depth, [m] Young's modulus, [psi] Force [N] Beam Width [m] Moment of Inertia about xaxis, [in4] Spring Constant for entire system [N/m] Beam length, [in] Proof Mass size [kg] Moment, [inlb] Position along length of beam, [in] Children Design Variables Design Variables, Parent Objective Function Angle of deflection, [rad] Random Value Density [kg/m3] Stress, [GPa]
Subscripts
a m x y Anchor Proof Mass xdirection ydirection
3
Background
MEMS accelerometers are used for a variety of applications, namely automobile airbag systems. Consumer products, such as computer games, cell phones, pagers, PDAs, advanced robotics, laptop computers, computer input devices, camcorders, digital cameras, and aftermarket SD card accessories are also common applications1. In each of these applications, size and accuracy are the most critical characteristics of the sensor's performance. As such, these factors will be considered throughout this investigation.
The MEMS accelerometer under study operates under the same principles of a springmass system, shown schematically in Figure 1. However, instead of springs, the accelerometer employs a double folded beam flexure system. The mass being displaced is the proof mass. To measure displacement, one capacitive sensor exists on each side of the proof mass. The sensitivity of these sensors is proportional to the size length of the mass. These sensors send back a voltage signal proportional to the displacement measured. By equating Hooke's Law to Newton's second law, kx = ma, the acceleration experienced by the mass can easily be determined.
Figure 1: Schemetic Diagram of SpringMass System, with sensor measuring displacement.
A representative diagram of the MEMS accelerometer is provided in Figure 2. The narrow beams of the assembly make up the equivalent spring for the entire system. The proof mass is shown in the middle and the anchors are restrained
Kionix. MEMS AccelerometersInertial Sensors. Accessed April 2, 2006 <http://www.kionix.com/Accelerometers/accelerometers.htm>
1
4
to the MEMS substrate. Ideally, the legs (beams) of the structure would be made as long as possible for softer springs and the mass as large as possible. This would result in the greatest displacement and thus maximize acceleration resolution. However, a maximum die area constrains the size of the
accelerometer, which we seek to minimize in order to produce more devices per wafer. As a result, these two factors fight each other when finding the optimal design.
Figure 2: MEMS Accelerometer Diagram
Theory
The governing differential equations for elastic beam bending serve as the basis for the theoretical analysis on the device. Since the focus of this project is
parametric optimization and not an indepth review of beam force analysis, all relevant equations and derivations are provided and shown below. They utilize the concepts of superposition to sum forces and moments at each beam's intersection. The results include the forces and moments at the points of 5
maximum stress (A and H) and the total effective spring constant for the entire system. These equations are later used to analyze the beams properties for a given set of design parameters.
Relevant Equations
Provided Derivations
6
7
Constraints & Engineering Goals
Coupled with the equations shown above are both equality and inequality constraints. These parameters are requirements determined by the
manufacturing process, cost constrictions, material properties, and design goals. They are given by the following equations:
8
Equality Constraints
Inequality Constraints
h1 = h2 = h3 = hb = 1.8 m b1 = b2 = b3 = bm = 1.8 m Lb = 150 m E = 160GPa g ρ = 2.33 3 cm y a = 100 m x a = Lb + 2h1 x m = Lb + 2 L2 + 2h1 + 2h3
L1 ≥ 20 m L2 ≥ 20 m L3 ≥ 20 m L3 ≥ h2 + y a + L1
With these constraints in mind, this investigation seeks to accomplish certain engineering goals. following: Minimum DC acceleration resolution: a r ≤ 0.0005 g Maximum DC acceleration survival: a max ≥ 2000 g Maximum die area: Ad = x m ( y m + 2 L3 + 2 h2 ) ≤ 90 ,000 m 2 Maximum stress in suspension: σ max ≤ 1.6 GPa Motion Resolution: x r The goals provided by the project descriptions are the
(10 =
1
m )
2
ym
After a thorough analysis of the various designs, it becomes clear that the maximum die area becomes the most significant limiting constraint. In fact, it is impossible to satisfy all of the constraints with a design that is less than 90,000m2 in area. As a result, this constraint is relaxed and the optimum design is found. Detailed descriptions of these findings are discussed in this report.
Genetic Algorithm Optimization Method
In many cases where parametric optimization is desired, methods such as the Lagrange Multiplier, Firstorder Necessary Conditions (FONC), and Secondorder Necessary Conditions (SONC) are preferred for analytical solutions. However, this design problem contains features too complicated to represent in a
9
system of differential equations.
Moreover, determining inactive constraints a result, a
through monotonicity is nearly impossible using logic tables. As genetic algorithm is developed and employed.
Although computationally
intensive, the basic concept behind the genetic algorithm is simple. This method also allows the freedom to change any of the constraints with relative ease.
The employment of the genetic algorithm follows a very simple iterative technique to minimize an objective function, given by ∏ . The details of this objective function will be explained later. The design variables are represented by ∧ .
∧ = {L1 , L2 , L3 , y m }
To apply the GA we, we take random values for the design variables within the following ranges,
20 m ≤ L1 ≤ 500 m 20 m ≤ L2 ≤ 100 m 100 m ≤ L3 ≤ 500 m 100 m ≤ y m ≤ 500 m
These ranges were chosen based on the minimum size constraints and maximum area constraints, in addition to general observation and intuition about the final design's optimal geometry.
The algorithm begins with populations of 100 random strings within the aforementioned ranges. evaluated. For each string, the following objective function is
∏ =A
where A represents the die area. Clearly, this method will find the beam lengths that minimize the total die area, which is only one of the design goals. However, it can be argued that die area is the most significant goal, since it directly relates to cost. By minimizing the die area we maximize the number of accelerometers that can be manufactured on the same silicon wafer. This means more
accelerometers for the same price. Also initial calculations have shown that this
10
constraint is the most difficult to satisfy. Therefore, this simple objective function is chosen. The other four design goals are evaluated for each design and thrown out if they do not meet the inequality constraints. This is done by giving ∏ a value of infinity. Once the design goals are checked, the values of ∏ are sorted to determine the top ten performing designs (10 smallest ∏ 's). These top ten 'parents' are then 'mated' to produce ten 'offspring' or 'children' design values. scheme is used to perform this 'mating' task, The following
λi = Φ ( I ) ∧i + (1 Φ ( I ) ) ∧i+1 λi+1 = Φ ( II ) ∧i + (1 Φ ( II ) ) ∧i+1
where
λi
are the children design values,
∧i
are the parent design values, and
Φ ( I ) , Φ ( II )
are random values between zero and one.
At this point, we take the ten parents and ten children and combine them with 70 new random design variables to create a second generation of 100 strings. With these 100 new strings, the objective function is calculated again, checked with the constraints, sorted, and so forth. This process is run for 50 generations (iterations). At the end of the 50 generations, the top ranked design variables give the values of L1, L2, L3, and ym. On top of performing 50 iterations of this genetic algorithm, we perform this entire scheme 1000 times, thus allowing 1000 different starting populations. As a result, 1000 different top performing design values are found for 1000 different starting populations. The top value
represents the design variables used for the optimum design, recommended by this report.
Code Verification
The genetic algorithm written in MATLAB is compared to hand derived calculations to prove accuracy within the code. The results for the most
11
significant accelerometer properties are shown in Table 1 for the optimum design variables. The values match perfectly, thus ensuring the program's precision.
Table 1: Code Verification Table.
Die Dimensions MATLAB Calculations Hand Calculations 197.2 m x 785.9 m 197.2 m x 785.9 m
Die Area
Proof Mass Size
Effective Spring Constant 0.32911 N/m 0.33 N/m
154,979.48 m2 0.23629 g 155,000 m2 0.2363 g
Summary of Results
The genetic algorithm developed finds the optimum design that minimizes die area while just satisfying all the other design constraints. The results are given in Table 2.
Table 2: Results of the Genetic Algorithm Optimization Method.
Optimum Design Property L1 L2 L3 ym Die Area Motion Resolution Acceleration Resolution Acceleration Survival Max Stress in Suspension Die Dimensions Proof Mass Size Effective Spring Constant
Value 146.0 m 20 m 248.3 m 285.7 m 154979.48 m2 3.5002 x 105 m 0.0049747g 20998.1619g 0.2166 GPa 197.2 m x 785.9 m 0.23629 g 0.32911 N/m
Initial calculations have shown that it is impossible to satisfy all of the constraints, as they are given in the project description. Specifically, the 90,000 m2
maximum die area and 0.0005g acceleration resolution are the limiting
12
constraints. As a result, the die area constraint is relaxed and minimized for values below 160,000 m2 while the minimum acceleration resolution is increased by a factor of ten. Although it is possible to satisfy the minimum
acceleration resolution given, it is found that changing this constraint can significantly minimize die area. The cost savings are chosen in favor of this engineering goal. Additionally, an acceleration resolution of 0.005g still satisfies most applications. The engineering goals and their corresponding constraints are shown below. The necessary modifications for the design are shown in the rightmost column.
Table 3: Design Goal Metric Table. The die area and acceleration resolution constraints are relaxed in order to minimize cost and satisfy the other constraints.
Design Goals L1 L2 L3 Die Area Motion Resolution Acceleration Resolution Acceleration Survival Max Stress in Suspension
Optimal Design 146.0 m 20 m 248.3 m 154979.48 m2 3.5002 x 105 m 0.0049747g 20998.1619g 0.2166 GPa
Initial Constraint > 20 m > 20 m > 20 m > h2 + ya + L1 < 90,000 m2 < 3.518 x 105 m < 0.0005g 2,000g 1.6 GPa
Modified Constraint < 160,000 m2 < 0.005g 
Theoretical Results
The theorybased genetic algorithm is designed using the MATLAB programming environment. The iterative genetic process minimizes the die area, represented by the objective function, ∏ , while satisfying all other design criteria. The best performing design is saved for each successive starting population to converge on the optimum values. Figure 3 illustrates this fact by displaying the minimum value of the objective function for the first 100 starting populations. Clearly, the genetic algorithm succeeds in progressively finding designs with smaller design
13
areas. Additionally, the genetic algorithm appears to converge to the best design asymptotically, as the starting population count increases.
Figure 3: Minimization of Objective Function for first 100 starting populations. The value of ∏ decreases as better performing designs are found.
After 1,000 starting populations of 50 generations (iterations) have been computed, the five best performing designs are output to the user. The final results are shown in Table 4. Note that L2 converges to 20 m, the minimum value possible. Although 20 m is never achieved exactly (due to the exclusive nature of the random generator function) it can be assumed that the optimum design has L2 = 20 m. Conversely, L1, L2, and ym do not appear to converge to a value. This must imply that there is a range of optimum values that can be used to achieve the best design. As a result, the optimum dimensions presented here are only one set of the possible values.
14
Table 4: Five Best Performing Designs computed using the Genetic Algorithm.
Rank 1 2 3 4 5
L1 (m) 146.0 142.6 128.8 117.7 159.7
L2 (m) 20.62 20.48 21.16 20.57 20.43
L3 (m) 248.3 253.0 235.7 232.8 267.1
ym (m) 285.7 281.6 312.3 324.2 259.8
∏ = A (m2) 154979.48639 156774.86310 157088.24299 157360.54186 157965.94329
The first ranked design is the one evaluated and recommended in this report. In order to analysis this design, three different operating conditions are considered. They include the sensor motion resolution (~0.005g), acceleration survival (2000g), and maximum displacement due to geometry constraints (20 m). For each condition, the displacement of the proof and mass and maximum stress is determined. The first test analyzes the accelerometer when it experiences a load just large enough to register a signal from the sensor. The second test considers the accelerometer's reaction to experiencing a load of 2000g, which is the minimum value for survival. The final test investigates the maximum stress in suspension of the accelerometer at the maximum displacement physically possible. Any displacement greater than 20 m would cause the beams to crash into the proof mass. The theoretical values for the displacement of the proof mass and maximum stress are shown in Table 5. Theses values are validated through finite element analysis.
Table 5: Theoretical Values of Proof Mass Displacement and Maximum Stress.
Case
Condition
Displacement (m) 3.5002 x 105 14.072 20
Stress (GPa) 3.81 x 107 0.1524 0.2166
1 2 3
Sensor Motion Resolution (~0.005g) Acceleration Survival (2000g) Maximum Displacement
15
FEM Results
The finite element analysis was performed within SolidWorks using the COSMOSWorks package. This program allows the user to select a custom material for the part by inputting userdefined values for the material properties. The material properties used for this analysis are given by the equality constraints. Additionally, a value of zero was assumed for the Poisson's ratio. Together, all of the material properties are
E = 160GPa kg ρ = 2330 3 m ν =0
A few simple techniques are applied to simulate the accelerometer under various loads (Figure 4). The two anchors are constrained to the substrate by placing fixed restraints on their bottom surfaces, shown by the green arrows. This
process eliminates all six degrees of freedom (DOFs). Simulating acceleration is an equally simple process: define the magnitude of gravity in the xdirection equal to the acceleration in question (shown schematically by the red arrow). Using these techniques, the three conditions are analyzed. The mesh size used for all case studies is equal to 6.6 m.
Figure 4: Accelerometer Model is SolidWorks COSMOS simulation study.
16
Case 1 – Sensor Motion Resolution
In the first case, the proof mass undergoes an acceleration of 0.005g, which roughly generates a displacement equal to the sensor's motion resolution. Figures 5 and 6 respectively show the maximum displacement and stress experienced by the device. The values computed by the COSMOS FEA are 3.546 x 105 m for the displacement and 3.766 x 107 GPa for the maximum stress.
Figure 5: Proof Mass Displacement at 0.005g.
Figure 6: Maximum Stress at 0.005g.
A more detailed view of the area under maximum stress is shown in Figure 7. Note that this area is the same area for which we assumed maximum stress would occur.
17
Figure 7: Detailed view of Maximum Stress at 0.005g.
18
Case 2 – Acceleration Survival
One of the engineering goals is to ensure that the device material will not fail for accelerations greater than 2000g. The theoretical calculations performed by
MATLAB state that failure occurs at an acceleration of about 21,000g, over 10 times greater than the minimum required value. Nonetheless, we wish to verify the theoretical calculations for stress and displacement at 2000g by performing FEA. The results from this analysis are shown in Figures 8 and 9.
Figure 8: Proof Mass Displacement at 2000g.
Figure 9: Maximum Stress at 2000g.
Figure 10 shows a detailed view of the area of maximum stress. It is evident from this test that the suggested accelerometer design meets the required acceleration survival constraint by over a factor of ten.
19
Figure 10: Detailed view of Maximum Stress at 2000g.
20
Case 3 – Maximum Displacement
As discussed previously, the proof mass may experience a displacement up to 20 m before the beams crash into the anchors. This limiting case provides useful information about the maximum possible stress. In order to model a
displacement of 20 m, MATLAB computations are performed to find the corresponding acceleration. This acceleration (2842.5279g) is input as the
magnitude of gravity in order to model a displacement of 20 m. Figure 11 exposes the approximate nature of this test case, since the FEA displacement result is slightly larger than 20 m at about 20.16 m. Note also how the beams come into contact with the anchors, thus confirming our theoretical computations.
Figure 11: Proof Mass Displacement at Maximum Possible Displacement.
Figure 12: Maximum Stress at Maximum Possible Displacement.
21
Shown in Figure 13 is a detailed view of the Von Mises near the proof mass corner. This stress is the maximum stress that the accelerometer can possibly experience, which easily satisfies the maximum stress constraint of 1.6 GPa by a factor of 7.5.
Figure 13: Detailed view of Maximum Stress at Maximum Possible Displacement.
22
Discussion & Conclusions
Theoretical Analysis vs. FEA
As described in the results, the MATLAB program agrees very well with the hand derived calculations. Both methods also agree with FEA for displacement and stress distributions. The results from both the theoretically driven genetic
algorithm and finite element analyses are provided in Table 6. Clearly, both methods provide very close results that help verify each method and satisfy the engineering constraints and goals.
It should be noted that the FEM analysis was performed with the default tetrahedral global mesh size of 0.6 m. This value is the maximum size able to accurately mesh the smallest feature of the device. As a result, controls meshes and/or finer meshes are unnecessary as the standard 0.6 m proves to be more than adequate.
Table 6: GA and FEM Results Comparison.
Analysis Method Displacement (m) Stress (GPa) Displacement (m) Stress (GPa) Displacement (m) Stress(GPa)
GA 3.518 x 105 3.81 x 107
FEM 3.546 x 105 3.766 x 107 14.18 0.1506 20.16 m 0.2141 GPa
Case 1 – Sensor Motion Resolution
Case 2 – Acceleration Survival 14.072 0.1524 Case 3 – Maximum Displacement 20 m 0.2166 GPa
Another notable difference between the two methods is that the theoretical method ignored the mass of the beams, whereas the FEA does account for this additional mass. Although this is a critical assumption for developing the EulerBeam equations, it appears that the proof mass dwarfs the mass of the beams. As a result, assuming that the beams' masses are negligible is verified.
23
For additional quantifiable verification of the error between the two methods, Table 7 is provided below. For all cases the error difference never exceeds
±1.2%, a remarkably accurate result.
Table 7: GA and FEM Results Comparison with Percentage Difference
Theoretical Value Percent Difference Analysis Method GA 3.518 x 105 m 3.81 x 107 GPa FEM
Case 1 – Sensor Motion Resolution Displacement Stress 0.7959 % 1.155 %
Case 2 – Acceleration Survival Displacement Stress 14.072 m 0.1524 GPa 0.7675 % 1.181 %
Case 3 – Maximum Displacement Displacement Stress 20 m 0.2166 GPa 0.8 % 1.154 %
Error due to Dimension Magnitude
Another minor source of error in both methods is roundoff and truncation error. The theoretical calculations were all performed in MKS units and therefore variable values existed between a range of 109 to 109. This broad range of numbers could easily lead to calculation errors. In contrast, COSMOS
calculations were performed in the MKS unit system. These calculations were characterized by a much smaller range of values, so roundoff and truncation error is less likely.
Anchor/Proof Mass Clearance
In an effort to minimize the die area and maximize the resolution, the genetic algorithm suggested legs as long as possible (softest springs), and a mass as large as possible. This fact leads to a minimal clearance between the anchor and mass. Although a structure characterized by these dimensions will produce
24
the most sensitive accelerometer possible, this small value may be an unacceptable distance for the clearance between a moving and stationary part. Although the analysis performed for this project only considered beam bending, the beam legs may also be susceptible to buckling, thus allowing the mass to crash into the anchor. This situation could occur if the device was dropped, thus creating damage. However, if the accelerometer is guaranteed to only encounter accelerations in onedirection, this is of no concern. The design proposed in this report has a clearance of only 2.3 m between the anchor and mass. However, if a greater clearance is required, the value of L1 can be decreased without a significant difference in resolution and no difference in die size.
Die Area and Acceleration Resolution
Initial optimization attempts revealed that satisfying all of the constraints and goals is nearly impossible with a maximum die size of only 90,000 m2. As a result, this constraint was relaxed and the optimization method attempted to minimize the area. Without modifying the acceleration resolution, it is possible to satisfy all of the constraints and engineering goals. However, the resulting In comparison, design will produce a die area of at least 240,000 m2.
increasing the acceleration resolution by a factor of ten produced a minimum die area of 154,979.48 m2. This is a decrease in die area by a factor of over 1.5, which is directly related to a decrease in cost by a factor of 1.5. Additionally, background research for most accelerometer applications proved that a resolution of 0.005g was more than adequate. Therefore, it appears preferable to relax the resolution constraint in favor of increasing cost savings, which should more than offset the increased resolution in terms of significance to the manufacturer/business.
ANSYS Complications
It was my original intention to provide a detailed analysis of the accelerometer design using ANSYS to augment the existing FEA performed with COSMOS. In fact, I spent a significant portion of my Spring Break in Etcheverry working with
25
ANSYS.
However, I was not successful in obtaining meaningful results that
approached values obtained theoretically and through COSMOS. I believe my problem had to do with the units input for the material properties and accelerations. First, I created a model in SolidWorks with the MKS units as suggested by the professor's email. After importing the IGES file into ANSYS and using MKS units for the material properties and accelerations, I received results that were off by 6 orders of magnitude. The GSI had told me that there were errors for exporting MKS models and suggested that I use MMKS. I did as she said and used MMKS units for the material properties and accelerations, but still obtain unsatisfactory results. The plots for displacement and stress for an acceleration of 2000g are shown in Figures 14 and 15, respectively. These
results are the best I could obtain with the limited ANSYS training we received. As a result, it is my recommendation that more lectures are dedicated to ANSYS training for future classes, as I obviously did not have the skills required to perform analysis on a MEMS scale device.
Figure 14: ANSYS Results for Displacement at 2000g.
26
Figure 15: ANSYS Results for Stress at 2000g.
Summary Conclusions
Using a theoretically derived genetic algorithm for optimization and FEM analysis, it is concluded that satisfying the die area requirement is impossible. Therefore, this constraint is relaxed and minimized by the software. In addition, the acceleration resolution is relaxed in favor of decreasing the die area, which improves the cost per unit ratio. As a result of the theoretical and FEM analyses, an accelerometer with the following dimensions is recommended: L1 = 146.0 m L2 = 20.0 m L3 = 248.3 m ym = 285.7 m This design meets the resolution criteria and exceeds the maximum stress criteria. In addition, the die area is exactly 154,979.48 m2, which is a 1.722 times increase above the maximum die area given. Both methods prove that this design is the most cost efficient while simultaneously satisfying all of the required engineering goals and constraints. 27
Source Code
C:\Documents and Settings\Scott Moura\My Documents\S...\gaopt.m April 5, 2006
%% Project #2  Parametric Design of a MEMS Accelerometer % Scott Moura % SID 15905638 % ME 128, Prof. Lin % Due Wed April 5, 2006 %% Genetic Algorithm clear %% Define Equality Constraints h_1 = 1.8e6; % Beam Width [m] h_2 = 1.8e6; % Minimum Feature Length of Device h_3 = 1.8e6; h_b = 1.8e6; b_1 b_2 b_3 b_m = = = = 1.8e6; 1.8e6; 1.8e6; 1.8e6; % Beam Thickness [m] % Minimum Feature Length of Device
Page 1 5:26:41 AM
L_b = 150e6; E = I = rho g = 160e9; h_1*b_1^3/12; = 2.33*100^3/1000; 9.8;
% Central Beam Width [m] % % % % Elastic Modulus [Pa] Area Moment of Inertia [m^4] Density [kg/m^3] (converted from gm/cm^3) Gravitational Acceleration [m/s^2]
y_a = 100e6; x_a = L_b + 2*h_1;
% Anchor Height [m] % Anchor Width [m]
A_max = 250000/(1e6)^2; % Maximum Die Area [m^2] %% Define Inequality Constraint Limits L_min = 20e6; % Minimum beam length a_r_min = 0.005*g; % Accerlation at minimum required resolution [m/s^2] a_stress = 2000*g; % Acceleration at maximum required stress [m/s^2] sigma_max = 1.6e9; % Maximum stress [Pa] %% Define L_1_min = L_1_max = L_2_min = L_2_max = L_3_min = L_3_max = y_m_min = y_m_max = boundary limits L_min; 500e6; L_min; 100e6; 100e6; 500e6; 100e6; 500e6;
% Initialize best1 & best2 start_pops = 1000; iters = 50; best1 = zeros(iters,start_pops); best2 = zeros(1,iters); % Loop for 1000 different starting populations for idx1 = 1:start_pops
Source Code
C:\Documents and Settings\Scott Moura\My Documents\S...\gaopt.m April 5, 2006
% Define Set of Random Strings d = 100; % Number of clear Lambda random = rand(d,4); Lambda(:,1) = L_1_min + random(:,1) * Lambda(:,2) = L_2_min + random(:,2) * Lambda(:,3) = L_3_min + random(:,3) * Lambda(:,4) = y_m_min + random(:,4) * strings
Page 2 5:26:41 AM
(L_1_max (L_2_max (L_3_max (y_m_max

L_1_min); L_2_min); L_3_min); y_m_min);
% % % %
L_1 L_2 L_3 y_m
% Repeat Genetic Algorithm for 50 iterations for idx2 = 1:iters L_1 L_2 L_3 y_m = = = = Lambda(:,1); Lambda(:,2); Lambda(:,3); Lambda(:,4);
% Calculate design variable dependent parameters x_m = L_b + 2*L_2 + 2*h_1 + 2*h_3; % Proof mass width [m] A = x_m .* (y_m + 2*L_3 + 2*h_2); % Die Area [m^2] m = rho * x_m .* y_m * b_m; % Mass of system (ignore beams) [kg] s_res = (1e7)^2 ./ y_m; % Motion resolution [m] % System Spring Constant [N/m] k = 24*E*I * (6*L_1.*L_2.^2 + L_b*L_2.^2 + 4*L_b*L_1.*L_2 + 24*L_1.*L_2.*L_3 + ... 12*L_b*L_1.*L_3 + 4*L_b*L_2.*L_3) ./ LOCALden(L_1,L_2,L_3,L_b); a_res = (k ./ m) .* s_res; s_stress = (m ./ k) .* a_stress; % Minimum resolvable acceleration [m/s^2] % Displacement at maximum survival accel [m]
% Moment at max survival accel [N*m] M = 6*s_stress*E*I .* (6*L_b*L_1.*L_3.^2 + 12*L_1.*L_2.*L_3.^2 + ... 2*L_b*L_2.*L_3.^2 + 6*L_1.*L_2.^2.*L_3 + L_b*L_2.^2.*L_3 + ... 4*L_b*L_1.*L_2.*L_3 + L_b*L_1.^2.*L_2) ./ LOCALden(L_1,L_2,L_3,L_b); % Forces at max survival accel [N] F_x = (k .* s_stress) / 4; F_y = 18*s_stress*E*I ./ L_2 .* (L_b*L_1.^2.*L_2 + 2*L_b*L_1.*L_3.^2  ... 2*L_b*L_1.^2.*L_3 + 6*L_1.*L_2.*L_3.^2 + L_b*L_2.*L_3.^2) ./ ... LOCALden(L_1,L_2,L_3,L_b); % Stresses at max survival accel [N/m^2] sigma_x = F_x / (h_1*b_1); sigma_y = F_y / (h_1*b_1); sigma_bend = M * (b_1/2) / I; % Calculate the Fitness of the Objective Function Pi = A*1e12; % Check Constraints for idx3 = 1:d % Anchor cannot petrude through proof mass if (L_3(idx3) < h_2 + y_a + L_1(idx3)) Pi(idx3) = inf; end
Source Code
C:\Documents and Settings\Scott Moura\My Documents\S...\gaopt.m April 5, 2006
% Check Die Area if (A(idx3) > A_max) Pi(idx3) = inf; end % Check minimum motion resolution if (a_res(idx3) > a_r_min) Pi(idx3) = inf; end % Check maximum stress sigma_norm = norm([sigma_x(idx3) sigma_y(idx3) sigma_bend(idx3)]); if (sigma_norm > sigma_max) Pi(idx3) = inf; end % Check if beams crash into anchor if (s_stress(idx3) > L_2(idx2)) %disp(s_stress(idx3)) Pi(idx3) = inf; end end % Rank the Genetic Strings output by Pi and save the best [sorted_Pi, order_Pi] = sort(Pi); if (sorted_Pi(1) == inf) %disp('INFINITY') end best1(idx2,idx1) = sorted_Pi(1); % Keep the Ten Best Parents nps = 10; % Number of Parents for i = 1:nps parent(i,:) = Lambda(order_Pi(i),:); end % Mate the Top Ten Parents to Generate Ten Children random2 = rand(nps,1); for j = 1:2:nps1 child(j,:) = random2(j,1) * parent(j,:) + (1  random2(j,1)) * parent(j+1,:); child(j+1,:) = random2(j+1,1) * parent(j,:) + (1  random2(j+1,1)) * parent(j +1,:); end % Generate New Random Strings to Combine with Parents and Children clear Lambda_new random3 = rand(d2*nps, 4); Lambda_newrand(:,1) Lambda_newrand(:,2) Lambda_newrand(:,3) Lambda_newrand(:,4) = = = = L_1_min L_2_min L_3_min y_m_min + + + + random3(:,1) random3(:,2) random3(:,3) random3(:,4) * * * * (L_1_max (L_2_max (L_3_max (y_m_max L_1_min); L_2_min); L_3_min); y_m_min); % % % % new new new new L_1 L_2 L_3 y_m
Page 3 5:26:41 AM
Lambda_new = vertcat(parent, child, Lambda_newrand); % Repeat Genetic Algorithm
Source Code
C:\Documents and Settings\Scott Moura\My Documents\S...\gaopt.m April 5, 2006
Lambda = Lambda_new; end % Save the best lambda values Lambda_final(idx1,:) = Lambda(1,:); end % Determine the runs that had the six best performing best2 = best1(iters,:); [sorted_best2, order_best2] = sort(best2); % Output the design variable values for top five for k = 1:5 best_Lambda(k,:) = Lambda_final(order_best2(k),:); fprintf('Rank: %2.0f L_1 = %6.3e L_2 = %6.3e L_3 = %6.3e y_m = %6.3e Pi = % 6.5f\n',... k,best_Lambda(k,1),best_Lambda(k,2),best_Lambda(k,3),best_Lambda(k,4),sorted_best 2(k)); end % Plot the minimum objective value as a function of generations for idx4 = 1:100 if (idx4 == 1) Pi_min(idx4) = best2(idx4); elseif (best2(idx4) <= Pi_min(idx41)) Pi_min(idx4) = best2(idx4); else Pi_min(idx4) = Pi_min(idx41); end end plot(1:100,Pi_min(1:100)) title(['\bfMinimization of Objective Function: \Pi_{\it{min}} = ' num2str(sorted_best2(1) )]) xlabel('Number of Starting Populations (First 100 Only)') ylabel('Minimum Objective Value Found')
Page 4 5:26:41 AM
C:\Documents and Settings\Scott Moura\My Documen...\enginprop.m April 5, 2006
function output = enginprop(L_1,L_2,L_3,y_m) %% Project #2  Parametric Design of a MEMS Accelerometer % Scott Moura % SID 15905638 % ME 128, Prof. Lin % Due Wed April 5, 2006 %% Calculate Engineering Properties %% Define Equality Constraints h_1 = 1.8e6; % Beam Width [m] h_2 = 1.8e6; % Minimum Feature Length of Device h_3 = 1.8e6; h_b = 1.8e6; b_1 b_2 b_3 b_m = = = = 1.8e6; 1.8e6; 1.8e6; 1.8e6; % Beam Thickness [m] % Minimum Feature Length of Device
Page 1 5:29:58 AM
L_b = 150e6; E = I = rho g = 160e9; h_1*b_1^3/12; = 2.33*100^3/1000; 9.8;
% Central Beam Width [m] % % % % Elastic Modulus [Pa] Area Moment of Inertia [m^4] Density [kg/m^3] (converted from gm/cm^3) Gravitational Acceleration [m/s^2]
y_a = 100e6; x_a = L_b + 2*h_1;
% Anchor Height [m] % Anchor Width [m]
%% Define Inequality Constraint Limits L_min = 20e6; % Minimum beam length a_r_min = 0.005*g; % Accerlation at minimum required resolution [m/s^2] a_stress = 2000*g; % Acceleration at maximum required stress [m/s^2] sigma_max = 1.6e9; % Maximum stress [Pa] % Calculate design variable dependent parameters x_m = L_b + 2*L_2 + 2*h_1 + 2*h_3; % Proof mass width [m] A = x_m .* (y_m + 2*L_3 + 2*h_2); % Die Area [m^2] zz = y_m + 2*L_3 + 2*h_2; % Die length [m] m = rho * x_m .* y_m * b_m; % Mass of system (ignore beams) [kg] s_res = (1e7)^2 ./ y_m; % Motion resolution [m] % System Spring Constant [N/m] k = 24*E*I * (6*L_1.*L_2.^2 + L_b*L_2.^2 + 4*L_b*L_1.*L_2 + 24*L_1.*L_2.*L_3 + ... 12*L_b*L_1.*L_3 + 4*L_b*L_2.*L_3) ./ LOCALden(L_1,L_2,L_3,L_b); a_res = (k ./ m) .* s_res; s_stress = (m ./ k) .* a_stress; % Minimum resolvable acceleration [m/s^2] % Displacement at maximum survival accel [m]
% Moment at max survival accel [N*m] M = 6*s_stress*E*I .* (6*L_b*L_1.*L_3.^2 + 12*L_1.*L_2.*L_3.^2 + ... 2*L_b*L_2.*L_3.^2 + 6*L_1.*L_2.^2.*L_3 + L_b*L_2.^2.*L_3 + ... 4*L_b*L_1.*L_2.*L_3 + L_b*L_1.^2.*L_2) ./ LOCALden(L_1,L_2,L_3,L_b); % Forces at max survival accel [N]
C:\Documents and Settings\Scott Moura\My Documen...\enginprop.m April 5, 2006
F_x = (k .* s_stress) / 4; F_y = 18*s_stress*E*I ./ L_2 .* (L_b*L_1.^2.*L_2 + 2*L_b*L_1.*L_3.^2  ... 2*L_b*L_1.^2.*L_3 + 6*L_1.*L_2.*L_3.^2 + L_b*L_2.*L_3.^2) ./ ... LOCALden(L_1,L_2,L_3,L_b); % Stresses at max survival accel [N/m^2] sigma_x = F_x / (h_1*b_1); sigma_y = F_y / (h_1*b_1); sigma_bend = M * (b_1/2) / I; sigma = norm([sigma_x sigma_y sigma_bend]); % Find maximum survivable acceleration (assume only bending stress) M_surv = sigma_max * I / (b_1/2); s_surv = M_surv / (6*E*I .* (6*L_b*L_1.*L_3.^2 + 12*L_1.*L_2.*L_3.^2 + ... 2*L_b*L_2.*L_3.^2 + 6*L_1.*L_2.^2.*L_3 + L_b*L_2.^2.*L_3 + ... 4*L_b*L_1.*L_2.*L_3 + L_b*L_1.^2.*L_2) ./ LOCALden(L_1,L_2,L_3,L_b)); a_surv = (k ./ m) .* s_surv; % Maximum survivable acceleration % Find displacements at 0.005g and 2000g s_sm = m/k * 0.005*g; s_lg = m/k * 2000*g; % Find stress at displacement of 20 um s_20 = 20e6; a_20 = (k ./ m) .* s_20; M_20 = 6*s_20*E*I .* (6*L_b*L_1.*L_3.^2 + 12*L_1.*L_2.*L_3.^2 + ... 2*L_b*L_2.*L_3.^2 + 6*L_1.*L_2.^2.*L_3 + L_b*L_2.^2.*L_3 + ... 4*L_b*L_1.*L_2.*L_3 + L_b*L_1.^2.*L_2) ./ LOCALden(L_1,L_2,L_3,L_b); F_x_20 = (k .* s_20) / 4; F_y_20 = 18*s_20*E*I ./ L_2 .* (L_b*L_1.^2.*L_2 + 2*L_b*L_1.*L_3.^2  ... 2*L_b*L_1.^2.*L_3 + 6*L_1.*L_2.*L_3.^2 + L_b*L_2.*L_3.^2) ./ ... LOCALden(L_1,L_2,L_3,L_b); sigma_x_20 = F_x_20 / (h_1*b_1); sigma_y_20 = F_y_20 / (h_1*b_1); sigma_bend_20 = M_20 * (b_1/2) / I; sigma_20 = norm([sigma_x_20 sigma_y_20 sigma_bend_20]); % Find stress at 0.005g M_sm = 6*s_sm*E*I .* (6*L_b*L_1.*L_3.^2 + 12*L_1.*L_2.*L_3.^2 + ... 2*L_b*L_2.*L_3.^2 + 6*L_1.*L_2.^2.*L_3 + L_b*L_2.^2.*L_3 + ... 4*L_b*L_1.*L_2.*L_3 + L_b*L_1.^2.*L_2) ./ LOCALden(L_1,L_2,L_3,L_b); F_x_sm = (k .* s_sm) / 4; F_y_sm = 18*s_sm*E*I ./ L_2 .* (L_b*L_1.^2.*L_2 + 2*L_b*L_1.*L_3.^2  ... 2*L_b*L_1.^2.*L_3 + 6*L_1.*L_2.*L_3.^2 + L_b*L_2.*L_3.^2) ./ ... LOCALden(L_1,L_2,L_3,L_b); sigma_x_sm = F_x_sm / (h_1*b_1); sigma_y_sm = F_y_sm / (h_1*b_1); sigma_bend_sm = M_sm * (b_1/2) / I; sigma_sm = norm([sigma_x_sm sigma_y_sm sigma_bend_sm]);
Page 2 5:29:58 AM
%% Display Engineering Properties disp(['L_1: ' num2str(L_1*1e6) ' um']) disp(['L_2: ' num2str(L_2*1e6) ' um']) disp(['L_3: ' num2str(L_3*1e6) ' um']) disp(['y_m: ' num2str(y_m*1e6) ' um']) disp(['Die Area: ' num2str(A*1e12) ' um^2']) disp(['Die Dimensions: ' num2str(x_m*1e6) ' um x ' num2str(zz*1e6) ' um'])
C:\Documents and Settings\Scott Moura\My Documen...\enginprop.m April 5, 2006
disp(['Proof Mass: ' num2str(m*1e9) ' ug']) disp(['Spring Constant: ' num2str(k) ' N/m']) disp(['Min Accel Resolution: ' num2str(a_res/g) 'g']) disp(['Motion Resolution: ' num2str(s_res*1e6) ' um']) disp(['Disp @ 0.005g: ' num2str(s_sm*1e6) ' um']) disp(['Stress @ 0.005g: ' num2str(sigma_sm/1e9) ' GPa']) disp(['Maximum Stress: ' num2str(sigma/1e9) ' GPA']) disp(['Disp @ 2000g: ' num2str(s_lg*1e6) ' um']) disp(['Max Survival Accel: ' num2str(a_surv/g) 'g']) disp(['Disp @ ' num2str(a_surv/g) 'g: ' num2str(s_surv*1e6) ' um']) disp(['Accel. @ 20 um disp: ' num2str(a_20/g) 'g']) disp(['Stress @ 20 um disp: ' num2str(sigma_20/1e9) ' GPa'])
Page 3 5:29:58 AM
C:\Documents and Settings\Scott Moura\My Document...\LOCALden.m April 5, 2006
Page 1 5:30:16 AM
function output = LOCALden(L_1,L_2,L_3,L_b) output = (3*L_1.^4.*L_2.^2 + 12*L_1.*L_2.^2.*L_3.^3 + ... 2*L_b*L_1.^3.*L_2.^2 + 2*L_b*L_1.^4.*L_2 + 12*L_1.^4.*L_2.*L_3 + ... 12*L_1.*L_2.*L_3.^4 + 6*L_b.*L_1.^4.*L_3 + 6*L_b*L_1.*L_3.^4 + ... 2*L_b*L_2.^2.*L_3.^3 + 6*L_b.*L_2.*L_1.^2.*L_3.^2 + 8*L_b*L_1.*L_2.*L_3.^3 + ... 2*L_b*L_2.*L_3.^4 + 8*L_b*L_1.^3.*L_2.*L_3);