Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0000831/00001
## Material Information- Title:
- Topology optimization under stress constraints
- Creator:
- Yu, Tae Joong (
*Dissertant*) Kumar, Ashok V. (*Thesis advisor*) - Place of Publication:
- Gainesville, Fla.
- Publisher:
- University of Florida
- Publication Date:
- 2003
- Copyright Date:
- 2003
- Language:
- English
## Subjects- Subjects / Keywords:
- Design engineering ( jstor )
Design optimization ( jstor ) Geometric shapes ( jstor ) Geometry ( jstor ) Material properties ( jstor ) Shape optimization ( jstor ) Stress concentration ( jstor ) Stress distribution ( jstor ) Structural design ( jstor ) Topology ( jstor ) Dissertations, Academic -- UF -- Mechanical Engineering Mechanical Engineering thesis, M.S City of Gainesville ( local ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Abstract:
- 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. ( , )
- Abstract:
- constraints, optimization, stress, topology
- General Note:
- Title from title page of source document.
- General Note:
- Includes vita.
- Thesis:
- Thesis (M.S.)--University of Florida, 2003.
- Bibliography:
- Includes bibliographical references.
- General Note:
- Text (Electronic thesis) in PDF format.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Yu, Tae Joong. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 9/9/1999
- Resource Identifier:
- 002903585 ( aleph )
53314675 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads:
yu_t ( .pdf )
yu_t_Page_01.txt yu_t_Page_57.txt yu_t_Page_35.txt yu_t_Page_58.txt yu_t_Page_37.txt yu_t_Page_30.txt yu_t_Page_51.txt yu_t_Page_28.txt yu_t_Page_36.txt yu_t_Page_09.txt yu_t_Page_20.txt yu_t_Page_02.txt yu_t_Page_40.txt yu_t_pdf.txt yu_t_Page_19.txt yu_t_Page_29.txt yu_t_Page_38.txt yu_t_Page_52.txt yu_t_Page_26.txt yu_t_Page_49.txt yu_t_Page_04.txt yu_t_Page_50.txt yu_t_Page_47.txt yu_t_Page_43.txt yu_t_Page_21.txt yu_t_Page_18.txt yu_t_Page_23.txt yu_t_Page_17.txt yu_t_Page_16.txt yu_t_Page_33.txt yu_t_Page_41.txt yu_t_Page_32.txt yu_t_Page_31.txt yu_t_Page_55.txt yu_t_Page_25.txt yu_t_Page_42.txt yu_t_Page_53.txt yu_t_Page_14.txt yu_t_Page_06.txt yu_t_Page_15.txt yu_t_Page_13.txt yu_t_Page_12.txt yu_t_Page_44.txt yu_t_Page_46.txt yu_t_Page_45.txt yu_t_Page_08.txt yu_t_Page_48.txt yu_t_Page_05.txt yu_t_Page_10.txt yu_t_Page_39.txt yu_t_Page_54.txt yu_t_Page_11.txt yu_t_Page_56.txt yu_t_Page_34.txt yu_t_Page_07.txt yu_t_Page_27.txt yu_t_Page_22.txt yu_t_Page_03.txt yu_t_Page_24.txt |

Full Text |

TOPOLOGY OPTIMIZATION UNDER STRESS CONSTRAINTS By TAE-JOONG 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 Tae-Joong 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: Young-Nam 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 I-D 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 2-1 Structural Optimization of sizing, shape and topology .............................................6 2-2 Homogrnization....................................................... ........8 2-3 Material Property as a Function of Density.................................... ......... ............9 3-1 Shape R presentation ...................... .. .......................... .. .. ...... ........... ... 13 3-2 Solution Flow chart ................................................ .. ............. ......... 19 5 1 E x a m p le 1 ............................................................................................................. 2 8 5-2 Example 1, When n=2, The Ratio of a/y is 332367......... ....... ........ ............... 29 5-3 Example 1, When n=3, The Ratio of a/y is 498550...............................................29 5 -4 E x am p le 5 .......................................................................... 3 0 5-5 Example 2, When n=2, The Ratio of a/y is 332367..................................................31 5-6 Example 2, When n=3, The Ratio of a/y is 498550..................................................31 5 -7 E x am p le 3 .......................................................................... 3 2 5-8 E xam ple 3, W hen n= 2 ....................................................................... ................... 33 5 -9 E x a m p le 4 ............................................................................................................. 3 3 5-10 Example 4, When the Height is 0.1m......................................................... 34 5-11 Example 4, When the Height is 0.2m .............. ............................. 34 5-12 Exam ple 4, W hen the H eight is 0.3m .................................... ........ ............... 35 5 -13 E x am p le 5 ......................................................................... 3 6 5-14 Example 5, When n=2 and the Safety Factor is 1 .................................................36 5-15 Example 5, When n=2 and the Safety Factor is 1.5 .............................................37 vi 5-16 A naly sis of E xam ple 5 ....................................................................... ..................38 5-17 Analysis of first modification of Example 5 ................................ ..... ..........38 5-18 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 Tae-Joong 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 trial-and-error process cannot be avoided if the designer really wants to select the minimum-weight 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, so-called 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, cross-sectional area, etc. It does not involve the redefinition of shapes of outer boundaries and inner holes of structure. For this reason, if a sub-optimal 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 user-defined 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 sub-optimal topology is chosen when formulating the optimization problem, the resulting structure will also be sub-optimal. 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 under-stressed 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 non-convex 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 all-or-nothing 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 2-1. 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 multiple-loading cases and three-dimensional 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 all-or-nothing approach). Areas void of material shows either holes, or areas outside of the structure's outermost boundary. During this optimization procedure, under-stressed 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 well-posed 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 2-2. 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 2-3. D 1-a2 Figure 2-3. Material Property as a Function of Density Figure 2-3 shows an example in which the hole within each unit cell is assumed to be square (a = b), so that 1-a2 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 low-density material. On the contrary, the hole size is decreased in areas where material has high stresses, resulting in high-density. 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 2-4(a) shows a microstructure consisting of unit cells with rectangular holes and Figure 2-4(b) illustrates a rank-2 microstructure. The rank-2 consists of alternating layers of stiff material and rank-1 material. Rank-1 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 rank-2 layer yields the stiffest composite. The optimal shape, which is solved for using rank-2 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 2-4. Microstructures. A)By Microscale Rectangular, B)By Rank-2 Papalambros et al. [Pap 90] used image-processing 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 two-step 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 (3-1) Figure 3-1 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 3-1. 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 4-node 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) (3-2) where, are the nodal density values and N, (s, t) are the isoparametric shape functions for the 4-node 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) (3-3) 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 3-1, 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 multi-objective problem may be written as Minimize : aM()) + 7L(x(o)) (3-4) where a and y are user specified control parameters. The mass, M(q!), can be expressed as M() = qd = {AT {(p (3-5) 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 (3-6) 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) (3-7) 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 3-4 through 3-7 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~ (3-8) 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 four-node isoparametric quadrilateral elements for the analysis. Therefore, the displacements are represented by piece-wise bilinear interpolation within each element. This representation yields the following strain-displacement relation for the element {e}= [B]{U} (3-9) 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 (3-7), for each element, we get Le(u) =f{&se} [D( e)]{g =}dQ {&}Q [J[B]Y [D(I)B]d[ e ]{4} (3-10) 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 (3-11) 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. (3-11) 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 2m-1. 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 3-2. 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 3-2 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 (4-1) 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 (4-2) v -v 1 where v is Poison's ratio, E is Young's modulus, and aT = {o- 2 o-3 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 v-on = 72V( -' 2 +2 (2 -3 2 +(3 -1 )2 (4-3) 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 -- (4-5) 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 (4-6) E 3 E Equation 4-6 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 (4-7) 2(1 + v) f2dA: 3E fTC 3E oa dA< a ciTCc ~ L(E) (4-8) O 2(1 + v), 2(1 + v)L() (4-8) 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. 4-8 can be simplified as 2 3E -max A < L(E) (4-9) 2(1 + v) where A is the area (2D) of Q. We obtain a lower bound on the mean compliance by -L() > 2 ) A (4-10) 3E Problem Formulation Our objective is minimizing the mass and the mean compliance with stress constraints. That is, minimize AM(q) + yL(x(q)) (4-11) 0< 0 1 (4-12) 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} (4-13) Using the strain-displacement relation, {(8= [B]{x},assumed in FEA model, the mean compliance equation Eq. 4-1 can be written as L(x) = {x [B]T [D][B]{x}d (4-14) 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 (4-15) L(x) = {x} [K]{x}, where [K] is the stiffness matrix. [K] [B]T [D][B]jd (4-16) 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 (4-17) where VM(O) = -M(O) = -{A}r {(p= {A} (4-18) L(x) = [{} [K]{ }]= x}T [K]{x}+ {x}T [K] }+ T [K]a(x} a, a0, Qa, Qa, (4-19) [K]{x}= {F}, therefore, a ([K]{x) =0 (4-20) (4-21) a [ {x}= -[K]- [K]{x a, a, Substituting, Eq. 4-21 into Eq. 4-19, we obtain VL(x) =-{x}T [K]{x} (4-22) Therefore, Eq. 4-17 can be written as aA, {X}T [C], X)}= 0 (4-23) where, [, [K] [BIT [D][B]d (4-24) awr, a[ , When we assume that stresses are uniform everywhere for the optimal solution and, then from Eq. 4-10 we have L(x) 2 2(1+v) (4-25) M(O) ).O 3E N) where, = Omax, Nis safety factor and a, is yield stress. N From Eq. 4-13 and Eq. 4-15 we got L(x) _{}T[K]{X} M() -A 1=1 At the optimal solution, from Eq. 4-23, we get A = {T [C] {X a n (4-26) (4-27) (4-28) S{xr([cl){x}= Y{x[E]{x} a a Hence, Eq. 4-25 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 1-v 0 0 0 0 1 0 and [Do]= E-2 v 1- 1 2 O (4-29) (4-30) (4-31) 0 0 1v 2_ Therefore, for this study we consider different polynomial material property density relations of the form E = E0" (4-32) 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 (4-33) aQ) aQ), aQa ()Qaaq! Then Eq. 4-24 can be written as [Ci]= [BT [D] N [B]dV (4-34) The stiffness matrix of each element is calculated by [K ]= [B]T [D][B]dV (4-35) V From Eq. 4-27, [E] matrix of each element can be computed by pe npe npe [Ee]= [Cik = [B]T[D][Bn dV = [B]T[D][B] ,NdV (4-36) z=1 z=1 V V where, npe is the number of nodes per element. npe Since -N, = 0, Eq. 4-35 and Eq. 4-36 can be combined to obtain 1=1 [Ee]= n[Ke] (4-37) We can rewrite Eq. 4-31 by this relation: ne ne {X}T [E]{X} = {Xe T [Ee ]{X } n {X} T [Ke ]{X} (4-38) e=l e=l where, ne is number of elements ne {x}T [K]{x} = X I { T[Ke { (4-39) e=l XT [E -n (4-40) {X}T [K]{X} 27 From Eq. 4-31 and Eq. 4-40, we get a 2(1 + v) Cr = n (4-41) 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. 4-41, 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 5-1, 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 X-direction 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 5-1. Example 1 4.a704E I 2.* 5I Si Figure 5-2. Example 1, When n=3, The Ratio of ca/// is 498550332367 SIbBJRES lIi SEt 2.DCOT57 Figure 5-2. Example 1, When n=2, The Ratio ofa a/y is 498550332367 Figure 5-2 and Figure 5-3 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 Y-direction 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 5-4. Example 5 Figure 5-5 and 5-6 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 5-5. Example 2, When n=2, The Ratio of a /7 is 332367 -In I. 1 H .- 1 EM -,, --,,0 .- i-i a10 Figure 5-5. Example 2, When n=2, The Ratio of al/y is 332367 II zt .1. . i .' Figure 5-6. 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 L-shape block shown in Fig. 5-7. 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. 5-7. 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 5-7. Example 3 In this L-shape 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 5-8. 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 5-9. 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 5-9. In Figure 5-10, h equals 0.1m, in Figure 5-11, h equals 0.2m and in Figure 5-12, h equals 0.3m. The ratio of a /y is 332367 IF eI Figure 5-10. Example 4, When the Height is 0. lm I PEl'SU lart .rt-I 1.73172E T kl*X 6 % -Tv; z 41' 'NE 00 0I~- Figure 5-11. Example 4, When the Height is 0.2m Figure 5-12. Example 4, When the Height is 0.3m When the height constraint is changed from 0.1m to 0.2m (Figure 5-10 and Figure 5-11), its maximum stress reduces from 7.586GPa to 2.4759GPa. In Figure 5-10, there are stress concentrations on each left side corner. In Figure 5-11, the height of the feasible region was doubled compared to that in Figure 5-10. 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 y-direction. 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 5-14, safety factor is 1 and the ratio of a y is 332367. In Figure 5-15, 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 5-13. Example 5 Figure 5-14. Example 5, When n=2 and the Safety Factor is 1 Average stress = 221711304 Max Stress = 810938231 0.5m , .:j" C I Figure 5-15. 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 I-DEAS Analysis In order to verify stress distribution in the optimal structure in example 5, I-DEAS software was used to compare the results. The optimal shape computed in Figure 5-15 is used for the comparison. After creating the solid model in I-DEAS that is almost identical to the result in Figure 5-15, FEA analysis was performed to obtain the stress distribution. Figure 5-16 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 5-16. Analysis of Example 5 Figure 5-17. Analysis of first modification of Example 5 To reduce the maximum stress, the designer has to modify the design further as shown in Fig. 5-17. In Figure 5-17, 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 5-18 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 5-18. Analysis of second modification of Example 5 ___________________________ II1 CHAPTER 6 CONCLUSION In sizing optimization, the design variables are the cross-sections 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 piece-wise 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 object-oriented 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 non-structural 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) { fem-feml; 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 size-1; double func, funcl 0.0; double func2 0.0; int tne fem->get tot num ele(); for(int i-l;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 **ele-new 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;l-nspn;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 *x-fem->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-(nspn-kl)]*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 *x-fem->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]*nspn-3]*dn[(k-1)*nspn+l][3]+ x[con[k]*nspn-2]*dn[(k-1)*nspn+2][3] + x[con[k]*nspn-l]*dn[(k-l)*nspn+3][3] + x[con[k]*nspn]*dn[(k-l)*nspn+4][3] + x[con[k]*nspn-3]*dn[(k-l)*nspn+l][4] + x[con[k]*nspn-2]*dn[(k-1)*nspn+2][4] + x[con[k]*nspn-l]*dn[(k-l)*nspn+3][4] + x[con[k]*nspn]*dn[(k-l)*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. 182-194, 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. 238-250, 1992. [Ben 91] Bendsoe, M.P. and Rodrigues, H.C., "Integrated topology and boundary shape optimization of 2-D solids," Computer Methods in Applied Mechanics and Engineering, vol. 87, pp. 15-34, 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, North-Holland Publishing Company, New York, 1978. [Bot 86] Botkin, M.E., Yang, R.J. and Bennett, J.A., "Shape optimization of three-dimensional 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. 2105-2118, 1969. [Dor 64] Dorn, W.S., Gormory, R.E. and Greenburg, H.J., "Automatic design of optimal structures," Journal de Mechanique, vol. 3, pp. 25-52, 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. 91-106, 1986. Kirsch, U., Optimum Structural Design: Concepts, Methods and Applications, McGraw-Hill 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. 113-137, 139-182, 333-350, 1986. Kumar, A.V., and Yu, Lichao, "An object-oriented modular framework for implementing the finite element method," Computers and Structures, vol. 79, no.9, pp. 919-928, 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. 735-742, 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. 73-96, 1990. Stadler, W., "Multicriteria optimization in mechanics." Applied Mechanics Review, vol. 37, pp. 277-286, 1984. Strang, Gilbert and Kohn, R.V., "Optimal design in elasticity and plasticity," International Journal for Numerical Methods in Engineering, vol.22, pp. 183-188, 1986. Suzuki, K. and Kikuchi, N., "A homogenization method for shape and topology optimization," Computer Methods in Applied Mechanics and Engineering, vol. 93, pp. 291-318, 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. 334-339, 1986. BIOGRAPHICAL SKETCH Tae-Joong 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 Kon-Kuk 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 Kon-Kuk 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. |

Full Text |

PAGE 1 TOPOLOGY OPTIMIZATION UNDER STRESS CONSTRAINTS By TAE-JOONG 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 PAGE 2 Copyright 2003 by Tae-Joong Yu PAGE 3 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: Young-Nam 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. iii PAGE 4 TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iii LIST OF FIGURES...........................................................................................................vi ABSTRACT.....................................................................................................................viii CHAPTER 1 INTRODUCTION............................................................................................................1 Structural Optimization...................................................................................................2 Objective and Motivation...............................................................................................3 2 STRUCTURAL OPTIMIZATION...................................................................................4 Sizing Optimization........................................................................................................4 Shape Optimization.........................................................................................................4 Topology Optimization...................................................................................................6 Homogenization Method................................................................................................7 3 TOPOLOGY OPTIMIZATION using SHAPE DENSITY FUNCTION.......................12 Shape Density Function................................................................................................12 Objective Function........................................................................................................15 Relationship between Material Properties and Density................................................16 Implementation by Finite Elements..............................................................................17 Solution Procedure........................................................................................................18 4 STRESS CONSTRAINTS..............................................................................................21 Mean Compliance.........................................................................................................21 Problem Formulation....................................................................................................23 iv PAGE 5 5 RESULTS.......................................................................................................................28 Effect of Material Property Density Relationship........................................................28 Example 1..............................................................................................................28 Example 2..............................................................................................................30 Example 3..............................................................................................................32 Effect of Design Space (Example 4).............................................................................33 Effect of Safety Factor (Example 5).............................................................................35 Comparison with I-DEAS Analysis..............................................................................37 6 CONCLUSION...............................................................................................................40 APPENDIX COMPUTER SOURCE CODE FROM OPTIMIZATION PROGRAM.....42 LIST OF REFERENCES...................................................................................................46 BIOGRAPHICAL SKETCH.............................................................................................49 v PAGE 6 LIST OF FIGURES Figure page 2-1 Structural Optimization of sizing, shape and topology................................................6 2-2 Homogrnization............................................................................................................8 2-3 Material Property as a Function of Density..................................................................9 3-1 Shape Representation.................................................................................................13 3-2 Solution Flowchart.....................................................................................................19 5-1 Example 1...................................................................................................................28 5-2 Example 1, When n=2, The Ratio of is 332367...................................................29 5-3 Example 1, When n=3, The Ratio of is 498550..................................................29 5-4 Example 5...................................................................................................................30 5-5 Example 2, When n=2, The Ratio of is 332367..................................................31 5-6 Example 2, When n=3, The Ratio of is 498550..................................................31 5-7 Example 3...................................................................................................................32 5-8 Example 3, When n=2................................................................................................33 5-9 Example 4...................................................................................................................33 5-10 Example 4, When the Height is 0.1m.......................................................................34 5-11 Example 4, When the Height is 0.2m.......................................................................34 5-12 Example 4, When the Height is 0.3m.......................................................................35 5-13 Example 5.................................................................................................................36 5-14 Example 5, When n=2 and the Safety Factor is 1....................................................36 5-15 Example 5, When n=2 and the Safety Factor is 1.5.................................................37 vi PAGE 7 5-16 Analysis of Example 5...............................................................................................38 5-17 Analysis of first modification of Example 5.............................................................38 5-18 Analysis of second modification of Example 5.........................................................39 vii PAGE 8 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 Tae-Joong 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 viii PAGE 9 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. ix PAGE 10 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 trial-and-error process cannot be avoided if the designer really wants to select the minimum-weight 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, 1 PAGE 11 2 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, so-called 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, cross-sectional area, etc. It does not involve the redefinition of shapes of outer boundaries and inner holes of structure. For this reason, if a sub-optimal 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 PAGE 12 3 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 user-defined 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. PAGE 13 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 sub-optimal topology is chosen when formulating the optimization problem, the resulting structure will also be sub-optimal. 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 4 PAGE 14 5 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. PAGE 15 6 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 under-stressed 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 non-convex 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 all-or-nothing dichotomy between material and holes. Bendsoe and Kikuchi [Ben 88] proposed the homogenization method to determine the properties of these composite materials. Figure 2-1. Structural Optimization of sizing, shape and topology PAGE 16 7 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 multiple-loading cases and three-dimensional 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 T (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 all-or-nothing approach). Areas void of material shows either holes, or areas outside of the structures outermost boundary. During this optimization procedure, under-stressed 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 well-posed and have better convergence performance. Kohn and Strang showed that relaxing the variational problem is identical to allowing composite materials for the PAGE 17 8 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. Figure 2-2. Homogrnization. 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 a i and b i where i is PAGE 18 9 the element number. These hole dimensions along with the angular orientation of the unit cells are the design variables. Initially, all elements are assumed to have the same hole dimensions and orientation giving the material a uniform porosity. i 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 2-3. Figure 2-3. Material Property as a Function of Density Figure 2-3 shows an example in which the hole within each unit cell is assumed to be square (a = b), so that 1-a 2 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 D ij A finite element analysis over the unit cell is required to calculate the material property coefficients for each value of hole size. PAGE 19 10 During the optimization procedure, the hole size is increased in areas where the stresses are low, which results in low-density material. On the contrary, the hole size is decreased in areas where material has high stresses, resulting in high-density. 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 2-4(a) shows a microstructure consisting of unit cells with rectangular holes and Figure 2-4(b) illustrates a rank-2 microstructure. The rank-2 consists of alternating layers of stiff material and rank-1 material. Rank-1 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 rank-2 layer yields the stiffest composite. The optimal shape, which is solved for using rank-2 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. PAGE 20 11 Figure 2-4. Microstructures. A)By Microscale Rectangular, B)By Rank-2 Papalambros et al. [Pap 90] used image-processing 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 two-step process has been used by Bendsoe [Ben 91] to obtain combined topology and shape optimization. PAGE 21 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 )(x which we refer to as the shape density function. The shape density function has values in the range from a threshold value T to 1. The contours of this function corresponding to the threshold value T of density express the boundaries of the geometry. The density field is defined over a feasible domain 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 0 0 x for which 1)(x T (3-1) Figure 3-1 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 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 12 PAGE 22 13 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 3-1. 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 4-node quadrilateral elements, the density function within each element can be expressed as, (3-2) where, i are the nodal density values and are the isoparametric shape functions for the 4-node quadrilateral element expressed in terms of the parametric coordinates s and t [Bathe 87]. The shape function can be expressed as ),(tsNi ),(),(),(),(),(44332211tsNtsNtsNtsNts PAGE 23 14 )1)(1(411tsN )1)(1(412tsN )1)(1(413tsN )1)(1(414tsN For an isoparametric element, the mapping between the parametric space and the real coordinates (x, y) is defined by 41),(iiitsNxx and (3-3) 41),(iiitsNyy 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 C 0 continuous shape density function ensures C 0 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 C 0 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 PAGE 24 15 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 3-1, 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 multi-objective problem may be written as Minimize : ))(()( xLM (3-4) 0 1 where and are user specified control parameters. The mass, ),( M can be expressed as ATdM)( (3-5) PAGE 25 16 where and dNA is nodal density values The mean compliance, is written by )(xL ddL)()())((xtxfx (3-6) 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, )()(00xDLdT (3-7) and are the strain and virtual strain in the structure caused by the displacement x and the virtual displacement x respectively. )( D 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 ( changes) the stiffness, the compliance and the deformation of the structure for the given load also change. The optimization problem given by the equation 3-4 through 3-7 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 PAGE 26 17 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 nEE0 (3-8) where, is the Youngs modulus of the material of the structure, 0E 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 four-node isoparametric quadrilateral elements for the analysis. Therefore, the displacements are represented by piece-wise bilinear interpolation within each element. This representation yields the following strain-displacement relation for the element eheuB (3-9) where, yNyNyNyNxNxNxNxN43214321B and ehu is the displacement vector corresponding to the quadrilateral finite element. Substituting this relation in (3-7), for each element, we get eheTTeheeTeeududLeeBDBDu (3-10) PAGE 27 18 eheTeheuuLKu where, eK is the element stiffness and ehu is the virtual displacement vector. The element stiffness matrix eK can be expressed as eTedeBDBK (3-11) 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. (3-11) 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 BDBT matrix must be determined. Since density is a linear function in the parametric coordinate s and t, the terms in the matrix 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 2m-1. BDBT 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 PAGE 28 19 that produce a structure having minimal compliance and minimal mass while simultaneously satisfying the stress constraint. Figure 3-2. 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 PAGE 29 20 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 3-2 illustrates a flowchart showing the solution procedure. PAGE 30 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 ddLCDTT)( (4-1) 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 1111EC (4-2) where is Poisons ratio, E is Youngs modulus, and 321T,,21 is the vector of the three principal stresses. In terms of the principal stresses and 3 the von Mises stress can be calculated as 21323222121von (4-3) A constraint on the allowable values of the von Mises stress in can be expressed as 2max2MTvon (4-4) where max is given maximum allowable stress and 21 PAGE 31 22 121212112121211M (4-5) In the case of plane stress, the following relation can be obtained [Ben 88] MCMTTT)1(213)1(21 EE (4-6) Equation 4-6 shows that the von Mises stress is bounded by a multiple of the square root of strain energy density at an arbitrary point in i.e., CT)1(23Evon or (4-7) )(123123CTLvEvEdA2von (4-8) 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. 4-8 can be simplified as )()1(23A2maxLE (4-9) where is the area (2D) of We obtain a lower bound on the mean compliance by A A3)1(2)(2max EL (4-10) PAGE 32 23 Problem Formulation Our objective is minimizing the mass and the mean compliance with stress constraints. That is, minimize ))(()( xLM (4-11) (4-12) 10 where and are user specified control parameters. ),( M the mass can be expressed as ATdM)( ,where dNA (4-13) Using the strain-displacement relation, xB ,assumed in FEA model, the mean compliance equation Eq. 4-1 can be written as dLxBDBxxTT)( (4-14) where, is the nodal displacement vector. x We can rewrite this equation as xBDBxxBDBxxTTTTddL)( (4-15) xKxxT)(L where is the stiffness matrix. K dBDBKT (4-16) In the objective function, when increases the optimal mass will be decreased and when decreases the optimal mass will be increased. By choosing appropriate and values, minimum mass and minimum compliance with fully stressed constraint can be obtained. Optimally, following conditions should be satisfied. 0)()( xLM (4-17) PAGE 33 24 where AATiiMM)()( (4-18) iiiiLxKxxKxxKxxKxxTTTT)( (4-19) FxK therefore, 0xKi (4-20) xKKxxKxKiiii10 (4-21) Substituting, Eq. 4-21 into Eq. 4-19, we obtain xKxxTiL)( (4-22) Therefore, Eq. 4-17 can be written as 0AXCXTii (4-23) where, diiiBDBKCT (4-24) When we assume that stresses are uniform everywhere for the optimal solution and, then from Eq. 4-10 we have 2*3)1(2)()(*NEMLyx (4-25) where, max Ny N is safety factor and y is yield stress. PAGE 34 25 From Eq. 4-13 and Eq. 4-15 we got niiiML1A)()(XKXxT (4-26) At the optimal solution, from Eq. 4-23, we get XCXTii A (4-27) XEXXCXTiiTniii1A (4-28) Hence, Eq. 4-25 can be written as 23)1(2)()(NEvMLyXEXXKXxTT (4-29) where, iniiiniKCEi11 (4-30) We obtain the ratio of and as 23)1(2NEvyXKXXEXTT (4-31) If where nEE0 n0DD 2100010112ED and 210001011200ED Therefore, for this study we consider different polynomial material property density relations of the form nEE0 (4-32) PAGE 35 26 where, is the Youngs modulus of the material and 0E is the density and n is an integer. Differentiating the [D] matrix with respect to the design variable iN iinininnn111DDDD00 D (4-33) Then Eq. 4-24 can be written as dVNniVBDBCTi (4-34) The stiffness matrix of each element is calculated by dVVBDBKTe (4-35) From Eq. 4-27, [E] matrix of each element can be computed by dVNnVdNnnpeiiiVnpeiiiVinpei111BDBBDBCETTie (4-36) where, npe is the number of nodes per element. Since Eq. 4-35 and Eq. 4-36 can be combined to obtain npeiiiN1 eeKEn (4-37) We can rewrite Eq. 4-31 by this relation: neeeeneeeeen11XKXXEXXEXeTTT (4-38) where, ne is number of elements neeee1XKXXKXeTT (4-39) nXKXXEXTT (4-40) PAGE 36 27 From Eq. 4-31 and Eq. 4-40, we get 23)1(2NEny (4-41) Therefore, when we know the yield stress of the material we can get the ratio of and To obtain and values, can be arbitrarily selected and is set equal to 1 in this thesis. By using the / ratio in Eq. 4-41, the stress in the optimal structure will be equal to if the stress distribution is uniform as assumed in the derivation. Ny PAGE 37 CHAPTER 5 RESULTS Effect of Material Property Density Relationship Example 1 The feasible region shown in Figure 5-1, 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 X-direction 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 Youngs modulus of the material used in the structure is 2.07GPa. The Poissons ratio is 0.29. Figure 5-1. Example 1 28 PAGE 38 29 Figure 5-2. Example 1, When n=2, The Ratio of / is 332367 Figure 5-3. Example 1, When n=3, The Ratio of / is 498550 PAGE 39 30 Figure 5-2 and Figure 5-3 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 Y-direction to six elements along the right side edge. The Youngs modulus of material used in the structure is 2.07GPa. The Poissons ratio is 0.29. The applied pressure is equal to the shear yield stress (=100MPa). Figure 5-4. Example 5 Figure 5-5 and 5-6 show the optimal geometries and its stress distributions. In these two pictures, we can see the stress concentration in each corner. Except for those two regions, the stress is almost uniform everywhere else. PAGE 40 31 Figure 5-5. Example 2, When n=2, The Ratio of / is 332367 Figure 5-6. Example 2, When n=3, The Ratio of / is 498550 PAGE 41 32 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 L-shape block shown in Fig. 5-7. 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. 5-7. The Youngs modulus of material was assumed to be 2.07GPa and the Poissons ratio equal to 0.29. The applied pressure is equal to the yield stress (=200MPa). Figure 5-7. Example 3 In this L-shape 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. PAGE 42 33 Figure 5-8. 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. Figure 5-9. Example 4 PAGE 43 34 The Youngs modulus of material used in the structure is 2.07GPa and the Poissons ratio is 0.29. The design domain is shown in Figure 5-9. In Figure 5-10, h equals 0.1m, in Figure 5-11, h equals 0.2m and in Figure 5-12, h equals 0.3m. The ratio of / is 332367 Figure 5-10. Example 4, When the Height is 0.1m Figure 5-11. Example 4, When the Height is 0.2m PAGE 44 35 Figure 5-12. Example 4, When the Height is 0.3m When the height constraint is changed from 0.1m to 0.2m (Figure 5-10 and Figure 5-11), its maximum stress reduces from 7.586GPa to 2.4759GPa. In Figure 5-10, there are stress concentrations on each left side corner. In Figure 5-11, the height of the feasible region was doubled compared to that in Figure 5-10. 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 y-direction. The Youngs modulus of the material used in the structure is 2.07GPa and the Poissons ratio is 0.29. The applied pressure is equal to the yield stress (=200MPa). In Figure 5-14, safety factor is 1 and the ratio of / is 332367. In Figure 5-15, safety factor is 1.5 and the ratio of / is 147718. PAGE 45 36 Figure 5-13. Example 5 Figure 5-14. Example 5, When n=2 and the Safety Factor is 1 Average stress = 221711304 Max Stress = 810938231 PAGE 46 37 Figure 5-15. 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 I-DEAS Analysis In order to verify stress distribution in the optimal structure in example 5, I-DEAS software was used to compare the results. The optimal shape computed in Figure 5-15 is used for the comparison. After creating the solid model in I-DEAS that is almost identical to the result in Figure 5-15, FEA analysis was performed to obtain the stress distribution. Figure 5-16 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. PAGE 47 38 Figure 5-16. Analysis of Example 5 Figure 5-17. Analysis of first modification of Example 5 PAGE 48 39 To reduce the maximum stress, the designer has to modify the design further as shown in Fig. 5-17. In Figure 5-17, 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 5-18 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 5-18. Analysis of second modification of Example 5 PAGE 49 CHAPTER 6 CONCLUSION In sizing optimization, the design variables are the cross-sections 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 piece-wise 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 40 PAGE 50 41 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 object-oriented 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 non-structural design requirements. PAGE 51 APPENDIX COMPUTER SOURCE CODE FROM OPTIMIZATION PROGRAM The following source code is for calculating the objectives and their gradients. #include PAGE 52 43 cout << "mass = << func1 << endl; if (BETA!=0.0) { for(i=1;i<=tne;i++){ c1_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 = "< PAGE 53 44 A_Type *atyp=ele[i]->get_atyp(); double t=atyp->get_thickness(i,0.5,0.5); // Thickness Interpolation *phi_intep=ele[i]->get_phi_interpolation(); int nnpe = ele[i]->get_nds_per_ele(); int tnsfs = nnpe*nspn; // num of shape functions per element if(g_mass==NULL) { g_mass = new double[tnsfs]; g_mass--; } for (int i1=1;i1<=tnsfs;i1++){ g_mass[i1]=0.0; } for (j=1;j<=ele[i]->get_num_ele();j++){ phi_intep->integration(fem,grad_object_fn,j,g_mass,tnsfs); for (k=1;k<=nnpe;k++){ for (int l=1;l<=nspn;l++) { int k1 = ele[i]->get_ele_con(j,k)*nspn nspn + l ; int k2 = k*nspn nspn + l ; g[k1]+=g_mass[k2]*t; } } } } double *grad_c; /** Note: The gradient of the function is stored as ggc in Opt **/ grad_c=fem->grad_compliance(); for(i=1;i<=nn;i++){ g[i] += GAMMA*grad_c[i]; } } // eval_grad_hess //====================object function===================== int *con; void object_fn1(FE_Model *fem,int en,int num_sfs_per_ele, double *nx,double **dn, double *value,int size){ value[1]=0.; Element *ele=fem->get_ele(1); int nds_per_ele=ele->get_nds_per_ele(); if(con==NULL) { con=new int[nds_per_ele];con--; } for(int j=1;j<=nds_per_ele;j++){ con[j]=ele->get_ele_con(en,j); } double *x=fem->get_phi(); Interpolation *c1_interp=ele->get_phi_interpolation(); int nspn=c1_interp->get_num_sfs_per_node(); // Add density to value[1] double ALPHA = C1topopt::get_alpha() ; for(int k=1;k<=nds_per_ele;k++){ for(int k1=1; k1<=nspn;k1++){ value[1]+=ALPHA*x[con[k]*nspn-(nspn-k1)]*nx[(k-1)*nspn+k1]; } } } //object_fn1 void object_fn2(FE_Model *fem,int en,int num_sfs_per_ele, PAGE 54 45 double *nx,double **dn, double *value,int size){ value[1]=0.; Element *ele=fem->get_ele(1); int nds_per_ele=ele->get_nds_per_ele(); if(con==NULL) { con=new int[nds_per_ele];con--; } for(int j=1;j<=nds_per_ele;j++){ con[j]=ele->get_ele_con(en,j); } double *x=fem->get_phi(); Interpolation *c1_interp=ele->get_phi_interpolation(); int nspn=c1_interp->get_num_sfs_per_node(); // Add laplacian of density to value[1] double BETA = C1topopt::get_beta(); for(int k=1;k<=nds_per_ele;k++){ value[1]+=BETA*(x[con[k]*nspn-3]*dn[(k-1)*nspn+1][3] + x[con[k]*nspn-2]*dn[(k-1)*nspn+2][3] + x[con[k]*nspn-1]*dn[(k-1)*nspn+3][3] + x[con[k]*nspn]*dn[(k-1)*nspn+4][3] + x[con[k]*nspn-3]*dn[(k-1)*nspn+1][4] + x[con[k]*nspn-2]*dn[(k-1)*nspn+2][4] + x[con[k]*nspn-1]*dn[(k-1)*nspn+3][4] + x[con[k]*nspn]*dn[(k-1)*nspn+4][4]); } } //object_fn2 //====================gradient object function===================== void grad_object_fn(FE_Model *fem,int en,int nsfns, double *nx,double **dn, double *shape,int size){ double ALPHA = C1topopt::get_alpha() ; double BETA = C1topopt::get_beta(); for(int i=1;i<=nsfns;i++){ shape[i]=ALPHA*nx[i]+BETA*(dn[i][3]+dn[i][4]); } } PAGE 55 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. 182-194, 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. 238-250, 1992. [Ben 91] Bendsoe, M.P. and Rodrigues, H.C., Integrated topology and boundary shape optimization of 2-D solids, Computer Methods in Applied Mechanics and Engineering, vol. 87, pp. 15-34, 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, North-Holland Publishing Company, New York, 1978. [Bot 86] Botkin, M.E., Yang, R.J. and Bennett, J.A., Shape optimization of three-dimensional 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. 2105-2118, 1969. [Dor 64] Dorn, W.S., Gormory, R.E. and Greenburg, H.J., Automatic design of optimal structures, Journal de Mechanique, vol. 3, pp. 25-52, 1964. [Esc 90] Eschenauer, H., Koski, J. and Osyczka, A., Multicriteria design optimization: Procedure and applications, Berlin: Springer, 1990. 46 PAGE 56 47 [Haf 92] Haftka, R.T. and Gurdal, Z., Elements of Structural Optimization, 3 rd Edition, Kluwer Academic Publisher, New York, 1992. [Haf 86] Haftka, R.T. and Grandhi, R.V., Structural shape optimization A survey, Computer Methods in Applied Mechanics and Engineering, vol. 57, pp. 91-106, 1986. [Kir 81] Kirsch, U., Optimum Structural Design: Concepts, Methods and Applications, McGraw-Hill Book Company, New York, 1981. [Koh 86] Kohn, R.V. and Strang, G., Optimal design and relaxation of variational problems, Communications in Pure and Applied Mathematics, vol. 39, pp. 113-137, 139-182, 333-350, 1986. [Kum 00] Kumar, A.V., and Yu, Lichao, An object-oriented modular framework for implementing the finite element method, Computers and Structures, vol. 79, no.9, pp. 919-928, 1998. [Kum 93a] 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. 735-742, m1993. [Kum 93b] Kumar, A.V., Shape and topology synthesis of structures using a sequential optimization algorithm, Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA., 1993. [Mor 82] Morris, V, Foundations of Structural Optimization: A Unified Approach, John Wiley and Sons, 1982. [Pap 90] Papalambros, Panos and Chirehdast, Hehran, An integrated environment for structural configuration design, Journal of Engineering Design, vol. 1, pp. 73-96, 1990. [Sta 84] Stadler, W., Multicriteria optimization in mechanics. Applied Mechanics Review, vol. 37, pp. 277-286, 1984. [Str 86] Strang, Gilbert and Kohn, R.V., Optimal design in elasticity and plasticity, International Journal for Numerical Methods in Engineering, vol.22, pp. 183-188, 1986. [Suz 91] Suzuki, K. and Kikuchi, N., A homogenization method for shape and topology optimization, Computer Methods in Applied Mechanics and Engineering, vol. 93, pp. 291-318, 1991. PAGE 57 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. 334-339, 1986. PAGE 58 BIOGRAPHICAL SKETCH Tae-Joong 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 Kon-Kuk 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 Kon-Kuk University and graduated in spring 1998. He began his masters 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. 49 |