UFDC Home  myUFDC Home  Help 



Full Text  
TOPOLOGY OPTIMIZATION UNDER STRESS CONSTRAINTS By TAEJOONG YU A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2003 Copyright 2003 by TaeJoong Yu ACKNOWLEDGMENTS This thesis would not have been completed without the help of several people to whom I wish to express my sincere appreciation. First, I wish to express my deepest gratitude to Dr. Ashok V. Kumar for his immense help. He is my best advisor and best mentor. He has provided excellent guidance all the time for my studies with infinite kindness and patience. I could never have imagined a better advisor in the world. Jongho Lee helped me to work out all problems and questions I had about programming in C++ and Java. He opened my eyes to the exciting world of object oriented programming. It became an exciting experience with his knowledge and endless help. I would also like to show my appreciation to other friendly and helpful lab mates: YoungNam Hwang and Shivakumar. They always gave me answers and academic stimulations. To my family in Korea, I would like to express my gratitude for their dedication, support, and love. We have been apart for three years, but my heart was always with them and their thoughtful help and care were driving forces for this achievement. I would like to acknowledge my dear Korean friends in Gainesville: Kyungha and Dongsik, Kyungmin and Sangwon, Sonho and Jawoo Ever since we met at the Atlanta airport on June 22, 1999, we have been through everything together in the U.S. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ................................................................................................. iii L IST O F FIG U R E S .... .............................. ....................... ........ .. ............... vi A B S T R A C T .......................................... .................................................. v iii CHAPTER 1 IN TR OD U CTION .............................................. .. ....... .... .............. . Stru ctu ral O ptim ization ............................................................. ................................ 2 O objective and M motivation .................................... ..................................................... 3 2 STRU CTURAL OPTIM IZATION .............................................................. ............... Sizing Optim ization .......................................................... ...... .......... 4 Shape Optimization....................... .. ............... 4 Topology Optimization ....................... .................... 6 H om og enization M eth od ................................................................................................ 7 3 TOPOLOGY OPTIMIZATION using SHAPE DENSITY FUNCTION....................12 Sh ap e D en sity F u n action ...................................................................... ................ .. 12 Objective Function................................................................. 15 Relationship between Material Properties and Density .............................................. 16 Implementation by Finite Elements ............... ............................................ 17 Solution Procedure ........... ...................... .......... ............... 18 4 STR E SS C O N STR A IN T S ...................................................................... ..................2 1 M ean Com pliance ...................................... ............ ................. 21 P problem F orm ulation ..... .. .. .. ................................................ ................. .. 23 5 R E S U L T S ............................................................................. 2 8 Effect of Material Property Density Relationship ..................................................... 28 E x am ple 1 ........................................ 2 8 E xam ple 2 ........................................ 30 E x am p le 3 ...................................................................... 3 2 Effect of D esign Space (Exam ple 4)................................................... .... ................ 33 Effect of Safety Factor (Example 5) ............... ....... ........................ 35 C om prison w ith ID EA S A nalysis.................................................... ... ................. 37 6 C O N C L U S IO N ................................................................................................. .. 4 0 APPENDIX COMPUTER SOURCE CODE FROM OPTIMIZATION PROGRAM.....42 L IST O F R E F E R E N C E S ........................................................................ .....................46 B IO G R A PH IC A L SK E TCH ..................................................................... ..................49 v LIST OF FIGURES Figure page 21 Structural Optimization of sizing, shape and topology .............................................6 22 Homogrnization....................................................... ........8 23 Material Property as a Function of Density.................................... ......... ............9 31 Shape R presentation ...................... .. .......................... .. .. ...... ........... ... 13 32 Solution Flow chart ................................................ .. ............. ......... 19 5 1 E x a m p le 1 ............................................................................................................. 2 8 52 Example 1, When n=2, The Ratio of a/y is 332367......... ....... ........ ............... 29 53 Example 1, When n=3, The Ratio of a/y is 498550...............................................29 5 4 E x am p le 5 .......................................................................... 3 0 55 Example 2, When n=2, The Ratio of a/y is 332367..................................................31 56 Example 2, When n=3, The Ratio of a/y is 498550..................................................31 5 7 E x am p le 3 .......................................................................... 3 2 58 E xam ple 3, W hen n= 2 ....................................................................... ................... 33 5 9 E x a m p le 4 ............................................................................................................. 3 3 510 Example 4, When the Height is 0.1m......................................................... 34 511 Example 4, When the Height is 0.2m .............. ............................. 34 512 Exam ple 4, W hen the H eight is 0.3m .................................... ........ ............... 35 5 13 E x am p le 5 ......................................................................... 3 6 514 Example 5, When n=2 and the Safety Factor is 1 .................................................36 515 Example 5, When n=2 and the Safety Factor is 1.5 .............................................37 vi 516 A naly sis of E xam ple 5 ....................................................................... ..................38 517 Analysis of first modification of Example 5 ................................ ..... ..........38 518 Analysis of second modification of Example 5 ...................................................39 vii Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science TOPOLOGY OPTIMIZATION UNDER STRESS CONSTRAINTS By TaeJoong Yu May 2003 Chair: Ashok V. Kumar Major Department: Mechanical and Aerospace Engineering In this thesis, shape and topology optimization of structures subject to stress constraints is studied. The geometry is represented using a shape density function interpolated over quadrilateral elements. This technique allows both the shape and the topology to be treated as a variable during the optimization process. A C++ object oriented finite element module was implemented to perform structural analysis and compute stiffness matrix. The optimization module iteratively modifies the geometry by calling the finite element module. To display the optimal shape and the stress distribution, a Java program was implemented. Multiple objectives are used to compute the optimal structure by simultaneously maximizing the stiffness and minimizing the mass. Maximizing stiffness is identical to minimizing the compliance of the structures. Therefore, the objective function to be minimized is defined as a weighted sum of the compliance and the mass of the structure. Assuming that the stress is more or less uniform in the optimal design, a relation between the compliance and mass can be derived for the optimal structure. Using this relation a ratio between the weighting factors in the objective function is computed such that the stress in the structure is equal to the design stress. It was found that by using this ratio between weight factors the optimal structures do have almost uniform stress distribution with stress values very close to the design stress everywhere except at stress concentrations. Due to constraints on available space for the structure, the stress concentration often occurs at the boundaries of the feasible domain. The stress can be reduced at these stress concentrations by using larger safety factors and also by increasing the size of the feasible domain when possible. CHAPTER 1 INTRODUCTION Structural topology optimization is the selection of the best configuration for the design of structures. In the last decade, the field of structural topology optimization has expanded significantly, successfully addressing many practical engineering problems. Structural optimization can be defined as maximizing a certain structural property subject to one or more structural constraints. The ability to design optimal structures allows the minimized material usage and maximized structural properties. This results in lower manufacturing costs, which can be passed on to the costumer as lower retail costs. Previous research on topology optimization focused primarily on global structural behavior such as stiffness and weight. However, to obtain a true optimum design of a structure, stresses must be considered. For example, to obtain the optimal geometry of a structure with the objective of minimizing compliance, usually a limit on its volume or mass is specified as the constraint. However, it is not clear what percentage of the material volume should be sufficient for supporting the applied loads. By using optimization method based on such a problem formulation, a trialanderror process cannot be avoided if the designer really wants to select the minimumweight optimal design for the given domain. Furthermore, it is possible that there are areas of stress concentration on the optimal structure. In practical optimization problems, there exists several design criteria that must be considered in the optimal design of structures. For example, reducing the weight and increasing the stiffness of structure are both typical design goals in a structural design, but these objectives conflict with each other. In order to accommodate many conflicting design goals, the sequential application of each single objective optimization can be considered. However, this method cannot produce an optimized solution because only a single objective is considered in each optimization process. As a consequence, socalled multiobjective optimization, which considers multiple objectives simultaneously, has recently been regarded as a methodology for solving optimization problems with several objective functions. A number of techniques and applications of multiobjective optimization have been developed over past few years. A comprehensive overview of the field of multiobjective optimization in mechanics was introduced by Stadler [Sta 84], and multiobjective design optimization was applied to engineering problems, including structural design problems, by Eschenauer et al.[Ech 90]. Structural Optimization Structural optimization can be classified into three categories: size, shape and topology optimization. Sizing optimization starts at a structure in which the configuration of the structure is already defined. It seeks the optimum combination of the element size such as length, thickness, crosssectional area, etc. It does not involve the redefinition of shapes of outer boundaries and inner holes of structure. For this reason, if a suboptimal initial topology is chosen the final optimization solution will be also be suboptimal. Shape optimization involves varying the domain of the structure. It seeks the geometric definition of the boundaries of outer circumference and inner holes of the structure. Design variables in shape optimization include often the location of key points that define a curve or some parameters that define a specific predetermined geometric shape such as edges or boundaries. During the optimization procedure, these boundaries may be manipulated; therefore new shape can be created. If the shape is changed significantly during the optimization procedure, the structural analysis package must recreate the mesh in order to ensure minimal element distortion. The primary disadvantage to shape optimization is similar to that of sizing optimization. If the shape optimization is based on a specific geometric configuration, any design that is not included in the set of the predetermined topology will not be created during the optimization procedure. Topology optimization aims to create the best configuration of a structure from a userdefined configuration set that comprised all possible configuration of a structure for the design problem. This means that new geometric boundaries, such as the boundary of an internal hole, may appear during the optimization procedure. This allows global shape and topology optimization so that the final optimal structure no longer depends on the initial design. Objective and Motivation This thesis presents the implementation of a structural topology optimization under stress constraint. The objective of optimization procedure is to minimize the compliance and mass of the structure, such that the stress is not over the yield stress at any point in the structure. The ideal optimal structure has uniform stress distribution with the stress lower than the yield stress. CHAPTER 2 STRUCTURAL OPTIMIZATION Structural optimization can be classified as sizing, shape, or topology optimization. They are classified by the design variables chosen to describe the geometry. They may also be classified based on the particular combination of objective function and constraints used to describe the optimization problem. The following sections give a brief description of each type of structural optimization techniques. Sizing Optimization Sizing optimization starts at a structure in which the configuration of structure is already defined, such as trusses, frames and plates [Haf 86, Haf 92, Kir 81, Mor 82]. It seeks the optimum combination of elemental size such as length, thickness, cross sectional area, etc. The geometry change is very small when varying these design variables, so it does not involve the redefinition of shapes of outer boundaries and inner holes of a structure. One of the primary disadvantages of sizing optimization is that the topology of the structure remains fixed throughout the optimization procedure. Therefore, if a suboptimal topology is chosen when formulating the optimization problem, the resulting structure will also be suboptimal. The optimal design of a sizing optimization is only best design that can come out of the predetermined structural geometric definition. Shape Optimization Shape optimization seeks the geometric definitions of the boundaries of outer circumference and inner holes of structure. Shape optimization requires the finite element model to change during the optimization procedure. The design variables used for shape optimization cause boundary variation and these are referred to as the "shape design variables". Compared to sizing optimization, the computational cost of shape optimization is higher due to the need to constantly update the finite element model. Shape optimization is divided into two categories; parametric variable variation and boundary variation. In parametric variable variation, the design variables are parameters defining certain features of the shape or important dimensions. For instance, the side of a square hole or the radius of a circular hole could be a design variable. Many examples of shape variation using parametric variables have been illustrated by Botkin et al. [Bot 86]. In boundary variation, parts of the boundary of the solid are treated as the design variables. For instance, the nodal coordinates of the nodes on the boundary of the shape could be the design variables. Yang [Yan 86] showed many examples of this type of shape optimization. A thorough review on shape optimization methods and applications was provided by Haftka and Grandhi [Haf 86]. The geometric configuration of a structure is required before the shape optimization can be performed. If the shape optimization is based on a specific geometric configuration, any design that is not included in the set of the predetermined geometric modeling will not be created during the optimization. Therefore, the shape optimization converges to different optimal shapes for different starting topologies. Another drawback to shape optimization is the finite element model must be continually updated to ensure that the mesh does not become highly distorted. This task is difficult because mesh refinement and recreation must be performed automatically during the optimization procedure. Topology Optimization From above sections, it was realized that to obtain optimal shapes the topology must be modified during the optimization procedure. Early works by Dorn et al. [Dor 64] and Dobbs and Felton [Dob 69] implemented topology optimization techniques in which simple numerical methods were used to find optimal layout and geometry for truss. Most topology optimization research in the next two decades was devoted to structural analysis using the finite element method followed by removal of understressed element. The final optimized shapes were highly dependent upon the density of the mesh used for the finite element analysis. Strang [Str 86] attributed this behavior to the nonconvex nature of the problem. Kohn and Strang [Koh 86], representing that the original problem was ill posed, proposed a relaxed variational problem, which allows composite or porous material instead of an allornothing dichotomy between material and holes. Bendsoe and Kikuchi [Ben 88] proposed the homogenization method to determine the properties of these composite materials. sizi nl shape optimization topology optimization ooo rsie of section III _ Figure 21. Structural Optimization of sizing, shape and topology The homogenization method is a mathematical tool for modeling the behavior of periodic structures, such as composites or porous materials, whose properties vary periodically at microscopic level [Ben 78]. The homogenization method has also been successfully applied to many structural design problems [Kikuchi 93, Ben 91] in which the objective function was the compliance of the structure with multipleloading cases and threedimensional applications. Homogenization Method The main idea of structural topology design using homogenization method is that a solid structural domain is modeled as a composite material with possibly perforated microstructures. Consequently, the homogenization method is utilized to analyze the composite structure. The behavior of periodic structure at a macroscopic level is predicted in the limit when the E, (the ratio of period of the structure to a typical macroscopic length scale) tends to zero. The early approaches to topology optimization allowed the material to either exist in a state of full density, or not exist at all (an allornothing approach). Areas void of material shows either holes, or areas outside of the structure's outermost boundary. During this optimization procedure, understressed element will be removed. This method may have difficult convergence problem. To overcome some of these difficulties, Kohn and Strang [Koh 86] proposed the use of homogenization as a means of allowing a uniform distribution of properties from 0 to 1. The difficult convergence problem suggested the need to relax the original variational problem statement so that it would be wellposed and have better convergence performance. Kohn and Strang showed that relaxing the variational problem is identical to allowing composite materials for the solution. By means of the homogenization method, the structural topology optimization problem is used as the optimal material distribution problem. By assuming the material to be porous and solving for the optimal distribution of porosity, Bendsoe and Kikuchi [Ben 88] have elaborated on the above finding. A design domain is defined as the space in which the structure must fit. This domain is divided into a rectangular mesh with respect to the given boundary conditions and loads. 1 1 A A B Figure 22. Homogmization. A) Unit Cell, B) Unit Cell Orientation in an Element The finite element method is used to analyze the structural behavior under the boundary conditions and the applied loads. Bendsoe and Kikuchi modeled the material as porous by assuming a specific microstructure. The material is represented as a porous medium containing infinitely many microscale cells. Suzuki and Kikuchi [Suzuki 91] have assumed a rectangular hole having edge dimension of 'a' and 'b' as shown in Figure 2 2(A). The dimensions of the hole within a unit cell determine the overall porosity or hole fraction of the material. Each finite element in the mesh is assumed to have a fixed porosity or hole fraction associated with its own hole dimensions 'ai' and 'bi' where 'i' is the element number. These hole dimensions along with the angular orientation of the unit cells ,0 are the design variables. Initially, all elements are assumed to have the same hole dimensions and orientation giving the material a uniform porosity. To solve for the optimal porosity distribution, Bendsoe and Suzuki [Ben 88, Suzuki 91] used an optimality criteria algorithm. The material properties are continuous function of the hole dimensions. For a given porosity or hole size, the material properties can be determined using the homogenization method. The relation between material properties and coefficients and hole size is given in Figure 23. D 1a2 Figure 23. Material Property as a Function of Density Figure 23 shows an example in which the hole within each unit cell is assumed to be square (a = b), so that 1a2 represents the density of the unit cell. If no hole exists in the unit cell, the material properties are the same as those of the fully dense material. When the hole size increases, the density of the material decreases and the material property coefficients decreases. For a given porosity, the homogenization method can be used to determine the material property coefficients Di. A finite element analysis over the unit cell is required to calculate the material property coefficients for each value of hole size. During the optimization procedure, the hole size is increased in areas where the stresses are low, which results in lowdensity material. On the contrary, the hole size is decreased in areas where material has high stresses, resulting in highdensity. When optimization procedure converges, some elements will have large hole sizes so that they have low density and other will have zero hole size (fully dense material). Elements with density value below a specific threshold value can be removed to reveal the final optimized topology. The relation between material property coefficients and hole size would be different for different microstructures. In the method adapted by Bendsoe and Kikuchi [Ben 88], a square or rectangular hole was assumed within a unit cell. Recently, Allaire and Kohn [All 92] have used rank 2 microstructures. Figure 24(a) shows a microstructure consisting of unit cells with rectangular holes and Figure 24(b) illustrates a rank2 microstructure. The rank2 consists of alternating layers of stiff material and rank1 material. Rank1 material in turn consists of alternating layers of stiff and soft material. Allaire and Kohn [All 92] have demonstrated that for a composite material with a fixed ratio of the components, the rank2 layer yields the stiffest composite. The optimal shape, which is solved for using rank2 layering, does not have clearly defined holes or boundaries. Instead, the optimal results obtained suggest that the truly optimal structure may have continuously varying density due to continuously varying ratio of two constituents. Since it is not economically or physically possible to control the density and ratio of constituent within a structure, such composites are difficult to manufacture. The geometry of the structure obtained by the above method does not yield smooth boundaries, especially when a coarse finite element mesh is used. 11 . ... .. . A' A B Figure 24. Microstructures. A)By Microscale Rectangular, B)By Rank2 Papalambros et al. [Pap 90] used imageprocessing techniques to extract a feasible initial shape from the topology generated by the homogenization method. They then used this shape as the starting shape for conventional shape optimization by boundary variation. A similar twostep process has been used by Bendsoe [Ben 91] to obtain combined topology and shape optimization. CHAPTER 3 TOPOLOGY OPTIMIZATION USING SHAPE DENSITY FUNCTION Shape Density Function In order to model the geometry as a variable, we define the boundaries of the geometry to be the contours of a density field O(x), which we refer to as the shape density function. The shape density function has values in the range from a threshold value 0, to 1. The contours of this function corresponding to the threshold value of density express the boundaries of the geometry. The density field is defined over a feasible domain Q0. The shape being expressed may be defined as the subset of this feasible domain within which the density value is greater than or equal to the threshold value. In other words, the geometry consists of the set of x e Qf for which 4 < (x) < 1 (31) Figure 31 shows geometry representation using a shape density function. The L shape region shown is the feasible domain within defined geometry. The final shape must fit within this region. The figure shows contours of the shape density function according to constant values of densities. The fully dense regions (where q = 1) are shown in black color and the white regions have density less than the threshold value. Therefore the white regions are void of material. All other density values are represented using a grayscale color distribution with darker shades representing higher density than lighter shades. Representation of shape using the shape density function allows the entire geometry to be treated as a variable, thus we can combine shape and topology optimization. In the shape optimization procedure, the internal hole shown in the figure would not be created if it were not included in the initial design provided by the designer. Figure 31. Shape Representation To define the shape density function within the feasible domain, the domain is divided into a quadrilateral mesh and the nodal values of density are specified. The density function distribution within each element is obtained by interpolating the nodal values. Using 4node quadrilateral elements, the density function within each element can be expressed as, q(s, t) = N, (s, t) + 02N2 (s, t) + 03N3(s, t) + 04N4 (s, t) (32) where, are the nodal density values and N, (s, t) are the isoparametric shape functions for the 4node quadrilateral element expressed in terms of the parametric coordinates s and t [Bathe 87]. The shape function can be expressed as 1 1 N, = (1+s)( +t), N2= (1 )( +t) 4 4 N3 = ( s)( t), N4= ( + S)( t) 4 4 For an isoparametric element, the mapping between the parametric space and the real coordinates (x, y) is defined by 4 4 x = x,N,(s,t) and y= yN, (s,t) (33) 1 1 The geometry can be graphically displayed by plotting gray scale border of the density function where different shades of gray are used for different ranges of density values. Since a bilinear interpolation is used to represent the density distribution within each element, the contours of the density are linear within each element. Contours are plotted by connecting points of equal density on the edges of an element. A Co continuous shape density function ensures Co continuous boundaries for the final shape, providing a means of combined shape and topology optimization. To have a fully dense final design, it is desirable that the density changes sharply from the threshold value to 1 at the boundary. However, since the shape density function is a Co continuous function, it cannot change discontinuously from the lowest value to the highest. The sharpest transition that can be represented is limited by the size of the quadrilateral elements used for the piece wise interpolation of this function. The density values can change from highest value at one end of the element to the lowest value at the other end. Therefore, the smaller the element size the sharper the density variation can be at the boundaries. From a purely shape representation point of view, if we define the threshold value to be equal to or nearly equal to 1, then within the boundaries of the geometry the material would be fully dense. However, when applying this shape representation scheme for topology optimization, a lower value of threshold is desirable so that optimal designs that are not fully dense can also be represented. The design problem definition consists of specifying the feasible domain, the loads the structure has to carry and the displacement boundary conditions that describe how the structure is supported. In Figure 31, the arrows at the top represent a uniformly distributed load supported by structure. The solid triangles represent nodes that are constrained to have zero displacement. Objective Function For the maximum stiffness topology design, minimizing the mean compliance of a structure is commonly used as the objective function and the constraint is imposed on a somewhat arbitrarily chosen material volume. However, it should be noted that the designer usually does not know what percentage of the material volume should be sufficient for supporting applied loads. Therefore, we should consider minimizing both the mean compliance and the mass as the objective. The design goal is to obtain the optimal topology of a structure, which maximizes the stiffness and minimizes the mass, which ensures that stress is lower than yield stress. Therefore, multiple objectives should be considered. We limited the study to planar problems such as plane stress and plane strain. The multiobjective problem may be written as Minimize : aM()) + 7L(x(o)) (34) where a and y are user specified control parameters. The mass, M(q!), can be expressed as M() = qd = {AT {(p (35) where {A} [ {N}/Q and {p} is nodal density values The mean compliance, L(x), is written by L(x(q)) =f f x())dQ + t x(()dF (36) Q F This mean compliance is twice the work done by the applied force (traction t and body force f) during the displacement x. To obtain the mean compliance, we can use additional equation, J {T [D(:)]{Jdo = L(sx) (37) Q0 {t} and {&c} are the strain and virtual strain in the structure caused by the displacement x and the virtual displacement Sx respectively. [D(O)] is the matrix of elastic constants that relate stresses and strains for a linear elastic material. Assuming that these elastic constants are functions of the shape density function, as the shape changes (0 changes) the stiffness, the compliance and the deformation of the structure for the given load also change. The optimization problem given by the equation 34 through 37 can be solved using a nonlinear programming algorithm. This work uses a variation of sequential programming, which is described in Kumar and Gossard [Kum 93a] and Kumar [Kum 93b]. Relationship between Material Properties and Density In order to affect the outcome of the finite element analysis with the shape density function, some relationship between material properties and density has to exist. The relation between the elastic modulus and the density has to be selected such that when material is removed, the structure should become more compliant. When the minimization of the compliance and the mass is subject to the stress constraint, the algorithm will distribute density such that material is available in locations where it contributes to minimizing compliance and mass. In this thesis, its relation between material properties and density is considered as E = E0~ (38) where, E, is the Young's modulus of the material of the structure, q is the density. Implementation by Finite Elements The deformation of the structure and its mean compliance were computed by finite element method. To construct the stiffness matrix, we integrate the principle of virtual work, over each element to obtain the element stiffness matrix and then assemble these together to obtain the global stiffness matrix. In this thesis, we have used fournode isoparametric quadrilateral elements for the analysis. Therefore, the displacements are represented by piecewise bilinear interpolation within each element. This representation yields the following straindisplacement relation for the element {e}= [B]{U} (39) aN, aN2 aN3 aN4, where, [B]= O x and t( ) is the displacement vector whr, ON2 ON3 ON4 8y 8y Oy Oy corresponding to the quadrilateral finite element. Substituting this relation in (37), for each element, we get Le(u) =f{&se} [D( e)]{g =}dQ {&}Q [J[B]Y [D(I)B]d[ e ]{4} (310) L'(Su) = {u& } [Ke ]' }, where, [K ] is the element stiffness and {(u } is the virtual displacement vector. The element stiffness matrix [K ] can be expressed as [K = J[B] [D(0)][BMd (311) Since the quadrilateral element has a bilinear interpolation for the displacement field within each element, [B] is not a constant matrix. The terms of matrix are linear functions of the parametric coordinates s and t. The elasticity matrix is a function of the parametric coordinates, since it depends on the density function. The integration in Eq. (311) can be calculated by using gaussian quadrature approach. Before using the gaussian quadrature approach, we need to determine the correct order of quadrature. To get this order, the degree of the polynomial terms of the [B] [D(b)][B] matrix must be determined. Since density is a linear function in the parametric coordinate s and t, the terms in the matrix [Bj [D][B] are polynomials of degree 2+n (where n is the degree of the polynomial in (3 8). Gaussian quadrature of order m (using m x m integration points) can integrate a polynomial of degree 2m1. Solution Procedure Two primary modules solve the topology optimization program; which are the optimization module and the finite element module. The optimization module controls the solution procedure, and calls the finite element module when the objective function is calculated. The first step of the optimization procedure is reading an input data file. Once the input data has been processed, the optimization module is called. The purpose of the optimization module is to calculate the nodal density values within the design domain that produce a structure having minimal compliance and minimal mass while simultaneously satisfying the stress constraint. Input Data Processed Optimization module starts the optimization Finite element module calculates the compliance and its gradient and also the mass and its gradient Optimization module modifies nodal densities. Convergence SNo Yes END Figure 32. Solution Flowchart To determine the amount to change the nodal density value, the optimization module has to calculate the objective function and the gradient of the objective function. This is evaluated by calling the finite element module. The finite element module uses nodal density values along with the material property density relation to execute a finite element analysis [Kum 00]. By the finite element analysis, the nodal displacements are calculated. The finite element module to calculate the compliance and its gradient values uses these nodal displacements. Then the optimization module uses these values to vary the nodal density values, resulting in the necessity to perform another finite element analysis. This iterative procedure continues until all convergence criteria are satisfied. With the final nodal displacement, the stress values are calculated. These stress values and the final density values are written to an output file to be used in the graphical display program. Figure 32 illustrates a flowchart showing the solution procedure. CHAPTER 4 STRESS CONSTRAINTS Mean Compliance The mean compliance is the same as twice of the total strain energy at equilibrium. Bendsoe [Ben 88] describes its relation. The mean compliance can be written as L(E) = JETDed = cTCad (41) D is the matrix of elasticity constants that relate stresses and strains for a linear elastic material. C is the compliance matrix. When we assume that the material is isotropic, C is given by C v 1 (42) v v 1 where v is Poison's ratio, E is Young's modulus, and aT = {o 2 o3 is the vector of the three principal stresses. In terms of the principal stresses a1, o2, and o3, the von Mises stress can be calculated as 2 )2 von = 72V( ' 2 +2 (2 3 2 +(3 1 )2 (43) A constraint on the allowable values of the von Mises stress in Q can be expressed as 2 ciTMi. where omax is given maximum allowable stress and 1 1 2 2 1 1 M= 1  (45) 2 2 1 1 2 2 In the case of plane stress, the following relation can be obtained [Ben 88] 1 2(1+v) TMl TC < 1 2(1 _v)aTMc (46) E 3 E Equation 46 shows that the von Mises stress is bounded by a multiple of the square root of strain energy density at an arbitrary point in Q, i.e., 3E vCn iT) or (47) 2(1 + v) f2dA: 3E fTC 3E oa dA< a ciTCc ~ L(E) (48) O 2(1 + v), 2(1 + v)L() (48) Using this result, the relation between the mean compliance problems and problems in which the stress is bounded over the design domain can be obtained. Assuming that the stress is nearly constant and the von Mises stress is almost the same as the maximum stress over the optimum layout, the mean compliance can be computed with the mean stress. Therefore, Eq. 48 can be simplified as 2 3E max A < L(E) (49) 2(1 + v) where A is the area (2D) of Q. We obtain a lower bound on the mean compliance by L() > 2 ) A (410) 3E Problem Formulation Our objective is minimizing the mass and the mean compliance with stress constraints. That is, minimize AM(q) + yL(x(q)) (411) 0< 0 1 (412) where a and y are user specified control parameters. M(0S), the mass can be expressed as M() = fJqdQI = {A}T{ ,where {A}= f {N} (413) Using the straindisplacement relation, {(8= [B]{x},assumed in FEA model, the mean compliance equation Eq. 41 can be written as L(x) = {x [B]T [D][B]{x}d (414) where, {x} is the nodal displacement vector. We can rewrite this equation as L(x) = ( x [B]T [D][B]Tx)dQ = (x' [ B]T [DI][BIIIdx (415) L(x) = {x} [K]{x}, where [K] is the stiffness matrix. [K] [B]T [D][B]jd (416) In the objective function, when a increases the optimal mass will be decreased and when a decreases the optimal mass will be increased. By choosing appropriate a and y values, minimum mass and minimum compliance with fully stressed constraint can be obtained. Optimally, following conditions should be satisfied. aVM() + 7yVL(x) = 0 (417) where VM(O) = M(O) = {A}r {(p= {A} (418) L(x) = [{} [K]{ }]= x}T [K]{x}+ {x}T [K] }+ T [K]a(x} a, a0, Qa, Qa, (419) [K]{x}= {F}, therefore, a ([K]{x) =0 (420) (421) a [ {x}= [K] [K]{x a, a, Substituting, Eq. 421 into Eq. 419, we obtain VL(x) ={x}T [K]{x} (422) Therefore, Eq. 417 can be written as aA, {X}T [C], X)}= 0 (423) where, [, [K] [BIT [D][B]d (424) awr, a[ , When we assume that stresses are uniform everywhere for the optimal solution and, then from Eq. 410 we have L(x) 2 2(1+v) (425) M(O) ).O 3E N) where, = Omax, Nis safety factor and a, is yield stress. N From Eq. 413 and Eq. 415 we got L(x) _{}T[K]{X} M() A 1=1 At the optimal solution, from Eq. 423, we get A = {T [C] {X a n (426) (427) (428) S{xr([cl){x}= Y{x[E]{x} a a Hence, Eq. 425 can be written as L(x) {x}T [K]{x 2(1+ v) XIT [E]X} 3E N Ca [KVk ,kiL^, where, [E]= [( 1=1 We obtain the ratio of a and y as a X}T EI(X) 2(1 +v) S(X)T [K]{X} 3E N If E = E0", [D] = [Dj0]" where 1 v [D]= v 1 1v 0 0 0 0 1 0 and [Do]= E2 v 1 1 2 O (429) (430) (431) 0 0 1v 2_ Therefore, for this study we consider different polynomial material property density relations of the form E = E0" (432) where, Eo is the Young's modulus of the material and q is the density and n is an integer. Differentiating the [D] matrix with respect to the design variable O[D] O[D]on 80 1 80 1 _= [Dl ,,Do' ,,=n[D]i = n [D] N (433) aQ) aQ), aQa ()Qaaq! Then Eq. 424 can be written as [Ci]= [BT [D] N [B]dV (434) The stiffness matrix of each element is calculated by [K ]= [B]T [D][B]dV (435) V From Eq. 427, [E] matrix of each element can be computed by pe npe npe [Ee]= [Cik = [B]T[D][Bn dV = [B]T[D][B] ,NdV (436) z=1 z=1 V V where, npe is the number of nodes per element. npe Since N, = 0, Eq. 435 and Eq. 436 can be combined to obtain 1=1 [Ee]= n[Ke] (437) We can rewrite Eq. 431 by this relation: ne ne {X}T [E]{X} = {Xe T [Ee ]{X } n {X} T [Ke ]{X} (438) e=l e=l where, ne is number of elements ne {x}T [K]{x} = X I { T[Ke { (439) e=l XT [E n (440) {X}T [K]{X} 27 From Eq. 431 and Eq. 440, we get a 2(1 + v) Cr = n (441) y 3E N Therefore, when we know the yield stress of the material we can get the ratio of a and y. To obtain a and y values, y can be arbitrarily selected and is set equal to 1 in this thesis. By using the aly ratio in Eq. 441, the stress in the optimal structure will be equal to if the stress distribution is uniform as assumed in the derivation. N CHAPTER 5 RESULTS Effect of Material Property Density Relationship Example 1 The feasible region shown in Figure 51, whose dimensions were 0.41x0.41, was divided into a uniform quadrilateral mesh containing 1764 nodes and 1681 elements. All nodes on the left side edge were fixed. Uniform pressure in the Xdirection was applied on six elements in the middle of the right side edge of the domain. The applied pressure is equal to the design stress (=200MPa) of the material. The Young's modulus of the material used in the structure is 2.07GPa. The Poisson's ratio is 0.29. Design 0.41m domain Pressure 0.41m Figure 51. Example 1 4.a704E I 2.* 5I Si Figure 52. Example 1, When n=3, The Ratio of ca/// is 498550332367 SIbBJRES lIi SEt 2.DCOT57 Figure 52. Example 1, When n=2, The Ratio ofa a/y is 498550332367 Figure 52 and Figure 53 show the optimal geometries and their stress distributions. The stresses are almost uniform everywhere, with values similar to those of the design (200MPa). Example 2 The rectangular feasible region, whose dimensions are 0.51(h)x0.26(w) was divided into a uniform quadrilateral mesh containing 1326 nodes and 1250 elements. All the nodes on the left side edge of the feasible were fixed. Uniform pressure was applied in the Ydirection to six elements along the right side edge. The Young's modulus of material used in the structure is 2.07GPa. The Poisson's ratio is 0.29. The applied pressure is equal to the shear yield stress (=100MPa). Design 0.51m domain Pressure domain 0.26m Figure 54. Example 5 Figure 55 and 56 show the optimal geometries and its stress distributions. In these two pictures, we can see the stress concentration in each corer. Except for those two regions, the stress is almost uniform everywhere else. Figure 55. Example 2, When n=2, The Ratio of a /7 is 332367 In I. 1 H . 1 EM ,, ,,0 . ii a10 Figure 55. Example 2, When n=2, The Ratio of al/y is 332367 II zt .1. . i .' Figure 56. Example 2, When n 3, The Ratio of aly is 498550 The range of the stresses is 1.6MPa 1.9MPa if the stress concentrations are ignored. The reason of the existence of stress concentration is that the optimal geometry must fit in the design domain. The high stress at the corners can be reduced if a larger feasible region was available. Example 3 The feasible region is Lshape block shown in Fig. 57. The feasible region was divided into a uniform quadrilateral mesh containing 781 nodes and 700 elements. All nodes on the right side edge were fixed. Uniform pressure was applied as shown in Fig. 57. The Young's modulus of material was assumed to be 2.07GPa and the Poisson's ratio equal to 0.29. The applied pressure is equal to the yield stress (=200MPa). Pressure O.m 0.3m T Design domain lm 0.3m m Figure 57. Example 3 In this Lshape structure, bending moment due to the applied load is too large to be supported by structure even if the geometry of the structure was identical to the feasible region. In addition, there is stress concentration at the inside corner of the shape. ,I Figure 58. Example 3, When n=2 Therefore, it is not possible to design a structure to carry this load that can fit in the feasible region. A valid design can be found only by enlarging the feasible domain. Effect of Design Space (Example 4) The design space in practical engineering design problems is often limited due to which the topology optimization process yields structures that have stress concentration at the boundary of the feasible region. This example investigates this by varying the height of the initial design domain.  O.1m Design Domain P 100MIPa 1 F Figure 59. Example 4 The Young's modulus of material used in the structure is 2.07GPa and the Poisson's ratio is 0.29. The design domain is shown in Figure 59. In Figure 510, h equals 0.1m, in Figure 511, h equals 0.2m and in Figure 512, h equals 0.3m. The ratio of a /y is 332367 IF eI Figure 510. Example 4, When the Height is 0. lm I PEl'SU lart .rtI 1.73172E T kl*X 6 % Tv; z 41' 'NE 00 0I~ Figure 511. Example 4, When the Height is 0.2m Figure 512. Example 4, When the Height is 0.3m When the height constraint is changed from 0.1m to 0.2m (Figure 510 and Figure 511), its maximum stress reduces from 7.586GPa to 2.4759GPa. In Figure 510, there are stress concentrations on each left side corner. In Figure 511, the height of the feasible region was doubled compared to that in Figure 510. This means that expanding the feasible region reduces the stress concentration. When the height increases from 0.2m to 0.3m, there is not significant difference. Effect of Safety Factor (Example 5) The rectangular feasible region, whose dimensions were 0.5(h)x0.3(w), was divided into a uniform quadrilateral mesh containing 1581 nodes and 1468 elements. All nodes on the left side edge were fixed. Uniform pressure was applied on the six elements at the bottom right side corner in the ydirection. The Young's modulus of the material used in the structure is 2.07GPa and the Poisson's ratio is 0.29. The applied pressure is equal to the yield stress (=200MPa). In Figure 514, safety factor is 1 and the ratio of a y is 332367. In Figure 515, safety factor is 1.5 and the ratio of a y is 147718. 11111 I Ii II II < 0.3m" Design Domain P=200MPa Figure 513. Example 5 Figure 514. Example 5, When n=2 and the Safety Factor is 1 Average stress = 221711304 Max Stress = 810938231 0.5m , .:j" C I Figure 515. Example 5, When n=2 and the Safety Factor is 1.5 Average stress = 158755368 Max Stress = 548166674 When the safety factor is increased from 1 to 1.5, average stress is reduced and its maximum stress is also reduced. Comparison with IDEAS Analysis In order to verify stress distribution in the optimal structure in example 5, IDEAS software was used to compare the results. The optimal shape computed in Figure 515 is used for the comparison. After creating the solid model in IDEAS that is almost identical to the result in Figure 515, FEA analysis was performed to obtain the stress distribution. Figure 516 shows the result of the analysis. There is a stress concentration on the left corner of the structure. The maximum stress (4.83e8 Pa) is greater than twice the yield stress (2.00e8 Pa). The topology optimization program is unable to eliminate this stress concentration due to the limited space available within the feasible domain. .. I Figure 516. Analysis of Example 5 Figure 517. Analysis of first modification of Example 5 To reduce the maximum stress, the designer has to modify the design further as shown in Fig. 517. In Figure 517, the maximum stress was reduced to 2.91e8 Pa by adding more material near the stress concentration. However, this maximum stress is still slightly above the yield stress and there is stress concentration also. Figure 518 shows the stress distribution after further modification. We can still see stress concentration on the right corner, but the value of stress of this area is almost same as the yield stress. Figure 518. Analysis of second modification of Example 5 ___________________________ II1 CHAPTER 6 CONCLUSION In sizing optimization, the design variables are the crosssections of the truss or the thickness of the plate etc. The geometry change is quite small when varying these design variables so it is not necessary to create a new analysis model each time the design variables are changed. This means that the topology of the structure remains fixed throughout the optimization procedure. For this reason, designers usually use sizing optimization after they have developed a detailed design that they wish to optimize. Shape and topology optimization is a powerful tool that can be used in the conceptual design phase. Since the entire geometry is considered as a variable, large structural changes are possible such as the creation of new holes and/or boundaries. Therefore, the optimal shape is independent of the initial guess. An implicit shape representation is used in this thesis where the contours of a shape function corresponding to a threshold value are treated as the boundaries of the shape. The shape density function is defined over a feasible region and is represented by a piecewise linear interpolation over quadrilateral finite element. The values of the density function at the nodes are used as the design variables. As the shape changes, the structural properties should change also. This means that the material properties should be related to the shape density function. In this thesis, the implementation of structural topology optimization under stress constraint is studied. A relation between the compliance and mass is developed for the optimal design assuming that the stress is uniform for the optimal structure. The design objective is to obtain the optimal topology of a structure, which maximizes the stiffness and minimizes the mass, such that the stress distribution over the design domain is uniform and equal to the desired design stress. An objectoriented topology optimization program was implemented using the C++ programming language. The program consists of two primary components, an optimization module and a finite element module. The optimization module was adapted from a previously developed program written by Dr. Ashok Kumar [Kum 93a]. To display the results of optimization program, a graphical display program was used that was originally developed by Wood [Woo 98]. This program was further modified to display the stress distribution for this thesis. The results in this thesis illustrate that the multiobjective optimization approach can be used to simultaneously minimize mass and maximize stiffness such that the stress constraints are satisfied everywhere except at stress concentration. This approach is good for conceptual design of structures. However, the geometry of the structure thus obtained has to be further modified by the designer to remove stress concentrations and to satisfy other nonstructural design requirements. APPENDIX COMPUTER SOURCE CODE FROM OPTIMIZATION PROGRAM The following source code is for calculating the objectives and their gradients. #include #include #include #include #include #include #include #include "cltopopt.h" #include "opt.h" #include "fe solver.h" #include "fe model.h" #include "interpolation.h" #include "element.h" #include "a type.h" #include "cl_4n_quad.h" double Cltopopt::ALPHA 0.; double Cltopopt::BETA 0.; double Cltopopt::GAMMA 0.; /* ========== Constructor for Topopt ============ */ Cltopopt::Cltopopt(int nnl,int necl,int nicl,double *xl,double *gl, double **hl,intnscl, double **scl, double *hcl,double **ghcl, double *gcl,double **ggcl,ofstream outl,FEModel *feml, Interpolation *clintegl): Opt(nnl,necl,nicl,xl,gl,hl,nscl,scl,hcl,ghc1,gc1,ggcl) { femfeml; clinteg clintegl; outstream outl; } /* =========== double eval function() =============== */ double Cltopopt::eval function() // Purpose : Returns the value of the objective function, for given x // Stores values of constraints in gc, for given x // Increment the counter num func calls++; double *value new double[l]; value; value[l] 0; int size1; double func, funcl 0.0; double func2 0.0; int tne fem>get tot num ele(); for(int il;i<tne;i++){ clinteg>integration(fem,object fnl,i,value,size); funcl+value[l]; } 43 cout << "mass = << func << endl; if(BETA!=0.0){ for(i=1;i<tne;i++){ cl integ>integration(fem,object fn2,i,value,size); func2+=value[1]; cout << "surface = << func2 << endl; //Compute compliance FE Solver *fes=fem>get fes(); fes>fill Kmax(); fes>assemble skylineK(); fes>assemble skylinef(); fes>skyline solver(1); double compliance = fes>compliance(); cout<"Current value of Compliance = " func = func + func2 + GAMMA*compliance; return func; }/* end eval function */ *==============void set final compliance()==================== * void Cltopopt::set abc(double a, double b, double c) ALPHA = a; BETA = b; GAMMA = c; } /*  double Cltopopt::get alpha() return ALPHA; } /* * double Cltopopt::get beta() return BETA; } /* ============void evalgrad hess()=======================* double *g mass; void Cltopopt::evalgrad hess() /* Purpose : evaluate the gradient & the Hessian of the /* of the objective fn. and gradient of the constraints /* Note: Call eval function before calling this function */ int ij,k; int nspn,tmp; nspn =0; int num eletyps=fem>get num eletyps(); Element **elenew Element *[numeletyps]; *ele[num eletyps]; for(i=l;i<=numeletyps;i++) ele[i]=fem>get ele(i); tmp= ele[i]>get_phi interpolation(>get num sfsper node(); if(tmp>nspn) nspn = tmp; int tot numnds=fem>get tot numnds(); int nn = nspn*totnumnds; //total number of design variables for (i=l;i<=nn;i++){ g[i]=0.0; } for (i=l;i<=num ele typs;i++) { 44 AType *atyp ele[i]>get atyp(); double t atyp>get thickness(i,0.5,0.5); //Thickness Interpolation *phiintep ele[i]>get_phiinterpolation(); int nnpe = ele[i]>get nds_perele(); int tnsfs = nnpe*nspn; //num of shape functions per element if(gmass =NULL) { gmass = new double[tnsfs]; g mass; for (int ill1;il<tnsfs;il++){ gmass[il] 0.0; } for (jl;j< ele[i]>get numele(;j++){ phiintep>integration(fem,gradobject fn,j,g mass,tnsfs); for (k 1;k< nnpe;k++){ for (int l 1;lnspn;l++) int kl = ele[i]>getelecon(j,k)*nspn nspn + 1; int k2 = k*nspn nspn + 1; g[kl]+=gmass[k2]*t; double *gradc; /** Note: The gradient of the function is stored as ggc in Opt ** gradc fem>gradcompliance(); for(il1;i< nn;i++){ g[i] += GAMMA*gradc[i]; } //eval_grad hess //===================object function============ int *con; void object fnl(FEModel *fem,int en,int numsfs_perele, double *nx,double **dn, double *value,int size){ value[l]=0.; Element *ele fem>getele(1); int nds_perele ele>get nds_per ele(); if(con =NULL) { con new int[nds_perele];con; for(intj 1;j< nds_per ele;j++){ con[j] ele>getele con(enj); } double *xfem>get_phi(); Interpolation *c linterp ele>get_phiinterpolation(); int nspn clinterp>get numsfs_per node(); //Add density to value[1] double ALPHA Cltopopt::get alpha(); for(int kl;k< nds_per ele;k++){ for(intkl 1; kl< nspn;kl++){ value[ ]+ ALPHA*x[con[k]*nspn(nspnkl)]*nx[(k 1)*nspn+kl]; S//object fnl void object fn2(FE Model *fem,int en,int num sfs_per ele, double *nx,double **dn, double *value,int size){ value[l]=0.; Element *ele fem>getele(1); int nds_perele ele>get nds_perele(); if(con =NULL) { con new int[nds_perele];con; for(intj 1;j< nds_per ele;j++){ con[j] ele>getele con(enj); } double *xfem>get_phi(); Interpolation *cl interp ele>get_phiinterpolation(); int nspn clinterp>get numsfs_per node(); //Add laplacian of density to value[1] double BETA Cltopopt::get beta(); for(int kl;k< nds_per ele;k++){ value[l]+ BETA*(x[con[k]*nspn3]*dn[(k1)*nspn+l][3]+ x[con[k]*nspn2]*dn[(k1)*nspn+2][3] + x[con[k]*nspnl]*dn[(kl)*nspn+3][3] + x[con[k]*nspn]*dn[(kl)*nspn+4][3] + x[con[k]*nspn3]*dn[(kl)*nspn+l][4] + x[con[k]*nspn2]*dn[(k1)*nspn+2][4] + x[con[k]*nspnl]*dn[(kl)*nspn+3][4] + x[con[k]*nspn]*dn[(kl)*nspn+4][4]); S//object fn2 /============= ======gradient object function============ void gradobject fn(FE Model *fem,int en,int nsfns, double *nx,double **dn, double *shape,int size){ double ALPHA Cltopopt::get alpha(); double BETA Cltopopt::get beta(); for(int i 1;i<nsfns;i++){ shape[i] ALPHA*nx[i]+BETA*(dn[i] [3]+dn[i] [4]); i LIST OF REFERENCES [All 92] Allaire, G. and Kohn, R.V., "Topology optimization and optimal shape design using homogenization," Topology Design of Structures, Ed. M. P. Bendsoe and Mota Soares, Kluwer Academic Publisher, New York, pp. 182194, 1992. [Ben 92] Bendsoe, M.P., Diaz, A. and Kikuchi, N., "Topology and generalized layout optimization of elastic structures," Topology Design of Structures, Ed. Bendsoe and Mota Soares, Kluwer Academic Publisher, New York, pp. 238250, 1992. [Ben 91] Bendsoe, M.P. and Rodrigues, H.C., "Integrated topology and boundary shape optimization of 2D solids," Computer Methods in Applied Mechanics and Engineering, vol. 87, pp. 1534, 1991. [Ben 88] Bendsoe, M.P. and Kikuchi, N., "Generating optimal topologies in structural design using a homogenization method," Computer Methods in Applied Mechanics and Engineering, vol. 71, pp. 197 224, 1988. [Ben 78] Bensoussan, A., Lions, J. and Papanicolaou, G., Asymptotic Analysis for Periodic Structures, NorthHolland Publishing Company, New York, 1978. [Bot 86] Botkin, M.E., Yang, R.J. and Bennett, J.A., "Shape optimization of threedimensional stamped and solid automotive components," The Optimum Shape, J.A. Bennett and M.E. Botkin Editors, Plenum Press, New York, 1986. [Dob 69] Dobbs, M.W. and Felton, L.P., "Optimization of truss geometry," Journal of Structural Engineering, vol. 95, pp. 21052118, 1969. [Dor 64] Dorn, W.S., Gormory, R.E. and Greenburg, H.J., "Automatic design of optimal structures," Journal de Mechanique, vol. 3, pp. 2552, 1964. [Esc 90] Eschenauer, H., Koski, J. and Osyczka, A., Multicriteria design optimization: Procedure and applications, Berlin: Springer, 1990. [Haf 92] [Haf 86] [Kir 81] [Koh 86] [Kum 00] Haftka, R.T. and Gurdal, Z., Elements of Structural Optimization, 3rd Edition, Kluwer Academic Publisher, New York, 1992. Haftka, R.T. and Grandhi, R.V., "Structural shape optimization  A survey," Computer Methods in Applied Mechanics and Engineering, vol. 57, pp. 91106, 1986. Kirsch, U., Optimum Structural Design: Concepts, Methods and Applications, McGrawHill Book Company, New York, 1981. Kohn, R.V. and Strang, G., "Optimal design and relaxation of variational problems," Communications in Pure and Applied Mathematics, vol. 39, pp. 113137, 139182, 333350, 1986. Kumar, A.V., and Yu, Lichao, "An objectoriented modular framework for implementing the finite element method," Computers and Structures, vol. 79, no.9, pp. 919928, 1998. Kumar, A.V. and Gossard, V, "Sequential approximation method for structural optimization using logarithmic barriers," Advances in Design Automation, American Society of Mechanical Engineering, Design Engineering Division (Publication) DE, vol. 65, pt 1, pp. 735742, m1993. Kumar, A.V., "Shape and topology synthesis of structures using a sequential optimization algorithm," Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA., 1993. Morris, V, "Foundations of Structural Optimization: A Unified Approach," John Wiley and Sons, 1982. Papalambros, Panos and Chirehdast, Hehran, "An integrated environment for structural configuration design," Journal of Engineering Design, vol. 1, pp. 7396, 1990. Stadler, W., "Multicriteria optimization in mechanics." Applied Mechanics Review, vol. 37, pp. 277286, 1984. Strang, Gilbert and Kohn, R.V., "Optimal design in elasticity and plasticity," International Journal for Numerical Methods in Engineering, vol.22, pp. 183188, 1986. Suzuki, K. and Kikuchi, N., "A homogenization method for shape and topology optimization," Computer Methods in Applied Mechanics and Engineering, vol. 93, pp. 291318, 1991. [Kum 93a] [Kum 93b] [Mor 82] [Pap 90] [Sta 84] [Str 86] [Suz 91] 48 [Woo 98] Wood, A. "Shape and topology optimization using a shape density function interpolated over quadrilateral elements," M.S. Thesis, University of Florida, Gainesville, FL, 1998. [Yan 86] Yang, R.J., Choi, K.K. and Haug, E.J., "Numerical considerations in structural component shape optimization," ASME Journal of Mechanics, Transmissions and Automation in Design, vol. 107, No. 3, pp. 334339, 1986. BIOGRAPHICAL SKETCH TaeJoong Yu was born on March 3, 1971, in Seoul, Korea. He grew up in an urban area. In spring 1991, he began an undergraduate program in the Mechanical Design and Production Engineering Department of KonKuk University. After two years of studying in the program, he enrolled in the military service from 1992 until 1993 to fulfill his mandatory service as a Korean citizen. After discharge, he returned to KonKuk University and graduated in spring 1998. He began his master's degree study in fall 1999 and from summer 2000 he studied with Dr. Ashok V. Kumar in the Mechanical Engineering Department of the University of Florida, Gainesville, Florida. 