.- I v . 7 :19. . .. 913.2 1 van. 3.. . wt. 3 «u. . . 5. :55... a. . .u Jrai.‘ 3!. .1; ¢ Ran. - \- ...uunm. ,. .I. :3? .: . . !. ‘31) um. :53... *mrznranuum 2 ‘. a. :1: “hang VP!!! (4.2.1. 5..., :.-._u.n....m. 1 .. z 9.. 9.3%.”. 5.... . My“... .4 5.2 .. kHz. 1 0‘ . qt. 3‘.“ 5.!) £1! ‘9 yr..lI I 2x?! $93.3. 712's it. 915.) , 2.4.. 9113.31 I~ .3 5.1.2. . tLHi: 5». 1:23: . .15 .. h§§s k 11499.13 \ T) «4‘ LIBRARY Michigan State University This is to certify that the thesis entitled TAILORING MATERIALS WITH PRESCRIBED CONSTITUTIVE PARAMETERS USING POLYGONAL CELLS presented by Weian Ou has been accepted towards fulfillment of the requirements for m degree in Mmhnn'nlfinsimm'ns ,lllzw Major professor Date ll bazaar 0-7 639 MS U is an Affirmative Action/Equal Opportunity Institution PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 6'01 cJCIRC/DateDue.p65-p.15 TAILORING MATERIALS WITH PRESCRIBED CONSTITUTIVE PARAMETERS USING POLYGONAL CELLS Weian On A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Mechanical Engineering 2001 ABSTRACT TAILORING MATERIALS WITH PRESCRIBED CONSTITUTIVE PARAMETERS USING POLYGONAL CELLS By Weian Ou The topic of this thesis is to use the topology optimization method to solve an inverse homogenization problem in which the micro-geometry of a composite (characterized by a polygonal cell) made of two constituents is designed to match a set of prescribed elastic properties. In this thesis, we also investigate whether using polygonal cells in the homogenization procedure can result in materials that reach extremal properties. The effect of different starting points is investigated. The results show that different starting points will lead to different answers with the same target tensor, which points at the non-unique of the solution. The question whether different shape parallelograms can lead to different answers matching the same target tensor is examined. The results show that one can achieve different microstructure that match the desired target tensor by changing the geometry of the base cell. We study the effect of changing of the amount of available material. The results show clearly that different volumes of strong material can generate different answers that match the same elastic tensor. The results obtained are close to the lower right comer of the Hashin-Shtrikman bounds. ACKNOWLEDGMENT I would like to thank my advisor Prof. Alejandro Diaz for guidance, along with my committee, Profs. Andre Benard, and Farhang Pourboghrat. Sudarsanam Chellappa provided me with the assistance in testing the program. I also would like to thank Nader Abedrabbo for his help on my programming work in Matlab. iii TABLE OF CONTENTS LIST OF FIGURES .................................................................................. v CHAPTER 1 INTRODUCTION .................................................................................... 1 Inverse Homogenization .......................................................................... l Topology Optimization ........................................................................... 2 Extreme Properties ................................................................................. 4 Periodic homogenization ......................................................................... 8 CHAPTER 2 HOMOGENIZATION .............................................................................. 10 CHAPTER 3 TILING OF 2-DIMENTIONAL DOMAINS USING POLYGON ........................... 19 CHAPTER 4 SOLVING THE INVERSE PROBLEM ......................................................... 26 CHAPTER 5 EXAMPLES ......................................................................................... 33 Examplel-effect of different starting point .................................................. 34 ExampleZ-effect of different base cell geometry ............................................ 38 Example3- effect of different strong material volume ...................................... 47 Example 4- material with extreme property ................................................. 56 CHAPTER 6 SUMMARY AND CONCLUSIONS ............................................................ 66 BIBLIOGRAPHY ................................................................................. 68 iv LIST OF FIGURES Figure 1.1. Plot of the theoretical bounds on effective bulk modulus and shear modulus for isotropic two-phase linear elastic well-ordered materials (Sigmund (1999)) ............ 5 Figure 2.1. Composite structure ................................................................... 10 Figure 2.2. Model of the base cell ................................................................. 11 Figure 2.3. Model of the three-test strain used in homogenization ........................... 16 Figure 3.1. Examples of various prototiles: (a) polyiamond (b)polyomino (c) polyhexe ( Benard and Diaz (2001)) ..................... 20 Figure 3.2. Fundamental domains associates with the tiling with polyiamond, polyhexe, and L-shape prototile ( Benard and Diaz (2001)) ............................................. 21 Figure 3.3. Fundamental domains associated with tiling of L-shape proto—tiles ............ 22 Figure 3.4. Periodic node numbering pattern of polyhexe proto-tile ......................... 23 Figure 3.5. Definition of the parallelogram cell ................................................ 24 Figure 4.1 Flowchart of the procedure to solve the inverse homogenization problem (Sigmund 1995) .................................................................................... 32 Figure 5.1. Examplela base cell .................................................................. 35 Figure 5.2. Example la solution .................................................................. 35 Figure 5.3. Example lb base cell .................................................................. 36 Figure 5.4. Example 1b solution .................................................................. 36 Figure 5.5. Example 1c base cell .................................................................. 37 Figure 5.6. Example 1c solution .................................................................. 37 Figure 5.7. Initial design of the parallelogram cell ............................................. 38 Figure 5.8. Example 2a Lleyzl cell and solution .............................................. 39 Figure 5.9. Example 2b Lx/Lyzl/Z cell and solution ........................................... 40 Figure 5.10. Example 2c Lx/Ly=2/3 cell and solution .......................................... 41 Figure 5.11. Example 2d Lx/Ly=2/5 cell and solution .......................................... 42 Figure 5.12. Example 26 (1:45 cell and solution ................................................ 43 Figure 5.13. Example 2f (1:60 cell and solution ................................................ 44 Figure 5.14. Example 2g (1:75 cell and solution ................................................ 45 Figure 5.15. Example 2h a=90 cell and solution ................................................ 46 Figure 5.16. Example 3a Lx/Lyzl a = 75° volume 0.2 cell and solution ................... 48 Figure 5.17. Example 3b Lx/Ly=l a = 75° volume 0.25 cell and solution .................. 49 Figure 5.18. Example 3c Lx/Lyzl a = 75° volume 0.3 cell and solution ................... 50 vi Figure 5.19. Example 3d Lx/Lyzl a = 75° volume 0.35 cell and solution .................. 51 Figure 5.20. Example 3e Lx/Lyzl a = 75° volume 0.4 cell and solution ................... 52 Figure 5.21. Example 3f Lx/Lyzl a = 75° volume 0.5 cell and solution ................... 53 Figure 5.22. Example 3g Lx/L,=l a = 90° volume 0.2 cell and solution ................... 54 Figure 5.23. Example 3h Lx/Lyzl a = 90° volume 0.3 cell and solution ................... 55 Figure 5.24. Convergence of the objective function with different strong material volume ................................................................................................ 56 Figure 5.25. Example 4a Lx/Ly=1 a = 900 cell and solution (k' = 0.3461 and u‘ = 0.0516.) ....................................................................................... 60 Figure 5.26. Example 4b Lx/Lyzl a = 90° cell and solution (k' = 0.3593 and u. = 0.0357 .) ....................................................................................... 61 Figure 5.27. Example 4c Lx/Ly=(3)"(3/2)/4 a = 79° cell and solution (k‘ = 0.3593 and u' = 0.0307) ....................................................................................... 62 Figure 5.28. Example 4c L,/L,=(3)"(3/2)/4 a = 79° cell and solution ( k. = 0.3593 and it. = 0.0287) ....................................................................................... 63 Figure 5.29. Example 4d Lx/Ly=(3)"(3/2)/4 a = 79° cell and solution (k' = 0.3593 and u' = 0.0282) ........................................................................................ 64 vii Figure 5.30. Bounds of effective bulk and shear moduli for two-phase isotropic composites. The crosses denote effective properties of the topology-optimized microstructures in Example 4. .................................................................... 65 viii Chapter 1 INTRODUCTION Composite materials are becoming increasingly important in industry. The weight and strength ratio of these materials makes them suitable for a large variety of applications. The application of the composite material determines what kind of property it should have. The purpose of this thesis is to use the topology optimization method to solve an inverse homogenization problem in which the micro-geometry of a composite made of two constituents is designed to match a set of prescribed elastic properties. In this thesis, we also investigate whether using polygonal cells in the homogenization procedure can result in materials that reach extremal properties. 1.1 Inverse Homogenization The goal of an inverse homogenization problem is to find a (periodic) microstructure in which the infinitly small-scale mixture of two constituents forms a composite material whose effective elastic tensor matches some arbitrarily prescribed values. The inverse homogenization problem here deals with linear elastic materials and small deformation kinematics in both macroscopic and the microscopic scales. It can be shown that after using the homogenization procedure, the constitutive parameters can be expressed as explicit functions of a finite number of parameters. The inverse homogenization problem can thus be formed as an optimization problem to possibly minimize the deviation from the target tensor with a prescribed amount of material. We can then use the topology optimization method to solve it. Sigmund (1994) solved the inverse homogenization by using square or rectangular base cells only. Although this is not a fundamental constraint, this restriction can be released by using polygonal cells without significant additional numerical complexity. It will increase the design space remarkably by using polygonal cells to solve the inverse homogenization problem. In Sigmund (1999), the inverse homogenization problem is solved to find a new class of extremal composites. Also in this work only the square cells are used and the bound of lower shear modulus is reached. The idea of using the polygonal cell is inspired by Benard and Diaz (2001). In that paper, features related to the discretization of problems characterized by simple periodic tiling using cells of various shapes are discussed. The paper shows that various cell geometries capable of tiling the plane can be replaced by an equivalent problem where the plane is tiled by polygonal base cells. The base cell will be discretized with finite elements. The result of the inverse homogenization will be presented as an optimized geometric distribution of the two constituents in the base cell. 1.2 Topology Optimization Topology optimization of plane structures involves the determination of features such as the number and location of holes and the connectivity of the domain. In the inverse homogenization problem, the goal is to allocate the prescribed volume fraction in the design domain to make its elastic tensor meet the arbitrarily prescribed elastic tensor, or target tensor. Composites can be designed to have many physical properties, such as Young’s modulus, Poison’s ratio, shear modulus, bulk modulus, conductivity, etc. Composite materials with extreme properties (e.g., the “stiffest” composite) are very important in material science and structural optimization. They can be used to establish bounds for structural performance and are widely used in the field of topology optimization. Bendsoe, Guedes, Haber, Pedersen, and Taylor (1994) answered in analytical form the question “what is the optimal configuration of material throughout a structure that is consistent with the objective for globally optimal design”. In that paper, they show that the measure of compliance predicted in their formulation bounds from below the compliance of similar structure, which presented in the paper, relative to all possible material configurations. Bendose, and Kikuchi (1988) use composite mixtures to produce a methodology for optimal shape design that avoids drawbacks that affect “traditional” shape optimizations. One such drawback is that the optimal shape design of structural elements based on boundary variations results in final designs that are topologically equivalent to the initial choice of design. Thus, the region of the design space explored by such method is limited. Another drawback is that general stable computational schemes for standard shape optimization often require some kind of remeshing of the finite element approximation of the analysis problem. 1.3 Extreme Properties The elastic tensor of an isotropic material, which is composed by mixing two isotropic constituents can be described by its effective bulk modulus k. and its shear modulus 11‘. Considering two isotropic constituents with bulk and shear moduli denoted by k1, k2 , [11 and #2 , respectively, if the constituent with larger bulk modulus also has the larger shear modulus ((k2 — kl Xuz — ”1).? 0) , the two constituents are called “well—ordered” . Otherwise they are called non-well-ordered ((k 2 — kl Xp 2 — M1 )< 0 ). There exist bounds of the effective bulk modulus and effective shear modulus that this kind of composite material can reach. The first bound was given by Hill(l952). Later, these bounds were improved by Hashin and Shtrikman (1963) for three-dimensional composites and well- ordered constituents. Hashin and Shtrikman (1963) applied variation principles, involving the elastic polarization tensor, to derive the upper and lower bounds for the effective elastic moduli of the materials described as mechanical mixtures of a number of different isotropic and homogeneous elastic phases. The bounds for the non-well-ordered case were then improved by Walpole (1966). Walpole presented general bounds for both well-ordered and non-well-ordered cases, as Hashin and Shtrikman (1963) only considered the well-order case. Subsequently, Milton and Phan-Thien (1982) and Berryman and Milton (1988) set up new bounds on the effective moduli of a two-constituent-composite material. Milton and Phan-Thien (1982) showed in detail how Hashin and Shtrikman’s (1963) bounds can be extended and how Walpole’s (1966) bounds can be improved using two inequalities on the two geometrical parameters that appear in the third-order bounds on the effective shear modulus. In Berryman and Milton (1988), examples of bounds on one effective material property from measurements of another are derived. These examples are somewhat more restrictive than the Hashin-Shtrikman bounds. Cherkaev and Gibiansky (1993) improved these bounds by the translation method. The set of the bulk modulus and the shear modulus pairs turns out to be bounded in the plane of the values of these moduli by straight lines and also by two fractional linear curves. Sketches of the Hashin-Shtrikman and Cherkaev-Gibiansky bounds for well-ordered case are shown in Figure 1.1. PuHS ii __ llashinShtrikhman bounds — Chm‘kaev“ Gibiansky bounds 111er .4 MW ll] HS ‘ K K 1H8 ITuHS Figure 1.1. Plot of the theoretical bounds on effective bulk modulus and shear modulus for isotropic two-phase linear elastic well-ordered materials (Sigmund (1999)). If there exist composites that attain the bounds, the effective properties of such composites are said to be optimal. An interesting question is whether one can identify the micro-geometry of composite materials that can attain the bounds. The first microstructure attaining the maximum bulk modulus bound was the spheres model composite presented by Hashin (1962). Francfort and Murat (1986) focused their attention on the mixture of two isotropic materials in prescribed volume fractions and strived to characterize all possible macroscopically isotropic composites. They produced a multi-layered composite with a finite number of element of layering directions. The results they gave proved that Hashin-Shtrikam bounds on the bulk and shear moduli can be achieved and suggested that the bounds can even be simultaneously achieved by multiple layering. Two-dimensional square symmetric rank-2 laminations attain the upper and lower bounds of the bulk modulus. Isotropic rank-3 laminations can attain the maximum bulk and shear modulus simultaneously (upper right comer of the Figure 1.1) and attain the lower-left comer of Figurel.l by inversing the two constituents. Lutie and Cherkaev (1984) developed a method for the calculation of the extreme conductivity. The method is based on the principle of the consecutive assembling of binary mixtures through the addition of infinitely small amounts of one of the initial compounds to the already-assembled isotropic composite. The process is assumed to produce an optimal isotropic binary mixture at each step, which is performed by the Hashin-Shtrikman procedure. They suggested that the same method can be used in the optimal design of elastic constructions. In other words, there exist finite rank laminates that can attain the bulk-shear muduli bounds. Milton (1992) produced examples of composites that can be rigorously proved to have low-bulk and high-shear moduli. He found a family of two- dimensional, two-phase, composite materials with hexagonal symmetry with Poison’s ratios close to -—1. He also showed that elastically isotropic two-dimensional composites with Poisson’s ratio approaching —1 can be generated simply by layering the component materials together in different directions on widely separated length scales. In Milton and Cherkaev (1995), for a two-phase composite comprised of a sufficiently compliant, isotropic phase and a sufficiently rigid, isotropic phase, the effective elasticity tensor can realize any given definite tensor satisfying the usual symmetries of elastic tensors. This means that the optimal bulk and shear moduli bounds can be attained by this kind of composite material. Composite with extremal properties such as rank-2 laminations are multi-scale mixtures. Instead of using the multi-length scale microstructures, Grabovsky and Kohn (1995) used the “Vigdergauz microstructure” to obtain the two-dimensional square symmetric composites with optimal shape of single inclusions. The composite is spatially periodic, consisting of properly shaped elastic inclusions embedded in an elastic matrix. Vigdergauz (1999) provided an alternative discussion of “Vigdergauz microstructure” and its properties. He also presented an analysis of various limits of the inclusion and extended the results of Grabovsky and Kohn to isotropic material. Si grnund (1993, 1994) used a numerical topology optimization approach to solve the material designing problem with extremal elastic properties. In his paper, examples of two-dimensional, two-phase microstructures with external bulk modulus are obtained by the inverse homogenization procedure. Ole Sigmund, (1999) tried to produce the isotropic material microstructure with maximum bulk modulus and minimum shear modulus. The results get very close to the lower bound but have not attained it. In this thesis we will try to find out whether a two-dimensional microstructure with a periodicity characterized by arbitrary polygonal cells can attain the lower shear bound. 1.4 Periodic homogenization In the 1970’s, much research focused on defining the equivalent mechanical properties of composite materials and on determining their dependence on the difference components. Some of the methods used in that research were based on engineering. They often show a good agreement with each other’s experimental data or with empirical methods. In Hashin’s (1970) paper, a survey in this approach is provided. At the same time, the counterpart of such engineering methods in mathematics appears under the name of homogenization theory. This theory has been the object of large amounts of research in the area of applied mathematics. A few examples are Lions (1981) and Oleinik (1984). A homogenization method as proposed by Guedes et al.(l990) can be used to calculate the effective constitutive parameters of complex materials. In homogenization theory, the composite material is assumed to be locally formed by the spatial repetition of microstructures infinitely small, when compared with the overall ‘macroscopic’ dimensions of the structure of interest. The microstructure is called ‘microscopic’ cells, or base cells. In other words, it is assumed that the material properties are periodic functions of the microscopic variable, where the period is infinitely small compared with the macroscopic variable. The assumption enables the computation of equivalent material properties by computation of the material property of the base cell. The homogenization procedure has been discussed several times in literature. One way to solve the homogenization problem is analytical calculation. For example, Aboudi and Jacob (1991) present the way to provide the overall behavior of composite materials by the micro-mechanical analysis from the known properties of the individual constituents (e. g. fibers and matrix).When more complicated microstructure need to be considered, analytical homogenization becomes too complex to be calculated. In this condition, numerical method based on the finite-element must be used. Guedes and Kikuchi (1990) presented an effective numerical homogenization procedure with adaptive finite elements methods. In this thesis, a modified version of this algorithm associated with element mutual energy will be used. This thesis is divided into 6 chapters. Chapter 2 will describe the homogenization method. Chapter 3 will discuss relevant features and concepts related to tiling two— dimensional domains using different polygonal base tiles. Chapter 4 will describe the inverse homogenization method. Chapter 5 will present examples of inverse homogenization problem solved and discuss whether it is possible to attain the lower shear modulus bound. Chapter 6 will present conclusions and summarize the work that has been done. Chapter 2 HOMOGENIZATION In this chapter, basic concepts in homogenization theory relevant in this work will be briefly reviewed. Consider a composite material formed by the spatial repetition of a base cell made of two constituents, as shown in Figure 2.1.The figure presents a two- dimensional case. We assume that the mixture is represented by a base cell that is very Constituentl (strong) Constituent2 (weak) =. Base cell Figure 2.1. Composite structure small (of order 8 ) compared with the dimensions of the body. As a result of the rapid variation of material properties, when the body is subject to some load, the resulting deformation and stresses vary rapidly from point to point. These quantities have two explicit spatial dependences: on the macroscopic scale variable x and on the small scale y = x/ E . Because of the periodic nature of the microstructure, the dependence of these functions on the microscopic variable y = x/ 8 is also periodic. However, discretization of the finite scale would be impossible, e. g. using finite element methods. It is therefore necessary to develop a method that averages the microscopic structure whenever the mechanical behavior of the macroscopic body is in question. The homogenization method is such a method. Let Q be an open subset of R2 with a smooth convex boundary F . Let Y be an open rectangle in R2 (Figure 2.2), defined as Y =10. y?rx10. ygr (2.01) Later in this thesis this will be generalized so that Y can be any polygonal shape that tiles the plane. Let :9 be an open sub set of Y with boundary 319 = S (2.02) and let w = Y \19 (2-03) Figure 2.2. Model of the base cell In our problemw is the area occupied by the strong material. The closure of 19 , denoted by 19 , is the area occupied by weak material, and Y represents the base cell of the 11 composite microstructure. The material properties vary inside Y and are defined using the following indicator function. 1 if y e w, 90’) = . (2-04) 0 If ye m, We can then define $28 = {ms 52 | 0(x/e) = 1} (2.05) as the subset of the domain 82 occupied by the strong material and Number of cells SE = USO, (2.06) a=1 as the “boundaries” between strong and weak materials. The following assumptions are made: 1. $28 is a connected domain. 2. In all cells the weak material part 19 has sufficient smooth boundary S. 3. None of the boundary S intersects with the boundary T‘ of 52. v8 =ire(Hl(S2£))2 |v|rd=0} (2.07) where v I111 is the value of v on the boundary I‘d . The problem of the deformation of a structure body S25 subjected to body force f and tractions t on the boundary 1", with the tractions p inside the boundary S E , and prescribed displacement on I‘ d , can be stated as 12 Find u‘eV‘, such that a“: 8‘)" 6 E 6 ax, 871.10- [9, fl. vidQ+IrltividF+Jst pivids We v (2.08) E EW 91' Here strain-stress relation is Psi=iEsiifgi (2.09) where the elastic tensor matrix is given in the standard engineering form, i.e., E 8 is written as E1111 E1122 E1112 [E]: E1122 E2222 E2212 (2-10) E1112 E2212 51212 E The conditions for existence of a unique solution u‘ to problem (2.08) is given in Necas and Hlavacek (1981), and require that the functions f , t, and p be sufficiently smooth, and that the boundaries 1" d , I“, , and S E be regular. The solution u‘ should depend on both x and y , i.e., u‘E =u(x,y), yzx/E (2.11) In the neighborhood of a fixed macroscopic point x it is assumed that there is a very large number of micro-scale cells, which are obtained as copies of a base cell. Dependence on y can be considered periodic or Y—periodic, at a fixed point x in the microscopic level. Furthermore, it is assumed that the form of the base cell varies in a smooth way as the macroscopic variable x changes. This means that at different points x the composite structure may vary, but a periodic pattern can always be found in a microscopic neighborhood of point x. 13 Since the solution as depends on both macroscopic and microscopic variables, it is reasonable to assume that u‘S can be expressed as an asymptotic expansion with respect to the parameters (the “scale”, or the ratio of microscopic/macroscopic dimensions), i.e., us =u0(x,y)+£ul(x,y)+82u2(x,y)+..., y=xl8 (2.12) where u1'(x, y) is defined on (x, y)e wa (2.13) u1'(x,y) is Y-periodic t In order to find the homogenized elastic constants such that the macroscopic equilibrium can be described by the same equation as (2. 08), let 1““) E V be the solution of L2,} av, av j Eiqu— ayq 3y ——-=dY [Y 5W 5;de VveV (2.14) v = {we (191(9))2 :v is Y — periodic} and let 1;! k e V be the solution of 3W a IYl‘:"1""_ay,—_keel—(”USP")(My Vvev (2.15) If only the first order terms of the us asymptotic expansion are considered, it can be shown (Sigmund 1994) that the global elastic constants of the material are given by akl x {Y EW— p dY (2.16) H l E Eiqu dyq ijk1:Y l4 Thrs equation represents the macroscopic equrlrbrrum, whrle Eijkl represents the homogenized elastic constants, i.e., 0 H auk 3"i Q U“ 3x, ax]- 1 1 all/k av,- “(V—IL f,dY)v,-dr2 + In i,v,-dr + [9(l—Y—ljy 15W F[211037de (2.17) VVE VQ As shown above, the homogenized elastic constants can be computed within the basic cell by solving the problem 2.14, and do not depend on the macroscopic deformation uo. If a composite material has a uniform cell structure in the whole domain £2 , the microscopic problem 2.14 needs to be solved only once. Following Sigmund (1994), equation (2.16) can be rewritten as H _ l _ *( kl) Eijkl ‘ Yiy [Eiqu Eiqugpq LY (2'18) where 62(qu = 81“} / Byq is still the solution to (2.14). Let now 825]“) represent one of three homogeneous tests strain-fields applied on the base cell, namely, two tensile strains and one shear strain (Figure 2.3), i.e., 1 0 0 {£01}: 0 {:02}: 1 {503}: 1 (2.19) 0 0 0 15 —-—_— ’— _,-. _.——. ........ 80‘ so' 80‘ Figure 2.3. Model of the three-test strain used in homogenization The solutions to problem (2.14) correspond to the deformation of the base cell, when subjected to the three independent cases of pre-strain {£0} and periodic boundary conditions, and can be obtained numerically using a finite-element method. Homogenized elastic constants can be obtained by substituting the result of the finite-element analysis *(kl) - 8P9 mto (2.18). In order to use a topology optimization method to find micro-geometries with prescribed effective constants, it is important to express the elastic constants E 5k! in (2. 18) in terms of element strain energies. This will make it possible to use the existing algorithms for the topology optimization. A homogenization method expressed in terms of element strain energies is described in the following. Rewriting equation (2.18) we obtain 0 H 0 _1 k0 * T 0 .. ]d giants/(1‘71), 517-517) Eijkl(8kj—8k1) Y (2.20) where a]?! and 83 are the test strains field and 8,;- and 8;] are the induced strains field coming from the in-homogeneities of the rrricrostructure . If the base cell is discretized by N finite elements, the element mutual energy Q,e in each element e, associated with test strain field {80 }, can then be written as 1 0 * T 0 * Qle =—e Ikeea) -Ee(1)) Eijkl (86(1) —8e(l))]1ye 6:1,...,N (2.21) Y Ye e—-—1—k°—*TE-°—*]ie-1N 222 92 - e1 8e(2) 8cm) W154» gem) y 6- (- ) Y Ye 1 O * T O * Q38 : —; J‘k88(3) —£e(3)) Eijkl (88(3) —Ee(3))Lye C:1,...,N (2.23) Y Ye 1 0 * T 0 * Q4e = '7 I18e(11‘8e(1)) 50111542) ‘84:»)er e=l,...,N (224) Y Ye 1 0 * T 0 * Qse “7 1 182(2) ‘5e(2)) Eijk1(5e(3) “512(3)er e=l,...,N (2-25) Y Ye 1 0 * T 0 * Q66 :3}? 111841) “540) Err/c1154» ’84:»)er e=l,...,N (226) Ye ' where test strain field {92 }= {£0 }; {15: } is the induced train field in element e . l7 Since we assume that the base cell is discretized by N finite elements, we can sum the element strain energies to obtain the expression of each entry in the homogenized elastic tensor E H . N E,” = 2Q," e=1,...,N ml,2,3,4,5,6 (2.27) e=l where we have used the following abbreviation llll—>l 2222—>2 1212—>3 1122—P4 1112—>5 2212—>6 Now we only need to calculate each element strain energy Q: from (2.21) to get the homogenized elastic constants of the composite material. 18 Chapter3 TILING OF 2-DIMENSIONAL DOMAINS USING POLYGON In the previous derivation (Chapter 2), it is usually enough to use square cells only (i.e., Y is a square or a rectangle). However, base cells with more general shapes could lead to other interesting results. In this thesis, we use the topology optimization method to solve an inverse, homogenization problem with polygonal cells. For this purpose, Y needs to be defined so that it can tile the plane periodically. If Y is defined as a simple square or rectangle, this requirement can be easily met. When we allow the base cell to take a polygonal shape, this restriction is more difficult to impose. We should review basic tiling concepts and introduce some basic results from Benard and Diaz (2001). We consider various shapes of cells capable of tiling the plane. Let T be a set of tiles in plane, which cover the plane without gaps or overlaps. T is a collection of tiles Ti. Here we only discuss periodic monohedral tilings which involve translations of a single polygonal proto-tile. All tiles have the same shape and size in monoheral tilings. In a periodic tiling, each set Tiis a tile. Every tile Tiin T has the form Ti = P + Mdl + "1612 where P is the proto-tile, representative tile in monohedral tiling, m and n are integers and d1 and (12 are two non-parallel vectors called the tiling vector. The dependency of Ton P, (11 and d2 will be expressed as T(P,d1,d2). 19 (a) (b) Figure 3.1. Examples of various prototiles: (a) polyiamond (b)polyomino (c) polyhexe ( Benard and Diaz (2001)) Figure 3.1 shows several interesting examples of different proto-tile. The proto-tiles used to tile the plane are polyiamond, L-shape, and polyhexe. These proto-tiles are made by simply connecting the triangles, squares, and hexagons respectively. We now introduce three important definitions leading to the concept of Y-periodicity when Y is a polygon. Definition: Lattice. 20 A lattice is a collection of translates of a single point p e P associated with a periodic tiling T along the tiling directions, i.e., the set of points L(d1,d2) = {q : q = p +md1+ nd2,m and n are Integers} Such lattice forms periodic parallelograms called fundamental domains, or fundamental parallelograms. Figure 3.2. Fundamental domains associates with the tiling with polyiamond, polyhexe, and L-shape prototile ( Benard and Diaz (2001)) Definition: Fundamental Domain. 21 A fundamental domain is any parallelogram F with its comers on the lattice L(d1,d2) , which tiles the plane with tiling vectors v1 and v2 and generate a lattice L(v1,v2)equal to L. There can be several fundamental domains associated with the same lattice. Figure 3.2 shows several examples of fundamental domain. In the figure polyiamond, polyhexe, and L-shape proto-tiles can be associated with different parallelogram individually. P“) d2 v 13(2) v2 Figure 3.3. Fundamental domains associated with tiling of L-shape proto-tiles Figure 3.3 shows a periodic tiling of L-shape proto-tiles. In the figure, two fundamental domains F (1) and F (2) associated with the same tiling, in which L(p,d1,d2)and L( p, v1,v2) are identical( the four points of the parallelograms F (I) and F (2) are the same points of the L-shape proto—tile). Definition: T -periodic function. 22 A function : R2 —> R is T(P) -periodic, if ( p) = ( p + mdl + mdz) for all pe P and integers m and n. A function (I) that is T(P) -periodic uses the base cell P in its definition. It can be very complicated to implement the numeric associated with the homogenization problem on such a base cell P. Figure 3.4 shows an example of a node numbering of a polyhexe base cell. In order to meet the periodic boundary conditions, the node numbers have to be Figure 3.4. Periodic node numbering pattern of polyhexe proto-tile assigned in the pattern shown in figure 3.4. If the shape of the base cell changes, the numbering of nodes will also need to be changed. This would be a very cumbersome procedure. 23 We can use a parallelogram as fundamental domain to avoid these difficulties. From the examples shown in Figure 3.2 and Figure 3.3, we can see that fundamental domains preserve the periodicity of T(P) -periodic function. Because proto-tile P and fundamental domain F are associated with the same lattice, a function that is T(P) - periodic is also T(F) -periodic. The material design problem then can use F instead of P as its base cell to solve the homogenization problem. So no mater how complex the polygonal base cell is, we can always find a fundamental domain to describe the periodicity of the problem. In others words, we can always use a parallelogram cell to solve the T(P) -periodic function in which P is a polygonal cell. This result will expand the design space of material design problems. Figure 3.5 shows the geometry of a fundamental domain. A more complete parameterization of design space of material tensors can be achieved by including in the analysis the sides Lx/Ly and the internal angle a . One may also find interesting solutions that cannot be reached when the base cell is restricted to be rectangular. Figure 3.5. Definition of the parallelogram cell 24 As in Figure 3.2, no mater which kind of polygonal proto-tile is used to tile the plane, we can always find a fundamental domain associated with it. After we transfer all these tiling problems into a basic parallelogram problem, solving the parallelogram cell-tiling problem can use a general routine. 25 Chapter 4 SOLVING THE INVERSE PROVLEM In this chapter, the numerical approach for material design or inverse homogenization method is briefly reviewed. The method is then applied to the design of two-phase E composites with mixtures characterized by a polygonal cell. As stated in Chapter 3, the material microstructure is assumed to be periodic and fully L described by its smallest repetitive base cell. In order to use the finite element method, the base cell is discretized by N finite elements. The effective properties are found by using the numerical homogenization method, which is stated in Chapter 2, using the finite element analysis to find the local fields and applying periodic boundary conditions. In approaching the problem of getting the effective average elastic properties E H of the base cell, the popular SIMP method, which is also used by Sigmund ( 1995), is used here. This method uses a variable thickness sheet formulation and penalized intermediate design values. The method assumes that the element elastic tensor E e is isotropic and expressed as E3 = prO p, 6 (0,1] (4.1) 26 Here p = {pl pN} are the design variables, and p is a real number, p >1, used to penalize the intermediate design values (i.e., values of p between 0 and 1). E 0 is an isotropic reference tensor, which is expressed in the standard (matrix) engineering form, E 1 v 0 [15°]: 2 v 1 0 (4.2) 1—1) 0 O (I ;‘U) where E is the modulus of elasticity and v is Poison’s ratio. In terms of shear modulus p and bulk modulus k , the Young’s modulus can be defined as = 9"" (4.3) 3k + u and Poison’s ratio can be defined as v:1 3"'2‘“) (4.4) 2 3k + u The isotropic reference tensor E 0 can also be expressed in terms of shear modulus u and bulk modulus k l 3k—2,u 9k” 6k+2u [E0]: 3k+u 3k—2,u l 0 (4.5) 1_(3k—2#)2 6k+2u 6k+2u 0 0 (3k+4u) L 12k+4u_ The inverse homogenization method is applied here to solve the following problem: 27 For a given amount of material, find the microstructure arrangement, which minimizes the deviation from a target tensor E i . In other words, in this material design problem, the goal is to find a set of design variables p = { p1 pN} that make the material effective tensor match the target tensor as close as possible. In addition to p , we also introduce an additional design variable, t, to scale the material effective tensor as a scaling parameter. The optimization problem is written as E Find {,01 , p2..., p N} and t that will : (4.6a) Min: ER (4.6b) N Subject to 2 Ac pe = VOA e=l (4.6c) O(pmin SpeSl where p, represents the design variables, in this case, the proportion of strong constituent in each element; N is the total number of elements that are used to discretize the base cell; ER is the deviation from the target tensor, which will be defined below; A, is the area of each element;v0 is a prescribed volume fraction, a number between 0 and l; A is the area of base cell; and pmin is a prescribed lower bound for p: , a small number (e.g., pm = 0.001 ). The objective function used in (4.6) is the weighted mean square deviation between E H , the homogenized elastic tensor, and the target tensor E i . It is expressed as 28 an As Opt the “'he Usin 6 ER =%2w1(tEIH —E;‘)2 (4.7) [=1 where compact notation E {I is used to denoted the six independent entries of the elastic tensor using abbreviations llll—bl l—Dl 2222—}22—52 1212—D33—D3 1122—>12—>4 1112—>l 3—>5 2212—b23—D6 and w I is a prescribed, non-negative weight factor (1:1,. . .,6); E * is the target tensor. As the function ER is an unconstrained, convex function oft for fixed p , we can get the optimal value of t ,t* , by solving the equation 95:15 = 0 to minimize ER. Using 6 955: E,”w,(rE,” -E}") (4.8) a: [:1 the minimizer is 1* = C—1 (4.9) C2 where 6 * 6 c1: ZE,w,E,” and c2 = zEfiwm,” (4.10) [:1 [=1 Using the t‘ in (4.5) and (4.6), the optimization problem is reduced to: Find {p1,p2...,pN} that will : (4.11a) 1 6 . H . 2 Min: ERz—Z-ZWIU E, —E,) (4.11b) [=1 29 We C an Defi and . L / /_ f N Subject to 2 Ac pe = VOA =1 (4.11c) 0 ' ~ 1 I.“ Figure 5.7. Initial design of the parallelogram cell Figure5.7). The results for fundamental domain with different length ratios and different internal angles are shown in Figure 5.8 to Figure 5.15. First four figures show the results of different length ratio, the length ratios are Lx/Ly =1, 1/2, 2/3, 2/5 respectively. The second four figure show the results of different angles, 45° , 60° , 75° and 90° 38 individually. The results show that one can achieve different microstructures that match the desired target tensor by change the geometry of the base cell. Figure 5.8. Example 2a Lx/Ly=l cell and solution 39 Figure 5.9. Example 2b Lx/Ly=l/2 cell and solution 40 o.‘ .1 E. r nan z... nan w {a t , -4 . _ 4%? Z a a a} Figure 5.10. Example 2c Lx/Ly=2/3 cell and solution 41 Figure 5.11. Example 2d Lx/Ly=2/5 cell and solution 42 Figure 5.12. Example 26 01:45 cell and solution 43 Figure 5.13. Example 2f 01:60 cell and solution ooooooooo 6 O c aaaaaaaaaaaa vvvvvvvvvvvv ooooooooo 000000000000 O C O O ' 74337674333333? 33333337633115 maeaaeaaaeaa Figure 5.14. Example 2g a=75 cell and solution 45 5.3 Example3- effect of different strong material volume In this example we will investigate that different volume v0 can lead to different answer with same target tensor. All the parameters are fixed, except the percentage of strong material. First, we set the shape of parallelogram as Lx/Ly=l and a = 75° . The volumes are 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, and 0.5. Since the figure for volume 0.45 has been shown in Figure 5.14, only the others will be shown here. Second, we set the shape of parallelogram as Lx/Ly=l and a = 900 . The volumes are 0.2, 0.3 and 0.45. For the same reason, the figure for volume 0.45 has been shown in Figure 5.15. The results show it clearly that different volume of strong material can generate different answers that match the same prescribed elastic tensor. When comparing these results; we can observe that small changes in volume induce small changes in the solution. Figure 5.24 shows the convergence curve of example 33, 3c and 3e. From the figure we can see that the volume has no effect on the pattern of convergence. Also we can see that the objective function value converges rapidly after an initial “flat” period. 47 33.3333 33,133.33. 3.33.33. 13.33 33.333321333333131. ”3.33;, 'o}§}‘o1'o1§;3* 9 ”£1933. 13.111. 311.31.111.13. 1:53:11: 9}: .“ 1.13.1'.1'.1'.1‘ 00000 I151“ Figure 5.17. Example 3b Lx/Ly=l a = 75° volume 0.25 cell and snlntinn 49 .. ... .33... .q volume 0.3 = 75° Figure 5.18. Example 3c Lx/Ly=l a cell and solution 50 Figure 5.19. Example 3d Lx/Ly=l a = 75° volume 0.35 cell and solution 51 Figure 5.20. Example 3e Lx/Lyzl a = 75° volume 0.4 cell and solution 52 3333333333. 3333333333: 3333333333~ 3333333333 :333333333 333333333. 3333333333 r3333333333 3333333333“ 0. :5: Figure 5.21. Example 3f Lx/Ly=l a = 75° volume 0.5 cell and solution 53 Figure 5.22. Example 3g Lx/Ly=l a = 90° volume 0.2 cell and solution 54 5.3.333: 55.35: $35.35: W 93.3.33 b, s 35.3.39: 33333 3.33333 u 3. e. . x a a .. a .w s 3 14 I I T T I l l I J -—— Volume 0.2 ---~ Volume 0.3 7 ~ - Volume 0.4 12i— H M 1 _ .4 0 2 9 8 0.8 — d ‘8 C 2 gm _ 25‘ O 0.4 - — 0.2 — — o I l I I L 1 0 10 20 30 40 50 60 70 80 90 100 Number of Iterations Figure 5.24. Convergence of the objective function with different strong material volume 5.4 Example 4- material with extreme property In this example, we will present several solutions which approach the lower right corner of Hashin-Shtrikman bounds for two-phase isotropic composites in well order case. In well order case, the upper Hashin-Shtrikman bound of bulk modulus is defined as 56 kuHS=k2+ 1 1"” p (5.2) + kl—kz k2+flz and the lower Hashin-Shtrikman bound of shear modulus is defined as ”(HS =u1+ p 1 1 (5.3) l +(l-lOXk +2M) #2 -u1 2u1(k‘+u‘) where k1, k 2, [11 and #2 denote the bulk moduli and shear moduli of two isotropic constituents , p is the given volume fraction of strong material. The improved lower Walpole bound of shear modulus is defined as u,“’”=u‘+ ’0 2 I (5.4) 1 +(l-p)(k +2u> #2 -M‘ 2u‘(k2 +u‘) Here our goal is set as upper Hashin-Shtrikman bound of bulk modulus and lower Walpole bound of shear modulus, which is the lower right corner of Figure 1.1. The solutions here are generated by fixing the bulk modulus at the upper bound and decreasing the shear modulus to approach the lower Walpole bound. We set p =0.45 and k2 =1.333,k1=0.01k2, #2 =1 and ul=0.01u2.From equations (5.2) and (5.4), the values of the target bulk modulus and shear modulus are k* = 0.3593 and if“ = 0.0257. 57 We first start with k' = 0.3461 and u. = 0.0516 . The geometry of the base cell is Lx/Lyzl a = 90° . The target tensor generated is 0.1807 0.0775 0 [15‘]: 0.0775 0.1807 0 (5.5) 0 0 0.0516 Figure 5.25 shows the result. Figure 5.26 shows the result for k‘ = 0.3593 and u. = 0.0357 . The geometry of the base cell also is Lx/Ly=l a = 90° . The target tensor generated is 0.1303 0.0589 0 [12‘]: 0.0589 0.1303 0 (5.6) 0 0 0.0357 Next step we set k‘ = 0.3593 and u' = 0.0307 . This shear modulus is middle point between it. = 0.0357 and 11* = 0.0257 (the lower Walpole shear modulus). The target tensor is 0.1133 0.0519 0 [5"]: 0.0519 0.1133 0 (5.7) 0 0 0.0307 The geometry of the base cell also is Ia/Ly=(3)"(3/2)/4 and a = 79° . The result is shown in Figure 5.27. Next step we set k‘ = 0.3593 and pi = 0.0287 . The geometry of the base cell also is Lx/L,=(3)"(3/2)/4 and a = 79° . The target tensor is 58 0.1064 0.0491 0 [5‘]: 0.0491 0.1064 0 (5.8) 0 0 0.0287 The result is shown in Figure 5.28. We then have tried k' = 0.3593 and ,u' = 0.0277. But the result can not converge to the target tensor. So we set the [f = 0.0282 , which is the middle point between it. = 0.0287 and u' = 0.0277 . The geometry of the base cell also is [alLy=(3)"(3/2)/4 and a = 79°. The target tensor is 0.1047 0.0483 0 [15"]: 0.0483 0.1047 0 (5.9) 0 0 0.0282 The result is shown in Figure 5.29. Comparison of these figures indicates that when the material properties approach the goal, the structure of the material is much more close to a hexagon. In Sigmund (1999), he uses a hexagon shape for the cell as a starting point to get a new class of extremal composites. Our routine generates the same trend in optimization. But still under the condition of the parameters we set, the optimization result cannot attain the theoretical value. The crosses in Figure 5.30 denote the effective properties of the topology-optimized microstructure we get in this example. 59 Figure 5.25. Example 4a Lx/Ly=l a = 90° cell and solution (k‘ = 0.3461and 11‘ = 0.0516 .) 6O or??? 90° cell and 0.0357 .) Figure 5.26. Example 4b Lx/Ly=l a 0.3593 and 11* solution (k ‘ 61 Figure 5.27. Example 40 Lx/Ly=(3)"(3/2)/4 a = 79° cell and solution (k. = 0.3593 and [1. = 0.0307) 62 79° cell Figure 5.28. Example 4c Lx/Ly=(3)"(3/2)/4 a and solution (k. =0.3593 and u‘ 0.0287 3 63 Infumfufnfm’ufuffln #WWWWWWWWWW Wtfiiflhmmfllmmm uuuuuuu In, quIquIqufqufu tttvttoWWW ;*;.;*;*;*;“;'4.;*4. ¢WWWWWW WWW mmmmmmmmmm 'WWWWWWWW¢¢ 3333:3333“: Ilffifl4fl1fl4fl4fi4flffl;fl Figure 5.29. Example 4d Lx/Ly=(3)"(3/2)/4 a = 79° cell and solution (k. = 0.3593 and ,u' = 0.0282) l o4 l—luHS {E HashimShtrikhman bounds —Cllerkaev~Gibiansky bounds _ _ 111W? #qu __ .1: K lilHS K' IHS K'ulls E Figure 5.30. Bounds of effective bulk and shear moduli for two-phase isotropic composites. The crosses denote effective properties of the topology-optimized microstructures in Example 4. 65 Chapter 6 SUMMARY AND CONCLUSIONS The inverse homogenization problem with polygonal cells has been successfully solved by topology optimization method by matching a target tensor. The effects of the initial material distribution, the shape of polygon cell and percentage of strong material on the solutions that match the same target tensor are studied. Besides these, material distributions approaching the lower right comer of Hashin-Shtrikman bound for well- order case have been found. A summary of the results is presented below: 0 Different starting points will lead to different answer with the same target tensor in the optimization procedure. Convergence to a particular solution depends on the initial material layout. 0 One can achieve different microstructures that match the desired target tensor by changing the geometry of the base cell (parallelogram). This will expand the design space when compared to using only a rectangular cell. 0 Different percentage of strong material can lead to different answer with same target tensor. Small changes in volume induce small changes in the solution. The volume has no effect on the pattern of convergence of the optimization problem. 0 The results, which are generated for approaching the lower right comer of Hashin-Shtrikman bound for well-ordered case, show that more and more layers 66 in the bars of the hexagonal cell result in lower and lower shear modulus (Sigmund 1999). 67 BIBLIOGRAPHY Aboudi , J. "Mechanics of Composite Material". Amsterdam: Elsevier Science Publishers B.V 1991. Benard, A., and A. R. Diaz. ”On the discretization of problems involving periodic planar tilings”. Commun.Numer.Meth.Engng (2001). Ashby, M.F., and D.R.H. Jones, Engineering Materials 2, Oxford: Pergamon (1986). Bendsoe, M.P., J .M. Guedes , R.B. Haber, P. Pedersen , and J .E. Taylor. “An Analytical Model to Predict Optimal Material Properties in the Contest of Optimal Structure Design”. J .Appl.Mech. 61(3), (1994),930-937 Bendose, M.P., N. Kikuchi. "Generating optimal topologies in optimal design using a homogenization method". Computational Methods in Applied Mechanics and Engineering. 71 (1988), 197-224 Berryman, J .G., and G.W. Milton. “Microgeometry of random composites and porous media “. Journal of Physics D: Applied Physics. 21 (1988), 87-94. Cherkaev, A.V., and LN. Gibiansky. "Coupled estimates for the bulk and shear moduli of a two dimentional isotropic elastic composite". Journal of the Mechanics and Physics of Solid. 41(5), (1993), 937-980. - Crolet, J .M., B. Aoubiza, and A. Meunier. “Compact bone: Numerical simulation of mechanical characteristics”. Journal of Biomechanics, v 26, n 6, Jun, 1993, p 677-687. Francfort, G., and F. Murat. ”Homogenization and optimal bounds in linear elasticity”. Archives of Rational Mechanics Analysis 94 (1986),307-334. 68 Grabovsky, Y., and RV. Kohn. ”Microstructures minimizing the energy of a two phase elastic composite in two space dimensions. 11: The Vigdergauz microstructure”. Journal of the Mechanics and Physics of solids 43(6), (1995), 949-972. Guedes, J.M., and N. Kikuchi. "Preprocessing and Postprocessing for Materials Based on the Homogenization Method with Adaptive Finite Element Method". Comp. Meth. Appl. Mech. Eng.Vol.83 (1990), 143-198. Hashin, Z. ”The elastic moduli of the heterogeneous materials”. ASME Journal of Applied Mechanics 29 (1962), 143-150. Hashin, Z. “ Theory of composite materials”. Mechanical of Composite Materials, F.W.Wend, H.Liebowitz and N.Perrone, eds. Oxford, Pergamon. (1970) Hashin, Z., and S. Shtrikman. "A variational approach to the theory of the elastic behaviour of multiphase materials". Journal of the Mechanics and Physics of Solids (1963), 127-140 Hill, R. ” The elastic behavior of the crystalline aggregate”. Proceedings of the Royal Society of London A65 (1952), 349-354 Lions, J .L. ”Some Methods in the Mathematical Analyses of System and their Control”. New York: Gordon and Breach 1981 Lutie, K.A. , and A.V. Cherkaev . ”The problem of formation of an optimal isotropic multi-component composite”. Technical report 895.A.F.Ioffe Physical Technical Institute, Acad. Sci. of the USSR, Leningrad (1984). Shorter version in J .opt.Th. Appl 46 (1985), 571-589. Milton, G.W. ”Composite materials with Poisson’s ratios close to -1”. Journal of the Mechanics and Physics of Solids 40(5) (1992),1105-1137. 69 Milton, G.W., and A.V. Cherkaev. ”Which elasticity tensors are realizable?”. Journal of Engineering Materials and Technology, Transactions of theASME 117(4)(1995), 483- 493. Milton, G.W., and N. Phan-Thien. "New bounds on the effective elastic moduli of two- component materials ". Proceedings of the Royal Society of London A380 (1982), 305- 33 1. Oleinik, O.A. “On homogenization problems”. Trends and Applications of Pure Mathematics to Mechanics, P.G. Ciarlet and M. Rouseau,eds., Berlin: Springer 1984. Sigmund, 0. ”Materials with Prescribed Constitutive Parameters; An Inverse Homogenization Problem”. DCAM Report 470, Technical University of Denmark 1993. Sigmund, O. ”Tailoring Materials with Prescribed elastic Properties”. DCAM Report 480, Technical University of Denmark 1994. Sigmund, O. “ A new class of extremal composites”. Journal of the Mechanics and Physics of Solids 48 (1999), 397-428. Svanverg, K. ”The method of moving asymptotes”. Int. J. Num. Meth. Eng. 24 (1987), 359-373 Vigdergauz, S.B. ”Energy-minimizing inclusions in a planar elastic structure with macro- isotropy”. Structure Optimization 17(203), (1999), 104-112. Walpole, L.J. “On bounds for the overall elastic moduli of inhomogeneous systems-I”. J. Mech. Phys. Solids. Vol.14 (1966), 151-162. 70