UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 EULERIAN SHAPE DESIGN SENSITIVITY ANALYSIS AND OPTIMIZATION FOR PLANE ELAS TICITY WITH FIXED GRID By YOUNG MIN CHANG A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2004 PAGE 2 Copyright 2004 by YOUNG MIN CHANG PAGE 3 ACKNOWLEDGMENTS I would like to thank my advisor, Dr. Nam Ho Kim, for his significant contributions to this thesis and his valuable guidance and support. I am also grateful to Dr. Raphael T. Haftka and Dr. Ashok V. Kumar for their helpful advice and for serving on my committee. I would also like to express my gratitude to friends in the design lab. I give special thanks to my family and friends. I especially appreciate their love and support. I would never have succeeded in my thesis without them. iii PAGE 4 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................iii LIST OF TABLES .............................................................................................................vi LIST OF FIGURES ..........................................................................................................vii ABSTRACT .......................................................................................................................ix CHAPTER 1 INTRODUCTION........................................................................................................1 1.1 Introduction.............................................................................................................1 1.2 Related Works........................................................................................................5 2 EULERIAN REPRESENTATION OF GEOMETRY.................................................7 3 FINITE ELEMENT ANALYSIS WITH PENALTY BOUNDARY CONDITIONS14 3.1 Finite Element Equation.......................................................................................14 3.2 Penalty Boundary Method....................................................................................17 4 DESIGN PARAMETERIZATION............................................................................23 4.1 Design Parameterization.......................................................................................23 4.2 Computational Geometry......................................................................................27 5 DESIGN SENSITIVITY ANALYSIS.......................................................................29 5.1 Direct Differentiation Method..............................................................................29 5.2 Adjoint Variable Method......................................................................................33 6 NUMERICAL EXAMPLEs.......................................................................................36 6.1 Torque Arm Model...............................................................................................36 6.2 Example of a bracket............................................................................................46 7 CONCLUSION...........................................................................................................53 iv PAGE 5 LIST OF REFERENCES...................................................................................................54 BIOGRAPHICAL SKETCH.............................................................................................56 v PAGE 6 LIST OF TABLES Table page 61 Design sensitivity results are compared with finite difference sensitivity results (perturbation size = 0.0001).....................................................................................40 62 Design variables at the initial and optimum designs................................................42 63 Comparison of the proposed method with the finite elelment analysis and meshfree method......................................................................................................................45 64 Design sensitivity results are compared with finite difference sensitivity results (perturbation size = 0.0001).....................................................................................50 65 Design variables at the initial and optimum designs................................................51 vi PAGE 7 LIST OF FIGURES Figure page 11 Algorithm for the design sensitivity analysis and optimization.................................4 21 Geometric representation methods.............................................................................8 22 Fixed grid approximation of a structural geometry....................................................9 23 Shape densities near the geometric boundary..........................................................10 24 Shape densities of elements in a row.......................................................................11 25 Approximation of a circle........................................................................................12 26 Algorithm for Eulerian representation of geometry.................................................13 31 Displacement boundary conditions on the boundary elements................................17 32 Algorithm for the finite element analysis.................................................................22 41 Design change in the fixed grid. Perturbed design occupies new region.................24 42 Shape design perturbation and corresponding change of shape density..................25 43 boundary at the intersection.....................................................................................27 44 Spline curve that represents the structural boundary................................................27 51 Algorithm for the design sensitivity analysis...........................................................35 61 The torque arm model with the initial shape and load.............................................37 62 Finite element analysis results of the torque arm.....................................................38 63 Finite element analysis results at the optimum design.............................................42 64 Finite element analysis results at the optimum design using MSC/NASTRAN......43 65 Design optimization history for normalized cost and constraint functions..............46 66 The bracket model with the initial shape and load...................................................47 vii PAGE 8 67 Finite element analysis result of the bracket............................................................47 68 Shape density and finite element analysis results at the optimum design................49 69 Design optimization history for (normalized) cost and constraint functions...........52 viii PAGE 9 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 EULERIAN SHAPE DESIGN SENSITIVITY ANALYSIS AND OPTIMIZATION FOR PLANE ELASTICITY WITH FIXED GRID By Young Min Chang May, 2004 Chair: Nam Ho Kim Major Department: Mechanical and Aerospace Engineering Conventional shape optimization based on the finite element method uses Lagrangian representation in which the finite element mesh moves according to the shape change, while modern topology optimization uses Eulerian representation. An approach to shape optimization is proposed using Eulerian representation such that the mesh distortion problem in the conventional approach can be resolved. A geometric model is defined on the fixed grid of finite elements. An active set of finite elements that defines the discrete domain is determined using a similar procedure as the topology optimization, in which each element has a unique value of shape density. The shape design parameter that is defined on the geometric model is transformed into the corresponding shape density variations of boundary elements. Using this transformation, it is shown that the shape design problem can be treated as a parameter or sizing design problem, which is much easier than the former. A detailed derivation of how the shape design velocity field can be converted into the shape density variation is presented along with sensitivity ix PAGE 10 calculation. The coefficients are calculated very efficiently by integrating only those elements that belong to the structural boundary. The sensitivity results obtained from the proposed approach are compared with those from the global finite difference method with excellent agreement. Two design optimization problems with plane stress are presented to show the feasibility of the proposed design approach. x PAGE 11 CHAPTER 1 INTRODUCTION 1.1 Introduction As CAD tools are rapidly developed, it is natural to couple CAD tools with shape optimization so that the design variables are chosen from CAD parameters, which makes a consistency between the design model and CAD model [1,2]. Many CAD tools now provide the shape design and optimization capability using the finite element method. A critical weakness of the geometrybased shape optimization is mesh distortion during the structural analysis process [3]. Regularly distributed mesh in the initial design is often distorted during the shape optimization. The distorted mesh deteriorates the solution accuracy of finite element analysis. In order to maintain the solution accuracy in a certain level, the adaptive mesh generation method in which the element properties are held constant while the mesh is iteratively improved has been studied [3,4]. However, work needs to be done for these methods to be effective. The conventional shape optimization is referred to as a Lagrangian shape representation method since both the geometry and the finite element mesh move together during the shape optimization process. In contrast to the Lagrangian method, the topology optimization method has recently been developed to determine the optimum shape of a structure without causing mesh distortion [5,6]. Although the material property of each element changes as a design variable changes, the initial geometry of the finite element mesh is unchanged during the design process. However, a great number of design variables make it difficult to find the optimum design. In addition, the optimum design often raises questions to the 1 PAGE 12 2 manufacturability of the final result. It is a nontrivial task to determine the structural boundary shape from the topology optimization result. This approach is referred to as an Eulerian method since the finite element mesh is fixed during the design process. The objective of this research is to develop a new shape optimization method based on Eulerian shape representation method. To enhance this approach, a fixed grid is used to represent the computational model, while a Lagrangian method is used to represent the geometric boundary accurately. For those elements on the boundary, the perturbed shape is interpreted as a shape density change. A new shape optimization method is proposed in this research that adopts advantages of both conventional shape and topology optimization methods. In other words, it uses advantages of accurately representing the geometric model in the Lagrangian method and that of maintaining the same mesh quality in the Eulerian method. In structural analysis, the geometric model is overlapped on the regularly meshed finite elements. The geometry of finite elements is fixed, while the geometry model is changed according to the shape design. The concept of shape density in the topology optimization is adapted in order to distinguish the structure and void. The finite elements that belong to the inside of the structure have a full magnitude of the shape density, while the finite elements that belong to the outside of the structure have a zero magnitude of the shape density (a void). Finally, the finite elements that are on the boundary of the structure have a shape density proportional to the area fraction between the material part and the void part. Thus, the finite elements on the boundary have shape densities between a full magnitude and a zero magnitude. This idea is similar to the homogenization method in PAGE 13 3 the topology optimization. Thus, it is referred to as boundary homogenization method in this research. Since the shape change of the geometric model causes the shape density change of the finite elements on the boundary in shape optimization, a new shape density should be calculated for those elements. In addition, some elements leave or enter the structure domain. Thus, accurate record keeping is an important part of the proposed method. The place where each element is located is identified by incrementally searching the boundary curve, and then the area fraction of each element on the boundary is calculated using Greens theorem. The finite elements that belong to the structural domain can be easily identified by counting the number of boundary elements in each row or column of the mesh. Unlike Lagrangian shape representation, this approach dose not require any mesh update process, and the solution accuracy can be maintained during the whole design process because the same size of the finite element is consistently used. A mathematical difficulty of the proposed method is how to represent the effect of shape change into the shape density change. Since the shape design variables are chosen from the geometric parameters, the explicit contribution of the boundary curve shape to the shape density of the boundary element is calculated based on the geometric relation. Accordingly, the boundary shape design velocity is related to the shape density of the boundary elements which is used in design sensitivity calculation. Thus, the complicated shape design sensitivity formulation can be converted to a simple, parametric design sensitivity formulation. In the conventional shape design sensitivity analysis and optimization, the shape design variable perturbs the boundary curve or surface [7]. Thus, theoretically it is PAGE 14 4 complete that the shape design sensitivity formulation can be expressed in terms of boundary functional. When the finite element method is used for the numerical approach, however, function evaluation on the boundary is not inherently accurate. Thus, the domain method has been developed in which the boundary perturbation induces the domain perturbation [8]. However, the mapping from the boundary perturbation to domain perturbation is not a onetoone relation. Thus, various methods have been developed to calculate the domain design velocity field [7]. However, the proposed method eliminates such an inconvenience because the formulation only affects those elements on the structural boundary. Numerical integration involved in the sensitivity calculation is limited for those elements on the boundary, which provides efficiency for the proposed approach. Figure 11. Algorithm for the design sensitivity analysis and optimization. PAGE 15 5 Figure 11 shows the whole process for the design sensitivity analysis and optimization. Each process will be explained in the following chapters in detail. 1.2 Related Works The structural shape design is an iterative process in which the geometry is modified and reanalyzed until the design objective function is minimized or maximized while constraints satisfy. In order to automate the design process, researchers have been interested in integrating the design and analysis procedures for decades. Recently, it has been shown that the fixed grid method can be used to integrate modeling and analysis in a single step of the design process. Fixed grid technique has been used in problems in which the geometry of the object or the physical properties of the body change with time [9]. A fixed grid representation of the finite element domain was used to solve an elasticity problem by Garca and Steven [10]. The 2dimensional elasticity problem was presented in order to evaluate the fixed grid finite element procedure. They showed that the accuracy of the fixed grid method is improved as the size of elements decreases. A least square method was used to approximate the stress on the boundary in order to reduce the error associated with the stepwise approximation of the boundary. They also used a fixed grid finite element method for structural design and optimization [11]. They achieved a fast reanalysis process by integrating the analysis and design process in single step which is similar to discrete semianalytical method in sensitivity analysis [12]. Recently, Kim et al. [13] applied the fixed grid approach to the evolutionary structural optimization (ESO) problem. An ESO process is started by generating a stiffness matrix of the given initial design. Once the matrix is defined, it is solved for displacement and stress values of each element. ESO then physically removes a small PAGE 16 6 percentage of elements that have low stress values. This completes one cycle of the ESO process. Repeating this process leads to the optimum design. Implementing the fixed grid methodology not only simplifies the mesh generation process, but also allows a significant reduction in the arithmetic calculation to update the stiffness matrix for the modified topology, instead of a full regeneration of the matrix. PAGE 17 CHAPTER 2 EULERIAN REPRESENTATION OF GEOMETRY The conventional geometric representation and shape optimization of a solid structure has been based on the Lagrangian approach in which the structural domain and boundary are changed according to the shape design parameters. Such geometric details as fillet surfaces or curvatures can be accurately represented in this approach. However when structural analysis of the structure is performed using the finite element method, mesh distortion occurs in the Lagrangian approach. It is difficult to create a good quality mesh from a complicated CAD geometry [see Figure 21(a)]. Besides, even if a regular mesh is initially created, the mesh quality deteriorates as the structural shape changes during the design optimization. In order to control the quality of finite element analysis, adaptive mesh refinement in shape optimization must be introduced [4]. Recently, many researchers started thinking about representing a structural domain using Eulerian approach in which the grid is fixed in space [5,6]. The region in which the material is occupied has a full shape density and the void part has a zero shape density. The shape change can be characterized using an analogy of fluid flow. The shape density in one region moves to the neighborhood region as the structural shape changes. After integrated with an optimization algorithm, this approach yields the modern form of a topology design [see Figure 21(b)]. Although the topology design approach can provide a creative conceptual design, it is difficult to extract geometric information in the case of complicated threedimensional structures. In addition, it is intricate to physically interpret those regions that have intermediate densities (gray area) between a full material and a 7 PAGE 18 8 void. However, the mesh distortion problem that exists in the Lagrangian shape design problem can be resolved because the mesh geometry is fixed during the design process. (a) Finite element Mesh (b) Topology design Figure 21. Geometric representation methods. As has been mentioned earlier, an advantage of the geometrybased shape parameterization is to accurately represent the structural domain while an advantage of Eulerian approach is to resolve the mesh distortion problem. By taking advantages from both methods, a geometrybased shape parameterization on a fixed grid is proposed in this research. A fixed grid is generated by superimposing a rectangular grid of equal sized elements on the given structure instead of generating a mesh to fit to the structure, which reduce the human effort in modeling the structure. In addition, the mathematical expression will be simplified when the grid is distributed parallel to the coordinate directions. In this research, a rectangular grid is used, which makes a fixed grid parallel to the x 1 and x 2 A solid geometry with domain reduces space and the boundary is defined on a regular, rectangular mesh. Figure 22 shows an example of a structure modeled by a fixed grid. In this way, elements are either inside, outside, or boundary of the structure. If an element belongs to the domain then it has a full shape density. If an element is out of the domain then it has a zero shape density. PAGE 19 9 Figure 22. Fixed grid approximation of a structural geometry. Although the approximation in Figure 22 seems straightforward, a technical difficulty exists for those elements that reside on the structural boundary. Part of the element belongs to the structural domain, while the other part is in a void. The idea of homogenization is used for the elements on the structural boundary. The participation of each element can be determined using the idea of shape density, which measures the amount of element area that belongs to the structural domain Let the area of an element m be A, and let the area that belongs to the domain m be a. The shape density of element m can be calculated by m 1,if0,if/,ifmmmmmm m A AuAaAA (2.1) where A=A means that element m is the inside of the domain(elements 2 and 3 in Figure 23), while A m m m = means that no portion of element m is located in the domain(elements 7, 8, and 9 in Figure 23). When an element is on the boundary PAGE 20 10 (elements 1, 4, 5, and 6 in Figure 23), the shape density u is a ratio of the area that belongs to the domain to A. m ma m Figure 23. Shape densities near the geometric boundary. In order to calculate the shape density, the domain integration is required for the elements on the geometric boundary. It is difficult to calculate the area a with the general domain integration procedure since the boundary curve cuts the element arbitrarily. Greens theorem is used to convert the domain integration into boundary integration. In general twodimensional problem area integration can be represented by m (2.2) 12mmACadx dx where x1 and x are coordinate directions, C is the curve that surrounds the area a, and the integration direction is counterclockwise. The curve C is comprised of straight element boundary lines and a geometric boundary curve. The evaluation of Eq.(2.2) is straightforward for the straight element boundary because either x 2 m 1 or x 2 is constant. This is because the fixed grid is initially defined to be parallel to the x 1 and x 2 In the case of a boundary curve, the geometry boundary is represented by parameter such that the expressions of x 1 () and x 2 () are available. Using the chain rule of differentiation, the PAGE 21 11 integral in Eq.(2.2) can easily be converted to an integral with respect to parameter After calculating a, the shape density can be obtained from Eq.(2.1). m After determining shape densities of boundary elements, the shape densities of other elements should be determined whether they are interior or exterior. It is assumed that the geometric boundary exists within the fixed grid. Starting from the leftmost element in a row, a value of zero is assigned to the element. When a boundary element is met, then a value of one is assigned. This process is repeated whenever a boundary element is encountered (see Figure 24). Boundary elements um = 1 um = 0 um = 0 Boundary curve Figure 24. Shape densities of elements in a row. Alternatively, if the surface geometry information is available in addition to the curve geometry, than that information can be used to identify those elements that belong to interior of the geometry. For example, when a parametric surface information x(,) is available, the interior elements can be found by incrementally searching parameters and During the structural analysis, the material property of each element is varied. It is calculated using the shape density by mmEuE (2.3) where E is the Youngs modulus of the material and E is the varied modulus. Since the Poissons ratio is related to the lateral contraction during the tensile deformation, it is m PAGE 22 12 fixed during this augmentation process. In the practical application, the shape density for the void part has a small value instead of zero in order to avoid numerical singularity during the finite element analysis process [6]. The approximation of domain in Figure 23 is different from the idea of pixel or voxel [15]. In which a continuum structure is divided by a number of squares. In order to approximate the boundary reasonably, a very fine pixel mesh is required. However, in the proposed method the effect of a continuous boundary is reflected on using the idea of boundary homogenization. As an example, a circle is approximated using pixel approximation and boundary homogenization in Figure 25. It is clear that the boundary homogenization method provides a smooth transition between the structural part and the void part. Indeed, the gray boundary of the topology design result in Figure 21(b) should be understood in the same context as boundary homogenization. However, in the proposed method the structural domain is still represented using boundary curves while Figure 25(b) is a mere approximation of the geometry (a) Pixel approximation (b) Boundary homogenization Figure 25. Approximation of a circle. Figure 26 shows that the algorithm for Eulerian representation of the geometry. First, the interior domain which includes inside elements and those elements on the PAGE 23 13 boundary is identified. Second, the finite elements on the boundary are identified. Then, the shape density on the boundary elements is calculated. In the next step, if the structural model is not optimized, the model is updated until it is optimized. Figure 26. Algorithm for Eulerian representation of geometry. PAGE 24 CHAPTER 3 FINITE ELEMENT ANALYSIS WITH PENALTY BOUNDARY CONDITIONS The proposed Eulerian shape representation method has an advantage in the viewpoint of finite element analysis. Since all elements have the identical shape in the proposed method, it is very efficient to construct one element stiffness matrix and use it repeatedly. Especially when the element is square, the element stiffness matrix can be calculated analytically [15]. 3.1 Finite Element Equation As a structure deforms, not only the internal force but also the structural energy increases. This stored energy is called the strain energy, defined as 1()()()2TU d zzCz (3.1) The strain energy U(z) is the energy required to produce the displacement z. For elastic problem, since U(z) dose not depend on the path chosen for deformation, it is a function of the configuration z. If force is applied to the structure and the structure deforms in the direction of the applied force, then work is done by the applied force. The work done by the applied load can be defined as ()sTTWd d zzfzT (3.2) The first integral in Eq. (3.2) represents the work done by the body force, while the second integral is the work is done by the surface traction load. 14 PAGE 25 15 The total potential energy of the structure is the difference between the strain energy and the work done by the applied loads, defined as ()()().UW zz z (3.3) A variational formulation of the structural problem can be written, using the first variation of (z), as (,)(,)(,)0,UW zzzzzz (3.4) for all z in where z is the variation of displacement z. If the kinematical boundary conditions are given in the continuum domain, the weak form of the structural problem can be written in the following form: (,)(),,a uuzzzz (3.5) where is the space of kinematically admissible displacements, and z means for all virtual displacements z that belong to Eq. (3.5) is a variational equation with displacement z as a solution. In Eq. (3.5), (,)()()Ta d uzzzCz (3.6) and ()sTTd d uzzfzT (3.7) are the structural bilinear and load linear forms, respectively. In Eq. (3.6), (z) is the engineering strain vector, and C is the linear elastic constitutive matrix. In Eq. (3.7), f is the body force and T is the surface traction on the traction boundary s The structural problem described in Eq. (3.5), with definitions in Eq. (3.6) and (3.7), is a standard form in the Lagrangian approach. In this case, represents the structural domain. PAGE 26 16 On the other hand, in the Eulerian approach is the whole domain, including both the structure and void. Let the domain be composed of NE subdomains (finite elements), and let each subdomain m have shape density u m Then, the structural bilinear and load forms can be written in the following forms: 1(,)()()mNETmmma ud uzzzCz (3.8) and 1(),smNETmmudd T uzzfz T (3.9) where in the definition of a u (,) and u (), the subscribed u is used to denote these forms dependence on the design variable vector u = [u 1 u 2 , u NE ] T Since u m is constant within the element, it can be taken outside the integral. It is assumed that the traction force is independent of the design. Even if the domain decomposition has been introduced in Eq. (3.8) and (3.9), all variables are still in the continuum level. Even if the proposed method can be applied for general anisotropic, nonsquare type element, we want to show the feasibility of the method through a simple isotropic planestress element. For twodimensional finite element with square shape, the stiffness matrix can be calculated analytically [15], as 313133113681281286813131313386868128123133113311286868128131313313868681281223113313131286868128[]1Ek 131331313812812868613313133168128128681331313138128128686. (3.10) PAGE 27 17 Note that the stiffness matrix [k] is not a function of geometry, but a function of material properties. In fact, it is independent of element size. All elements have the same [k] matrix with different shape density. In order to calculate the stiffness matrix at each element, the stiffness matrix is multiplied by the shape density at that element. That is, a small value for the outside element of the boundary, one for the inside element of the boundary, and the shape density for the element on the boundary are multiplied by the stiffness matrix. The element stiffness matrix can be calculated by []mmuk .k (3.11) 3.2 Penalty Boundary Method In the Lagrangian approach, there exists a discrete set of nodes along the geometric boundary. Thus, displacement boundary condition can be applied to those nodes on the boundary. In the Eulerian approach since the geometry moves around within a fixed set of finite elements, it is better to apply the displacement boundary condition on the geometric curve or point. Boundary curve h Boundary elements Figure 31. Displacement boundary conditions on the boundary elements. However, the geometric boundary is often located in the interior of the boundary elements. Thus, it is not trivial to apply the displacement boundary conditions along the geometric curve. As an approximation, one can fix all elements that intersect with the PAGE 28 18 displacement boundary curve (see Figure 31). However, this method overestimates the effect of boundary conditions. In this research, the penalty boundary method (PBM) is used to solve the boundary value problem, which is proposed by Clark and Anderson [16]. The PBM employs the penalty method to apply boundary conditions on a simple, regular mesh that completely overlaps the geometry boundary. Traditional methods for applying boundary conditions in finite element analysis require the mesh to be conformed with the geometry boundaries. However, The PBM dose not require the mesh to be conformed with the geometry boundaries. In the PBM, traditional finite element approximation methods are used to generate a system of linear equations for a regular mesh. The boundary conditions are treated as constraints on the system and are incorporated early with the finite element formulation using penalty methods. Let h be the essential boundary of the structure in which the displacement is prescribed. In the penalty boundary method, the prescribed boundary condition is approximated by using the penalty function, as 1()()()2hTP ,d zzgzg (3.12) where g is the prescribed displacement (usually zero for linear elastic problems), and is the penalty parameter. If the displacement z on the boundary h is different from the given value, then Eq. (3.12) penalizes the total potential energy. In order to incorporate Eq. (3.12) with the weak form of the structural problem, the variation of the penalty function needs to be obtained, as (,)(),hTP d zzzzg (3.13) PAGE 29 19 where the superposed  denotes the variation of the function. This variation needs to be added to the weak form in Eq. (3.5). When the finite element method is used, a discrete version of Eq. (3.13) needs to be developed. Let the boundary h intersects with element m. Then, the approximation of Eq. (3.13) becomes (,),hhmmTTTTPd d zzdNNddNg (3.14) where N is the matrix of shape functions using Lagrange interpolation, d is the vector of nodal displacements, and d is the vector of virtual displacements. Note that the integration is only performed along hm The same approach can be applied to the natural boundary condition in which the traction force is applied along the boundary curve. However, in such a case the derivation of the penalty term involves more mathematical elaborations. The penalty function for the traction boundary condition can be stated, as 1()[()][()],2sTQ d zznTznT (3.15) where n is the unit outward normal vector to the boundary and (z) is the stress matrix. In the small deformation problem, the normal vector is calculated based on the initial geometry. Thus, stress is the only function of displacement in Eq. (3.15). The penalty parameter in Eq. (3.15) can be different from that in Eq. (3.12). Same as the displacement penalty function, the variation of Q(z) can be taken, as (,)[()][()].sTQ d zzznznT (3.16) This variation needs to be added to the weak form in Eq. (3.5). In conjunction with the finite element method, the variation in Eq. (3.16) can be approximated by PAGE 30 20 (,),ssmmTTTTTTTTQdzzdBCSSCBddBCST d (3.17) where B is the displacementstrain relation matrix, C is the elasticity matrix, and S is the matrix of normal vectors. Different types of finite elements may have different expressions for these matrices. In twodimensional plane stress problem with the size of aa square finite element using four nodes, these three matrices are defined as 00010000myayayy 0, x axxAxayaxyaxyxayB xa (3.18) 2121010,100E C (3.19) and 0,0xyxnnnn y S (3.20) Note that the integration in Eq. (3.15) is only performed along sm As mentioned in the previous section, the stiffness matrix can be calculated analytically for twodimensional finite element with square shape. In order to apply the boundary condition, the first parts on the righthand sides in Eq. (3.12) and (3.15) are added to the element stiffness matrix [k m ] as follows (3.21) []hsmmTTTTmmudkkNNBCSSCB ,d where the second part on the righthand side is the contribution from the penalty displacement boundary method. It is only applied to those elements that reside on the essential boundaries. The third part is the contribution from the traction force. The applied force vector for the element can also be calculated by PAGE 31 21 ,hsmmmTTTTmmuddfNfNgBCS TdT (3.22) where the second part on the righthand side is the contribution from the penalty boundary method. In most linear static problems, the prescribed displacement is zero, i.e., g = 0. In such a case, the penalty contribution vanishes. The element stiffness matrix and force vector are assembled to construct the global system of matrix equations, as []{}{} KDF (323) The theoretical aspect of the penalty boundary method is that it is unnecessary the virtual displacement to belong to the space of kinematically admissible displacements. The numerical aspect of the method is that the coefficient matrix can be illconditioned as the magnitude of the penalty parameter increases. Thus, there exists a possible difficulty when an iterative matrix solver is employed. Even if the proposed method has many attractive features in design and simulation points of view, it requires a numerically intensive procedure due to the excessive number of finite elements in high resolution. For example, the torque arm structure in Section 6.1 has about 37,000 degreesoffreedom even if it is a simple twodimensional example. It would be very expensive to store the global stiffness matrix in the computer memory. Thus, a multifrontal sparse matrix solver is employed to store only nonzero components of the global stiffness matrix and to solve the finite element matrix equations [17]. Figure 32 shows the algorithm for finite element analysis. The shape density DEN is already available from Figure 27. The element stiffness matrix [k] can also be calculated analytically as explained above. For all elements in the grid, the element position in the global matrix is calculated. Then, if the shape density is equal to zero, the PAGE 32 22 small number is multiplied to the stiffness matrix. Otherwise, the shape density is multiplied to the stiffness matrix. Figure 32. Algorithm for the finite element analysis. PAGE 33 CHAPTER 4 DESIGN PARAMETERIZATION 4.1 Design Parameterization A major difference between the proposed method and the topology design method exists in design parameterization. A design engineer does not have any freedom to control the design direction in the topology optimization. The optimum shape (or topology) of the structure is determined by finding the shape density of each element, which dose not guarantee any continuity or smoothness of the boundary. In the proposed method, the design parameterization is the same as the conventional shape design method in which the structural boundary changes according to the design velocity field that designates the direction of shape change. As will be shown later, it is unnecessary to define the domain design velocity field in the proposed method. The boundary design velocity field is enough to calculate design sensitivity information. In structural modeling, the physical problem is represented by mathematical expressions, which contain parameters for defining the problem. These parameters define the system and are called design variables. Especially in the shape design problem, the parameters that determine the boundary curve are chosen as design variables. For example, when spline curves are used to represent the boundary, the location of control points, which define the spline curves, can be chosen as design variables. As a design variable changes, the structural boundary and domain are changed continuously. Let the initial boundary and domain are changed to the perturbed boundary and domain 23 PAGE 34 24 respectively. Figure 41 shows the perturbed design after the initial design is changed within the fixed grid. Perturbed DesignInitial Design Fixed Grid Figure 41. Design change in the fixed grid. Perturbed design occupies new region. Such a shape perturbation process has an analogy to the dynamic process, having play the role of time [18]. At the initial time = 0, the structural domain is and the boundary is When the firstorder perturbation is used, the material point x can be denoted by (),, xxVxx (4.1) where V(x) is called the design velocity field that designates the direction of shape change, and is a scalar parameter that controls the amount of shape change. Eq. (4.1) describes the shape perturbation of the continuum model. If a discrete model follows the same perturbation with Eq. (4.1), then it is called the Lagrangian approach. As with the continuum model, the initial shape of each finite element geometry changes according to the design velocity field, which frequently causes the mesh distortion problem. However, in the Eulerian approach the discrete finite element model is fixed in the space as illustrated in Figure 42, and each element has a shape density PAGE 35 25 value between zero and one. The effect of shape change appears through the shape density change. This effect only appears in those elements on the structural boundary. Initial Boundary V(x) Perturbed Boundary Figure 42. Shape design perturbation and corresponding change of shape density. An important theoretical issue is how the shape perturbation can be interpreted as a shape density change on the boundary. Note that the shape perturbation is given as a vector quantity (design velocity field), but the shape density variation is a scalar quantity. As the structural shape changes in the direction of the design velocity field, u can be denoted by m ,mmmuuu (4.2) where mu is the design variation [17]. For those elements that reside within the structural domain, mu is zero. Therefore, perturbation in Eq. (4.2) is only applied for those elements on the structural boundary. When the boundary curve is perturbed in the direction of design velocity V(x), the shape density u also changes as shown in Figure 42. For the purpose of explanation, only element m on the boundary curve is considered. The shape density at the perturbed design can be defined as m PAGE 36 26 11,mmmAAmmudAA d J (4.3) where J is the Jacobian matrix of shape perturbation in Eq. (4.1), defined as xVJIxx (4.4) The material derivative formulas for the Jacobian can be found in Choi and Haug [18]. For example, the material derivative of the Jacobian becomes 0,ddivdJ V (4.5) where divV is the divergence of the design velocity. If the shape density in Eq. (4.3) is differentiated by and using the formula in Eq. (4.5), the relation between V(x) and mu can be obtained from the following relation 11,mTmACmmudivdAAVV dn (4.6) where n is the outward unit normal vector to the boundary and C is the boundary of area in the counterclockwise direction. The second equality in the above equation can be obtained from the divergence theorem. It is interesting and important to note that only the normal component of the boundary velocity appears in Eq. (4.6) because the tangential component does not contribute to the shape change. ma When the boundary passes the intersection point P as in Figure 43, the shape densities for elements 1 and 4 are not differentiable because their change is different when the boundary curve moves upward or downward. In this research, however, this case is not considered partly because the contribution of one element is usually small and the sensitivity is summation of the contributions from all elements on the boundary curve. PAGE 37 27 Figure 43. boundary at the intersection. 4.2 Computational Geometry After design parameterization is finished, the corresponding design velocity V(x) is calculated on the boundary curve. For those elements on the boundary, the variation of the shape density can be calculated by integrating the normal component of the design velocity along the boundary curve. There are many available methods for representing the boundary geometry of the structure. For examples, for handling 2dimensional shape design problems, six basic curves exist: algebraic, geometric, fourpoint, Bezier, spline, and Bspline. Here the geometric curve in Figure 44 is employed for the purpose of explanation. The geometric curve is represented by the position vectors and tangent vectors at its two endpoints. xy upu0up1 p0p1 Figure 44. Spline curve that represents the structural boundary. To parameterize the geometric curve, eight geometric coefficients in matrix [G] in Eq. (4.7) can be defined as shape design variables. One important aspect of the geometric curve for shape design purposes is that continuities of the adjacent curves can be maintained at the joining point of two curves by linking shape design variables. In the PAGE 38 28 geometric curve, the coordinates of the curve are represented using a parameter [0,1] as 32010123012300()[,,,]121011001[][][] xppppGM (4.7) where and 00[,]Txyppp 1[,]Txypp 1 p are locations of two end points, respectively, and 00[,]Txydpddpdp and 1[,Txydpddpd 1] p are tangent vectors at both points, respectively. Design variables for the shape problem can be defined by choosing a component of the matrix [G]. For example, when a design variable moves point p 0 in the xdirection, we can define a unit perturbation matrix by 1000[].0000 B (4.8) Then, the design velocity vector can be defined by replacing matrix [G] with matrix [B], as ()[][][], VBM (4.9) where the matrix [B] is a unit perturbation of the matrix [G] in the direction of design variable. The tangent vector and normal vector to the boundary can be calculated by differentiating Eq. (47) with respect to parameter PAGE 39 CHAPTER 5 DESIGN SENSITIVITY ANALYSIS Design sensitivity analysis is used to develop relationships between a variation in shape and resulting variations in functionals that arise in shape design problems. For demonstration purposes, a linear elastic problem is considered in the following sensitivity development. However, a general nonlinear problem can also be taken into account using a similar approach. 5.1 Direct Differentiation Method In this chapter, design parameterization is utilized to derive the shape sensitivity expression in terms of mu Displacement z in Eq. (3.5) implicitly depends on the design through the structural problem in Eq. (3.5), which must be calculated from design sensitivity analysis, as explained below. An important component of design sensitivity analysis is calculating the variation of the state variable by differentiating Eq. (35) with respect to the design, or equivalently, The variation of the state variable can be defined as 00(;)dd zzzxuuu u (5.1) Note that z depends on the design u, where the variation is evaluated, and on the direction u of the design variation. In the direct differentiation method, z is calculated first, and then the chain rule of differentiation is used to calculate the sensitivity of performance functions. 29 PAGE 40 30 Similar to Eq. (5.1), the structural bilinear and load linear forms can be differentiated with respect to the design. Although the design vector and its variation contain NE components, only boundary elements need to be considered in the calculation of mu because it is zero for those elements inside or outside of the structural domain. Let M be the number of elements that belong to the structural boundary. The variation of the structural bilinear form can be obtained using the chain rule of differentiation, as 0(;),(,)(,),daad a uuuuzxuuzzzzz (5.2) where 1(,)()()mMTmma du uzzzCz (5.3) is the bilinear forms dependence on the design. If the structural problem in Eq. (3.5) is solved for z and the design variation mu in Eq. (4.6) is available as a result of design parameterization, then (,)auzz can be readily calculated by following the same integration procedure used in finite element analysis. The second term on the righthand side of Eq. (5.2) is the same as the bilinear form in Eq. (3.8) if displacement z is replaced by z, which will be solved. The variation of the load linear form can also be obtained by following a similar procedure, as 10()().mMTmmddud uuuzzzf (5.4) When the traction boundary is changed according to the design, a careful treatment is required regarding boundary homogenization, which is not developed in this research. PAGE 41 31 When a concentrated load is applied to the structure, the variation of the load linear from in Eq. (5.4) vanishes because the load is independent of the design. After differenting Eq. (3.5) at the perturbed design and using the formulas in Eq. (5.2) and (5.4), the following design sensitivity equation can be obtained (,)()(,),,aa uuuzzzzzz (5.5) where the solution z is desired. If the righthand side is considered to be an applied load, Eq. (5.5) is similar to the structural problem in Eq. (3.5) with a different load, which is called the fictitious load. For a given design variable, the variation of shape density can be calculated from Eq. (4.6). The righthand side Eq. (5.5) can then be calculated using mu and z. By following the same discretization as the finite element method, the matrix equation for the design sensitivity analysis in Eq. (5.5) can be obtained, as (5.6) []{}{},ficKDF where {D} is the sensitivity of the nodal displacement vector and {F fic } is the fictitious load vector, defined as (5.7) 11{}mmMMficTTmmmududFNfBCBd .m All information in Eq. (5.7) is already available from the structural analysis. Thus, calculating the integrals in Eq. (5.7) is relatively convenient with the design variation mu Consider a general performance function defined in the form of integral over the domain as ,(),(),b d uzuuzu (5.8) PAGE 42 32 The performance can be a pointwise function if Diracdelta measure is used as an integrand. If the structural volume or area is a performance function, then the shape density mu can be an integrand with summation over all boundary elements. If stress at element m is concerned, then the integrand is the stress function and integrated over subdomain m The sensitivity of in Eq. (5.8) can be obtained by taking variation with respect to the design u, as bbd uzuz (5.9) Equation (59) can be readily evaluated using the solution of Eq. (5.5) and mu The gradient information that is necessary for design optimization is equivalent to the coefficient of u. Thus, the coefficient of u in Eq. (5.9) is called the sensitivity coefficient. Compared to the shape design sensitivity formulation in the Lagrangian approach, the expressions in Eq. (53) and (54) provide significantly simple computational methods, since their expressions also appear during regular finite element analysis. In geometrybased shape optimization, domain integration is involved in Eq. (53) and (54). However, only the boundary integral is sufficient for the proposed method. In the computational viewpoint, the lefthand side of Eq. (5.6) is the same as that of Eq. (3.23) if {D} is replaced by {D}. Thus, solving the sensitivity equation becomes very efficient when a direct matrix solver is used. For example, the coefficient matrix [K] in Eq. (3.23) is factorized during finite element analysis. In design sensitivity analysis, the factorized coefficient matrix can be reused for the calculation of {D}. Thus, the PAGE 43 33 major computational effort involved in sensitivity analysis is to construct the fictitious load in Eq. (5.7) and then, forwardand backwardsubstitutions to solve for {D}. 5.2 Adjoint Variable Method The main idea of the adjoint variable method is to avoid the direct calculation of z in Eq. (5.5). Since the sensitivity expression in Eq. (5.9) requires the calculation of z, the adjoint problem is defined, as using the coefficient of z in Eq. (59) as a load term. Accordingly, the adjoint problem is defined, as (,),,bad uz (5.10) where the adjoint solution is unknown and its variation plays the same role as z in Eq. (3.5). By replacing with z in Eq. (5.10) and by replacing z with in Eq. (5.5), it can be shown that the second integrand of Eq (5.9) can be represented by ()(,).bdauuzzz (5.11) In deriving the above equation, the symmetric property of the energy bilinear form has been used. The physical meaning of Eq. (5.11) is that the implicit dependence of the performance function has been eliminated using the adjoint solution. Thus, the sensitivity of the performance function can be expressed in terms of the structural solution z and the adjoint solution as (,)au ()(,).bdauuuzu (5.12) It would be beneficial to compare the efficiency of the direct differentiation and adjoint variable methods in computational viewpoint. It is interesting to note that the adjoint Eq. (5.10) is independent of the design. In fact, each performance function has PAGE 44 34 different adjoint load on the righthand side. Thus, Eq. (5.10) needs to be solved per each performance function, whereas the sensitivity equation needs to be solved per each design variable in the direct differentiation method. Thus, the adjoint method is more efficient when the number of performance measures is smaller than the number of design variables, which is the case for most optimization problems. In the numerical approach, the adjoint problem in Eq. (5.10) needs to be discretized using the same method as structural analysis. The matrix equation for the adjoint problem becomes (5.13) []{}{},adjKF where {} is the nodal solution of the global adjoint vector and {F adj } is the adjoint load vector, defined by {}Tadjbd Fz (5.14) The calculation of sensitivity in Eq. (5.12) involves a similar procedure as the fictitious load in Eq. (5.7). This process must repeat for each design variable, as the fictitious load depends on the design. Figure 51 shows the algorithm for design sensitivity analysis. The adjoint load is calculated first to solve the adjoint problem. The design variation mu is calculated using Eq. (4.6) as explained in the Chapter 4 and then, the sensitivity coefficient is also calculated for the performances. PAGE 45 35 Calculate adjoint load Solve adjoint problem Calculate ummu Calculate sensitivity coefficient I = 0 I = nperf nperf : the number of performance functionsNo Figure 51. Algorithm for the design sensitivity analysis. PAGE 46 CHAPTER 6 NUMERICAL EXAMPLES Two example problems will be shown to demonstrate Eulerian optimization with fixed grid. These examples are compared with the shape optimization results in the literature [3,15]. 6.1 Torque Arm Model The torque arm model has been used to demonstrate the shape optimization by many different methods. Bennett and Botkin [3] used the Lagrangian method with parameteric boundary geometry. Kim et al. [19] used the meshfree method for shape optimization. Jang et al. [20] performed the design optimization of the torque arm model using the multiscale wavelet approach. The torque arm model is considered as a benchmark problem in the shape optimization. The initial shape and load applied of the torque arm is shown in Figure 61(a). The torque arm consists of 32 points, 28 curves, and 16 surfaces. The rectangular domain is established with lowerleft corner being (7, 8) and upperright corner being (49, 8), which covers the whole structure. The size of each element is 0.21cm 0.21cm. The fixed grid consisting of 267 77 elements is used to solve the problem. Figure 61(b) shows the structural domain that is identified using the boundary homogenization method. The black region which is inside to the structural boundary has a full shape density (u = 1.0), while the gray boundary which is on the structural boundary represents intermediate shape density (0 < u <1) calculated using Eq. (2.1). All void areas are represented using a 36 PAGE 47 37 white color. For material properties, the following values are used: Youngs modulus = 207. 4 GPa, Poissons ratio = 0.3, and thickness = 0.3 cm. (a) (b) b1b2 b3b4b5 b6b7b85066N 2789N b1b2 b3b4 b8 u6Fixed Figure 61. The torque arm model with the initial shape and load. (a) Design parameterization. (b) Boundary homogenization of the torque arm with pixel size = 0.21cm. In finite element analysis, the left circle is fixed and the horizontal and vertical forces are applied at the center of the right circle. In order to apply for the displacement boundary conditions, the boundary curves that correspond to the left circle are identified first. It is easy to retrieve boundary element information corresponding to the displacement boundary curves. Then, the penalty boundary method explained in Section 3 is used to impose the displacement boundary condition. Thus, in this approach displacement boundary conditions are applied in a layer of elements. The penalty parameter for applying the boundary conditions is hundred times greater than Youngs modulus. The force boundary condition can also be applied in the same manner. Figure 62 shows the initial stress distribution for the torque arm. The sparse matrix equation is solved using UMFPACK 4.1 [17]. Von Mises stress is calculated at the center of each element. Maximum stress of about 248 MPa appears at the top and bottom surface of the torque arm (see Figure 62). This result is expected because the applied PAGE 48 38 force is superposition of compressive and bending loads. In addition, relatively high stress concentration is observed at the end of interior slot, which is caused by distortion at the small radius region. Figure 62. Finite element analysis results of the torque arm. In the mathematical point of view, the pixelbased geometric representation may cause singularity at the nonsmooth boundary, which is inevitable when inclined boundary is approximated by xand ydirectional squares. However, the proposed approach reduces such singularity by gradually reducing the shape density at the boundary. However, different material properties between interior and boundary elements cause stress discontinuity. Smoothening algorithm in stress may help to reduce the discontinuity [10,20]. Since design parameters are defined on the geometric model, the horizontal and vertical movements of geometric points can be selected as design parameters. Eight design parameters are chosen (see Figure 61). They are the control points of the spline curves that describe the shape of the structure. Design parameters are linked in order to maintain symmetric geometry. As design parameters are changed, new shape density for each finite element is calculated from which the material constants are changed as depicted in Eq. (2.3). PAGE 49 39 Since there is no analytical solution, it is difficult to verify the accuracy of the sensitivity results. The finite difference method would be the only choice available for verification purpose. This verification by no means guarantees the accuracy of sensitivity results. It assures the consistency with the numerical method employed. The finite difference method perturbs the design variable with a small amount () and solves the structural problem again. The sensitivity can then be approximated by ()(,u )u (6.1) where u is the current design value and u + is the perturbed design value. The process must be repeated for each design variable. Design sensitivity results are compared with finite difference sensitivity results in Table 61. Design variables are perturbed with a small perturbation of = 10 4 Since there are eight design variables, the finite difference method needs to perturb the design eight times and solves the structural problem repeatedly. Three different types of performance functions are considered: structural area, maximum von Mises stress, and ydirectional displacement at the location (10, 0). In Table 61, the first column is the design variable, the second column is the performance type and its value, and the third column is the performance change calculated from the finite difference method, the fourth column is the performance change estimated from the proposed sensitivity results, and the last column is ratio between the third and fourth columns. It shows an excellent agreement between two methods. Thus, the proposed sensitivity calculation method can be used for the gradient calculation during the optimization. The biggest advantage of the proposed method is its computational efficiency. The cost of sensitivity calculation is less than 5% of the structural analysis cost per design variable. PAGE 50 40 Table 61. Design sensitivity results are compared with finite difference sensitivity results (perturbation size = 0.0001) Design / 100 (%) Area3.749200E+021.036225E041.036146E04100.01MAX2.484748E+026.082722E046.084701E0499.97zy9.215100E021.237382E071.237448E0799.99Area3.749200E+023.566629E033.566612E03100.00MAX2.484748E+022.092790E022.094751E0299.91zy9.215100E024.259725E064.259550E06100.00Area3.749200E+021.011942E041.011766E04100.02MAX2.484748E+022.957883E052.937141E05100.71zy9.215100E021.231901E081.231175E08100.06Area3.749200E+023.482701E033.482692E03100.00MAX2.484748E+021.010746E031.011070E0399.97zy9.215100E024.238100E074.237938E07100.00Area3.749200E+021.999743E041.999999E0499.99MAX2.484748E+023.757007E043.758011E0499.97zy9.215100E023.233982E073.234727E0799.98Area3.749200E+021.654841E031.654802E03100.00MAX2.484748E+023.386682E033.386067E03100.02zy9.215100E027.485533E077.484632E07100.01Area3.749200E+022.000052E041.999999E04100.00MAX2.484748E+025.511106E045.510200E04100.02zy9.215100E023.767109E073.766835E07100.01Area3.749200E+021.654854E031.654803E03100.00MAX2.484748E+022.847409E032.847083E03100.01zy9.215100E021.359770E061.359566E06100.01b8b4b5b6b7b1b2b3 A simple design optimization problem is formulated to minimize the area of the structure with the maximum stress constraint. In order to induce large shape change, a constraint limit that is far away from the initial value is deliberately provided. Thus, the design optimization problem can be stated that 0max0areaMinimizeSubjectto10A (6.2) PAGE 51 41 In Eq. (6.2), the cost function and constraint are normalized such that the cost function is one at the initial design and the constraint is zero at the optimum design. The lower and upper bounds of the design parameters are selected in order to maintain the topology structure. For this particular example, A 0 = 374.9 cm 2 and 0 = 800 MPa are used. Even if the maximum stress is not a differentiable function, we calculate the sensitivity of the element that has a maximum stress at the current design, assuming that it is a differentiable with respect to the design. This assumption can cause difficulty during design optimization when the maximum stress is evenly distributed throughout the structure. The design optimization problem is solved using the sequential quadratic programming method implemented in DOT [21]. The optimization procedure is to find a search direction and then change design in the search direction. The general optimization procedure follows the next steps. The first step in finding the search directions is to determine which constraints, if any, are active or violated. Here active constraint is defined as one with a value between a small negative number and a small positive number. On the other hand, any constraint is more than a small positive number, it is defined as violated. After identifying active and violated constraints, the gradients of the objective function and all active and violated constraints are calculated. Then, a search direction is found using the gradients. Three methods are available in DOT for solving the constrained minimization problem. These include the modified feasible direction method, sequential linear programming method, and sequential quadratic programming method. The torque arm model is solved using the modified feasible direction method. The advantage of this method is that it only moves within the feasible design so that PAGE 52 42 every intermediate design satisfies the constraint. Function values and sensitivity information are provided to the gradientbased optimization algorithm. The design optimization problem is converged after the fifth iteration. Table 62 shows the design variables values at the initial and optimum designs. All initial design variables start from zero so that the design variables value represents the relative change of the dimension from the initial design. Two design variables are on the boundary, five design variables are very close to the boundary, and u 8 is in the middle of design domain. The lower and upper bounds are selected such that the topology of the structure maintains. Table 62. Design variables at the initial and optimum designs DesignLower BoundInitial DesignUpper BoundOptimum Design b 1 3.00000.01.0002.999660 b 2 0.50000.01.0000.500000 b 3 1.00000.01.0000.999722 b 4 2.70000.01.0002.697690 b 5 5.50000.01.0005.499770 b 6 0.50000.02.0002.000000 b 7 1.00000.06.0005.999490 b 8 0.50000.00.0000.238489 Figure 63. Finite element analysis results at the optimum design. (Unit: MPa) PAGE 53 43 Figure 63 shows the shape density and stress contour plots at the optimum design. The optimization algorithm chose the geometry such that the maximum stress is evenly distributed along the upper and lower regions of the structure. The optimum design conforms to engineering sense because in such a beamlike structure the moment of inertia needs to be increased as the moment arm increases. (a) (b) Figure 64. Finite element analysis results at the optimum design using MSC/NASTRAN. (a) 459 elements. (b) 1608 elements. (c) 3589 elements. (d) 5358 elements. (Unit: MPa) PAGE 54 44 (c) (d) Figure 64. Continued. In order to demonstrate the accuracy of optimum result, the Lagrangian finite element analyses at the optimum model were performed by using MSC/NASTRAN with a triangular grid of gradually increasing the number of elements because of the mesh distortion problem when quadrilateral elements were used. Figure 64 shows the four finite element analysis results with the different number of elements, which show the PAGE 55 45 stress contour plot and the maximum stress. Each model consists of 459 elements, 1608 elements, 3589 elements and 5358 elements, respectively. For all cases, the maximum stress appears at the same location: above the inside slot. The maximum stress with the lowest number of elements in Figure 64(a) was a little higher. However, the maximum stress became a little lower than that of the proposed method as the number of elements was increased. The stress distributions also show consistent trends. The optimum result also can be compared with literature. The same torque arm model was solved using meshfree method by Kim et al. [19]. In that paper, the analysis result at optimum design shows the similar shape and stress distribution with a little lower maximum stress. However, they calculated the stress value at the center of the integration cell. Table 63 shows the comparison of the proposed method with the finite element analyses using MSC/NASTRAN and meshfree method. In Table 63, the second row is the number of elements at the Lagrangian finite element analysis, the third row is the maximum stress, and the last row is the ratio of maximum stress between the proposed method and the other methods. Table 63. Comparison of the proposed method with the finite elelment analysis and meshfree method Eulerian Meshfreewith FG methodNumber of element459160835895358Maximum stress (MPa)800810778779791762Ratio (%)100.0101.397.397.498.995.3Lagrangian FEA using a triangular grid Figure 65 shows the scaled history plot of cost and constraint functions. The optimum design reduces more than 48% of the structural area. The maximum stress at the optimum design appears to be 800MPa. Even if the design converges in the fifth iteration, PAGE 56 46 the structural analysis has been performed fifty four times. This is because the modified feasible direction method requires many line searches. 0.80.60.40.20.00.20.40.60.81.0012345Iter Cost Constraint Figure 65. Design optimization history for normalized cost and constraint functions. 6.2 Example of a Bracket The second example is the bracket model. The initial shape of the bracket is shown in Figure 66. The bracket model has been used to demonstrate the shape optimization by different methods [3,20]. Two holes in the bottom are fixed in space, while a horizontal force of 15,000N is applied at the upper hole. The geometry of the bracket consists of 52 points, 44 curves, and 27 surfaces. The rectangular domain is set with lowerleft corner being (1, 8) and the upperright (8, 20) which covers the whole structure. The size of each element is 0.11cm 0.11cm. The fixed grid consisting of 146 191 elements is used to solve the problem. The Youngs modulus of the material is 207.4 GPa and the Poissons ratio is 0.3. The thickness is 0.3cm. Von Mises stress is calculated at the center of each element. PAGE 57 47 The finite element analysis result is shown in Figure 67. The maximum stress of 410 MPa appears in the upper circle. b1b1b2b3b4b5b4b5b8b8b9b9 b6b6b7b7b10b11b10b11b12 15 kN (a) (b) Figure 66. The bracket model with the initial shape and load. (a) Design parameterization. (b) Boundary homogenization of the bracket with pixel size = 0.11cm. Figure 67. Finite element analysis result of the bracket. (Unit : MPa) PAGE 58 48 In the bracket model, twelve design parameters have been selected to change the inner and outer boundary of the bracket, while maintaining symmetry. Three performance functions have been chosen: structural area, maximum von Mises stress, and xdirectional displacement at the location (0.37, 14.15). Since the number of design variables is greater than the number of performance, it is clear that the adjoint variable method would be more efficient than the direct differentiation method as explained at Chapter 5. The sensitivity coefficients for three performances are calculated using the adjoint variable method. It took 10.3 seconds for design sensitivity analysis using 2.0 GHz desktop computer. On the other hand, the structural analysis took 11.9 seconds for solving the matrix equation. If the finite difference method is employed, it would take 11.9 seconds 12 design variables = 142.8 seconds for calculating sensitivity information without mentioning design perturbation effort. Thus, the proposed sensitivity calculation method is quite efficient. Table 64 shows the sensitivity results compared to the finite difference results. For the finite difference method, a perturbation size of = 10 4 is employed. The first column is the design variable, the second column is the performance type and its value, and the third column is the performance change calculated from the finite difference method, the fourth column is the performance change estimated from the proposed sensitivity results, and the last column is ratio between the third and fourth columns. Two methods agree very well, as all ratios are very close to 100%. The same design optimization problem as with Eq. (62) is used to minimize the area of the bracket structure. For this example, the sequential quadratic programming method is employed. The basic concept is quite simple. First, a Talyer series PAGE 59 49 approximation of the objective and constraint functions is created. However, instead of minimizing the linearized objective, we create a quadratic approximating of the objective function. A linearized constraint with the quadratic objective function is used to create a search direction. The objective function and the constraint are normalized by A 0 = 145 cm 2 and 0 = 800 MPa, respectely. Figure 68 shows the shape density plot and the structural analysis results at the optimum design. In the optimum design, the size of the inner triangle is bigger than the initial design. The maximum stress increases to the constraint boundary so that the constraint becomes active. This result can be compared with literature. The same bracket model was solved using meshfree method by Kim et al. [19]. In that paper, the optimum result shows a similar shape and stress distribution with a little lower maximum stress. Figure 68. Shape density and finite element analysis results at the optimum design. (Unit: MPa) PAGE 60 50 Table 64. Design sensitivity results are compared with finite difference sensitivity results (perturbation size = 0.0001) Design/100(%)Area1.452729E+026.048413E046.048400E04100.00MAX4.100496E+027.123357E077.122662E07100.01zx1.731155E021.677348E081.677140E08100.01Area1.452729E+022.916115E042.916118E04100.00MAX4.100496E+027.314450E077.313131E07100.02zx1.731155E022.852526E102.851279E10100.03Area1.452729E+024.499358E044.499347E04100.00MAX4.100496E+021.712231E071.712719E0799.97zx1.731155E021.829742E101.829556E10100.01Area1.452729E+021.068117E031.068115E03100.00MAX4.100496E+023.249783E063.247563E06100.07zx1.731155E023.614318E073.614374E07100.00Area1.452729E+023.555796E043.555831E04100.00MAX4.100496E+021.083351E061.081315E06100.19zx1.731155E021.203238E071.203252E07100.00Area1.452729E+021.170124E031.170127E03100.00MAX4.100496E+021.666535E031.666363E03100.01zx1.731155E022.674115E072.674292E0799.99Area1.452729E+023.895439E043.895433E04100.00MAX4.100496E+025.547770E045.547374E04100.01zx1.731155E028.902625E088.902916E08100.00Area1.452729E+024.320476E044.320467E04100.00MAX4.100496E+021.557322E021.557651E0299.98zx1.731155E021.129704E081.129773E0899.99Area1.452729E+021.438321E041.438314E04100.00MAX4.100496E+025.185300E035.185555E03100.00zx1.731155E023.761005E093.761081E09100.00Area1.452729E+020.000000E+001.127218E200.00MAX4.100496E+020.000000E+006.600613E160.00zx1.731155E020.000000E+003.648030E200.00Area1.452729E+023.999396E043.999419E04100.00MAX4.100496E+026.968617E086.967609E08100.01zx1.731155E026.246176E116.245705E11100.01Area1.452729E+022.916109E042.916118E04100.00MAX4.100496E+023.804684E063.801079E06100.09zx1.731155E025.256567E085.256648E08100.00b1b2b3b4b5b6b7b12b8b9b10b11 PAGE 61 51 Table 65. Design variables at the initial and optimum designs DesignLower BoundInitial DesignUpper BoundOptimum Designb10.9000.0001.6001.600b21.0000.0001.5001.000b31.0000.0001.8001.800b41.0000.0000.5000.908b50.2000.0000.2000.172b60.8000.0000.6000.800b70.2000.0000.2000.200b81.0000.0000.5000.936b90.2000.0000.2000.196b100.8000.0000.0000.500b111.0000.0001.8001.800b121.0000.0005.5005.150 Table 65 shows the values of design variables at the initial and optimum designs. The lower and upper bounds are selected such that the structure maintains its topology during design optimization. All design varialbes start from zero, which means the values of design variables are relative change of its coordinates. Five variables are changed up to their lower or upper bounds. Figure 69 shows the design optimization history of cost and constraint functions. There are no active or violated constraints during the first six iterations. Thus, the cost function is significantly reduced until the constraint becomes activated. Next nine iterations are carried out to further reduce, while maintaining the constraint active. An attempt to further reducing the cost function at the sixteenth iteration caused the constraint violation. The rest of iterations are consumed to overcome the violated constraint without reducing the cost function. The optimization problem coverged at the twenty second iteration. The cost function is reduced by 48% from the original value. The increased height of the inner triangle (b 12 ) produces the most significant reduction in PAGE 62 52 structural area, until the stress constraint is activated at the side frame. The maximum value of stress, located at the top of inner triangle, moves to the side frame at 800 MPa. A total of thirty two response analyses and twenty two sensitivity analyses were carried out during twenty two optimization iterations. If we consider the cost of three sensitivity analyses is equivalent to one structural analysis, the total optimization cost is less than fourty structural analyses, which is a significant reduction compared to three hundred structural analyses when the finite difference sensitivity is used. When the Lagrangian approach is used with the remeshing process, the optimization process converges at thirty four iterations with seven remeshing processes [3]. When the meshfree method [19] is employed, the optimization problem converged with seventeen iterations and thirty seven function evaluations. Thus, the Eulerian approach performs more design iterations, but the number of structural analyses are slightly less than the meshfree method. The optimum cost function of the Eulerian approach is about 1.5% smaller than that of the meshfree method. 0.60.40.20.00.20.40.60.81.005101520Iter Cost Constraint Figure 69. Design optimization history for (normalized) cost and constraint functions. PAGE 63 CHAPTER 7 CONCLUSION The conventional shape optimization has an advantage of accurately representing the geometric model. However, mesh distortion problem happens because the finite element mesh moves according to the shape change during the optimization. On the other hand, mesh distortion problem is resolved in the topology optimization method. But it is difficult to obtain the geometry information. Thus, we used the advantages from both methods. A new domain approximation method using an Eulerian description is developed for the structural optimization problem. For that end, the fixed grid is used to integrate modeling and analysis in a single step of the design process. Boundary homogenization provides a unique approximation of the structural domain and boundary on the fixed grid of finite elements. Design parameterization on the geometric model provides accurate representation of design intent. For elements on the boundary, the variation of the shape density can be calculated by integrating the design velocity along the boundary curve, which plays a key role in making this approach possible. The main contribution of the current work is that we formulate the shape design sensitivity analysis which involves in geometry change, into the parameter sensitivity analysis. In order to be a practical engineering tool, the proposed approach needs to be extended to threedimensional structures, which involves in boundary homogenization of a volume. As systems degreesoffreedom increase significantly, an iterative matrix solver may need to be incorporated. 53 PAGE 64 LIST OF REFERENCES 1. E. Hardee, K. H. Chang, J. Tu, K. K. Choi, I. Grindeanu, X. Yu, A CADbased Design Parameterization for Shape Optimization of Elastic Solids, Advances in Engineering Software, Vol. 30, pp.185199, 1999. 2. N. Olhoff, M. P. Bendse, J. Rasamussen, On CADIntegrated Structural Topology and Design Optimization, Computer Methods in Applied Mechanics and Engineering, Vol. 89, pp. 259279, 1991. 3. J. A. Bennette, M. E. Botkin, Structural Shape Optimization with Geometric Description and Adaptive Mesh Refinement, AIAA Journal, Vol. 23, pp. 458464, 1985. 4. J. A. Bennett, M. E. Botkin (Editors), The Optimum Shape Automated Srtuctural Design, Plenum Press, New York, 1986. 5. K. Suzuki, N. Kikuchi, A Homogenization Method for Shape and Topology Optimization, Computer Methods in Applied Mechanics and Engineering, Vol. 93, pp. 291318, 1991. 6. M. P. Bendse, Optimization of Structural Topology, Shape, and Material, SpringerVerlag, Berlin, Heidelberg, 1995. 7. K. K. Choi, K. H. Chang, A Study of Design Velocity Field Computation for Shape Optimal Design, Finite Elements in Analysis and Design, Vol. 15, pp. 317341, 1994. 8. K. K. Choi, H. G. Seong, A Domain Method for Shape Design Sensitivity Analysis of BuiltUp Structures, Computer Methods in Applied Mechanics and Engineering, Vol. 57, pp. 115, 1986 9. V. R. Voller, C. R. Swaminathan, B. G. Thomas, Fixed Grid Techniques for Phase Change Problems : a review, International Journal for Numerical methods in Engineering, Vol. 130, pp. 17191731, 1990 10. M. J. Garca, G. P. Steven, Fixed Grid Finite Elements in Elasticity Problem, Engineering Computations, Vol. 16, pp. 145164, 1999. 54 PAGE 65 55 11. M. J. Garca, G. P. Steven, Fixed Grid Finite Element Analysis in Structural Design and Optimization, Proc. 2 nd ASMO/AIAA Internet Conference Approximations and Fast Reanalysis Engineering Optimization, http://wwwtm.wbmt.tudelft.nl/~wbtmavk/2aro_conf 2000, Last accessed December 15, 2003. 12. K. Horimatsu, N. Kikuchi, A Shape Optimization Method Based on Fixed Grid Analysis, IACM Third World Congress in Computational Mechanics, Chiba, Japan, Vol. II, pp. 10701071 13. H. Kim, O. M. Querin, G. P. Steven, Y. M. Xie, Improving Efficiency of Evolutionary structural Optimization by Implementing Fixed Grid Mesh, Structural and Multidisciplinary Optimization, Vol. 24, pp. 441448, 2003. 14. K. Terada, T. Miura, N. Kikuchi, Digital ImageBased Modeling Applied to The Homogenization Analysis of Composite Material, Computational Mechanics, Vol.20, pp. 1997 15. O. Sigmund, A 99 Line Topology Optimization Code with the Matlab, Structural and Multidisciplinary Optimization, Vol. 21, pp. 120127, 2001 16. B. W. Clark, D. C. Anderson, Penalty Boundary Method, Finite Elements in Analysis and Design, Vol. 39. pp. 387401, 2003 17. T. A. Davis, UMFPACK 4.1 User Guide, Univ. of Florida, 2003 18. K. K. Choi, E. J. Haug, Shape Design Sensitivity Analysis of Elastic Structures, Journal of Structural Mechanics, Vol. 11, pp. 231269, 1983. 19. N. H. Kim, K. K. Choi, M. E. Botkin, Numerical Method for Shape Optimization Using Meshfree Method, Structural Multidisciplinary Optimization, Vol. 24, pp. 418429, 2003 20. G. W. Jang, Y. Y. Kim, K. K. Choi, Remeshfree Shape Optimization Using The WaveletGalerkin Method, International Journal for Solids and Structures, submitted, 2003 21. G. N. Vanderplaats, DOT Users Manual, VMA Corp., Colorado Springs, Colorado, 1997 PAGE 66 BIOGRAPHICAL SKETCH Young Min Chang was born in Seoul, Korea. He received his bachelors degree in naval architecture and ocean engineering in 1997 from Inha University, Incheon, Korea. In August of 2001, He began to study as a masters student in mechanical engineering at the University of Florida. He will graduate in May of 2004 with his masters degree in mechanical and aerospace engineering. 56 