Lecture Notes of the CISM Course " Shape and Layout Optimization in Structural Design ". ed. by G. Rozvany Structural Optimization of a Linearly Elastic Structure using the Homogenization Method Report No. UM-MEAM-91-03 by Noboru Kikuchi and Katsuyuki Suzuki Computational Mechanics Laboratory College of Engineering The University of Michigan Ann Arbor, MI 48109 U.S.A. July 1990 Revised December 1990

h 1r ut q P -b3 9

Structural Optimization of a Linearly Elastic Structure using the Homogenization Method Noboru Kikuchi and Katsuyuki Suzuki The University of Michigan, Ann Arbor, MI 48109, USA 1. Introduction There are three major structural optimization problems of a linearly elastic structure; namely, 1) sizing, 2) shape, and 3) layout(topology) optimization problems. The characteristics of these problems can be summarized as follows: Sizing Problem A typical setting of the problem is to find the optimal thickness distribution of a linearly elastic plate that is supported on its boundary and is subject to a given loading condition on the plate or its boundary. The optimal thickness of the plate is obtained so as to, e.g., minimize ( or maximize ) a certain physical quantity such as the mean compliance, while the state (i.e. equilibrium) equation of the structure is maintained as well as various constraints on the state and design variables. Here, the state variable can be the deflection of the plate, and the thickness of the plate becomes the design variable. If a beam is subject to axial torsional and bending forces, the sizes of the cross section of the beam such as the radius, height, and width would be optimized to construct the strongest beam structure within a given design restriction. That is, the design variable of the sizing problem is a physical dimension of the structure. The main feature of the sizing problem is that the domain of the design and state variables is a priori known and is fixed in the optimization process. For example, the 1

optimum thickness distribution is considered over the whole plate, the domain of the thickness and deflection functions is the whole middle surface of the plate that is not altered by design optimization. Shape Problem On the other hand, the shape problem is defined on a domain which is unknown a priori. For certain restricted problems, the shape of the domain may be described parametrically using a finite number of parameters with appropriate basis functions such as Bezier splines, but the state equation is defined on a "variable" domain. The optimum shape of the domain is obtained so as to minimize e.g. the mean compliance of the structure, while the equilibrium( i.e. state ) equation is satisfied on such a domain to be determined. Thus, if the finite element method is applied to solve the equilibrium equation, its discrete model must be developed at each iteration of optimization process according to the change of the domain. Layout ( Topology ) Problem Ambiguity of the layout ( topology ) problem is substantial in its mathematical description. In this case, only the known quantities are applied loads, desired support conditions, various design restrictions, and possibly the volume of a structure constructed. The purpose of optimization is to find the optimal layout of a structure so that given applied loads are transmitted to desired supports in a specified region using a given amount of material while equilibrium and design constraints are satisfied. In this problem, not only the physical size but also the shape of the structure are unknown, and further we cannot define these geometrical quantities with appropriate parametric representation! In order to define a form of layout, topology of the domain ( structure ) must be specified, but it is not known a priori. A formal description of the structural optimization problem may be defined by minimizing ( or maximizing ) a functional with respect to the design variables consisting of [size, shape, topology) subject to (the state equation, constraints on the state and design variables}. A typical mathematical setting is given by Minimize f(d,u) d subject to ueV: a(u,v)=L(v), VveV gd (u)< gugdax gd (d)<gd (1a) (1

Here,f(d,u) is the "cost" function, d and u are the design and state variables, respectively, "a(u,v)=L(v) Vv " is a generalized ( weak ) form of the state equation, and gu and gd represent the constraints on the state and design variables, respectively. As a concrete example, we shall define a shape problem for a linearly elastic plane structure to determine the optimal shape h(x) which is parametrically represented by a set of basis functions { hl(x), h2(x),..., hn(x) }, i.e., n h(x) = Xdihi(X) i=l In this case the design variable is defined by d = dl,....., dn }. If the constitutive relation of the material of the structure is given by O ll1 D2 D13 E1 E01 ] (F 1T2 =1 2 D22 D23 ~2 - ~02 D(E-0) T12J D13 D23 D33 7 Y12 1012J, (2) where a and ~ are the stress and strain vectors in the contracted notation, respectively, and go is a specified initial strain, the generalized form of the state ( equilibrium ) equation is given by the principle of virtual displacement: uEV: a(u,v)= L(v) VveV (3) Here a(u,v) = E(v)T D~(u)dM + fvTkudl (4) L(v) (= e(v)TD~od + VD + bd T (kg + t)dF (5) V={v={vv2}:vi E H1()} (6)

*(vy = {e1(v) ~2(v) 712(v)}= { ~_ __ 2 MVI 1 )aX1 aX2 aX1 - aX1 (7) b is the applied body force in the domain Q, H1(Q) is the Sobolev space defined on Q, and r is the boundary of the domain Q. Here the generalized boundary condition 1{ 01(u) _ nj 0.2 2(u:_ kl k,21(2u1l +{12 0 n2 nl] c2(U)j kl2 k2l({u2} {g2} {2} [12(u)J i.e. Na(u) = -k(u - g) + t (8) is assumed on the boundary r, where n = { ni, n2 } is the unit vector outward normal to the boundary, k is the elastic tensor of spring distributed, g is a given specified displacement vector, and t is the applied traction on the boundary. If the "cost" functionf is defined as the mean compliance of the elastic structure, this is given by f(d,u) = L(u) (9) If the shape optimization is restricted by the upper bound Do of the volume of a structure, the constraint on the design variables d is specified so as to 2 < o0. However, the explicit representation of this volume constraint is a rather difficult task in the parametric definition of the shape of a structure. Furthermore, the domain K2 and its boundary r are functions of the design variables but they are written in very indirect manner. Thus, finding the derivative of the cost function and the constraints with respect to d is not obvious. In other words, existence of the optimum shape design may be shown by using compactness argument, since only a finite fixed number of design variables, but computing the optimum shape could be very difficult. In addition, a different choice of the basis functions hi(x) may yield a significantly different optimum shape. There is no assurance that a unique optimum shape can be obtained by different choice of the basis functions. 2. Difficulties in Structural Optimization 4

If the structural optimization problem can be defined using a finite number of parameters together with well-behaved appropriately chosen basisfunctions, we do not have much difficulties in showing existence of a solution by standard compactness argument. Furthermore, if the cost function and the constraints are explicit functions of such a finite number of parameters, i.e., the design variables, computation of their sensitivity is straightforward, and then application of appropriate optimization algorithm yields the optimum without much difficulty in computation. However, there are so many optimum design problems outside of this category. For example, let a typical sizing problem be considered to determine the optimum thickness of a linearly elastic so as to minimize the mean compliance, or that maximizes the first eigenvalue of the free vibration problem of a plate. If the thickness is represented by a fixed finite number of C"- basis functions, the optimum can be easily obtained, see, e.g., Banichuk[1]. However, as Cheng and Olhoff[2] shown, a different selection of a finite number of basis functions yields the different optimum, and the true optimum does contains many discrete ribs in various size. We also see many discrete rib reinforcement in plates and shells in engineering practice. It is clear that infinitely many basis functions are required to represent such discrete ribs in various size and location. Thus, it is unrealistic to define the thickness optimization problem using a fixed finite number of parameters and basis functions. From the examination of the principle of virtual displacement that is a weak form of the equilibrium equation, it can be also determined that convergence of parametric representation of the thickness hn = spank hj,..., hnI is expected at most in L""(2) as n goes to co. This is implied by the fact that the coefficient tensor of the differential operator need be merely in L'(92) in the weak form, and then the assumption of the isotropic plate could be insufficient, since convergence of the thickness of a plate 2(i - 1) 2iha1 2(i-1) < x<2ij-1 h" (xy)=i 2n 2n ad 2n 2n hb if 2i-1 < x < 2i and < < 2n 2n 2n 2n where ij = 1,......, n, as n goes to oo, yields an orthotropic plate. That is, we have to state the equilibrium equation in a "relaxed" form using the assumption of an orthotropic plate that possibly possesses variety of microstructures instead of restricting our attention only for too well-behaved isotropic plate that assumes smooth variation of thickness as well as Young's 5s

modulus and Poisson's ratio without having many microstructure. Lurie, Fedorov, and Cherkaev[3] studied this nature mathematically using the theory of G-closure, and provided a clear answer how this problem should be stated and be solved, see also Bends0e[4]. In other words, even in the sizing problem, if the design variable is distributed and is expected in Ln space, the state equation must be extended to the relaxed form that can allow existence of microstructure. In the shape problem difficulty we encounter is not only for the convergence of parametric representation of the shape of a structure, but also for representation itself. Furthermore, there is no assurance that a singly connected domain is optimal. In fact, we know there should be holes in a structure for a certain case from our experience. But, we do not know the number of holes and their location as well as shapes a priori. Can we explicitly represent this situation parametrically with a finite number of preassigned basis functions? The answer is negative. Thus, the shape problem can only be solved with many additional restriction such as a specified number of holes, i.e., topology. In this case only sub-optimal solutions can be obtained, and we do not have any guarantee that change of the number of holes would not yield significant difference on the optimum shape obtained. Even with this restriction to a fixed topology, there still exists considerable difficulty in the shape problem when it is solved by using the finite element method. Because of a finite number of parametric representation of the shape, the sensitivity analysis for defining a sequential linearized problem is rather straightforward, but solving the state equation using the finite element method requires development of a discrete finite element model of the structure. If the initial shape of the domain is considerably different from the "current" one, development of a discrete model of the structure at each optimization step becomes a significant task, since it must be done automatically without any manual help nor terminating computation for optimization. Since the design variable is defined on a fixed domain, the sizing problem is free from this modeling difficulty. An automatic finite element mesh generation must be introduced to solve this into an optimization system. It is needless to say that what kind of difficulties we can find in the layout (topology) problem, since it is a complex combination of the sizing and shape problems if the topology of a structure is fixed with a finite number of parameters. Again it is not known whether or not the optimal topology can be represented by a finite number of parameters. Only available approach is to set up an initial set of truss joints and to generate additional ones within the limited total volume of the structure using, e.g. dynamic programing methods, see Palmer[5]. In this case, the type of structure must be predetermined, for example, as a truss 6

structure that cannot transmit the bending moment at joints. But, this choice need not be the best again. Furthermore, we can deal with only a fixed finite number of joints. There is, again, no assurance that a convergent optimal structure can be obtained as the number of joints is increased. 3. A New Approach Based on the Homogenization As shown in above, there are many difficulties in structural optimization. A common difficulty in the sizing, shape, and layout problems is convergence of a finite number of parametric representation. This means that once a method is introduced to solve the convergence problem for one of these problems, it may be applicable for the others, and that there is a possibility of existence of a method to solve these three problems at once. To solve the shape problem, a common approach is to define the shape using certain functions and to vary it to achieve the optimality. Variation of the shape generates substantial difficulty in computation. Thus, we should abandon this approach. That is, the shape should not be defined by a set of functions, and the problem must be solved in a fixed domain as for the sizing problem so that finite element(difference) model need not be renewed during optimization. Furthermore, since the topology is unknown a priori, it should not be specified. From the very beginning, an "infinite" number of holes should be prepared at "everywhere" to solve both the shape and topology problems. Does there exist such a method? To answer to this, we shall recall that convergence in the sizing problem is solved by applying the concept of the homogenization, i.e., Gclosure by Lurie, Fedorov, and Cherkaev[3], and their method is applied to solve the thickness optimization problem by e.g. Bends0e[4]. Now let us also recall the paper by Murat and Tartar[6] presented in 1983 on optimality conditions and homogenization. In this paper, Murat and Tartar stated "Where the variable is a domain, the computing of variations was done as early as 1950 by Hadamard, by pushing the boundary along the normal and then computing the induced variation of the functional. Turning this idea into theorems is not easy, and moreover already one can see another defect of this method: a given domain is compared only with domains of the same shape; it is impossible to make a hole inside the domain or to add a few small pieces far away by means of this technique...... The real difficulty lies in the fact that the set of domains, i.e., characteristic functions, does not possesses natural paths from one domain to another: there is no manifold structure that enables us to use classical derivatives...... In the problem where the variable appears as a domain and some partial differential equation is involved, there is another phenomenon that 7

we discovered ten years ago ( it was later called homogenization ): generalized domains appear which are the analogue of a mixture of two different materials and the effective properties of these mixtures have to be understood ( they are not obtained by averaging certain quantities in more than one dimension )." The very similar but more directly related idea for structural optimization is given by Kohn and Strang[7]: "The need for relaxation is reflected in the design problem by the possibility that there may be no optimal design. Though initially surprizing, this phenomenon is easily understood: it is sometimes advantageous to perforated parts of fQ by many fine holes with a suitably chosen geometry. If the optimal characteristics can be realized only in the limit as the scale of the perforation tends to zero, then the design problem has no geometrical solution - rather, a solution exists only in a suitably generalized class of designs, which allows the use of composites obtained by perforation at some points of Q2...... We first extend the design problem, allowing for the use of composite materials. Then we choose the best composite at each point." It should be now clear that there is a method to solve quite wide range structural optimization problems using the notion of homogenization after extending the design problem by allowing possibly perforated structures. Since perforation is expected in the microstructure, and since the degree of perforation as well as location is not specified a priori, arbitrary shape and topology of the structure may be well represented. It is certain that we may end up a perforated structure as the optimum because of the relaxation ( i.e. extension ) of the original design problem that is restricted to the use of only a solid nonperforated material. Because of the convergence property in the sense of homogenization, a precise mathematical theory should be able to furnish to the relaxed i.e. extended problem of structural optimization. In short, application of the theory of homogenization is an answer to the optimal design problem in structural analysis. 4. Relaxed Optimal Design Problem We shall describe a relaxed design problem using microscale rectangular holes to perforate a structure, see Bends0e and Kikuchi[1]. Suppose that the total volume of microscale holes is specified in a given design domain Q, that is, the volume Qs of "solid" material distributed in the design domain is specified. For simplicity, the design domain is plane so that plane stress analysis is sufficient to compute displacements and stresses, while the shape of microscale voids is assumed to be rectangular as shown in Figure 1.

e a - I b 1 1 1 Figure 1. A Unit Cell Describing the Microstructure with a Rectangular Hole Rectangular holes are chosen because they can realize the complete void ( a=b=l) and solid (a=b=O) as well as generalized perforated medium ( O<a<1, O<b<l ). If circular holes are assumed, they cannot reach to the complete void, and then are not appropriate to our purpose. It is noted that there are other choices to represent microscale holes such as using a generalize ellipse defined by (yl/a)n+(y2/b)n = 1, where a and b are the principal radii of the ellipse and n is the power defining the shape. However, in order to develop the complete void in the unit cell both a and b must be 1 as well as n goes to infinity for the generalize ellipse, while only a and b need be 1 if the rectangular hole is applied. Thus, using rectangular holes is simpler. Since holes are rectangular in the unit cell that characterize the microstructure of a perforated medium for the design problem, their orientation is important in the macroscopic problem for stress analysis. Indeed the anisotropic elasticity tensor in the macroscopic problem strongly depends on the orientation of microscale holes. Thus, the sizes a and b and the orientation 0 of microscale rectangular holes must be the design variables of the relaxed problem. Suppose that a, b, and 0 are functions of the position x of an arbitrary point of a macroscale domain of a linearly elastic perforated structure Q in the two-dimensional Euclidean space RI2: a=a(x), b=b(x), and =0-(x). Functions a, b, and 0 may not be so smooth, but we shall assume they are sufficiently smooth, for example, a,b,0E H (e ). Assuming that a periodic microstructure characterized by a(x), b(x), and 0(x) exists in a small neighborhood of an arbitrary point x in Q, and assuming that such microstructure at x 9

need not be the same with the one at a different point x*, see Figure 2 which shows a schematic setting of varying microstructures, a homogenized elasticity tensor EH(x) is computed in order to solve a macroscopic stress analysis problem of a perforated structure. Figure 2. Assumption of "Continuous" Change of Microstructures The homogenized elasticity tensor is computed by solving the problem defined in the unit cell in which a rectangular hole is placed: Find the characteristic deformations x(kOl){=(kl) kx(kl) Vy satisfying E __y dY,VVEVy ijmn dzjkl)yjfwy dY,VVE, (10) i,j,m,n=l dyn dyJ, I,' where Vy = (v=(vl,v21: vi EH1(Y), vi is Y-periodic in the unit cell Y, i=1,2}, and the unit cell is defined by Y=(-1/2,1/2)x(-1/2,1/2). The elasticity tensor E is chosen either for the plane stress or for plane strain problem depending on the structure to be designed. The elasticity tensor E is zero if y is located in the hole, and coincides with the one of the "solid" material that is utilized to form a structure if y is outside of the hole. It is noted that Young's modulus does not affect to the optimum perforation of the relaxed design problem while Poisson's ratio may imply change to the optimum. It is also noted that the material of the 10

"solid" portion need not be isotropic. For example, if the layout of a fiber-reinforced laminate is considered, the elasticity tensor of the "solid" portion must be the one for the laminate, and is anisotropic. After obtaining the characteristic deformations X(kl), the homogenized elasticity tensor IEH is computed by kl= E |(jkl 3imn y (11) Since the sizes {a,b) of rectangular holes are functions of the position x, the homogenized elasticity tensor EH varies in Q. This means that the characteristic deformations must be obtained everywhere in the design domain U2. Solving the unit cell problem (10) at everywhere is unrealistic. Thus, we shall solve (10) for several sampling points {ai, bj: ij=l,....,n) of the sizes {a,b} of rectangular holes, where 0 < ai < 1 and 0 < bj < 1, and we shall form a function EH = EH(a,b) by an appropriate interpolation. The last step of obtaining the elasticity tensor for stress analysis of the macroscopic perforated structure is rotation of EH by the angle 0. Defining the rotation matrix R by cos 0 -sin 0 R()= sin0 cosO the elasticity tensor EG for stress analysis is computed by 2 EiG= EHKL(a,b)R;i(O)Rjj(O)Rk(O)RlL(0) (12) I,J,K,L=1 for ij,k,l = 1, and 2, at arbitrary point x in U. It is clear that EG is a function of the sizes {a,b) and the rotation 0 of microscale rectangular holes, i.e. the design variables. Now these EG define the D matrix in the principle of virtual work (4) for stress analysis of the perforated structure. For the optimization problem, d={a,b,)}, and the upper bound Q2s is specified that is the total volume of "solid" material forming the perforated structure, and it defines a constraint on the design variable: 11

(l - ab)dQ < gts < Q J(1-a,<2 (13) where Q is the total volume of the design domain in which an optimum perforated structure is placed, that is, the layout of a structure is given. The main feature of this setting of the optimum layout of a structure is that size, shape, and topology are represented by the distributed design variables d= a,b, 01, and that the problem is defined on a fixed design domain Q. In other words, the layout problem is writen as a sizing optimization problem. Thus, a geometric model of the finite element analysis defined at the beginning need not be modified during the optimization process, and if the design domain is sufficiently simple, its discrete model by finite elements can be developed easily without introducing any special spline functions to define the shape of a structure designed in the design domain. Only the change in optimization is the design variables d=(a,b,0), that is, the D matrix in (4), and the domain Q is independent of the design variables. Thus, sensitivity of the cost function and the constraints, and then the optimality condition of the minimization problem (1) can be obtained without any difficulty. Once sensitivity and the optimality condition are obtained, introduction of an optimization method is straightforward, despite that the number of discrete design variables obtained from d=fa,b,0} would be fairly large. For example, if the design domain gQ is decomposed into N finite elements, and if d=(a,b,01 are represented by constants in each finite element, we generate 3N discrete design variables. Details of computational procedure including an optimization method to solve the design problem (1) can be found in Suzuki and Kikuchi[9], and will not be discussed here. Since the design variables a, b, and 0 are distributed functions on the design domain U, they must be approximated by the discrete ones to solve the optimization problem (1). To this end, these are approximated by constant functions in each finite elements of the discrete model of the design domain. This yields that the microstructure is assumed to be constant in an element, and then the homogenized elasticity tensor is also constant in each finite element. From the experience we have so far, the approximation by piecewise constant functions can provide satisfactory results. Because of this choice, the optimality condition of the discrete approximation of the constrained minimization problem (1) can be obtained easily, and then an iteration scheme based on the optimality criteria method can be derived from this. In the following examples, we shall assume that Young's modulus and Poisson's ratio of the isotropic "solid" material in the unit cell are given as E=100GPa and v=0.3, 12

respectively. The homogenized elasticity tensor EH is computed for 6 x 6 sampling points in the design variables a and b, respectively, and is interpolated by the Legendre polynomials. The homogenization problem (10) is solved using 12 x 12 finite elements in the unit cell as shown in Figure 3, where a rectangular hole with a=1/3 and b=2/3 is assumed. The components EH1111, EH1122, and EH1212 of the homogenized elasticity tensor EH are computed as shown in Figures 4, 5, and 6, respectively. IL I 2/3..>>.. l......... H 1.,............i................. _.....,... Figure 3 A finite element model of the unit cell with a rectangular hole 100 o8 25 0.2: )o. 4 ~~~~~~0.4~01 0.8

Figure 4 Homogenized elasticity constant EH1 1l 1 with respect to the size a and b 30 2C 10,:;.cs4fPBwre0.8 0 6. 0 0.2 4w0.4 0.4 - 0.6 ]_Fn~p0.2 0.8 Figure 5 Homogenized elasticity constant EH1122 with respect to the size a and b 3 0 2 0 0 6 4 10 0J ~L::~JRIA[B~~~~.6 0 0.2 ~ e "PRRC~Ps~abTT).4q 0.4 0.6 ~ ]D.2 0.8 Figure 6 Homogenized elasticity constant EH1212 with respect to the size a and b 5. A Verification of the Relaxed Design Problem We shall examine the present method to solve a simple structural design problem, the exact solution of which can be obtained analytically using a simple structural model such as a truss or a beam. To this end, a two-bar framed structure shown in Figure 7 is considered. 14

H - HL Figure 7. A Two-Bar Truss r the fixed values of the applied load P and the horizontal length L of the frame, the:imal height H is determined by minimizing the mean compliance. If the cross section of L1It truss-bars is unchanged and is rectangular with the unit width, the optimal height H is obtained as H=2L under the constant volume constraint. We shall now obtain the optimum layout by minimizing the mean compliance with a single load P applied at the center of the right edge of the design domain, shown in Figure 8, that is larger than the size Lx2L, where L=lOcm, so that it can contain the optimum layout obtained by a two bar truss structure inside of the design domain. The relaxed design problem (1) is solved for the discrete finite element model of the design domain consisting of 40 x 96 uniform rectangular finite elements. As shown in Figure 9, the present homogenization method forms a two-bar framed structure if the amount of the "solid" material is specified to be rather small without assuming any type of structures such as trusses, beams, and frames. It automatically identifies the optimum structure and its configuration. Furthermore, if details are examined, it is clear that the size of the bar is also determined. As the point is close to the end portion where a load P is applied, the height of the cross section of the bar becomes is small. 15

24 2 Design Domain Figure8 DesignDomain andthe Figure 9. Optimal Relaxed Design Support/Loading Condition (Volume =*40 cm2) Figure 10 The "optimum" configuration with constraints in the design variables........din Codto' Vouegn0 m FIgur 10Te"pImu"cnIgrtowihcntatsntedsgnvabs

Now we shall solve the same problem by assuming that the design variables { a, b, 0 1 are restricted to be a=b and 0=0, that is, the microscale hole is square and is not rotated during the optimization. Figure 10 shows the configuration obtained under this restriction. Its configuration is similar to the one shown in Figure 9, but it cannot reproduce the two-bar framed structure, and its compliance is also much higher than the case of no constraints on the design variables. The main reason of this difference is caused by the restriction on the angle 0 of the microscale hole. Since holes cannot be oriented to the direction of the principal stresses if the angle 0 is restricted, the mean compliance is not minimized. This suggests that the simplistic idea of removing finite elements whose stress is low does not work well to obtain the optimum layout of a structure, because the angle 0 is not involved in the design problem. 6. Layout of a Beam Structure We shall solve a layout problem of a beam structure using the above formulation and make a comparison to the result of the beam height optimization that can be found in standard textbooks and monographs of structural optimization. If the problem is regarded as a sizing problem that finds the optimum thickness of the rectangular cross section of a beam with a fixed width, see Figure 11, this can be solved easily by applying a standard technique of structural optimization. Figure 11 Thickness optimization of a clamped beam subject to a point load P Minimizing mean compliance subject to equilibrium equation and volume constraint is equivalent to 17

max min Ebh3-P (14) subject to Jbhdx - V O (15) Here w is vertical displacement, h is height of beam that is to be designed, and b is width of beam that is fixed, E is Young's modulus, P is applied point load, and V is upper bound on volume. Using the Lagrange multiplier method, necessary condition for optimum can be derived as h2( )= constant (O < x < L) ~~\~d~~~X)2~~~~~ ~(16) It follows from the boundary conditions that the optimal thickness is given by h(x) = 3- 1- x 2bL (L / 4) ( < x L / 2) (17) Symmetry is applied for the rest of the beam: L/2 < x < L, and the optimum thickness is shown in Figure 12. Figure 12. Analytical solution of beam thickness design The optimal thickness involves two hinges at x=L/4 and 3L/4, i.e., h = O at x=L/4 and 3L/4, while the maximum thickness is obtained at x=O, L/2, and L. The gradient h' of the thickness is oo at a=L/4 and 3L/4. 18

_ Load t1 eDesign Domain 4 20 Figure 13. A Design Problem of a Beam Structure Now, if this problem is solved by the present method by specifying different volume of the solid material, see Figure 13, while the design domain is fixed, the optimal layout of a structure is obtained as shown in Figure 14. It is clear that hinges appear at x=5 and 15, and is the same to the case of the standard thickness optimization. Since it is possible to generate internal holes inside the structure, holes are "naturally" formed to increase the stiffness against bending moment. i.e. forms sandwich beam. In most of thickness optimization published so far, internal holes are not presented. Thus, the method introduced in this paper can provide far more sophisticated optimal structures. If the volume of solid material is reduced, very truss like structures are formed. Because of the restriction on the design domain, reinforcement can be placed only inside the rectangular domain. Large reinforcement is observed in the vicinity of the fixed end points as well as the center at which the point force is applied. It is clear that the present method can provide the shape and topology of the overall structure as well as the optimal sizes of members of a "truss" structure, and that there have been no other methods which have comparable capability to the present one. Volume Constraint 10 Volume Constraint 15 19

Volume Constraint 18 Figure 14. Optimal Layout of a Beam like Structure 7. Convergence of the Finite Element Approximation Next issue to be discussed is whether the shape and topology, i.e., the layout of the structure obtained as the optimum in the relaxed design problem, converges to the unique design as finite element meshes are uniformly refined, while other conditions are fixed. Importance of this test can be recognized by the observation in Cheng and Olhoff that shows if the relaxed approach is not applied, the optimal rib distribution is highly mesh dependent. That is, refinement of the discrete model may lead completely different optimal solutions. 16 Design domain 10 Figure 15. Design Domain for Bending of a Short Cantilever To check convergence property, let us solve the relaxed optimization problem (1) for a short cantilever subject to the vertical force at the free end, see Figure 15. For the volume 60 of the "solid" material, (1) is solved by using 32x20, 48x30, 64x40, and 80x50 equal size uniformly divided rectangular finite elements covering a rectangular design domain Q. Applying the same homogenized elasticity tensor to the first example, the optimal configurations are obtained as shown in Figure 16. It is clear that the optimal configurations are convergent as the size of finite elements is reduced. Even every coarse meshes can 20

provide sufficient idea of the topology and shape of the optimal structure. A very truss-like framed structure is built if the amount of the solid material is sufficiently smaller than that of the design domain. Each member is even straight. If the amount of the solid material becomes large, the optimal design may not be a truss-like structure. Curved frames can be generated, and more continuum-like shapes are formed. In this problem, despite of the relaxation that allows perforated composite at everywhere of the design domain, the optimal structure is not perforated at all. In other words, the relaxed problem provides the "classical" optimal solution to the design problem using only solid structural members. Mesh 32x20 meAsh 48x30 Mesh 64x40 Mesh 80x50 Volume 60 oohune 60 Figure 16. Convergence of the Optimal Configuration ( fs = 60 ) Figure 17 presents convergence of the mean compliance as the size of finite element goes to zero. Monotonic decreasing of the mean compliance is clearly observed. We shall now examine the stress distribution of the optimum layout of a structure obtained by minimizing the mean compliance. Since the mean compliance is a global quantity related to the strain energy of a structure, it cannot control the stress locally. In other words, the stress cannot be bounded by a given value at every point of a structure. 21

Figure 18 shows the optimum layout and the distribution of the Mises stress for the three different volume of the "solid" material. Here the design domain and the loading/support Convergence Test 80 78 [ 1/MC c. 77 L. 76 75, 0.0 0.1 0.2 0.3 Squares of the Mesh Size Figure 17. Convergence of the Mean Compliance condition is the same to the one in Figure 15. It is clear that the fully stressed condition is achieved in most of the identified structure except the neighborhood of the support/loading points, especially if the volume of the "solid" material becomes small. In this sense, minimizing the mean compliance tends to yield a fully stressed structure in the present homogenization approach because not only the shape and size but also the topology of a structure can be arbitrarily changed in the present formulation. However, if the optimum layout is not a "structure" such as a truss, beam, and frame, then the stress distribution is not unidirectional, and is constant. Thus, the fully stressed condition cannot be related to the mean compliance. 8. Layout of a Triangular Membrane So far only rectangular finite elements are applied to solve the optimal design problem (1) and have obtained rather discrete truss-like structures. Now let us examine a case that non-rectangular finite elements must be applied to form a discrete model, and that the present approach can also provide smooth shape of a structure if the volume of the "solid" portion is large enough in the design domain. To do this, let a "triangular" design domain is considered as shown in Figure 19. 22

..... Figure 18 Optimum layout and Mises stress distribution of a structure for the three different volume of the "solid"' material 23

Figure 19 A Triangular Membrane Subject to Tension Forces Applying symmetry condition, only 1/6 part of the design domain is modeled by the finite element method. The boundary condition and physical dimensions are specified in Figure 20. The domain is discretized by 40 x 40 meshes which are not orthogonal. Using this example, we shall examine transition from the shape to the topology optimization. The volume of the "solid" material is varied in a wide range from 50 cm2 to 180 cm2, while the area of the 1/6 of the original domain is fixed to be 317.5 cm2. 30" Figure 20 Design Domain for the Optimization(1/6 portion of the original domain) As shown in Figure 21, if there are sufficient amount of "solid" material, the optimal structure is simply connected without having holes inside of the domain. The obtained shape 24

is also very smooth that is almost same quality to the one obtained by the boundary variation method which is traditionally applied to solve the shape optimization problem. When the amount of "solid" is reduced, holes starts appearing interior of the domain. It is also noted that the loading portion becomes porous when the amount of "solid" material is really small in order to diverse the applied traction uniformly in all the directions equally, while the three compressed link bars in 30,150, and 270 degree directions, are also porous. Thus, if these porous bars are replaced by solid bars maintaining the "solid" material volume, the size of bars becomes very small. Thin solid bars which are in tension, are allocated in 90, 210, and 330 degree directions in the optimal structure. The outer frames are completely solid. Transition from the simply connected domain to the multi-connected domain occurs at 9s=130 cm2 - 140 cm2. It is also noted that from our experience the perforated composite appears only in the vicinity of the portion where distributed traction or body forces are applied. If a finite number of discrete "point" forces is applied, the optimal structure by the relaxed problem tends to generate a very discrete solid structure rather than a perforated composite structure. In the present example, a distributed load is applied on the three end surfaces, the optimal structure is perforated in the vicinity of distributed loads while most of the rest of the part of the structure forms a discrete solid structure, i.e. a truss like structure without any perforation. Instead of a distributed load applied along a portion of the boundary, if the displacement is specified over there, perforation does not also occur in the vicinity of such boundary. In this sense, despite of possibility of perforation and extensive composite type microstructure, the optimum structure is rather discrete, and it tends to avoid perforation. Voimi 70 25

Volume 180 Volm 120 Figure 21 Optimal Configurations for Various Volumes of "solid" Material 9. Optimal Layout of a Structure Subject to Multi-Loads In practice, a structure must often be designed under the multi-loading condition. That is, a structure must be safe for a set of various loadings. How can we extend the present homogenization method to solve such multi-loaging optimization problem? Here we shall introduce a "non-standard" way to define the problem, although most of multi-loading problems are defined by introducing an appropriate linear combination of the cost functions for each loading. To this end, if the mean compliance is regarded as the cost function, we have the following relation: L(u)=a(u,u). (18) Thus, the optimization problem (1) for a single load can be written as Minimize a(u,u) d subject to ueV: a(u,v)=L(v), VveV f,"(1-ab)dQ<Q2 (19) Now, let the principle of virtual work be represented by 26

ui E V: a(ui,v) = Li(v), Vv E V, i = 1,.....,m (20) for the i-th loading case Li(v)= e(V)TDeoid + vTbidQ + VT (kgi + ti)d (21 (21) where i = 1,....., m, and then let a special functional am(u, u) be defined by am(u,u) = fMax Ce(ui)T D(ui))}dd + IMax {i Tkui}dF Ji=1.....,m i=1.......m (22) Using this special functional, we shall define the optimization problem Minimize am(u,u) d subject to ui EV: a(ui,v)=L(v), VveV,i=l,...,m Ia(1-ab)dn-<Qo (23) This formulation yields a "conservative" optimum design in the sense that the resulted structure is safe for all the loads applied independently, and yields a sub-optimum solution of the standard multi-lodaing optimization problem: Minimize Max. a(ui,ui) d i=1.......m subject to i eV: a(ui,v)=L(v), VveV,i=l,...,m I(1-ab)dn<Qn (24) Despite of disadvantage of the possibility that only a sub-optimum solution may be obtained from formulation (23), this problem can be solved very similarly to the single load case with minor modification, while solving the standard multi-loading optimization problem (24) is rather difficult. As an example of multi-loading, we shall solve the design problem shown in Figure 22 in which a structure pin-supported at two different size circular holes is subject to three loadings, tension P1, bending downward P2, and bending upward P3. Furthermore, we 27

shall restrict design to the case that parts of the constructed structure cannot pass through the center portion of the design domain. Design Domain P3=P rl=2R r2=R H=4R...... ________ _ L P1 =P No Material Distribution P2=P. - L=10R O Figure 22. Design Domain, Restriction on Design, and Multiple Loadings The relaxed design problem for the multi-loading (23) is solved using the similar homogenization method for the single load case (1), and the optimum layout is obtained for two different volume of the solid material, as shown in Figure 23. Figure 24 shows the optimum layout for the case that the design restriction is ignored, i.e., that structure can be placed along the center horizontal line connecting two pins. Significant design 28

Figure 23 Optimal Designs for the Multi-Loading Problem (with the Design Restriction ) Figure 24 Optimal Design of the Multi-Loading Problem ( without the Design Restriction ) change is observed in this problem with/without design restriction. If design restriction is abandoned, the optimum structure contains elements of Michell truss, i.e., is formed by nearly orthogonal truss/beam network. Design restriction yields rather stiff two bars in most of the portion, and they are supported by two small scale Michell truss-like substructures. It is natural that the "mean compliance" of the case without design restriction is lower than that with design restriction. 10. Extension to Plates Subject to Bending Since the homogenized D matrix for plane stress has been obtained in above, this can be applied to formulate plate bending problems by replacing the D matrix for isotropic plate. However, it must be realized that this simple replacement may not yield an orthotropic plate obtained by the homogenization process by assuming microstructure due to infinitely refined 29

ribs as studied by Cheng and Olhoff, and Bendsoe[4]. The plate model obtained by simple replacement of the D matrix computed for plane stress problems is an intermediate model of an isotropic plate and an orthotropic plate with microscale ribs, and it yields rather discrete "finite" size ribs in the optimum layout of a plate-like structure as shown in the following example. Since microscale ribs might be only interested in theory, the present formulation due to simple replacement may have much practical meaning for design of plate structures. As an example of the optimum layout of a plate subject to transverse loads, let a square design domain be considered subject to a concentrated force at the center and a uniformly distributed pressure force. Suppose that a plate is simply supported along the boundary of the square design domain, and suppose that there are no other design constraints except the total volume of the "solid" material. The optimum layouts are obtained as shown in Figure 25. Majority of the solid portion forms a cross-like structure, and is reinforced by thin beams along the diagonals.in both loading cases. It is also noted that hinges are formed at the quarter points of the diagonals where the slip lines are passing through in the limit analysis of a plate. As shown in Figure 25, the optimum structures are not perforated in most of the design domain.'V Figure 25 Optimum Layout of Plate-Like Structure Subject to a Transverse Point Force at the Center and a Uniformly Distributed Pressure 11. Extension to a Shell-like Structure 30

Extension of the present method can be also made for shell-like three-dimensional structures similarly to plate-like structures. As an example, we shall consider a layout of reinforcement of a thin shell supported by three hinges and is subject to a uniformly distributed pressure. As shown in Figure 26, three edges are straight, while a curved edge is defined by a sin function. Interior portion of the shell is interpolated linearly. hinge support hinge support Figure 26 Shell Structure Reinforced The design problem is to find the optimal reinforcement by adding material on the already buildup thin shell structure, say 1mm thick shell. As shown in Figure 27, the optimum reinforcement is very discrete, i.e., there are little portion where perforated composite is assigned to minimize the mean compliance. Very clear rib reinforcement is obtained. Basically two main ribs are generated along the lines connecting the three support points, while two rather narrow ribs are formed along the flat boundary edges to increase rigidity. If the shell is curved, the applied force is decomposed into the membrane and bending ones, and then the curved edge can resist to transverse loads more than the flat edge. Thus reinforcement is required along the flat edges of the shell considered as shown in Figure 27. 12. Extension to Three-Dimensional Solid Structures The same concept is applicable to the layout problem of three-dimensional solid structures. In this case, we shall make perforation by microscale rectangular parallelepipeds with rotation about the three coordinate axes. Thus the design variables become six at an 31

arbitrary point x: d = { a, b, c, 01, 02, 03 1., i.e. three sizes and three rotations of rectangular parallelepipeds at each point of a design domain. The rest of the part is exactly the same to the case of plane problems. As an example, let us solve an optimum layout H ing Hinge Hinge Hinge Figure 27 Top View of the Rib Reinforcement problem of a structure subject to multi-loads, Px=Py=3Pz=1, at the free end of a rectangular parallelepiped design domain as shown in Figure 28. Figure 29 shows the optimum layout obtained by the present method. It is clear that there are many difficulty to find such a layout by applying the standard boundary variation method for the shape optimization of threedimensional structures by specifying the shape using appropriate spline surface representation. Even in this simple design domain and applied multi-loads, the optimum layout is rather complicated in the sense that it requires many patches of spline surfaces to develop its solid modeling, and then many control points must be introduced for shape optimization to vary the shape of a structure in optimization process. Although the number of finite elements and discrete design variables becomes fairly large in the homogenization method for three-dimensional solid structures, we are free from defining 32

Figure 28 Design Domain of a Three-Dimensional Solid Structure Figure 29 Three-dimensional description of the optimum layout of a three-dimensional solid under the multi-loading condition 33

,,..~ 1% "" r Y ~~~~~~~~~~~~~~~~~~~~~2 Y,~~r

.<~~~~~~~~~~L.~~~~~~~~~ wSK ok ark-.. Figure 32 Optimum layout in a cross-section perpendicular to the y axis the solid model of a structure and automatic decomposition into finite elements. Topology and shape are obtained as the result of optimization by accumulating "pixels" for the solid portion. The result in Figure 29 is obtained by a 34x14x8 uniform mesh, and then it involves 22,848 discrete design variables and 14,175 degrees of freedom in the finite element model for stress analysis. Figures 30,31, and 32 shows the layouts of the structure in its three orthogonal cross sections. It is clear that a hole is generated inside of the structure. 13. Construction of an Integrated System As shown in above the homogenization method can provide the optimum layout of an elastic structure subject to a single load as well as multiloads, and is a fixed domain method that makes it possible to solve the problem using a fixed finite element model during optimization. As the result obtained, the topology and shape of a structure are identified by the distribution of the microscale holes of perforation. Since these are represented by "gray scale of pixels" of a finite element model of the design domain, we have to recognize them by defining appropriate spline surfaces based on such gray scale map if function representation

of the shape must be obtained. Furthermore, it should note that the homogenization method introduced here is defined by minimizing the mean compliance only with volume constraint. However, practice of design optimization of elastic structures requires much more: the "cost" function minimized can be the maximum Mises stress with various other constraints on displacement, strain, stress, and requirement to be stable structures as well as constraints due to manufacturability and others. Thus, the method presented here cannot answer to all of the requirement in practice of design optimization. In order to fulfill such requirement, we may construct an integrated system of a topology-generating into the design optimization process with the following general scheme: Phase I: Generate information about the optimum topology for the structure. Phase II: Process and interpret the topology information. Phase III: Set up a model for detailed design optimization and optimize the design applying standard structural optimization techniques. This scheme is the basis for an integrated structural optimization system which is rudimentarily drafted in Figure 33. Phase I is accomplished by the homogenization method describe in above. As demonstrated, the method requires only loading and boundary conditions, a definition of the design domain (i.e., the domain where the designer wants material to be distributed in an optimal way), the objective function ( i.e. the mean compliance ) and the volume constraint. What we achieve here is an analytically derived approximate configuration, that corresponds to results of a "concept evolution" stage of the design process, based on satisfying the primary design requirement ( based on the stiffness of a structure ) from the structural viewpoint. For a practical design tool, the ability to refine the design and account for other requirements is necessary. The method-produces an optimum material distribution described by density, 1 - ab, data which can be represented as images. The density information is generally "noisy" and has to be processed and smoothed with computer vision tools, producing higher level representations that are easier to be interpreted than the original ones. Based on these smoothed pictures the designer can now generate an "initial design", choosing geometric elements (like straight lines, circles, and appropriate splines) that meet the manufacturing requirements. Any commercial CAD program featuring standard data 36

DTOPOLOG GENERATION MODULE Volume ConstraintsVaibeItalDmn Global Specifications n=1 Boundary Conditions HOMOGENIZATION Density Data integratediii~i~ ~jiii~i~i~iiiijjii s c lnoisy"......... iitii~jiiiiiiiiii...............ii IMAGE INTERPRETER MIODULE j interfaces may be used for thiis step. For example, if acomponent.is intended.for.automated assembly, certain.......i...i...d ehasm IMAGE PROCESSING ij i C~~~........... INTERPRETATION:;-:~~~~~:::::ii~iiiiii~ M o d ify V o lu m e...... i iiiiiii Constraint......~.~.~...~.~.~........~i~i~ii~ii3 set n-n +1......~l%~t~ Resul.......iiiiiiiiiiiiiiiiii t::::::i::::i:::::::::::i:::::::iclear & ~iiiiii~iiiiii~~'~ti M anuufactu ringg iiiiiiiiiiiiiiiri~iiiii~iisiiii~i....................Xii~iiiii Requi re ments 1ii:~::~::~;.~.:~.' ~ ~:~~................:~::~:: ~ on Topology ~i..........5ii~i:iiitiiiti Boundary D.......... tructral Mmbers DETAILED DESIGN MODULE DETAILED DESIGN MODEL DETAILED DESIGN COMPUTATION No ign Satisfactor Yes FIAL DESIGN Figure 33 Bas'i' flow cha-rt for the three-Phase integrated stmuturalo pti za on st m, interfaces mnv be used for this step. For example, if a component is intended for automated assembly, certain symmetries will substantially reduce the assembly cost, such as adding a non-functional hole to achieve easy part orientation during assembly. If the topological

change is sufficiently drastic, a new topology optimization may be necessary with a modified initial domain specification for the homogenization. At the same time a mechanical interpretation has to be done, i.e., structural elements (like trusses, beams or solids) have to be recognized. The steps "image processing" and "interpretation" constitute Phase II of the design process and are represented in the center part of the diagram shown in Figure 33. The final result of Phase II is an initial design represented by boundary data and information about structural members. Starting from the Phase II design, a detailed design model can be set up in order to perform Phase III operation: standard sizing and shape optimization using structural optimization methods. As the number of design variables and constraints is usually moderate in this step, general mathematical programming algorithms can be utilized to locate the solution, which makes the solution procedure adjustable concerning versatile objective and constraint specifications. This means that the designer may add specific model definitions that have been neglected in Phase I. The structural optimization procedure consists of a finite element program for stress and sensitivity analyses, several different optimization algorithm and standard pre- and post-processors that allow for a very flexible definition of design variables, objective functions and constraints. Since the topology of a structure is fixed to the one obtained by the homogenization method ( Phase I), and since its geometry is represented by a set of appropriate spline functions, Phase II involves only standard optimization procedures for sizing and shape optimization. As stated in Introduction, if the topology of a structure is fixed, the shape optimization problem can be solved by introducing an automatic mesh generation scheme to the system of optimization consisting of finite element analysis, design sensitivity analysis, and optimization methods, that is, one of the key features of the integrated system is the implementation of an automatic mesh generation method inside of the optimization algorithm. For the shape optimization problem, the domain discretized into finite elements may be largely deformed during the optimization process despite that the topology is fixed. If the initial finite element connectivity and the total numbers of nodes and elements are kept throughout the optimization, no matter how sophisticated mesh moving schemes are introduced, many finite elements would be excessively distorted. It is natural that this implies unnecessary approximation error resulting in non-reliable stress analysis for the optimization. In order to avoid this difficulty, implementation of an automatic mesh generation method to regenerate finite element discretization at each optimization step, is indispensable as pointed out by Bennett and Botkin [10] 38

Various algorithms of automatic mesh generation have been introduced in the last two decades after Fukuda and Suhara [11] first derived an algorithm for two-dimensional domain using triangular elements. This work is extended by Cavendish [12], and then the concept of automatic mesh generation becomes widely recognized among engineers as well as researchers. A significant conceptual development is made by Lo [13]. After this, automatic mesh generation becomes a commonly applied tool at least for two-dimensional and three-dimensional shell-like structures. A program for automatic mesh generation can be developed relatively easily, if the speed of computation is ignored. It is also noted that algorithms applicable both for triangular and quadrilateral elements are developed based on the QUADTREE concept by Shephard [14], and extensive success of their application in practical fields is reported. This is extended even to three-dimensional solid structures. In the present work, an algorithm developed by Tezuka [15] which is a variation of Lo's method, is applied. We shall briefly describe the algorithm of automatic mesh generation applied in the present study: (1) A given domain is divided into several subdomains in order to introduce different mesh "density", i.e., refinement. If a uniform discretization is required, only one single domain needs to be introduced. (2) Nodal points are placed on the boundary of each subdomain in order to define their geometry. Boundary segments of subdomains are represented by Besizer splines with given control points. If the total number of nodes on each segment is specified, nodes are automatically or manually generated on the segment. Some of segments become the design boundary of the shape optimization. (3) In each subdomain, compute the average mesh size from the nodes on the boundary, and then span a rectangular grid over the subdomain. (4) Identify grid points which are interior of the subdomain. Note that grid points being too close to the boundary should be excluded to generate finite elements which are not excessively distorted. 39

(5) Apply a triangulation algorithm from the nodes on the boundary using the interior nodes identified in the previous step. Guideline for triangulation is that a triangle having the maximum of the minimum inner angle should be chosen. (6) After triangulation, apply a smoothing scheme to have smooth gradual refinement as well as almost regular shape triangular elements. Step 4 and Step 5 could be time consuming, since search algorithms are extensively used. Especially, if a large number of nodes is introduced, these search algorithms are executed over "all" of the nodes. Sophistication must be introduced to reduce the necessary computing time for these searches. The mesh generator described in above is included into the pre-processor of a model of finite element stress analysis, and the coordinates of control points of the design boundary is modified by an optimization algorithm based on the result of sensitivity analysis. Only the fixed information is the control points of the design boundary segments during the optimization, while nodal points on such segments are changed at each design step, sensitivity analysis must be performed for the coordinate of the control points. Since the total number of the control points of the design boundary is rather small, say at most 10 - 20 for plane problems, and since its explicit relation between to the nodal coordinates of the finite element model generated can be expressed by the coordinate transformation of all the elements, sensitivity may be computed by the finite difference approximation. If semianalytical or analytical method is applied to compute sensitivity, we have to evaluate derivatives of the stiffness matrix for all the nodes in the model either approximately or analytically, and this is unrealistic. The best way in this setting is application of the finite difference method to evaluate design sensitivity with respect to the location of the control points. It is also noted that representation of design constraints on geometry such as the total volume of a structure and geometric restriction to define its configuration is rather difficult if spline representation is applied for the design boundary. Figure 34 shows a bracket design problem in plane. Phase I, the homogenization method for the optimum layout results the configuration shown in Figure 35, and its result is processed through image processing and interpretation of the geometry ( Phase II ), and implies the finite element model generated by the automatic mesh generation method. Applying SAPOP system, [16,17] the optimum shape is obtained as shown in Figure 36. 40

Details of the integrated system and an example described here can be found in Brenicker, Chirehdast, Kikuchi, and Papalambros[18]. 120 4.2 10 Thickness: 18mm Supports 70 Square Holes: 10mm x 10 mm F = 1780N B y / M = 45.25 Nm 4.2 4.2 4.2 Figure 34 Boundary specifications and the initial design domain Figure 35 Optimum layout by the homogenization method ( Phase I) 41

Figure 36 The initial model and the optimum shape of the structure ( Phase I ) Initial design model and its optimum Acknowledgement During the work, the authors are supported by NASA Lewis Research Center, NAG 3-1160, Office of Naval Research, ONR N-00014-88-K-0637, DHHS-PHSG-2-R01-AR34399-04, and RTB Corporation.,Ann Arbor. The authors also express their sincere appreciation to Professors Martin Bendsoe, Technical University of Denmark, Alejandro Diaz, Michigan State University, and Panos Papalambros, The University of Michigan, for their valuable discussions on this subject, and to Dr. Mike Bremicker and Mr. Y.K. Park for their help to complete some of examples described in this article. 42

References [1] Banichuk, N.V., Problems and methods of optimal structural design, Plenum Press, New York (1983) [2] Cheng, K.T., and Olhoff, N., An investigation concerning optimal design of solid elastic plates, Int. J. Solids and Structures 17 (1981) 305-323 [3] Lurie, K.A., Fedorov, A.V., and Cherkaev, A.V., Regularization of optimal design problems for bars and plates, Parts I and II, J. Optim. Theory Appl. 37(4) (1982) 4999-521, 523-543 [4] Bendsoe, M.P., Generalized plate models and optimal design, in J.L. Eriksen, D. Kinderlehrer, R. Kohn and J.L. Lions, eds., Homogenization and effect moduli of materials and media, The IMA Volumes in Mathematics and Its Applications, Spriger-Verlag, Berlin, 1986, 1-26 [5] Palmer, A.C., Dynamic Programing and Structural Optimization, in R.H. Gallagher and O.C. Zienkiewicz, eds., Optimum Structural Design, John Wiley & Sons, Chichester (1973) 179-200 [6] Murat, F., and Tartar, L., Optimality conditions and homogenization, in A. Marino, L. Modica, S. Spagnolo, and M. Degiovanni, eds., Nonlinear variational problems, Pitman Advanced Publishing Program, Boston, 1985, 1-8 [7] Kohn, R., and Strang, G., Optimal Design and relaxation of variational problems, Parts I, II, and III, Communications on Pure and Applied Mathematics, XXXIX (1986) 113-137, 139-182, 353-378 [8] Bendsoe, M.P., and Kikuchi, N., Generating optimal topologies in structural design using a homogenization method, Comput. Mechs. Appl. Mech. Engrg., 71 (1988) 197-224 [9] Suzuki, K., and Kikuchi, N., Shape and topology optimization using the homogenization method, Comput. Mechs. Appl. Mech. Engrg. too appear (1991) [10] Bennett, J. A. and Botkin, M.E., Structural shape optimization with geometric problem description and adaptive mesh refinement, AIAA J., 23(3), 458-464 (1985) [11] Fukuda, J., and Suhara, J., Automatic Mesh Generation for Finite Element Analysis, in Advance, in: Computational Methods in Structural Mechanics and Design, ( Ed. Oden, J.T., and Yamamoto, Y. ), UAH Press, Huntsville, Alabama, U.S. (1972) [12] Cavendish, J.C., Automatic Triangulation of Arbitrary Planar Domains for the Finite Element Method, International Journalfor Numerical Methods in Engineering, 8, 679-696 (1974) [13] Lo, S.H., A New Mesh Generation Scheme for Arbitrary Plannar Domains, International Journalfor Numerical Methods in Engineering, 21, 1403-1426 43

(1985) [14] Shephard, M.S., and Yerry, M.A., An approach to automatic finite element mesh generation, in: Computers in Engineering 1982, 3., edited by Hulbert, L.E., The American Society of Mechanical Engineers, New York, 21-28 (1982) [15] Tezuka, A, A Development of Automatic Mesh Generator with Arbitrary Geometry-Based Input Description, MS Thesis, Department of Mechanical Engineering and Applied Mechanics, The University of Michigan, Ann Arbor, MI, U.S., (1988) [16] Eschenauer, H., Post, P.U. and Bremicker. M., Einsatz der Optimierungsprozedur SAPOP zur Auslegung von Bauteilkomponenten. Bauingenieur 63, 515-526 (1988) [17] Bremicker,M., Eschenauer, H., Post, P., Optimization Procedure SAPOP - A General Tool for Multicriteria Structural Design. In: Eschenauer, H., Koski, J., Osyczka, A., Multicriteria Design Optimization. Berlin, Springer Verlag (to appear May 1990) [18] Bremicker, M., Chirehdast, M., Kikuchi, N., and Papalambros, P.Y., Integrated Topology and Shape Optimization in Structural Design, Journal of Mechanics of Structures and Machines, ( in Review ) 1989 44

UNIVERSITY OF MICHIGAN 3 9015 02527 7883