Integrated topology and shape optimization ... - Semantic Scholar
Advances in Engineering Software 39 (2008) 1–14 www.elsevier.com/locate/advengsoft
Integrated topology and shape optimization software for compliant MEMS mechanism design Gang-Won Jang a, Kyung Joo Kim b, Yoon Young Kim b
b,*
a School of Mechanical Engineering, Kunsan National University, Kunsan, Chonbuk 573-701, Republic of Korea School of Mechanical and Aerospace Engineering and National Creative Research Initiatives Center for Multiscale Design, Seoul National University, Seoul 151-742, Republic of Korea
Received 31 March 2006; received in revised form 22 November 2006; accepted 4 December 2006 Available online 20 February 2007
Abstract The objective of this work is to develop a specialized integrated design tool for MEMS compliant mechanisms from topology to shape optimization. We use the topology optimization technology to obtain an optimal mechanism layout achieving the best overall system performance. This design process is followed by shape optimization to improve stress-based local performance around hinges of the mechanism. Numerical instabilities in topology optimization such as checkerboard layout are overcome by use of nonconforming finite elements. To avoid the cumbersome remeshing process required to handle boundary shape changes, we employ the remesh-free shape optimization method based on the wavelet-Galerkin analysis. A few examples are considered to illustrate the integrated design procedure from topology to shape optimization and to demonstrate the effectiveness of the present approach. 2007 Elsevier Ltd. All rights reserved. Keywords: Compliant mechanism; Topology optimization; Nonconforming element; Wavelet-Galerkin method; Shape optimization; Integrated design; MEMS
1. Introduction Recently, the topology optimization of the compliant mechanism has been successfully applied to the design of the MEMS actuators [1,2]. Unlike the size or shape optimization, the topology optimization does not require a baseline design but uses a simplified design domain, e.g., a rectangle for a two-dimensional case. The design domain is discretized by mesh of the finite element method, and the density of each element in the mesh becomes a design variable. After optimization, the resulting topology is retrieved by connecting elements with high densities. Creative micro-scale design, not based on existing macro-scale counterpart, can be obtained through topology optimization. The resulting image of the topology optimization shows, however, rough boundaries. Furthermore, relatively high *
Corresponding author. Tel.: +82 2 880 7154; fax: +82 2 883 1513. E-mail address: [email protected] (Y.Y. Kim).
0965-9978/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2006.12.003
stresses are observed on specific parts of the boundaries, especially around the hinges, even after executing the boundary smoothing procedure. Moreover, careless boundary smoothing causes the extreme drop of the system performance. The control of local physical quantities such as stress or magnetic flux density in the stage of topology optimization is difficult. Duysinx and Bensøe [3] solved the compliance minimization problem with stress constraints. They suggested e-constraint relaxation of the stress constraints to deal with singularity and introduced as many design constraints as the number of elements so that the numerical expense of the optimization was extremely high. Thus, the final design obtained from the image of topology optimization needs an additional optimization procedure considering local design constraints, which has mainly been conducted based on experience of the designers. Such disadvantage of topology optimization can be overcome to some extent by applying shape optimization to the result of topology optimization. While topology
2
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
optimization deals with global system performances such as compliance or natural frequencies, shape optimization can control local performances by locally adjusting the boundaries of the structure. Olhoff et al. [4] proposed a design method to integrate topology and shape optimization. Tang and Chang [5] approximated the boundaries of the resulting image after topology optimization by using B-spline curves and performed shape optimization by treating control points of the B-spline curves as design parameters. In this paper, using the approach of Tang and Chang, we set the stress of a topologically optimized actuator as the design constraint of shape optimization to prevent the fatigue failure of the structure. During the shape optimization procedure, adopting the standard finite element method as a solver often leads to breakdown of the optimization by excessive mesh distortion [6,7]. Therefore, remeshing and setting reasonable criteria for remeshing are essential to calculate accurate local stress in the finite element method. For the problems undergoing large shape changes or structures with complex boundaries, such remeshing is not so simple task that sometimes it requires the manual intervention of the designer. To resolve this problem, some remesh-free approaches utilizing meshfree [8] or fixed-grid-based method using finite elements [9] have been recently proposed. The fixed-grid-based method has generally low accuracy of stress on a boundary due to its poor boundary approximation. Accurate stress prediction requires very high grid density, which leads to an extraordinarily high computational cost when the fixedgrid-based method is employed. The present software adopts the adaptive multiscale wavelet-Galerkin method by Jang et al. [10] to deal with the mesh distortion problem. In their wavelet-Galerkin method, the fictitious domain approach as in the fixedgrid-based method [9] is employed to avoid remeshing, and the multiscale wavelet bases facilitate the adaptive analysis to reduce the numerical cost for the system analysis. Since the wavelet coefficients represent the difference between solutions of different resolution levels, they can be used as error indicators for the adaptive analysis. To accurately approximate a boundary, they suggested the approximation with piecewise oblique lines connecting intersection points between fixed mesh and the original boundaries. Thus, three main procedures are necessary for the integrated design of MEMS actuators; topology optimization, boundary approximation with B-spline curves, and shape optimization. For topology optimization, the nonconforming finite elements [11] are used to minimize numerical instability such as checkerboard pattern [12]. Then, we smooth the boundaries of the topologically optimized structure and obtain the parametric model of the boundaries by using least square fitting with B-spline curves. The stresses around hinges are controlled through shape optimization. While the shape optimization proceeds, the multiscale adaptive wavelet-Galerkin solver requires the intersection points between fixed-mesh and original boundaries, which are provided using the commercial library [13].
The efficiency of the software will be shown by examining designs of a micro force inverter and a micro gripper. 2. Integrated design concept for the design of MEMS actuators The whole integrated design procedure for the present software is illustrated in Fig. 1. First, the overall system performance and the global constraint are to be defined. In the paper, we consider the deformed displacement at the output point of an actuator as the overall system performance, and the maximization of the output displacement will be treated as a design objective. In topology optimization, the total mass of an actuator is usually constrained to avoid distorted topology by unnecessary mass. The topologically optimized image has zigzag-shaped boundaries. High-order polynomials or B-spline curves (including Bezier curves) are mainly used to approximate the rough boundaries. Because the approximation with high-order polynomials may fluctuate in many cases, Bspline curves are used to represent boundaries in the research, and the control points of the B-spline curves are treated as design variables for shape optimization. For shape optimization, local stress constraints are considered as design constraints, and the objective is to increase or at least maintain the system performance which was obtained in the topology optimization procedure. In this section, we briefly explain the theoretical issues for three main functions of the integrated design scheme: topology optimization, boundary approximation with Bspline curves, and shape optimization. The nonconforming finite element and the multiscale adaptive wavelet-Galerkin method are explained as well because they are important in relieving numerical instabilities for the whole optimization procedures. 2.1. Topology optimization Structural topology optimization can be thought of as a method to determine the optimal spatial material distribution. To obtain an optimized objective function, i.e., minimum compliance, maximum natural frequency or maximum output displacement, we must solve the problem of material redistribution, in other words, how to determine spatial point-wise material or no material condition. By means of the finite element discretization, the topology optimization problem can be generally written as Minimize f ðqÞ ð1aÞ q
subject to KU ¼ F; Ne Z X H ðqÞ ¼ qi dX M 0 6 0; i¼1 T
ð1bÞ ð1cÞ
Xi
q ¼ fqi g ; e 6 qi 6 1; 0 < e 1; i ¼ 1; 2; . . . ; N e ; ð1dÞ Ci ¼ C0 qni ðn P 2Þ:
ð1eÞ
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
3
Problem definition for topology optimization - Set design domain - Set the design objective and constraints
Topology optimization
Data transfer to ANSYS - Examine high stress regions from the ANSYS result
MMA optimizer
Boundary extraction and smoothing Boundary approximation using the B-Spline curve
Problem definition for shape optimization - Select points of stress constraints
Set DOT parameters Wavelet-Galerkin method
Shape optimiaztion
Final optimized structure Fig. 1. Integrated design process.
Eq. (1b) is the discretized system equation, and the total mass of the structure is constrained below M0 by (1c). In (1d), note that the design variable qi, the density of the ith element, has continuous value by adopting the artificial material method (or SIMP method). In Eq. (1e), Ci denotes the matrix form of the elasticity tensor of the ith element, and C0 is that of the solid element (the element with q = 1). In Eq. (1e), the power form of the density is used to converge the density to 1 or e (n = 3 is used in this paper). The finite element method, especially using low order approximation, often suffers from some numerical instabilities during topology optimization [12]. Among them, checkerboard layout of a material distribution is a wellknown problem. For compliance minimization problems, if bilinear finite elements are used, checkerboard patterns are evaluated to be overly stiff and thus employed as optimum layout, which is physically incorrect. Checkerboard patterns are not completely removed even with higher order elements. For compliant mechanisms, an additional problem of one-point hinges [14] also occurs. We use the DSSY nonconforming element [11] to avoid checkerboards or one-point hinges. Fig. 2 illustrates the checkerboard patch composed of 2 strong elements with
þ þ Fig. 2. Four-element checkerboard patch ðC þ 1111 ¼ 1; C 1122 ¼ 0:3; C 1212 ¼ þ 1 0:35 and C ¼ C Þ. ijkl 1000 ijkl
the elasticity tensor C þ ijkl and 2 weak elements with C ijkl . Jang et al. [11] calculated the homogenized elasticity tensor C Hijkl of the checkerboard patch in Fig. 2. They showed that the homogenized stiffness of the checkerboard patch is correctly evaluated with the DSSY nonconforming elements, so the optimized topology by using the nonconforming elements is free from the checkerboard layout. On the other
4
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
hand, the homogenized stiffness of the checkerboard patch consisting of 4-node bilinear conforming elements is evaluated too high, which is physically incorrect. In Jang et al. [11], for example, they showed that C H1111 of 4-node bilinear elements is 250 times larger than that of DSSY noncon 1 forming elements when C þ ijkl ¼ 1000 C ijkl . The basis functions for the DSSY nonconforming element are given by [11]:
where Pj is the (j + 1)th point on a boundary of the result of topology optimization. The second and rth points are averaged with neighboring 3 points, and the first and last points are not averaged. The B-spline curve with the (k 1)th order polynomial blending functions is expressed as [15] n X xðuÞ ¼ Bi N i;k ðuÞ; ð0 6 u 6 1Þ; ð3Þ
1 1 hðxÞ hðyÞ þ xþ ; 4 2 4hð1Þ 1 1 hðyÞ hðxÞ _ u2 ðx; yÞ ¼ þ y þ ; 4 2 4hð1Þ 1 1 hðxÞ hðyÞ _ u3 ðx; yÞ ¼ x þ ; 4 2 4hð1Þ 1 1 hðyÞ hðxÞ _ u4 ðx; yÞ ¼ y þ ; 4 2 4hð1Þ
where the blending function Ni,k is recursively defined as
_
u1 ðx; yÞ ¼
i¼0
N i;k ðuÞ ¼ N i;1 ðuÞ ¼
ðu ti ÞN i;k1 ðuÞ ðtiþk uÞN iþ1;k1 ðuÞ þ ; tiþk1 ti tiþk tiþ1 1 ti 6 u 6 tiþ1 ; 0 otherwise:
ð4bÞ
In Eq. (3), Bi means the control point, and ti in Eq. (4) is the knot value. The approximated boundaries with B-spline curves for shape optimization are obtained by minimizing r X Pj xðuj Þ2 f ¼ ð5Þ
and 5 hðtÞ ¼ t2 t4 : 3
j¼0
Note that nodes of the DSSY element are located at the midpoints of four edges, not on the vertices. For instance, _ u1 ðx; yÞ takes the value of 1 at (1, 0) and vanishes at (0, 1), (1, 0) and (0, 1). Therefore, the continuity along the element edges is not guaranteed except on the midpoints of the edges. In topology optimization, the overall system performance such as mean compliance, or the displacement at a specific point can be effectively controlled. However, topologically optimized design considering local stress constraint is difficult to obtain. Because more than one stress constraints are required for each element, it results in too many design constraints and thus severely raises the cost of the optimization. Even with the optimized topology considering stress constraints, a post-process procedure might change the location of the maximum stress and its magnitude. In this research, only global performance measure is considered for topology optimization. The stress control to prevent fatigue failure of the structure is considered during the following shape optimization procedure.
where Pj denotes the smoothed boundary point calculated in Eq. (2). The values of uj are determined by Pj Pj1 u0 ¼ 0 and uj ¼ uj1 þ ; j ¼ 1; 2; . . . ; r; d ð6aÞ r X Pj Pj1 : d¼ ð6bÞ j¼1
To minimize Eq. (5), the derivatives of f with respect to n + 1 control points are set to be zero, which leads to the following matrix equation: NT NB ¼ NT P;
After topology optimization, the boundaries of the optimized image are parameterized by using B-spline curves. To reduce boundary wiggles inherited from the finite element mesh, smoothing procedure is performed before boundary approximation [5]. New smoothed boundary points are obtained by applying a simple averaging filter to the boundary nodes: j ¼ 2; . . . ; r 2; ð2Þ
ð7Þ
where B and P denote the vectors of control points and smoothed boundary points, respectively, and the matrix N is defined as 3 2 N 0;k ðu0 Þ N 1;k ðu0 Þ N n;k ðu0 Þ 6 N 0;k ðu1 Þ N 1;k ðu1 Þ N n;k ðu1 Þ 7 7 6 7 N¼6 : ð8Þ . . .. .. 7 6 .. 5 4. . . . N 0;k ður Þ
2.2. Boundary parameterization with B-spline curves
1 Pj ¼ ðPj2 þ Pj1 þ Pj þ Pjþ1 þ Pjþ2 Þ; 5
ð4aÞ
N 1;k ður Þ
N n;k ður Þ
ðrþ1Þðnþ1Þ
Some of the control points obtained from Eq. (7) will be selected as design variables for shape optimization. 2.3. Remesh-free shape optimization using the wavelet-Galerkin method The interest on the shape optimization based on various fixed-grid-type analyses has been recently increased because the fixed-grid analyses offer simple ways to avoid cumbersome remeshing processes which cannot be avoided if the standard finite element approach is employed. Another advantage of the fixed-grid-based shape optimization
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
method is that it requires only the boundary velocity field for design updates while the standard finite-element-based shape optimization method generally requires design velocities for all nodes. Since the analysis nodes are independent of the shape of a given structure in the fixed-grid method, only the boundary velocity field is needed for shape optimization. Most fixed-grid-based analyses adopt the fictitious domain method to deal with domains with general shapes. In Fig. 3, the example of the fictitious domain is illustrated. The real domain (or package domain) x is embedded into the simplified fictitious domain X. Because the fixed mesh of the fictitious domain is independent of the shape of the package domain which changes during shape optimization, there is no mesh distortion problem. The most important issue of the fictitious domain approach lies in the boundary approximation. Poorlyapproximated boundaries will lead to inaccurate stress prediction and spoil the stability of optimization process. The original boundary ox consisting of B-spline curves needs to be approximated for stiffness calculation of the fixed elements. In Fig. 4, various approximation methods for the fictitious domain method are illustrated. Compared to the zigzag approximation in Fig. 4b and the area-fractionbased approximation in Fig. 4c, the approximation with piecewise oblique lines by Jang et al. [10] in Fig. 4d gives accurate stress results on the boundaries. The piecewise boundary lines are approximated by connecting the intersection points between the fixed mesh and the original boundary ox. The stiffness of an element in the fixed mesh depends on its location, which is given by 8 < C0 if element x; Ci ¼ cC0 else if element X n x; ð9Þ : Ci else if element is on bounary;
Fig. 4. The approximation of a circle: (a) the original domain of the circle, and its approximated domains using (b) zigzag, (c) area-fraction-based method, and (d) present approach with piecewise oblique lines.
element B has C0. In Eq. (9), Ci is piecewise-constant over the element. In Fig. 5, the element C in Fig. 3b is illustrated with its original boundary ox (dotted line), approximated boundary oxA (thick solid line), and local element coordinates n, g. The element stiffness matrix for the element is given by Kei ¼ ¼
Z
Z
1
1
Z Z
1
1
BT ðn; gÞCi ðn; gÞBðn; gÞjJjdndg;
BT ðn; gÞC0 Bðn;gÞjJjdndg; xAi 1
Z
1
where c = 0.001 is used in this paper. In Fig. 3b, for example, the element A has the elasticity tensor of cC0 and the
a
5
1
ð10aÞ ð10bÞ
BT ðnðr; sÞ;gðr; sÞÞC0 Bðnðr; sÞ; gðr;sÞÞjJjj b Jðr; sÞjdrds;
1
ð10cÞ
b
A
B
C
D
ω Ω
Fig. 3. (a) The simplified fictitious domain X enclosing the package domain x, and (b) its fixed grid mesh.
6
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14 0
1
2
3
4
5
6
8
7
0
1
0,0
1,0
2,0
2,1
0,1
1,1
2,2
2,3
0,2
Fig. 5. Zoomed view of the element C in Fig. 4b.
0
where B and jJj mean the strain–displacement matrix and the determinant of the Jacobian matrix, respectively, and xAi denotes the region of the approximated x within the element. In Eq. (10c), we introduce additional coordinates: the coordinates r, s which map xAi or Xi n xAi onto the normalized rectangle. Note that the other Jacobian matrix b J is introduced to concern the mapping between the coordinates r, s and the local element coordinates n, g: (see [10] for more details) " # og on or or b¼ J : ð11Þ og os
on os
If Xi n xAi is a triangle as in the element D of Fig. 3b, the stiffness matrix for the element is obtained by subtracting the stiffness of Xi n xAi from the stiffness of an inside element, i.e., Kei
¼
Z
1
Z
1
Z
1
BT ðn; gÞC0 Bðn; gÞjJj dn dg
1
Fig. 6. Equivalent vector space spanned by (a) the single scale basis functions, and (b) multiscale hat interpolation wavelets.
b J , respectively, the system equation for the multiscale and F wavelet-Galerkin method becomes KJ WJ ¼ FJ ;
ð13aÞ
where the subscript J denotes the highest uniform-grid density or resolution level, and WJ is the wavelet coefficients vector for nodal displacements. The relation between bJ; F b J Þ and (KJ, FJ) can be stated by using a transformaðK tion matrix TJ as (see [16]) b J TJ KJ ¼ TTJ K
and
bJ : F J ¼ TTJ F
ð13bÞ
Fig. 7 illustrates the flowchart for the multiscale adaptive analysis where the analysis level increases from a certain resolution level j0 to the highest resolution level J. The adaptive strategy is implemented by utilizing the difference-checking nature of wavelets for error estimation [16].
1
BT ðn; gÞC0 Bðn; gÞjJj dn dg:
ð12Þ
Begin analysis
Xi nxAi
The above element stiffness calculations are necessary only for those elements lying on the boundaries. For elements inside or outside the boundaries, their stiffness matrices are all the same, the first term of the right hand side of Eq. (12) or its c times. Thus, the advantages of fixedgrid-based methods, the fast mesh and fast formulation of the stiffness matrix, are still preserved at a small additional expense required to calculate intersection points and element stiffness matrices on the boundaries Fig. 6. When an adaptive analysis is incorporated into the fixed-grid analysis strategy, the shape optimization can become more efficient. To this end, we will use the adaptive multiscale wavelet-Galerkin method employing hat interpolation wavelets. The details on the method may be found in [16], so we only emphasize that the system matrix (in Eq. (13a)) of the multiscale wavelet-Galerkin method can be easily transformed form the single-scale system matrix, i.e., the standard system matrix resulting from the bilinear finite element formulation. If the single-scale stiffness bJ matrix and the single-scale load vector are denoted by K
Resolution j0 analysis Calculate intersecting points
Build single-scale system equations
^ ^ ^ U K J J = FJ
j=j0
Adaptive scheme j=j+1 Resolution j analysis
No
Wavelet transform
j=J? Multiscale equations
K JWJ = FJ
Yes End analysis
Fig. 7. Flowchart for multiscale adaptive analysis.
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
7
3. Software details and design examples In order to implement the idea of the integrated design for compliant MEMS’s into a user-friendly software, additional considerations besides the analysis modules for topology and shape optimization are necessary as follows. (1) To show the progress of the topology optimization, element densities are expressed as a black–white image. (2) The boundary points of the topologically optimized image are extracted, smoothed, and parameterized with B-spline curves. The commercial software NLib [13] is employed to manipulate the B-spline curves. (3) As a shape optimizer, the DOT library [17] is loaded, and a user interface is graphically expressed to set up the parameters for the optimizer. (4) To determine the design variables and stress constraints for shape optimization, the control points and boundary points estimated to have high stress are graphically selected with a mouse. (5) To calculate the element stiffness matrices of the wavelet-Galerkin analysis, the intersection points between B-spline curves and fixed mesh are calculated using the NLib and are handed over to the solver. (6) The geometry of a structure can be converted to the ANSYS file, and the performance of the optimized structure is confirmed though the commercial analyzer. The overall procedure of the software, including the above processes, follows the integrated design scheme in Fig. 1. We explain the details of the software step-by-step through a couple of examples.
Fig. 8. The design domain for topology optimization of a micro force inverter. (64 · 32 design variables for the upper half of the structure.)
Table 1 The parameters for topology optimization of the micro force inverters defined in Figs. 8 and 17 Young’s modulus Poisson’s ratio Fin kin kout Mesh Number of design variables Thickness M0
1.7 · 105 lN/lm2 0.26 5000 lN 200 lN/lm 2 lN/lm 64 · 32, DSSY 2048 6 lm 30%
3.1. Design of micro force inverters Using the integrated design software, the two micro force inverters are designed. The micro force inverter exerts a force on a target in the opposite direction to an input force while the energy loss by elastic deformation of compliant mechanism is minimized. The topology of a force inverter having maximum performance is determined through topology optimization, and then the local stresses of the structure around hinges are controlled below the allowed level during shape optimization. The problem definition for topology optimization of the micro force inverter is given as [18] Maximize subject to
uout ðqÞ H ðqÞ ¼
designed. The MMA optimizer [19], which is efficient in dealing with many design variables, is used for this problem. The optimized topology of the force inverter is illustrated in Fig. 9. A total of 34 analyses are performed during the optimization. Fig. 10 shows the example of boundary smoothing process for region A in Fig. 9. The dark points denote the extracted boundary from the result of topology optimization, and the bright points mean the smoothed boundary points obtained by using Eq. (2). Before approximating
ð14aÞ Ne X i¼1 T
Z
qi dX M 0 6 0;
ð14bÞ
Xi
q ¼ fqi g ; e 6 qi 6 1; i ¼ 1; 2; . . . ; N e :
0 < e 1;
The objective is to maximize the displacement at the output point, and the total mass of the structure is constrained. Fig. 8 shows the design domain for topology optimization where the dimensions, the input point A, and the output point B are denoted. The parameters for the problem are listed in Table 1. Considering the symmetry of the structure, only upper half of the domain is analyzed and
Fig. 9. The topology optimization result by the proposed software.
8
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
Fig. 10. The smoothed boundary points of the region A in Fig. 10. (Dark circles: original boundary points, white circles: smoothed boundary points.)
the boundaries with use of B-spline curves, the user should determine the order of the curves, number of control points for each curve, and start and end tangential angles of the curves via the dialog box in Fig. 11a. Excessive number of control points often leads to an oscillatory boundary while too small number of control points produces a poorly approximated curve quite different from the original boundary. Fig. 11b shows approximated boundaries of the force inverter by using 11 B-spline curves with second or third orders. The straight lines in the figure denote control polygons of B-spline curves. Note that the hinges are
Fig. 11. The boundary approximation with B-spline curves: (a) the dialog box for B-spline curve information, and (b) the approximated boundaries of the force inverter.
approximated by separate curves (curve number 0, 2, and 4 in Fig. 11b) to increase the degree of design changes. The control points around the hinges are mainly selected as the design parameters for shape optimization. Fig. 12 shows the design parameters around three hinge regions. The proposed integrated design software uses the DOT [17] as mathematical optimization algorithms for shape optimization. To set up various parameters for the optimizer, the graphical user interface in Fig. 13 is provided. Through the dialog box in Fig. 13, the use can determine file names for various purposes, number of design variables, number of stress constraints, upper and lower bound limits for design variables, initial value for design variables, optimization strategy, and other detail parameters for the optimization strategy. Other information required by DOT is • Design variables: Using a mouse, the user can select or unselect them on the screen among the control points of B-spline curves. • Stress constraints: The boundary points expected to have high stress are examined. The parametric value u in Eq. (3) for each point under stress constraint is stored to keep track of the point during the optimization. • Sensitivity: The sensitivities of the objective and constraints are obtained by simply using the finite difference method. The offset value for the finite difference method is determined through the dialog box in Fig. 13. The default for the offset value is set to be equal to the length of the fixed elements. The user can determine the upper and lower bounds for design variables uniformly by using the dialog box in Fig. 13, or can impose different bounds for each variable by picking it on the screen. The design variables lying on domain edges have only one-directional degree of freedom. For example, only horizontal movement of the control point lying on the bottom edge in Fig. 12 can be selected as a design variable. To constrain the stress of the micro force inverter during shape optimization, the boundary regions with high stress need to be located. The proposed integrated design software provides the data exporting function to ANSYS that the geometry of a structure can be converted to the ANSYS input file. After examining the analysis result of ANSYS, the points of high stress are selected, and their von Mises stresses are set as design constraints. Fig. 14a shows the stress contour plot by ANSYS. Based on the result of ANSYS, the points under stress constraints are selected as in Fig. 14b. Note that the points under stress constraints are concentrated around elastic hinge regions. Other information for the shape optimization of the micro force inverter is listed in Table 2. During shape optimization, the wavelet-Galerkin solver requires the intersection points between B-spline boundaries and the fixed mesh. Fig. 15 shows the fixed mesh generation by the software. The intersection points between
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
9
Fig. 12. The design parameters (blue circles). (For interpretation of colour representation in this figure legend the reader is referred to the web version of this article.)
Fig. 13. The user interface for the DOT parameters.
the curves and the fixed mesh are automatically calculated and are used to construct stiffness matrices of the boundary elements according to Eqs. (10) and (12). The optimized result with its deformed shape is illustrated in Fig. 16. Noticeable change can be seen around the central hinge region. The system performance improves 81.4% compared to the structure before shape optimization. Thus the hinge region is regarded to have much higher sensitivity to the system performance than other regions of the micro structure. In Fig. 16, the stress distribution around the hinge verified by using ANSYS is also shown. The maximum stress around the hinge is 982.67 lN/lm2 which satisfies the design constraint. Table 3 lists the output displacements and the maximum stresses of the optimized result.
Fig. 14. Stress constraint imposition: (a) the stress contour (by ANSYS), and (b) the points under stress constraints.
Because a force inverter experiences large displacement, the performance of the optimized result may have to be checked by nonlinear analysis. Thus, ANSYS was used to carry out geometrically-nonlinear analysis of the force converter shown in Fig. 16. The comparison with the nonlinear and linear analyses showed that the output displacement and the von Mises stress by the nonlinear analysis decrease
10
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
Table 2 Parameters used to obtain the results in Figs. 16 and 19 by shape optimization Objective
Output displacement uout
Number of design variables Number of stress constraints Maximum allowable stress Side constraints Optimizer Mesh
42 23 1000 lN/lm2 6 lm 6 qi 6 6 lm SQP 256 · 128, wavelet-Galerkin method
Fig. 15. The fixed mesh of the wavelet-Galerkin analysis (a) for the whole domain and (b) around the hinge A.
by 3.4% and 2.8% , respectively in comparison with the calculations by the linear analysis. So, for the present problem in consideration, the proposed linear wavelet-Galerkin method appears to predict the actual structural behavior quite accurately. Next, a micro force inverter having two output points is designed. Fig. 17 shows the design domain for topology optimization procedure. The design objective is to maximize the output displacements at B and C in the figure. Other conditions are given in Tables 1 and 2. Fig. 18 shows the optimized topological layout and its boundary approximation by B-spline curves. Before extracting curves from the optimized result in Fig. 18a, manual modifications of the resulting topology are necessary. Typical results by topology optimization may include several holes, some of which can be removed before applying shape optimization because their contributions to mass constraints and structural performance may be insignificant. For instance, a few small holes in Fig. 18a are eliminated during the modification and the modified layout is illustrated in Fig. 18b. Nevertheless, the system performance should be confirmed by additional analysis on the modified topology. To facilitate the modification process, the proposed integrated design software has an editing function. For shape optimization, 7 control points around hinges are selected for design parameters, and 48 points around high stress regions, mostly around hinges, are selected for design constraints. The final shape of the force inverter and its deformed shape are illustrated in Fig. 19. The input force is inverted at the output points via three elastic hinges. The stress distributions of the hinges are also illustrated in Fig. 19. Now we have presented two design examples, le us explain the date structure and communication flow of the proposed software: see Fig. 20. The input data for topology optimization is provided in a text file that contains design domain information, boundary conditions and material
Fig. 16. The deformation shape of the optimized micro force inverter and the von Mises stress distribution around the hinge (verified by using ANSYS).
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
11
Table 3 Output displacements and stresses after shape optimization Output displacement [lm]
Force inverter in Fig. 6 Force inverter in Fig. 19 Gripper in Fig. 21
Initial value
Optimized value
Performance increase [%]
1.893 0.648 0.744
3.435 2.563 2.260
81.4 295.5 203.8
Maximum von Mises stress of the optimized structure [lN/lm2] by wavelet-Galerkin method (by ANSYS) 953.50 (982.67) 974.16 (975.63) 983.29 (1021)
Fig. 17. The design domain for topology optimization of a micro force inverter with two output points (64 · 32 design variables for the upper half of the structure).
properties. The optimized densities q are transferred to the boundary handling program as a form of a data array, from which the boundary information data including smoothed boundary points and B-spline curves are
Fig. 18. (a) The topology optimization result of the force inverter, and (b) its boundary approximation.
Fig. 19. The optimized shape of the micro force inverter with two output points and von Mises stress distribution around hinges (verified by using ANSYS).
12
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14 Topology optimization module Problem configuraion file (text file) 1. Design domain dimensions 2. Spring stiffness 3. Input and output points 4. Boundary conditions 5. Material properties
Scalar data: design domain size Array data: rho.dat (design variables, text file)
Boundary handling module (Extraction, smoothing and parameterization) Wavelet-Galerkin solver
Scalar data: number of boundaries Array data: 1. bnd.dat (boundary points, text file) 2. spline.dat (B-spline curve information, text file)
Array data: int.dat (text file) intersection points
Scalar data: objective value Array data: 1. cnst.dat (constraint values, text file) 2. objgrd.dat, cnstgrd.dat (text files) sensititvities for objective and constraints
DOT optimizer (dot.dll)
Shape optimization module
Array data: Updated design parameters
Fig. 20. Data structure and communication flow chart of the integrated design software developed in this study.
extracted. The shape optimization program manages three data groups: one for boundary manipulation, another for optimization algorithm, and the other for data communication with the analyzer. 3.2. Design of a micro gripper A micro gripper grasping and moving a micro-scale object is designed in this example. The design domain for topology optimization and its dimensions are shown in
Fig. 21. The design domain for topology optimization of a micro gripper (kin = 400 lN/lm, kout = 20 lN/lm and Fin = 4000 lN).
Fig. 21. The objective is to maximize the displacement at the output point B. To increase the gripping force of the structure, the stiffness of the output spring 10 times large than that used in the previous examples is chosen. 10 times larger stiffness for the output spring than in the previous examples is used. The result by the topology optimization is shown in Fig. 22a. After boundary extraction and smoothing, the approximated boundaries by B-spline curves are illustrated in Fig. 22b. A total of 8 B-spline curves are employed, and the orders of the curves vary from linear to cubic splines. Note that 6 small holes in Fig. 22a are removed through the topology modification procedure. Next, the regions with high stresses are located by exporting the boundary data to ANSYS and the points subject to stress constraints are selected. The control points around hinge regions are selected as design variables, which are denoted in Fig. 22b. A total of 11 design parameters and 54 stress constraints are employed for shape optimization. Fig. 23 shows the final design and its deformation shape. Table 3 shows that due to high stiffness of the output spring, the output displacement of the gripper is relatively small. Fig. 24 shows the distribution of the adaptivelyadded wavelet bases of the wavelet-Galerkin method. Note that the support sizes of the multiscale wavelets are different according to their resolution levels. The dots in the figure denote the centers of the wavelets. In Fig. 24, it is observed that while few wavelets are located on the domain Xnx having weak material, most of wavelets exist within
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
13
Fig. 22. (a) The topology optimization result, and (b) boundary approximation (blue circles denote design parameters). (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 23. The optimized shape of the micro gripper and the von Mises stress distribution around the hinge (verified by using ANSYS).
the package domain, and high stress regions require dense distribution of wavelets. 4. Conclusions The integrated design software for complaint MEMS actuators was developed. Using topology optimization, a rough image for an actuator having high overall system performance was obtained. Nonconforming finite elements were used to overcome the numerical instability of checkerboards or one-point hinges. After boundary extraction, smoothing, and approximation with the B-spline curves, remesh-free shape optimization based on the multiscale wavelet-Galerkin method was performed to control local stresses around hinges. Through design examples, the validity and effectiveness of the proposed integrated design software were verified.
Acknowledgements
Fig. 24. The distribution of adaptively-added wavelets (upper: x-directional wavelets, lower: y-directional wavelets).
This research was supported by the National Creative Research Initiatives Program (Korea Science and Technology Foundation Grant No. 2006-033) contracted through
14
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
the Institute of Advanced Machinery and Design at Seoul National University. References [1] Sigmund O. Design of multiphysics actuators using topology optimization – part I: one material structures. Comput Methods Appl Mech Engng 2001;190:6577–604. [2] Chen BC, Silva CN, Kikuchi N. Advances in computational design and optimization with application to MEMS. Int J Numer Meth Engng 2001;52:23–62. [3] Duysinx P, Bendsøe MP. Topology optimization of continuum structures with local stress constraints. Int J Numer Meth Engng 1998;43:1453–78. [4] Olhoff N, Bendsøe MP, Rasmussen J. On CAD-integrated structural topology and design optimization. Comput Meth Appl Mech Eng 1992;89:259–79. [5] Tang PS, Chang KH. Integration of topology and shape optimization for design of structural components. Struct Multidiscip Opt 2000; 22:65–82. [6] Bennett JA, Botkin ME. Structural shape optimization with geometric description and adaptive mesh refinement. AIAA J 1985; 23(3):458–64. [7] Yao TM, Choi KK. 3-D shape optimal design and automatic finite element regridding. Int J Numer Meth Engng 1989;28:369–84. [8] Kim NH, Choi KK, Botkin ME. Numerical method for shape optimization using meshfree method. Struct Multidis Opt 2003; 24:418–29.
[9] Garcia MJ, Steven GP. Fixed grid finite elements in elasticity problems. Eng Comput 1998;16(2):154–64. [10] Jang GW, Kim YY, Choi KK. Remesh-free shape optimization using the wavelet-Galerkin method. Int J Sol Struct 2004;41:6465–83. [11] Jang GW, Jeong JH, Kim YY, Sheen D, Park C, Kim MN. Checkerboard-free topology optimization using nonconforming finite elements. Int J Numer Meth Engng 2003;57:1717–35. [12] Sigmund O, Petersson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, meshindependencies and local minima. Struct Opt 1998;16:68–75. [13] Solid Modeling SolutionsTM. Available from: http://www.smlib.com/ nlib.html. [14] Yoon GH, Kim YY, Bendsøe MP, Sigmund O. Hinge-free optimization with embedded translation-invariant differentiable wavelet shrinkage. Struct Multidis Opt 2004;27:139–50. [15] Lee K. Principles of CAD/CAM/CAE systems. Addison Wesley; 1999. [16] Jang GW, Kim JE, Kim YY. Multiscale Galerkin method using interpolation wavelets for two-dimensional elliptic problems in general domains. Int J Numer Meth Engng 2004;59:225–53. [17] Vanderplaats Research & Development, Inc. DOT Users Manual, 1999. [18] Sigmund O. On the design of compliant mechanisms using topology optimization. Mech Struct Mach 1997;25(4):495–526. [19] Svanberg K. The method of moving asymptotes – a new method for structural optimization. Int J Numer Meth Engng 1987;24:359–73.
Integrated topology and shape optimization software for compliant MEMS mechanism design Gang-Won Jang a, Kyung Joo Kim b, Yoon Young Kim b
b,*
a School of Mechanical Engineering, Kunsan National University, Kunsan, Chonbuk 573-701, Republic of Korea School of Mechanical and Aerospace Engineering and National Creative Research Initiatives Center for Multiscale Design, Seoul National University, Seoul 151-742, Republic of Korea
Received 31 March 2006; received in revised form 22 November 2006; accepted 4 December 2006 Available online 20 February 2007
Abstract The objective of this work is to develop a specialized integrated design tool for MEMS compliant mechanisms from topology to shape optimization. We use the topology optimization technology to obtain an optimal mechanism layout achieving the best overall system performance. This design process is followed by shape optimization to improve stress-based local performance around hinges of the mechanism. Numerical instabilities in topology optimization such as checkerboard layout are overcome by use of nonconforming finite elements. To avoid the cumbersome remeshing process required to handle boundary shape changes, we employ the remesh-free shape optimization method based on the wavelet-Galerkin analysis. A few examples are considered to illustrate the integrated design procedure from topology to shape optimization and to demonstrate the effectiveness of the present approach. 2007 Elsevier Ltd. All rights reserved. Keywords: Compliant mechanism; Topology optimization; Nonconforming element; Wavelet-Galerkin method; Shape optimization; Integrated design; MEMS
1. Introduction Recently, the topology optimization of the compliant mechanism has been successfully applied to the design of the MEMS actuators [1,2]. Unlike the size or shape optimization, the topology optimization does not require a baseline design but uses a simplified design domain, e.g., a rectangle for a two-dimensional case. The design domain is discretized by mesh of the finite element method, and the density of each element in the mesh becomes a design variable. After optimization, the resulting topology is retrieved by connecting elements with high densities. Creative micro-scale design, not based on existing macro-scale counterpart, can be obtained through topology optimization. The resulting image of the topology optimization shows, however, rough boundaries. Furthermore, relatively high *
Corresponding author. Tel.: +82 2 880 7154; fax: +82 2 883 1513. E-mail address: [email protected] (Y.Y. Kim).
0965-9978/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2006.12.003
stresses are observed on specific parts of the boundaries, especially around the hinges, even after executing the boundary smoothing procedure. Moreover, careless boundary smoothing causes the extreme drop of the system performance. The control of local physical quantities such as stress or magnetic flux density in the stage of topology optimization is difficult. Duysinx and Bensøe [3] solved the compliance minimization problem with stress constraints. They suggested e-constraint relaxation of the stress constraints to deal with singularity and introduced as many design constraints as the number of elements so that the numerical expense of the optimization was extremely high. Thus, the final design obtained from the image of topology optimization needs an additional optimization procedure considering local design constraints, which has mainly been conducted based on experience of the designers. Such disadvantage of topology optimization can be overcome to some extent by applying shape optimization to the result of topology optimization. While topology
2
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
optimization deals with global system performances such as compliance or natural frequencies, shape optimization can control local performances by locally adjusting the boundaries of the structure. Olhoff et al. [4] proposed a design method to integrate topology and shape optimization. Tang and Chang [5] approximated the boundaries of the resulting image after topology optimization by using B-spline curves and performed shape optimization by treating control points of the B-spline curves as design parameters. In this paper, using the approach of Tang and Chang, we set the stress of a topologically optimized actuator as the design constraint of shape optimization to prevent the fatigue failure of the structure. During the shape optimization procedure, adopting the standard finite element method as a solver often leads to breakdown of the optimization by excessive mesh distortion [6,7]. Therefore, remeshing and setting reasonable criteria for remeshing are essential to calculate accurate local stress in the finite element method. For the problems undergoing large shape changes or structures with complex boundaries, such remeshing is not so simple task that sometimes it requires the manual intervention of the designer. To resolve this problem, some remesh-free approaches utilizing meshfree [8] or fixed-grid-based method using finite elements [9] have been recently proposed. The fixed-grid-based method has generally low accuracy of stress on a boundary due to its poor boundary approximation. Accurate stress prediction requires very high grid density, which leads to an extraordinarily high computational cost when the fixedgrid-based method is employed. The present software adopts the adaptive multiscale wavelet-Galerkin method by Jang et al. [10] to deal with the mesh distortion problem. In their wavelet-Galerkin method, the fictitious domain approach as in the fixedgrid-based method [9] is employed to avoid remeshing, and the multiscale wavelet bases facilitate the adaptive analysis to reduce the numerical cost for the system analysis. Since the wavelet coefficients represent the difference between solutions of different resolution levels, they can be used as error indicators for the adaptive analysis. To accurately approximate a boundary, they suggested the approximation with piecewise oblique lines connecting intersection points between fixed mesh and the original boundaries. Thus, three main procedures are necessary for the integrated design of MEMS actuators; topology optimization, boundary approximation with B-spline curves, and shape optimization. For topology optimization, the nonconforming finite elements [11] are used to minimize numerical instability such as checkerboard pattern [12]. Then, we smooth the boundaries of the topologically optimized structure and obtain the parametric model of the boundaries by using least square fitting with B-spline curves. The stresses around hinges are controlled through shape optimization. While the shape optimization proceeds, the multiscale adaptive wavelet-Galerkin solver requires the intersection points between fixed-mesh and original boundaries, which are provided using the commercial library [13].
The efficiency of the software will be shown by examining designs of a micro force inverter and a micro gripper. 2. Integrated design concept for the design of MEMS actuators The whole integrated design procedure for the present software is illustrated in Fig. 1. First, the overall system performance and the global constraint are to be defined. In the paper, we consider the deformed displacement at the output point of an actuator as the overall system performance, and the maximization of the output displacement will be treated as a design objective. In topology optimization, the total mass of an actuator is usually constrained to avoid distorted topology by unnecessary mass. The topologically optimized image has zigzag-shaped boundaries. High-order polynomials or B-spline curves (including Bezier curves) are mainly used to approximate the rough boundaries. Because the approximation with high-order polynomials may fluctuate in many cases, Bspline curves are used to represent boundaries in the research, and the control points of the B-spline curves are treated as design variables for shape optimization. For shape optimization, local stress constraints are considered as design constraints, and the objective is to increase or at least maintain the system performance which was obtained in the topology optimization procedure. In this section, we briefly explain the theoretical issues for three main functions of the integrated design scheme: topology optimization, boundary approximation with Bspline curves, and shape optimization. The nonconforming finite element and the multiscale adaptive wavelet-Galerkin method are explained as well because they are important in relieving numerical instabilities for the whole optimization procedures. 2.1. Topology optimization Structural topology optimization can be thought of as a method to determine the optimal spatial material distribution. To obtain an optimized objective function, i.e., minimum compliance, maximum natural frequency or maximum output displacement, we must solve the problem of material redistribution, in other words, how to determine spatial point-wise material or no material condition. By means of the finite element discretization, the topology optimization problem can be generally written as Minimize f ðqÞ ð1aÞ q
subject to KU ¼ F; Ne Z X H ðqÞ ¼ qi dX M 0 6 0; i¼1 T
ð1bÞ ð1cÞ
Xi
q ¼ fqi g ; e 6 qi 6 1; 0 < e 1; i ¼ 1; 2; . . . ; N e ; ð1dÞ Ci ¼ C0 qni ðn P 2Þ:
ð1eÞ
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
3
Problem definition for topology optimization - Set design domain - Set the design objective and constraints
Topology optimization
Data transfer to ANSYS - Examine high stress regions from the ANSYS result
MMA optimizer
Boundary extraction and smoothing Boundary approximation using the B-Spline curve
Problem definition for shape optimization - Select points of stress constraints
Set DOT parameters Wavelet-Galerkin method
Shape optimiaztion
Final optimized structure Fig. 1. Integrated design process.
Eq. (1b) is the discretized system equation, and the total mass of the structure is constrained below M0 by (1c). In (1d), note that the design variable qi, the density of the ith element, has continuous value by adopting the artificial material method (or SIMP method). In Eq. (1e), Ci denotes the matrix form of the elasticity tensor of the ith element, and C0 is that of the solid element (the element with q = 1). In Eq. (1e), the power form of the density is used to converge the density to 1 or e (n = 3 is used in this paper). The finite element method, especially using low order approximation, often suffers from some numerical instabilities during topology optimization [12]. Among them, checkerboard layout of a material distribution is a wellknown problem. For compliance minimization problems, if bilinear finite elements are used, checkerboard patterns are evaluated to be overly stiff and thus employed as optimum layout, which is physically incorrect. Checkerboard patterns are not completely removed even with higher order elements. For compliant mechanisms, an additional problem of one-point hinges [14] also occurs. We use the DSSY nonconforming element [11] to avoid checkerboards or one-point hinges. Fig. 2 illustrates the checkerboard patch composed of 2 strong elements with
þ þ Fig. 2. Four-element checkerboard patch ðC þ 1111 ¼ 1; C 1122 ¼ 0:3; C 1212 ¼ þ 1 0:35 and C ¼ C Þ. ijkl 1000 ijkl
the elasticity tensor C þ ijkl and 2 weak elements with C ijkl . Jang et al. [11] calculated the homogenized elasticity tensor C Hijkl of the checkerboard patch in Fig. 2. They showed that the homogenized stiffness of the checkerboard patch is correctly evaluated with the DSSY nonconforming elements, so the optimized topology by using the nonconforming elements is free from the checkerboard layout. On the other
4
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
hand, the homogenized stiffness of the checkerboard patch consisting of 4-node bilinear conforming elements is evaluated too high, which is physically incorrect. In Jang et al. [11], for example, they showed that C H1111 of 4-node bilinear elements is 250 times larger than that of DSSY noncon 1 forming elements when C þ ijkl ¼ 1000 C ijkl . The basis functions for the DSSY nonconforming element are given by [11]:
where Pj is the (j + 1)th point on a boundary of the result of topology optimization. The second and rth points are averaged with neighboring 3 points, and the first and last points are not averaged. The B-spline curve with the (k 1)th order polynomial blending functions is expressed as [15] n X xðuÞ ¼ Bi N i;k ðuÞ; ð0 6 u 6 1Þ; ð3Þ
1 1 hðxÞ hðyÞ þ xþ ; 4 2 4hð1Þ 1 1 hðyÞ hðxÞ _ u2 ðx; yÞ ¼ þ y þ ; 4 2 4hð1Þ 1 1 hðxÞ hðyÞ _ u3 ðx; yÞ ¼ x þ ; 4 2 4hð1Þ 1 1 hðyÞ hðxÞ _ u4 ðx; yÞ ¼ y þ ; 4 2 4hð1Þ
where the blending function Ni,k is recursively defined as
_
u1 ðx; yÞ ¼
i¼0
N i;k ðuÞ ¼ N i;1 ðuÞ ¼
ðu ti ÞN i;k1 ðuÞ ðtiþk uÞN iþ1;k1 ðuÞ þ ; tiþk1 ti tiþk tiþ1 1 ti 6 u 6 tiþ1 ; 0 otherwise:
ð4bÞ
In Eq. (3), Bi means the control point, and ti in Eq. (4) is the knot value. The approximated boundaries with B-spline curves for shape optimization are obtained by minimizing r X Pj xðuj Þ2 f ¼ ð5Þ
and 5 hðtÞ ¼ t2 t4 : 3
j¼0
Note that nodes of the DSSY element are located at the midpoints of four edges, not on the vertices. For instance, _ u1 ðx; yÞ takes the value of 1 at (1, 0) and vanishes at (0, 1), (1, 0) and (0, 1). Therefore, the continuity along the element edges is not guaranteed except on the midpoints of the edges. In topology optimization, the overall system performance such as mean compliance, or the displacement at a specific point can be effectively controlled. However, topologically optimized design considering local stress constraint is difficult to obtain. Because more than one stress constraints are required for each element, it results in too many design constraints and thus severely raises the cost of the optimization. Even with the optimized topology considering stress constraints, a post-process procedure might change the location of the maximum stress and its magnitude. In this research, only global performance measure is considered for topology optimization. The stress control to prevent fatigue failure of the structure is considered during the following shape optimization procedure.
where Pj denotes the smoothed boundary point calculated in Eq. (2). The values of uj are determined by Pj Pj1 u0 ¼ 0 and uj ¼ uj1 þ ; j ¼ 1; 2; . . . ; r; d ð6aÞ r X Pj Pj1 : d¼ ð6bÞ j¼1
To minimize Eq. (5), the derivatives of f with respect to n + 1 control points are set to be zero, which leads to the following matrix equation: NT NB ¼ NT P;
After topology optimization, the boundaries of the optimized image are parameterized by using B-spline curves. To reduce boundary wiggles inherited from the finite element mesh, smoothing procedure is performed before boundary approximation [5]. New smoothed boundary points are obtained by applying a simple averaging filter to the boundary nodes: j ¼ 2; . . . ; r 2; ð2Þ
ð7Þ
where B and P denote the vectors of control points and smoothed boundary points, respectively, and the matrix N is defined as 3 2 N 0;k ðu0 Þ N 1;k ðu0 Þ N n;k ðu0 Þ 6 N 0;k ðu1 Þ N 1;k ðu1 Þ N n;k ðu1 Þ 7 7 6 7 N¼6 : ð8Þ . . .. .. 7 6 .. 5 4. . . . N 0;k ður Þ
2.2. Boundary parameterization with B-spline curves
1 Pj ¼ ðPj2 þ Pj1 þ Pj þ Pjþ1 þ Pjþ2 Þ; 5
ð4aÞ
N 1;k ður Þ
N n;k ður Þ
ðrþ1Þðnþ1Þ
Some of the control points obtained from Eq. (7) will be selected as design variables for shape optimization. 2.3. Remesh-free shape optimization using the wavelet-Galerkin method The interest on the shape optimization based on various fixed-grid-type analyses has been recently increased because the fixed-grid analyses offer simple ways to avoid cumbersome remeshing processes which cannot be avoided if the standard finite element approach is employed. Another advantage of the fixed-grid-based shape optimization
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
method is that it requires only the boundary velocity field for design updates while the standard finite-element-based shape optimization method generally requires design velocities for all nodes. Since the analysis nodes are independent of the shape of a given structure in the fixed-grid method, only the boundary velocity field is needed for shape optimization. Most fixed-grid-based analyses adopt the fictitious domain method to deal with domains with general shapes. In Fig. 3, the example of the fictitious domain is illustrated. The real domain (or package domain) x is embedded into the simplified fictitious domain X. Because the fixed mesh of the fictitious domain is independent of the shape of the package domain which changes during shape optimization, there is no mesh distortion problem. The most important issue of the fictitious domain approach lies in the boundary approximation. Poorlyapproximated boundaries will lead to inaccurate stress prediction and spoil the stability of optimization process. The original boundary ox consisting of B-spline curves needs to be approximated for stiffness calculation of the fixed elements. In Fig. 4, various approximation methods for the fictitious domain method are illustrated. Compared to the zigzag approximation in Fig. 4b and the area-fractionbased approximation in Fig. 4c, the approximation with piecewise oblique lines by Jang et al. [10] in Fig. 4d gives accurate stress results on the boundaries. The piecewise boundary lines are approximated by connecting the intersection points between the fixed mesh and the original boundary ox. The stiffness of an element in the fixed mesh depends on its location, which is given by 8 < C0 if element x; Ci ¼ cC0 else if element X n x; ð9Þ : Ci else if element is on bounary;
Fig. 4. The approximation of a circle: (a) the original domain of the circle, and its approximated domains using (b) zigzag, (c) area-fraction-based method, and (d) present approach with piecewise oblique lines.
element B has C0. In Eq. (9), Ci is piecewise-constant over the element. In Fig. 5, the element C in Fig. 3b is illustrated with its original boundary ox (dotted line), approximated boundary oxA (thick solid line), and local element coordinates n, g. The element stiffness matrix for the element is given by Kei ¼ ¼
Z
Z
1
1
Z Z
1
1
BT ðn; gÞCi ðn; gÞBðn; gÞjJjdndg;
BT ðn; gÞC0 Bðn;gÞjJjdndg; xAi 1
Z
1
where c = 0.001 is used in this paper. In Fig. 3b, for example, the element A has the elasticity tensor of cC0 and the
a
5
1
ð10aÞ ð10bÞ
BT ðnðr; sÞ;gðr; sÞÞC0 Bðnðr; sÞ; gðr;sÞÞjJjj b Jðr; sÞjdrds;
1
ð10cÞ
b
A
B
C
D
ω Ω
Fig. 3. (a) The simplified fictitious domain X enclosing the package domain x, and (b) its fixed grid mesh.
6
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14 0
1
2
3
4
5
6
8
7
0
1
0,0
1,0
2,0
2,1
0,1
1,1
2,2
2,3
0,2
Fig. 5. Zoomed view of the element C in Fig. 4b.
0
where B and jJj mean the strain–displacement matrix and the determinant of the Jacobian matrix, respectively, and xAi denotes the region of the approximated x within the element. In Eq. (10c), we introduce additional coordinates: the coordinates r, s which map xAi or Xi n xAi onto the normalized rectangle. Note that the other Jacobian matrix b J is introduced to concern the mapping between the coordinates r, s and the local element coordinates n, g: (see [10] for more details) " # og on or or b¼ J : ð11Þ og os
on os
If Xi n xAi is a triangle as in the element D of Fig. 3b, the stiffness matrix for the element is obtained by subtracting the stiffness of Xi n xAi from the stiffness of an inside element, i.e., Kei
¼
Z
1
Z
1
Z
1
BT ðn; gÞC0 Bðn; gÞjJj dn dg
1
Fig. 6. Equivalent vector space spanned by (a) the single scale basis functions, and (b) multiscale hat interpolation wavelets.
b J , respectively, the system equation for the multiscale and F wavelet-Galerkin method becomes KJ WJ ¼ FJ ;
ð13aÞ
where the subscript J denotes the highest uniform-grid density or resolution level, and WJ is the wavelet coefficients vector for nodal displacements. The relation between bJ; F b J Þ and (KJ, FJ) can be stated by using a transformaðK tion matrix TJ as (see [16]) b J TJ KJ ¼ TTJ K
and
bJ : F J ¼ TTJ F
ð13bÞ
Fig. 7 illustrates the flowchart for the multiscale adaptive analysis where the analysis level increases from a certain resolution level j0 to the highest resolution level J. The adaptive strategy is implemented by utilizing the difference-checking nature of wavelets for error estimation [16].
1
BT ðn; gÞC0 Bðn; gÞjJj dn dg:
ð12Þ
Begin analysis
Xi nxAi
The above element stiffness calculations are necessary only for those elements lying on the boundaries. For elements inside or outside the boundaries, their stiffness matrices are all the same, the first term of the right hand side of Eq. (12) or its c times. Thus, the advantages of fixedgrid-based methods, the fast mesh and fast formulation of the stiffness matrix, are still preserved at a small additional expense required to calculate intersection points and element stiffness matrices on the boundaries Fig. 6. When an adaptive analysis is incorporated into the fixed-grid analysis strategy, the shape optimization can become more efficient. To this end, we will use the adaptive multiscale wavelet-Galerkin method employing hat interpolation wavelets. The details on the method may be found in [16], so we only emphasize that the system matrix (in Eq. (13a)) of the multiscale wavelet-Galerkin method can be easily transformed form the single-scale system matrix, i.e., the standard system matrix resulting from the bilinear finite element formulation. If the single-scale stiffness bJ matrix and the single-scale load vector are denoted by K
Resolution j0 analysis Calculate intersecting points
Build single-scale system equations
^ ^ ^ U K J J = FJ
j=j0
Adaptive scheme j=j+1 Resolution j analysis
No
Wavelet transform
j=J? Multiscale equations
K JWJ = FJ
Yes End analysis
Fig. 7. Flowchart for multiscale adaptive analysis.
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
7
3. Software details and design examples In order to implement the idea of the integrated design for compliant MEMS’s into a user-friendly software, additional considerations besides the analysis modules for topology and shape optimization are necessary as follows. (1) To show the progress of the topology optimization, element densities are expressed as a black–white image. (2) The boundary points of the topologically optimized image are extracted, smoothed, and parameterized with B-spline curves. The commercial software NLib [13] is employed to manipulate the B-spline curves. (3) As a shape optimizer, the DOT library [17] is loaded, and a user interface is graphically expressed to set up the parameters for the optimizer. (4) To determine the design variables and stress constraints for shape optimization, the control points and boundary points estimated to have high stress are graphically selected with a mouse. (5) To calculate the element stiffness matrices of the wavelet-Galerkin analysis, the intersection points between B-spline curves and fixed mesh are calculated using the NLib and are handed over to the solver. (6) The geometry of a structure can be converted to the ANSYS file, and the performance of the optimized structure is confirmed though the commercial analyzer. The overall procedure of the software, including the above processes, follows the integrated design scheme in Fig. 1. We explain the details of the software step-by-step through a couple of examples.
Fig. 8. The design domain for topology optimization of a micro force inverter. (64 · 32 design variables for the upper half of the structure.)
Table 1 The parameters for topology optimization of the micro force inverters defined in Figs. 8 and 17 Young’s modulus Poisson’s ratio Fin kin kout Mesh Number of design variables Thickness M0
1.7 · 105 lN/lm2 0.26 5000 lN 200 lN/lm 2 lN/lm 64 · 32, DSSY 2048 6 lm 30%
3.1. Design of micro force inverters Using the integrated design software, the two micro force inverters are designed. The micro force inverter exerts a force on a target in the opposite direction to an input force while the energy loss by elastic deformation of compliant mechanism is minimized. The topology of a force inverter having maximum performance is determined through topology optimization, and then the local stresses of the structure around hinges are controlled below the allowed level during shape optimization. The problem definition for topology optimization of the micro force inverter is given as [18] Maximize subject to
uout ðqÞ H ðqÞ ¼
designed. The MMA optimizer [19], which is efficient in dealing with many design variables, is used for this problem. The optimized topology of the force inverter is illustrated in Fig. 9. A total of 34 analyses are performed during the optimization. Fig. 10 shows the example of boundary smoothing process for region A in Fig. 9. The dark points denote the extracted boundary from the result of topology optimization, and the bright points mean the smoothed boundary points obtained by using Eq. (2). Before approximating
ð14aÞ Ne X i¼1 T
Z
qi dX M 0 6 0;
ð14bÞ
Xi
q ¼ fqi g ; e 6 qi 6 1; i ¼ 1; 2; . . . ; N e :
0 < e 1;
The objective is to maximize the displacement at the output point, and the total mass of the structure is constrained. Fig. 8 shows the design domain for topology optimization where the dimensions, the input point A, and the output point B are denoted. The parameters for the problem are listed in Table 1. Considering the symmetry of the structure, only upper half of the domain is analyzed and
Fig. 9. The topology optimization result by the proposed software.
8
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
Fig. 10. The smoothed boundary points of the region A in Fig. 10. (Dark circles: original boundary points, white circles: smoothed boundary points.)
the boundaries with use of B-spline curves, the user should determine the order of the curves, number of control points for each curve, and start and end tangential angles of the curves via the dialog box in Fig. 11a. Excessive number of control points often leads to an oscillatory boundary while too small number of control points produces a poorly approximated curve quite different from the original boundary. Fig. 11b shows approximated boundaries of the force inverter by using 11 B-spline curves with second or third orders. The straight lines in the figure denote control polygons of B-spline curves. Note that the hinges are
Fig. 11. The boundary approximation with B-spline curves: (a) the dialog box for B-spline curve information, and (b) the approximated boundaries of the force inverter.
approximated by separate curves (curve number 0, 2, and 4 in Fig. 11b) to increase the degree of design changes. The control points around the hinges are mainly selected as the design parameters for shape optimization. Fig. 12 shows the design parameters around three hinge regions. The proposed integrated design software uses the DOT [17] as mathematical optimization algorithms for shape optimization. To set up various parameters for the optimizer, the graphical user interface in Fig. 13 is provided. Through the dialog box in Fig. 13, the use can determine file names for various purposes, number of design variables, number of stress constraints, upper and lower bound limits for design variables, initial value for design variables, optimization strategy, and other detail parameters for the optimization strategy. Other information required by DOT is • Design variables: Using a mouse, the user can select or unselect them on the screen among the control points of B-spline curves. • Stress constraints: The boundary points expected to have high stress are examined. The parametric value u in Eq. (3) for each point under stress constraint is stored to keep track of the point during the optimization. • Sensitivity: The sensitivities of the objective and constraints are obtained by simply using the finite difference method. The offset value for the finite difference method is determined through the dialog box in Fig. 13. The default for the offset value is set to be equal to the length of the fixed elements. The user can determine the upper and lower bounds for design variables uniformly by using the dialog box in Fig. 13, or can impose different bounds for each variable by picking it on the screen. The design variables lying on domain edges have only one-directional degree of freedom. For example, only horizontal movement of the control point lying on the bottom edge in Fig. 12 can be selected as a design variable. To constrain the stress of the micro force inverter during shape optimization, the boundary regions with high stress need to be located. The proposed integrated design software provides the data exporting function to ANSYS that the geometry of a structure can be converted to the ANSYS input file. After examining the analysis result of ANSYS, the points of high stress are selected, and their von Mises stresses are set as design constraints. Fig. 14a shows the stress contour plot by ANSYS. Based on the result of ANSYS, the points under stress constraints are selected as in Fig. 14b. Note that the points under stress constraints are concentrated around elastic hinge regions. Other information for the shape optimization of the micro force inverter is listed in Table 2. During shape optimization, the wavelet-Galerkin solver requires the intersection points between B-spline boundaries and the fixed mesh. Fig. 15 shows the fixed mesh generation by the software. The intersection points between
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
9
Fig. 12. The design parameters (blue circles). (For interpretation of colour representation in this figure legend the reader is referred to the web version of this article.)
Fig. 13. The user interface for the DOT parameters.
the curves and the fixed mesh are automatically calculated and are used to construct stiffness matrices of the boundary elements according to Eqs. (10) and (12). The optimized result with its deformed shape is illustrated in Fig. 16. Noticeable change can be seen around the central hinge region. The system performance improves 81.4% compared to the structure before shape optimization. Thus the hinge region is regarded to have much higher sensitivity to the system performance than other regions of the micro structure. In Fig. 16, the stress distribution around the hinge verified by using ANSYS is also shown. The maximum stress around the hinge is 982.67 lN/lm2 which satisfies the design constraint. Table 3 lists the output displacements and the maximum stresses of the optimized result.
Fig. 14. Stress constraint imposition: (a) the stress contour (by ANSYS), and (b) the points under stress constraints.
Because a force inverter experiences large displacement, the performance of the optimized result may have to be checked by nonlinear analysis. Thus, ANSYS was used to carry out geometrically-nonlinear analysis of the force converter shown in Fig. 16. The comparison with the nonlinear and linear analyses showed that the output displacement and the von Mises stress by the nonlinear analysis decrease
10
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
Table 2 Parameters used to obtain the results in Figs. 16 and 19 by shape optimization Objective
Output displacement uout
Number of design variables Number of stress constraints Maximum allowable stress Side constraints Optimizer Mesh
42 23 1000 lN/lm2 6 lm 6 qi 6 6 lm SQP 256 · 128, wavelet-Galerkin method
Fig. 15. The fixed mesh of the wavelet-Galerkin analysis (a) for the whole domain and (b) around the hinge A.
by 3.4% and 2.8% , respectively in comparison with the calculations by the linear analysis. So, for the present problem in consideration, the proposed linear wavelet-Galerkin method appears to predict the actual structural behavior quite accurately. Next, a micro force inverter having two output points is designed. Fig. 17 shows the design domain for topology optimization procedure. The design objective is to maximize the output displacements at B and C in the figure. Other conditions are given in Tables 1 and 2. Fig. 18 shows the optimized topological layout and its boundary approximation by B-spline curves. Before extracting curves from the optimized result in Fig. 18a, manual modifications of the resulting topology are necessary. Typical results by topology optimization may include several holes, some of which can be removed before applying shape optimization because their contributions to mass constraints and structural performance may be insignificant. For instance, a few small holes in Fig. 18a are eliminated during the modification and the modified layout is illustrated in Fig. 18b. Nevertheless, the system performance should be confirmed by additional analysis on the modified topology. To facilitate the modification process, the proposed integrated design software has an editing function. For shape optimization, 7 control points around hinges are selected for design parameters, and 48 points around high stress regions, mostly around hinges, are selected for design constraints. The final shape of the force inverter and its deformed shape are illustrated in Fig. 19. The input force is inverted at the output points via three elastic hinges. The stress distributions of the hinges are also illustrated in Fig. 19. Now we have presented two design examples, le us explain the date structure and communication flow of the proposed software: see Fig. 20. The input data for topology optimization is provided in a text file that contains design domain information, boundary conditions and material
Fig. 16. The deformation shape of the optimized micro force inverter and the von Mises stress distribution around the hinge (verified by using ANSYS).
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
11
Table 3 Output displacements and stresses after shape optimization Output displacement [lm]
Force inverter in Fig. 6 Force inverter in Fig. 19 Gripper in Fig. 21
Initial value
Optimized value
Performance increase [%]
1.893 0.648 0.744
3.435 2.563 2.260
81.4 295.5 203.8
Maximum von Mises stress of the optimized structure [lN/lm2] by wavelet-Galerkin method (by ANSYS) 953.50 (982.67) 974.16 (975.63) 983.29 (1021)
Fig. 17. The design domain for topology optimization of a micro force inverter with two output points (64 · 32 design variables for the upper half of the structure).
properties. The optimized densities q are transferred to the boundary handling program as a form of a data array, from which the boundary information data including smoothed boundary points and B-spline curves are
Fig. 18. (a) The topology optimization result of the force inverter, and (b) its boundary approximation.
Fig. 19. The optimized shape of the micro force inverter with two output points and von Mises stress distribution around hinges (verified by using ANSYS).
12
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14 Topology optimization module Problem configuraion file (text file) 1. Design domain dimensions 2. Spring stiffness 3. Input and output points 4. Boundary conditions 5. Material properties
Scalar data: design domain size Array data: rho.dat (design variables, text file)
Boundary handling module (Extraction, smoothing and parameterization) Wavelet-Galerkin solver
Scalar data: number of boundaries Array data: 1. bnd.dat (boundary points, text file) 2. spline.dat (B-spline curve information, text file)
Array data: int.dat (text file) intersection points
Scalar data: objective value Array data: 1. cnst.dat (constraint values, text file) 2. objgrd.dat, cnstgrd.dat (text files) sensititvities for objective and constraints
DOT optimizer (dot.dll)
Shape optimization module
Array data: Updated design parameters
Fig. 20. Data structure and communication flow chart of the integrated design software developed in this study.
extracted. The shape optimization program manages three data groups: one for boundary manipulation, another for optimization algorithm, and the other for data communication with the analyzer. 3.2. Design of a micro gripper A micro gripper grasping and moving a micro-scale object is designed in this example. The design domain for topology optimization and its dimensions are shown in
Fig. 21. The design domain for topology optimization of a micro gripper (kin = 400 lN/lm, kout = 20 lN/lm and Fin = 4000 lN).
Fig. 21. The objective is to maximize the displacement at the output point B. To increase the gripping force of the structure, the stiffness of the output spring 10 times large than that used in the previous examples is chosen. 10 times larger stiffness for the output spring than in the previous examples is used. The result by the topology optimization is shown in Fig. 22a. After boundary extraction and smoothing, the approximated boundaries by B-spline curves are illustrated in Fig. 22b. A total of 8 B-spline curves are employed, and the orders of the curves vary from linear to cubic splines. Note that 6 small holes in Fig. 22a are removed through the topology modification procedure. Next, the regions with high stresses are located by exporting the boundary data to ANSYS and the points subject to stress constraints are selected. The control points around hinge regions are selected as design variables, which are denoted in Fig. 22b. A total of 11 design parameters and 54 stress constraints are employed for shape optimization. Fig. 23 shows the final design and its deformation shape. Table 3 shows that due to high stiffness of the output spring, the output displacement of the gripper is relatively small. Fig. 24 shows the distribution of the adaptivelyadded wavelet bases of the wavelet-Galerkin method. Note that the support sizes of the multiscale wavelets are different according to their resolution levels. The dots in the figure denote the centers of the wavelets. In Fig. 24, it is observed that while few wavelets are located on the domain Xnx having weak material, most of wavelets exist within
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
13
Fig. 22. (a) The topology optimization result, and (b) boundary approximation (blue circles denote design parameters). (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 23. The optimized shape of the micro gripper and the von Mises stress distribution around the hinge (verified by using ANSYS).
the package domain, and high stress regions require dense distribution of wavelets. 4. Conclusions The integrated design software for complaint MEMS actuators was developed. Using topology optimization, a rough image for an actuator having high overall system performance was obtained. Nonconforming finite elements were used to overcome the numerical instability of checkerboards or one-point hinges. After boundary extraction, smoothing, and approximation with the B-spline curves, remesh-free shape optimization based on the multiscale wavelet-Galerkin method was performed to control local stresses around hinges. Through design examples, the validity and effectiveness of the proposed integrated design software were verified.
Acknowledgements
Fig. 24. The distribution of adaptively-added wavelets (upper: x-directional wavelets, lower: y-directional wavelets).
This research was supported by the National Creative Research Initiatives Program (Korea Science and Technology Foundation Grant No. 2006-033) contracted through
14
G.-W. Jang et al. / Advances in Engineering Software 39 (2008) 1–14
the Institute of Advanced Machinery and Design at Seoul National University. References [1] Sigmund O. Design of multiphysics actuators using topology optimization – part I: one material structures. Comput Methods Appl Mech Engng 2001;190:6577–604. [2] Chen BC, Silva CN, Kikuchi N. Advances in computational design and optimization with application to MEMS. Int J Numer Meth Engng 2001;52:23–62. [3] Duysinx P, Bendsøe MP. Topology optimization of continuum structures with local stress constraints. Int J Numer Meth Engng 1998;43:1453–78. [4] Olhoff N, Bendsøe MP, Rasmussen J. On CAD-integrated structural topology and design optimization. Comput Meth Appl Mech Eng 1992;89:259–79. [5] Tang PS, Chang KH. Integration of topology and shape optimization for design of structural components. Struct Multidiscip Opt 2000; 22:65–82. [6] Bennett JA, Botkin ME. Structural shape optimization with geometric description and adaptive mesh refinement. AIAA J 1985; 23(3):458–64. [7] Yao TM, Choi KK. 3-D shape optimal design and automatic finite element regridding. Int J Numer Meth Engng 1989;28:369–84. [8] Kim NH, Choi KK, Botkin ME. Numerical method for shape optimization using meshfree method. Struct Multidis Opt 2003; 24:418–29.
[9] Garcia MJ, Steven GP. Fixed grid finite elements in elasticity problems. Eng Comput 1998;16(2):154–64. [10] Jang GW, Kim YY, Choi KK. Remesh-free shape optimization using the wavelet-Galerkin method. Int J Sol Struct 2004;41:6465–83. [11] Jang GW, Jeong JH, Kim YY, Sheen D, Park C, Kim MN. Checkerboard-free topology optimization using nonconforming finite elements. Int J Numer Meth Engng 2003;57:1717–35. [12] Sigmund O, Petersson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, meshindependencies and local minima. Struct Opt 1998;16:68–75. [13] Solid Modeling SolutionsTM. Available from: http://www.smlib.com/ nlib.html. [14] Yoon GH, Kim YY, Bendsøe MP, Sigmund O. Hinge-free optimization with embedded translation-invariant differentiable wavelet shrinkage. Struct Multidis Opt 2004;27:139–50. [15] Lee K. Principles of CAD/CAM/CAE systems. Addison Wesley; 1999. [16] Jang GW, Kim JE, Kim YY. Multiscale Galerkin method using interpolation wavelets for two-dimensional elliptic problems in general domains. Int J Numer Meth Engng 2004;59:225–53. [17] Vanderplaats Research & Development, Inc. DOT Users Manual, 1999. [18] Sigmund O. On the design of compliant mechanisms using topology optimization. Mech Struct Mach 1997;25(4):495–526. [19] Svanberg K. The method of moving asymptotes – a new method for structural optimization. Int J Numer Meth Engng 1987;24:359–73.
