- Email: [email protected]
Topology optimization with design-dependent loads
Finel=812=Ravi=Venkatachala=BG
Finite Elements in Analysis and Design 37 (2001) 57}70
Topology optimization with design-dependent loads Bing-Chung Chen*, Noboru Kikuchi Department of Mechanical Engineering and Applied Mechanics, The University of Michigan, Ann Arbor, MI, 48109-2125, USA
Abstract This paper presents a generalization of topology optimization of a linearly elastic structure subject to design-dependent loads. The direction and location of the design-dependent loads may change as the shape of the structure changes. A new approach based on applying a "ctitious thermal loading to simulate this dependent load is proposed. The optimization of the structural topology is transformed to a three-phase material distribution problem within the design domain in which the solid, void, and hydrostatic #uid phases are optimally distributed. Examples of both internally and externally pressurized structures are presented. 2001 Elsevier Science B.V. All rights reserved. Keywords: Topology optimization; Design-dependent loads; Pressure
1. Introduction In the last decade, the "eld of structure topology optimization has expanded signi"cantly, successfully addressing many practical engineering problems. The explosive expansion of this "eld can be traced back to the pioneer work by Kohn and Strang [1,2], Bends+e and Kikuchi [3]. This approach has been applied to optimize structures subject to the single static load [4] and multiple loads [5] using the homogenization method. Bends+e et al. [6,7] generalized the homogenization method using the most general form of elastic tensors as the material model to simultaneously optimize the local material properties and overall structure for minimum compliance. The problem of tailoring structures for optimal vibration response has also been addressed [8,9]. For a comprehensive review on the structural topology optimization, readers are referred to the monograph by Bends+e [10] and the references therein.
* Corresponding author. E-mail addresses: [email protected] (B.-C. Chen), [email protected] (N. Kikuchi). 0168-874X/01/$ - see front matter 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 0 ) 0 0 0 2 1 - 4
Finel=812=Ravi=VVC
58
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
Fig. 1. The structural topology optimization subject to the "xed and design-dependent loads. (a) shows the "xed extended design domain with the boundary condition and the "xed load. (b) is a possible optimal topology for a "xed loading (c) is a possible optimal design where the pressure is allowed to change the direction and location as the loaded surface changes.
Nonetheless, existing methods have not solved the class of problems in which the applied loads depend on the design of the structure itself. The design-dependent loads acting on a structure may alter their direction and location of applications as the shape of the structure changes. This class of problems manifests its importance in many engineering disciplines. The hydrostatic pressure acting on the dam or the underwater container, the weight of snow acting on the roof are typical examples of such design-dependent loads. Generally speaking, all structures involving solid and #uid interaction including ducts, pipes, and airfoils carry such design-dependent loads. The reason for not addressing this problem adequately may be attributed to the way the topology optimization is formulated. The underlying scheme of topology optimization [3] is based on a speci"c choice of a "xed reference region upon which the optimal material distribution is to be found. This "xed region, the extended design domain, on which the loads and constraints are speci"ed should be properly chosen to accommodate the loads and constraints as shown in Fig. 1(a). This paradigm has been regarded as insu$cient to address the problems involving designdependent loads such as the hydrostatic pressure [10]. With the existing approach, one can "nd only the optimal structure subject to the "xed load as illustrated in Fig. 1(b) while fails to design a structure subject to the design-dependent load as the loaded surface changes as shown in Fig. 1(c). The strategy of simultaneous design of the shape of the reference domain and the topology has been proposed since then [10]. The implication of such setting requires very involved sensitivity analysis [11] and incorporating the re-meshing capacity into the analysis code. More recently, Hammer and Olho! [12] proposed a new approach to solve this problem. Their formulation retains the "xed grid scheme to avoid the complexity inherited from the domain shape change. In their work, the pressure acts on a smooth surface extracted from the iso-volumetric density curve and is represented by the spline functions whose control points are determined by the density distribution. The end points or curves of this pressure loaded surface are de"ned a priori. The pressure acting on this smooth surface is then transformed into nodal forces in the "nite element analysis model. As a result, the sensitivity analysis becomes more involved under such settings. Since the pressure acts on this parameterized surface, the sensitivity of the pressure load with respect to the change of the element density needs to go through the relation of iso-volumetric density and the control points of the spline functions.
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
59
In this work, we propose a new approach to simulate the design-dependent loads while staying within the context of the classical topology optimization paradigm. The "xed mesh scheme is retained to reduce the computational complexity associated with the domain shape change. The design-dependent loads are simulated by the "ctitious thermal loads rather than constructing a parameterized surface for the pressure to act on. The structural topology will be optimized for minimum compliance with the constraint on the structural volume. The problem is transformed to a three-phase material distribution problem within the design domain in which the solid, void, and hydrostatic #uid phases are optimally distributed. The underlying material model is based on an extension of the solid isotropic material with penalization (SIMP) [13] to simulate the three-phase material distribution. An algorithm is proposed to track the interface between the #uid and solid regions. Examples of both the internally and externally pressurized structures will be presented. It will also be demonstrated that the proposed approach can be used to simulate the uni-directional loading in addition to the hydrostatic pressure loads.
2. General problem formulation and the design-dependent loads Let X be the design domain composed of three distinctive regions, the solid X , the void X and the #uid X . Each region is non-overlapping, i.e., X"X 6X 6X , X 5X ", X 5X ", X 5X ". The material properties are homogeneous within each region. The #uid exerts the hydrostatic pressure p on the solid region along the interface C as shown in Fig. 2(a). One can N introduce the energy bilinear form
a(*, u)"
X
v C u dX" G H GHIJ I J
X
e (*)C e (u) dX GH GHIJ IJ
with the linearized strain
1 *u *u G# H e (u)" GH 2 *x *x H G
Fig. 2. Schematics of the pressure loading on a structure. (a): Three-phase material distribution. (b): Three-phase decomposition. The domain is composed of #uid and non-#uid region, exerting pressure at the interface. The non-#uid region is further decomposed into solid and void as in the classical two-phase topology optimization.
Finel=812=Ravi=VVC
60
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
Fig. 3. The loads acting on the solid and #uid region when subject to the unit temperature rise. (a): Two homogeneous regions with the common interface C and the free boundary C , C . (b): Thermal forces acting on the free boundaries N @ @ and the common interface.
and the load linear form due to the hydrostatic pressure p
¸(*)"
C
N
v (!pn ) dS. G G
(1)
The equilibrium equation written in the variational form is a(*, u)"¸(*) ∀*3K (X), (2) where the elastic tensor C has the major and minor symmetries, i.e., C "C "C ; the GHIJ GHIJ IJGH GHJI space K (X) is the space of kinematically admissible displacement "elds; n is the unit outward G normal vector of the interface C and (!pn ) is the hydrostatic pressure force acting perpendicular N G to the surface C . The goal of the present work is to simulate this pressure force with a "ctitious N thermal load. The hydrostatic pressure force is simulated by the thermal load due to the mismatch of thermal expansion coe$cients of two materials along the material interface. Let the domain X be composed of two regions, the region 1, X , and the region 2, X , made of homogeneous material phase 1 and phase 2, respectively. Region 1 joins region 2 on the common interface C with the unit outward N normal vector n. The phase 1 and phase 2 regions have the free boundary C, C, respectively, as @ @ shown in Fig. 3(a). The thermal stress tensors band b are homogeneous in each region where GH GH b"C a, b"C a and C , a, C , a are elastic and thermal expansion coe$cGH GHIJ IJ GH GHIJ IJ GHIJ IJ GHIJ IJ ient tensors of each phase, respectively. The thermal virtual work associated with the uniform unit temperature rise can be expressed as
X
v b dX" G H GH
X
v b dX# G H GH
X
v b dX. G H GH
(3)
Applying the Gauss divergence theorem to each region individually, one can rewrite the thermal virtual work on the homogeneous region 1 as
X
v b dX" G H GH
C @
v bn dS# G GH H
C
N
v bn dS! G GH H
X
v b dX. G GH H
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
61
The last term vanishes due to homogeneity. Similarly, one can transform the virtual work on the region 2. Then Eq. (3) can be simpli"ed by the relation that the outward normal n of the interface C in phase 2 region point to the opposite direction as in region 1: N
X
v b dX" G H GH
C @
v bn dS# G GH H
C @
v bn dS# G GH H
C
N
v (b!b)n dS. G GH GH H
(4)
The "rst two terms are due to the thermal force acting on the free boundary and the last term derives from the mismatch of thermal stress tensors along the interface C as shown in Fig. 3(b). By N manipulating terms in Eq. (4), one can de"ne the `pseudoa-linear form ¸ (*) as the virtual work due N to thermal load subtracting the virtual work done by the forces acting on the free boundary as
¸ (*)" N
X
"
C
v b dX! G H GH
N
C @
v bn dS# G GH H
C @
v bn dS G GH H
v (b!b)n dS. G GH GH H
For a special case, when both the b and b are isotropic, namely, b"bd , b"bd , and GH GH GH GH GH GH b and b di!er by the amount of p, then the `pseudoa-linear form ¸ (*) is identical with the N linear form ¸(*) due to the hydrostatic pressure p as in Eq. (1):
¸ (*)" N
C
N
v (b!b)n dS" G G
C N
v (!pn ) dS"¸(*). G G
In the actual implementation, it is straightforward to single out the forces acting on the free boundary. One can construct the load vectors in the "nite element calculation of a homogeneous solid subject to the constant temperature change under the free expansion condition and identify the non-vanishing forces. These are the forces acting on the free boundary and are intended to be eliminated in the `pseudoa-linear form ¸ (*). N It should be noted that this formulation is versatile enough to simulate many di!erent loading conditions in addition to the hydrostatic pressure. By manipulating the thermal stress tensors b and b, one can simulate the interface force with directional preference, in particular the uni-directional force. If the thermal stress tensor has only one non-zero component, for example, b"be e , b"be e , the interface force acts in the x -direction only. 3. The material model for the three-phases design problem When extending this two-phase result to the three-phase system as in Fig. 2(a), some special techniques need to be applied. The hydrostatic pressure forces act only at the interface between the solid and #uid regions and no interface forces should exist between the solid and void regions. In other words, it is free of interface force within the non-#uid region. Consequently, one needs to set the thermal stress tensor of the non-#uid area constant throughout the region regardless of the density distribution.
Finel=812=Ravi=VVC
62
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
In the case of optimal topology design, one seeks the optimal elastic tensor C distribution over GHIJ the domain of the body, treating C (x) as a function of the spatial variable x3X. Bends+e and GHIJ Kikuchi [3] used the homogenization method to compute the local elastic tensor by assuming that the structure is a porous composite with voids at the micro level. At each local point, this composite is then modeled by a microscopic square cell with a rectangular hole. The dimensions of the hole and orientation of the micro cell are subject to change during the optimization process. The material property of such composite is orthotropic in nature. In this work, however, to simulate the hydrostatic pressure, the thermal stress tensor needs to be isotropic. To realize such material characteristics, an arti"cial mixing assumption for the local elastic CC , thermal expansion GHIJ coe$cient aC and thermal stress tensors bC in element e of the "nite element model are stated as GH GH CC (o )"(o )NC , GHIJ C C GHIJ 1 a, aC (o )" GH C (o )N GH C
bC (o )"C a"b, GH C GHIJ IJ GH
(5)
where C , a, b are elastic, thermal expansion coe$cient and thermal stress tensors of the base GHIJ IJ GH solid material. The local density o 3[o , 1] is assumed to be constant for element e. The element C
is regarded as a full solid material with o "1 while it denotes void for o "o . The small C C
number o is used to prevent singularity in the "nite element calculation. This local material
mixing assumption can be regarded as a direct extension of the SIMP material model, provided the base material properties C , a are isotropic. The power p is a penalty factor to discourage the GHIJ IJ intermediate density in the "nal solution. It has been shown [13] that for an isotropic material with the Poisson's ratio v", this penalty factor must satisfy p*3 to stay within the Hash in}Shtrikman bounds [14], thus making it physically realizable. It should be noted that the a may GH not necessarily be isotropic in this work. To simulate the interface force with directional preference, one needs to furnish the a with peculiar properties other than isotropy. If the desired thermal GH stress tensor b has only one non-zero component, b "b "0, b "b, for example, while the GH elastic tensor C remains isotropic with the Poisson's ratio v, the components of a can be set as GHIJ GH a "a "(!v/(1#v))a . In the #uid region, we still use the classical SIMP to model the local elastic properties while a modi"ed SIMP material law is employed to simulate the local thermal stress tensor as CC (o )"(o )NC , GHIJ C C GHIJ bC (o )"[h#(1!h)oN]b. GH C GH
(6)
The constant h3[0,1] is a `lift-upa coe$cient to modify the power-law curve of the classical SIMP model. This `lift-upa SIMP material model serves two important purposes. First, this coe$cient places the bC within a speci"c range of values for a wider spectrum of density values. This GH characteristic is vital to prevent excessive dispersion of the interface force throughout the #uid region if the classical SIMP is used instead. The classical SIMP corresponds to the curve for h"0 as illustrated in Fig. 4 where the material property spans a wide range of values, 0)bC (o ))b. GH C GH For the `lift-upa curve, on the other hand, the span of allowable material properties is more con"ned as hb)bC (o ))b. It can be readily shown that under such a scheme, the pressure GH GH C GH
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
63
Fig. 4. The `lift-upa SIMP curve of b /bQ , vs o for p"3. The classical SIMP curve corresponds to h"0.
acts only on a small range of density contour curves. As a result, we are relieved from the cumbersome task of constructing the iso-volumetric density curve for the pressure to act on since the process of extracting a speci"c range of density contour is automatically imbedded in this formulation. The second advantage of using this material model is that it discredits the deployment of the intermediate density within the design domain, especially along the #uid and solid interface. Elements with intermediate density result in similar pressure forces along the solid and #uid interface compared with the void elements, while the intermediate density elements consume more volumetric resources than the void elements. If the total volume ratio of the solid phase is imposed as a constraint, the intermediate density is less e$cient in terms of volumetric resource allocation, thus making it less favorable for the optimizer. Similarly, this material model also has the e!ect of eliminating the isolated intermediate density elements in the #uid region. An additional parameter is necessary to distinguish between the #uid and non-#uid region. The `dryness coe$cienta m for each "nite element is introduced for this purpose. Those elements free C from #uid entrenchment are designated with m "1, regardless of their local density o while the C C #uid elements are distinguished with m "0. It should be emphasized that this coe$cient is used C primarily to track the #uid and solid interface and is not a design variable in the optimization process. The coe$cient m for each "nite element can be regarded as an implicit function of its own C local density o and all the coe$cients m of its adjacent elements. An algorithm to determine this C C coe$cient will be described in detail later. Incorporating the material law in the #uid and non-#uid regions into Eqs. (6) and (5), the local material properties in element e can be written as functions of the coe$cients o and m : C C CC (o , m )"(o )NC , GHIJ C C C GHIJ bC (o , m )"m b#(1!m )b GH C C C GH C GH "+m #(1!m )[h#(1!h)oN],b, C C GH where the thermal stress tensor bC is the arithmetic mean of the #uid b and the solid b. The GH GH GH material properties and the interface force for this three-phase system are summarized in Table 1 where the diagonal items denote the elastic and thermal stress tensors of each phase and the o!-diagonal items depict the interface force between the constituent phases.
Finel=812=Ravi=VVC
64
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
Table 1 Material properties of the three-phase material system Regions
Solid (m"1, o"1)
Void/gray (m"1)
Fluid (m"0)
Solid
C b"C : a
No interface force b!bT"0
Interface force existed b!b"(1!h)b
CT"oNC bT"b
Interface force existed b!bT"[oN#(1!oN)h!1]b
Void/gray Fluid
C"oNC b"[oN#(1!oN)h]b
4. The optimization procedure Since the purpose of the present paper is to demonstrate the extension of classical topology optimization with design-dependent load, we shall follow to a large extent the standard procedure of the classical mean-compliance-based topology optimization [4]. The optimization problem consists of minimizing the total potential energy of the structure subject to the volumetric constraints P(*)"a(*, *)!¸ (*) N
max min MC *Z) s.t.
X o(x) dX)< ,
0(o )o )1.
C The equilibrium equation, Eq. (2), is expressed in terms of the principle of minimum potential energy in which the displacement "eld u is a minimizer of the total potential energy functional P(*): min P(*)"P(u) 0 a(*, u)"¸ (*) ∀*3K (X) N * Z) and at equilibrium, the total potential energy becomes
(7)
P(u)"a(u, u)!¸ (u)"!¸ (u)"!a(u, u). N N Then the optimization problem is equivalent to seeking the minimum of the mean compliance ¸ (u): N
¸ (u)" N
C N
u (b!b)n dX" G GH GH H
X
e (u)b dX# GH GH
C @
u (!b n ) dS, G GH H
where C is the free boundary of the design domain. The mean compliance can be written in terms @ of the design variable o and the displacement and strain "eld u, e(u):
¸ (u)" N
X
g(e, o) dX#
C @
h(u, o) dS.
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
65
When g, h are continuous and di!erentiable functions of their arguments, the "rst variational of ¸ (u) equals N *g *h *g *h d¸ " do# de dX# do# du dS. N *e *u *o X *o C @
"
X
e
*b GH do#b de dX# GH *o GH GH
C @
*b u ! GH n do#(!b n )du dS. G GH H G *o H
(8)
The explicit terms, (*g/*o)do and (*h/*o)do, depend on the change in design o directly while the implicit parts, (*g/*e)de and (*h/*u)du, depend on the change through the strain and displacement "elds. To evaluate the implicit parts, the procedure of adjoint variable sensitivity analysis is employed [15,16]. As a result, the "rst variation of the mean compliance is
d¸ (u)"! N
*C e (u) GHIJ doe (u) dX#2 IJ *o X GH
#2
*b e (u) GH do dX *o X GH
*b u ! GH don dS. G H *o C @
For ease of implementation, we use the sequential linear programming (SLP) method to solve this optimization problem, which consists of solving a linear approximation of the original problem sequentially. In each iteration step, the optimization problem is linearized around the current design point. The optimal design change in each iteration is found by solving the linear programming problem with bounds on the design variables. The design variables are the local densities of each element o . It should be noted that the coe$cient m is not a design variable itself and is used C C to track the interface of the #uid region. The determination of the coe$cient m depends on its own C local density o and the coe$cients m of its surrounding elements. The updating algorithm of m is C C C somewhat like the process of #uid #ooding over the design domain. The initial design domain is always composed of a "xed area from which the #uid originates, and is not subject to change during the subsequent optimization iterations. The #uid will propagate to all the neighboring
Fig. 5. The schematic of the #uid #ooding algorithm to track the interface of solid and #uid regions. (a): The #uid originates from the "xed #uid region and will propogate to the adjacent neighbors. This propogating process is stopped at the solid interface. (b): After the propagating process, the design domain is separated into three regions, the solid, void, and #uid.
Finel=812=Ravi=VVC
66
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
elements unless stopped by the solid elements which the #uid cannot penetrate. The #uid propagating process is illustrated in Fig. 5. By using this algorithm, one can keep track of the interface of the #uid region. The pressure force is kept "xed during the "rst 25 iterations and results in a constant distributed load acting on the boundary of the "xed #uid region within the design domain. In other words, it is not until the 26th iteration that the `#ood-overa algorithm will be applied to track the #uid interface. The primary reason for this scheme is to prevent the #uid region from #ooding all over the domain, should the #ood-over algorithm be applied in the "rst place when there is no signi"cant structure acting as a containment for the #uid region.
5. Numerical examples In the following numerical examples, we apply the aforementioned method to design internally and externally pressurized structures. It will also be demonstrated that the proposed scheme can be used to simulate the uni-directional loading besides the hydrostatic pressure load. The "lter algorithm, as described in the Petersson and Sigmund paper [17], is applied to prevent the checkerboard and mesh-dependent problems in the optimized topology. 5.1. Externally pressurized structure design In this example, an underwater structure subject to the hydrostatic pressure is optimized. The initial setup of the problem is shown in Fig. 6 where the gray area denotes the "xed #uid region. The "nite element analysis model consists of 1800 quad elements. The structure is "xed near the two end corners at the bottom. The hydrostatic pressure forces come from the top, left and right sides. The thermal stress tensor of the solid region is chosen to be isotropic, i.e., b"!bd , to GH GH simulate the hydrostatic pressure acting in the direction of inward normal of the solid-phase interface. The `lift-upa coe$cient h in the modi"ed SIMP material model is chosen to be 0.5 to prevent intermediate density elements. The solid region shall occupy no more than 30% of the total volume in the design domain. The optimized structure has a near semi-circular topology as shown in Fig. 7(a) and conforms to the shape of underwater structures widely adopted in the engineering practice. Fig. 7(b) illustrates
Fig. 6. The design domain with the "xed #uid region located on top, left and right sides.
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
67
Fig. 7. The underwater structure subject to the hydrostatic pressure. (a): The optimal structural topology. (b): The region of non-#uid phase and the hydrostatic pressure force vectors acting on the structure.
Fig. 8. The civil structure subject to the uni-directional snow weight. (a): The optimal structure with an arch-like topology (b): The region of non-#uid phase and the uni-directional force vectors acting on the structure.
the non-#uid phase region and the corresponding hydrostatic pressure force vectors acting at the interface. It is noted that Hammer and Olho! [12] have also designed a similar structure but with a less circular shape. In the following, we will use the identical problem setting of the previous example (Fig. 6) to design a civil structure subject to the uni-directional load due to the snow's weight by simply modifying the thermal stress tensor b of the solid phase. Instead of choosing the b to be GH GH isotropic, one can set the b to have only one non-zero component, namely, b "0, b "!b, GH b "b "0 to imitate the downward snow weight acting on a house roof in the !x direction. The optimized topology is an arch-like structure as shown in Fig. 8(a). The non-#uid phase region and the downward force vectors are illustrated in Fig. 8(b). It is interesting to note that the structure is similar to the arch computed by Rozvany and Prager [18]. 5.2. Internally pressurized structure design In this example, we will design a cover-like structure subject to the internal pressure acting on the top. The problem setup is shown in Fig. 9(a) with the "xed #uid region in the rectangular shape on the top. The structure is clamped in both endpoints at the top corners. The "nite element analysis model consists of 2400 quad elements. The solid region is to use no more than 30% of the total
Finel=812=Ravi=VVC
68
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
Fig. 9. The internally pressurized container design. (a): The design domain and the initial pressure distribution. (b): The optimized topology.
Fig. 10. The internally pressurized container design with enlarged "xed #uid region. (a): The design domain and the initial pressure distribution. (b): The optimized topology.
volume. Again, the thermal stress tensor of the solid region is chosen to be isotropic to simulate the hydrostatic pressure. The optimized topology is shown in Fig. 9(b). Unlike the cover-like structure obtained in the Hammer and Olho! work [12], the optimized structure bears a distant resemblance to the MBB beam [19] and has a relative small #uid-phase area. The structure is basically composed of two overlapping arch trusses. The larger inverted arch truss spans the two supports and the smaller one is located adjacent to the "xed #uid region to resist the pressure force acting on the top. These two arch trusses are further reinforced by two additional supporting trusses connecting them. To form a cover-like structure, we enlarge the "xed #uid region to form an inverted triangle shape as shown in Fig. 10(a). The optimized topology is a cover-like structure as shown in Fig. 10(b). By extending the "xed #uid-phase further into the design domain, we essentially rule out this occupied area as a candidate region for the structure to reside. It should be noted that the mean compliance of the double-arch structure of Fig. 9(b) is 20.26% smaller than that of the cover-like structure of Fig. 10(b). In other words, the double-arch structure is sti!er when subjected to the hydrostatic pressure.
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
69
Fig. 11. The #uid chamber design. (a): The design domain and the "xed #uid and "xed solid regions. (b): The optimal chamber pro"le with inlet from a channel to a bigger pressure chamber.
5.3. Flow chamber design The last example models the internally pressurized structure with inlet from a channel to a larger chamber as shown in Fig. 11(a). The inlet is designated as the solid region and is not subject to change in the optimization process as shown by the black color, while the gray color denotes the "xed #uid region. The "nite element analysis model consists of 2280 quad elements and the solid region occupies no more than 25% of the total volume. The thermal stress tensor of the solid region is isotropic. The optimized topology is shown in Fig. 11(b). The optimized structure is similar to the result in the Hammer and Olho! work [12]. 6. Conclusion A generalization of the classical topology optimization problem of designing structures subject to design-dependent loads has been presented. A "ctitious thermal force is used to simulate this design-dependent load such as the hydrostatic pressure force. The formulation stays within the context of the classical topology optimization paradigm and uses the "xed "nite element mesh in the optimization process. The major advantage is to circumvent the computational complexity associated with the domain shape variation if the shape of the analysis domain changes. A modi"ed SIMP material model has been proposed to con"ne the pressure force acting on a small range of density contour curves without explicitly constructing the iso-volumetric density curve. Both internally and externally pressurized structures have been designed. It has also been demonstrated that by manipulating the components of the thermal stress tensor, one can simulate a wide variety of interface force loading conditions in addition to the hydrostatic pressure. This work serves as the starting point to design structures involving solid and #uid interaction in which the magnitude of interaction force is a design-dependent quantity itself. Acknowledgements The "nancial support of DARPA for this research under Grant No. DARPA/NIST F000765 is gratefully acknowledged.
Finel=812=Ravi=VVC
70
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
References [1] R.V. Kohn, G. Strang, Optimal design in elasticity and plasticity, Numer. Methods Eng. 22 (1986) 183}188. [2] R.V. Kohn, G. Strang, Optimal design and relaxation of variational problems, Comm. Pure Appl. Math (New York). 39 (1986) 1}25 (Part I), 139}182 (Part II), 353}357 (Part III). [3] M.P. Bends+e, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Eng. 71 (2) (1988) 197}224. [4] K. Suzuki, N. Kikuchi, A homogenization method for shape and topology optimization, Comput. Methods Appl. Mech. Eng. 93 (1991) 291}318. [5] A. Diaz, M.P. Bends+e, Shape optimization of structures for multiple loading situations using a homogenization method, Struct. Optim. 4 (1992) 17}22. [6] M.P. Bends+e, J.M. Guedes, R.B. Haber, P. Pedersen, J.E. Taylor, An analytical model to predict optimal material properties in the context of optimal structural design, Trans. ASME J. Appl. Mech. 61 (4) (1994) 930}937. [7] M.P. Bends+e, A.R. DmH az, R. Lipton, J.E. Taylor, Optimal design of material properties and material distribution for multiple loading conditions, Int. J. Numer. Methods Eng. 38 (7) (1995) 1149}1170. [8] A. Diaz, N. Kikuchi, Solution to shape and topology eigenvalue optimization problems using a homogenization method, Int. J. Numer. Methods Eng. 35 (1992) 1487}1502. [9] Z.-D. Ma, N. Kikuchi, I.H Hagiwara, Structural topology and shape optimization for a frequency response problem, Comput. Mech. 13 (3) (1993) 157}174. [10] M.P. Bends+e, Optimization of Structural Topology, Shape, and Material, Springer, London, UK, 1995. [11] M.P. Bends+e, J. Soko"owski, Shape sensitivity analysis of optimal compliance functionals, Mech. Struct. Mach. 23 (1) (1995) 35}58. [12] V.B. Hammer, N. Olho!. Topology optimization with design dependent loads, WCSMO-3 Bu!alo, 1999, pp. 629}631. [13] M.P. Bends+e, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. 69 (1999) 635}654. [14] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behavior of multiphase materials, J. Mech. Phys. Solids 11 (1963) 127}140. [15] K. Dems, Z. MroH z, Variational approach by means of adjoint systems to structural optimization and sensitivity analysis * I, Int. J. Solids Struct. 19 (8) (1983) 677}692. [16] K.K. Choi, E.J. Haug, V. Komkov, Design Sensitivity Analysis of Structural Systems, Academic Press, London, 1986. [17] J. Petersson, O. Sigmund, Slope constrained topology optimization, Int. J. Numer. Methods Eng. 41 (8) (1998) 1417}1434. [18] G.I.N. Rozvany, W. Prager, New class of structural optimization problems: optimal archgrids, Comput. Methods Appl. Mech. Eng. 19 (1) (1979) 127}150. [19] N. Olho!, M.P. Bends+e, J. Rasmussen, On CAD-integrated structural topology and design optimization, Comput. Methods Appl. Mech. Eng. 89 (1992) 259}279.
Finite Elements in Analysis and Design 37 (2001) 57}70
Topology optimization with design-dependent loads Bing-Chung Chen*, Noboru Kikuchi Department of Mechanical Engineering and Applied Mechanics, The University of Michigan, Ann Arbor, MI, 48109-2125, USA
Abstract This paper presents a generalization of topology optimization of a linearly elastic structure subject to design-dependent loads. The direction and location of the design-dependent loads may change as the shape of the structure changes. A new approach based on applying a "ctitious thermal loading to simulate this dependent load is proposed. The optimization of the structural topology is transformed to a three-phase material distribution problem within the design domain in which the solid, void, and hydrostatic #uid phases are optimally distributed. Examples of both internally and externally pressurized structures are presented. 2001 Elsevier Science B.V. All rights reserved. Keywords: Topology optimization; Design-dependent loads; Pressure
1. Introduction In the last decade, the "eld of structure topology optimization has expanded signi"cantly, successfully addressing many practical engineering problems. The explosive expansion of this "eld can be traced back to the pioneer work by Kohn and Strang [1,2], Bends+e and Kikuchi [3]. This approach has been applied to optimize structures subject to the single static load [4] and multiple loads [5] using the homogenization method. Bends+e et al. [6,7] generalized the homogenization method using the most general form of elastic tensors as the material model to simultaneously optimize the local material properties and overall structure for minimum compliance. The problem of tailoring structures for optimal vibration response has also been addressed [8,9]. For a comprehensive review on the structural topology optimization, readers are referred to the monograph by Bends+e [10] and the references therein.
* Corresponding author. E-mail addresses: [email protected] (B.-C. Chen), [email protected] (N. Kikuchi). 0168-874X/01/$ - see front matter 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 0 ) 0 0 0 2 1 - 4
Finel=812=Ravi=VVC
58
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
Fig. 1. The structural topology optimization subject to the "xed and design-dependent loads. (a) shows the "xed extended design domain with the boundary condition and the "xed load. (b) is a possible optimal topology for a "xed loading (c) is a possible optimal design where the pressure is allowed to change the direction and location as the loaded surface changes.
Nonetheless, existing methods have not solved the class of problems in which the applied loads depend on the design of the structure itself. The design-dependent loads acting on a structure may alter their direction and location of applications as the shape of the structure changes. This class of problems manifests its importance in many engineering disciplines. The hydrostatic pressure acting on the dam or the underwater container, the weight of snow acting on the roof are typical examples of such design-dependent loads. Generally speaking, all structures involving solid and #uid interaction including ducts, pipes, and airfoils carry such design-dependent loads. The reason for not addressing this problem adequately may be attributed to the way the topology optimization is formulated. The underlying scheme of topology optimization [3] is based on a speci"c choice of a "xed reference region upon which the optimal material distribution is to be found. This "xed region, the extended design domain, on which the loads and constraints are speci"ed should be properly chosen to accommodate the loads and constraints as shown in Fig. 1(a). This paradigm has been regarded as insu$cient to address the problems involving designdependent loads such as the hydrostatic pressure [10]. With the existing approach, one can "nd only the optimal structure subject to the "xed load as illustrated in Fig. 1(b) while fails to design a structure subject to the design-dependent load as the loaded surface changes as shown in Fig. 1(c). The strategy of simultaneous design of the shape of the reference domain and the topology has been proposed since then [10]. The implication of such setting requires very involved sensitivity analysis [11] and incorporating the re-meshing capacity into the analysis code. More recently, Hammer and Olho! [12] proposed a new approach to solve this problem. Their formulation retains the "xed grid scheme to avoid the complexity inherited from the domain shape change. In their work, the pressure acts on a smooth surface extracted from the iso-volumetric density curve and is represented by the spline functions whose control points are determined by the density distribution. The end points or curves of this pressure loaded surface are de"ned a priori. The pressure acting on this smooth surface is then transformed into nodal forces in the "nite element analysis model. As a result, the sensitivity analysis becomes more involved under such settings. Since the pressure acts on this parameterized surface, the sensitivity of the pressure load with respect to the change of the element density needs to go through the relation of iso-volumetric density and the control points of the spline functions.
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
59
In this work, we propose a new approach to simulate the design-dependent loads while staying within the context of the classical topology optimization paradigm. The "xed mesh scheme is retained to reduce the computational complexity associated with the domain shape change. The design-dependent loads are simulated by the "ctitious thermal loads rather than constructing a parameterized surface for the pressure to act on. The structural topology will be optimized for minimum compliance with the constraint on the structural volume. The problem is transformed to a three-phase material distribution problem within the design domain in which the solid, void, and hydrostatic #uid phases are optimally distributed. The underlying material model is based on an extension of the solid isotropic material with penalization (SIMP) [13] to simulate the three-phase material distribution. An algorithm is proposed to track the interface between the #uid and solid regions. Examples of both the internally and externally pressurized structures will be presented. It will also be demonstrated that the proposed approach can be used to simulate the uni-directional loading in addition to the hydrostatic pressure loads.
2. General problem formulation and the design-dependent loads Let X be the design domain composed of three distinctive regions, the solid X , the void X and the #uid X . Each region is non-overlapping, i.e., X"X 6X 6X , X 5X ", X 5X ", X 5X ". The material properties are homogeneous within each region. The #uid exerts the hydrostatic pressure p on the solid region along the interface C as shown in Fig. 2(a). One can N introduce the energy bilinear form
a(*, u)"
X
v C u dX" G H GHIJ I J
X
e (*)C e (u) dX GH GHIJ IJ
with the linearized strain
1 *u *u G# H e (u)" GH 2 *x *x H G
Fig. 2. Schematics of the pressure loading on a structure. (a): Three-phase material distribution. (b): Three-phase decomposition. The domain is composed of #uid and non-#uid region, exerting pressure at the interface. The non-#uid region is further decomposed into solid and void as in the classical two-phase topology optimization.
Finel=812=Ravi=VVC
60
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
Fig. 3. The loads acting on the solid and #uid region when subject to the unit temperature rise. (a): Two homogeneous regions with the common interface C and the free boundary C , C . (b): Thermal forces acting on the free boundaries N @ @ and the common interface.
and the load linear form due to the hydrostatic pressure p
¸(*)"
C
N
v (!pn ) dS. G G
(1)
The equilibrium equation written in the variational form is a(*, u)"¸(*) ∀*3K (X), (2) where the elastic tensor C has the major and minor symmetries, i.e., C "C "C ; the GHIJ GHIJ IJGH GHJI space K (X) is the space of kinematically admissible displacement "elds; n is the unit outward G normal vector of the interface C and (!pn ) is the hydrostatic pressure force acting perpendicular N G to the surface C . The goal of the present work is to simulate this pressure force with a "ctitious N thermal load. The hydrostatic pressure force is simulated by the thermal load due to the mismatch of thermal expansion coe$cients of two materials along the material interface. Let the domain X be composed of two regions, the region 1, X , and the region 2, X , made of homogeneous material phase 1 and phase 2, respectively. Region 1 joins region 2 on the common interface C with the unit outward N normal vector n. The phase 1 and phase 2 regions have the free boundary C, C, respectively, as @ @ shown in Fig. 3(a). The thermal stress tensors band b are homogeneous in each region where GH GH b"C a, b"C a and C , a, C , a are elastic and thermal expansion coe$cGH GHIJ IJ GH GHIJ IJ GHIJ IJ GHIJ IJ ient tensors of each phase, respectively. The thermal virtual work associated with the uniform unit temperature rise can be expressed as
X
v b dX" G H GH
X
v b dX# G H GH
X
v b dX. G H GH
(3)
Applying the Gauss divergence theorem to each region individually, one can rewrite the thermal virtual work on the homogeneous region 1 as
X
v b dX" G H GH
C @
v bn dS# G GH H
C
N
v bn dS! G GH H
X
v b dX. G GH H
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
61
The last term vanishes due to homogeneity. Similarly, one can transform the virtual work on the region 2. Then Eq. (3) can be simpli"ed by the relation that the outward normal n of the interface C in phase 2 region point to the opposite direction as in region 1: N
X
v b dX" G H GH
C @
v bn dS# G GH H
C @
v bn dS# G GH H
C
N
v (b!b)n dS. G GH GH H
(4)
The "rst two terms are due to the thermal force acting on the free boundary and the last term derives from the mismatch of thermal stress tensors along the interface C as shown in Fig. 3(b). By N manipulating terms in Eq. (4), one can de"ne the `pseudoa-linear form ¸ (*) as the virtual work due N to thermal load subtracting the virtual work done by the forces acting on the free boundary as
¸ (*)" N
X
"
C
v b dX! G H GH
N
C @
v bn dS# G GH H
C @
v bn dS G GH H
v (b!b)n dS. G GH GH H
For a special case, when both the b and b are isotropic, namely, b"bd , b"bd , and GH GH GH GH GH GH b and b di!er by the amount of p, then the `pseudoa-linear form ¸ (*) is identical with the N linear form ¸(*) due to the hydrostatic pressure p as in Eq. (1):
¸ (*)" N
C
N
v (b!b)n dS" G G
C N
v (!pn ) dS"¸(*). G G
In the actual implementation, it is straightforward to single out the forces acting on the free boundary. One can construct the load vectors in the "nite element calculation of a homogeneous solid subject to the constant temperature change under the free expansion condition and identify the non-vanishing forces. These are the forces acting on the free boundary and are intended to be eliminated in the `pseudoa-linear form ¸ (*). N It should be noted that this formulation is versatile enough to simulate many di!erent loading conditions in addition to the hydrostatic pressure. By manipulating the thermal stress tensors b and b, one can simulate the interface force with directional preference, in particular the uni-directional force. If the thermal stress tensor has only one non-zero component, for example, b"be e , b"be e , the interface force acts in the x -direction only. 3. The material model for the three-phases design problem When extending this two-phase result to the three-phase system as in Fig. 2(a), some special techniques need to be applied. The hydrostatic pressure forces act only at the interface between the solid and #uid regions and no interface forces should exist between the solid and void regions. In other words, it is free of interface force within the non-#uid region. Consequently, one needs to set the thermal stress tensor of the non-#uid area constant throughout the region regardless of the density distribution.
Finel=812=Ravi=VVC
62
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
In the case of optimal topology design, one seeks the optimal elastic tensor C distribution over GHIJ the domain of the body, treating C (x) as a function of the spatial variable x3X. Bends+e and GHIJ Kikuchi [3] used the homogenization method to compute the local elastic tensor by assuming that the structure is a porous composite with voids at the micro level. At each local point, this composite is then modeled by a microscopic square cell with a rectangular hole. The dimensions of the hole and orientation of the micro cell are subject to change during the optimization process. The material property of such composite is orthotropic in nature. In this work, however, to simulate the hydrostatic pressure, the thermal stress tensor needs to be isotropic. To realize such material characteristics, an arti"cial mixing assumption for the local elastic CC , thermal expansion GHIJ coe$cient aC and thermal stress tensors bC in element e of the "nite element model are stated as GH GH CC (o )"(o )NC , GHIJ C C GHIJ 1 a, aC (o )" GH C (o )N GH C
bC (o )"C a"b, GH C GHIJ IJ GH
(5)
where C , a, b are elastic, thermal expansion coe$cient and thermal stress tensors of the base GHIJ IJ GH solid material. The local density o 3[o , 1] is assumed to be constant for element e. The element C
is regarded as a full solid material with o "1 while it denotes void for o "o . The small C C
number o is used to prevent singularity in the "nite element calculation. This local material
mixing assumption can be regarded as a direct extension of the SIMP material model, provided the base material properties C , a are isotropic. The power p is a penalty factor to discourage the GHIJ IJ intermediate density in the "nal solution. It has been shown [13] that for an isotropic material with the Poisson's ratio v", this penalty factor must satisfy p*3 to stay within the Hash in}Shtrikman bounds [14], thus making it physically realizable. It should be noted that the a may GH not necessarily be isotropic in this work. To simulate the interface force with directional preference, one needs to furnish the a with peculiar properties other than isotropy. If the desired thermal GH stress tensor b has only one non-zero component, b "b "0, b "b, for example, while the GH elastic tensor C remains isotropic with the Poisson's ratio v, the components of a can be set as GHIJ GH a "a "(!v/(1#v))a . In the #uid region, we still use the classical SIMP to model the local elastic properties while a modi"ed SIMP material law is employed to simulate the local thermal stress tensor as CC (o )"(o )NC , GHIJ C C GHIJ bC (o )"[h#(1!h)oN]b. GH C GH
(6)
The constant h3[0,1] is a `lift-upa coe$cient to modify the power-law curve of the classical SIMP model. This `lift-upa SIMP material model serves two important purposes. First, this coe$cient places the bC within a speci"c range of values for a wider spectrum of density values. This GH characteristic is vital to prevent excessive dispersion of the interface force throughout the #uid region if the classical SIMP is used instead. The classical SIMP corresponds to the curve for h"0 as illustrated in Fig. 4 where the material property spans a wide range of values, 0)bC (o ))b. GH C GH For the `lift-upa curve, on the other hand, the span of allowable material properties is more con"ned as hb)bC (o ))b. It can be readily shown that under such a scheme, the pressure GH GH C GH
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
63
Fig. 4. The `lift-upa SIMP curve of b /bQ , vs o for p"3. The classical SIMP curve corresponds to h"0.
acts only on a small range of density contour curves. As a result, we are relieved from the cumbersome task of constructing the iso-volumetric density curve for the pressure to act on since the process of extracting a speci"c range of density contour is automatically imbedded in this formulation. The second advantage of using this material model is that it discredits the deployment of the intermediate density within the design domain, especially along the #uid and solid interface. Elements with intermediate density result in similar pressure forces along the solid and #uid interface compared with the void elements, while the intermediate density elements consume more volumetric resources than the void elements. If the total volume ratio of the solid phase is imposed as a constraint, the intermediate density is less e$cient in terms of volumetric resource allocation, thus making it less favorable for the optimizer. Similarly, this material model also has the e!ect of eliminating the isolated intermediate density elements in the #uid region. An additional parameter is necessary to distinguish between the #uid and non-#uid region. The `dryness coe$cienta m for each "nite element is introduced for this purpose. Those elements free C from #uid entrenchment are designated with m "1, regardless of their local density o while the C C #uid elements are distinguished with m "0. It should be emphasized that this coe$cient is used C primarily to track the #uid and solid interface and is not a design variable in the optimization process. The coe$cient m for each "nite element can be regarded as an implicit function of its own C local density o and all the coe$cients m of its adjacent elements. An algorithm to determine this C C coe$cient will be described in detail later. Incorporating the material law in the #uid and non-#uid regions into Eqs. (6) and (5), the local material properties in element e can be written as functions of the coe$cients o and m : C C CC (o , m )"(o )NC , GHIJ C C C GHIJ bC (o , m )"m b#(1!m )b GH C C C GH C GH "+m #(1!m )[h#(1!h)oN],b, C C GH where the thermal stress tensor bC is the arithmetic mean of the #uid b and the solid b. The GH GH GH material properties and the interface force for this three-phase system are summarized in Table 1 where the diagonal items denote the elastic and thermal stress tensors of each phase and the o!-diagonal items depict the interface force between the constituent phases.
Finel=812=Ravi=VVC
64
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
Table 1 Material properties of the three-phase material system Regions
Solid (m"1, o"1)
Void/gray (m"1)
Fluid (m"0)
Solid
C b"C : a
No interface force b!bT"0
Interface force existed b!b"(1!h)b
CT"oNC bT"b
Interface force existed b!bT"[oN#(1!oN)h!1]b
Void/gray Fluid
C"oNC b"[oN#(1!oN)h]b
4. The optimization procedure Since the purpose of the present paper is to demonstrate the extension of classical topology optimization with design-dependent load, we shall follow to a large extent the standard procedure of the classical mean-compliance-based topology optimization [4]. The optimization problem consists of minimizing the total potential energy of the structure subject to the volumetric constraints P(*)"a(*, *)!¸ (*) N
max min MC *Z) s.t.
X o(x) dX)< ,
0(o )o )1.
C The equilibrium equation, Eq. (2), is expressed in terms of the principle of minimum potential energy in which the displacement "eld u is a minimizer of the total potential energy functional P(*): min P(*)"P(u) 0 a(*, u)"¸ (*) ∀*3K (X) N * Z) and at equilibrium, the total potential energy becomes
(7)
P(u)"a(u, u)!¸ (u)"!¸ (u)"!a(u, u). N N Then the optimization problem is equivalent to seeking the minimum of the mean compliance ¸ (u): N
¸ (u)" N
C N
u (b!b)n dX" G GH GH H
X
e (u)b dX# GH GH
C @
u (!b n ) dS, G GH H
where C is the free boundary of the design domain. The mean compliance can be written in terms @ of the design variable o and the displacement and strain "eld u, e(u):
¸ (u)" N
X
g(e, o) dX#
C @
h(u, o) dS.
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
65
When g, h are continuous and di!erentiable functions of their arguments, the "rst variational of ¸ (u) equals N *g *h *g *h d¸ " do# de dX# do# du dS. N *e *u *o X *o C @
"
X
e
*b GH do#b de dX# GH *o GH GH
C @
*b u ! GH n do#(!b n )du dS. G GH H G *o H
(8)
The explicit terms, (*g/*o)do and (*h/*o)do, depend on the change in design o directly while the implicit parts, (*g/*e)de and (*h/*u)du, depend on the change through the strain and displacement "elds. To evaluate the implicit parts, the procedure of adjoint variable sensitivity analysis is employed [15,16]. As a result, the "rst variation of the mean compliance is
d¸ (u)"! N
*C e (u) GHIJ doe (u) dX#2 IJ *o X GH
#2
*b e (u) GH do dX *o X GH
*b u ! GH don dS. G H *o C @
For ease of implementation, we use the sequential linear programming (SLP) method to solve this optimization problem, which consists of solving a linear approximation of the original problem sequentially. In each iteration step, the optimization problem is linearized around the current design point. The optimal design change in each iteration is found by solving the linear programming problem with bounds on the design variables. The design variables are the local densities of each element o . It should be noted that the coe$cient m is not a design variable itself and is used C C to track the interface of the #uid region. The determination of the coe$cient m depends on its own C local density o and the coe$cients m of its surrounding elements. The updating algorithm of m is C C C somewhat like the process of #uid #ooding over the design domain. The initial design domain is always composed of a "xed area from which the #uid originates, and is not subject to change during the subsequent optimization iterations. The #uid will propagate to all the neighboring
Fig. 5. The schematic of the #uid #ooding algorithm to track the interface of solid and #uid regions. (a): The #uid originates from the "xed #uid region and will propogate to the adjacent neighbors. This propogating process is stopped at the solid interface. (b): After the propagating process, the design domain is separated into three regions, the solid, void, and #uid.
Finel=812=Ravi=VVC
66
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
elements unless stopped by the solid elements which the #uid cannot penetrate. The #uid propagating process is illustrated in Fig. 5. By using this algorithm, one can keep track of the interface of the #uid region. The pressure force is kept "xed during the "rst 25 iterations and results in a constant distributed load acting on the boundary of the "xed #uid region within the design domain. In other words, it is not until the 26th iteration that the `#ood-overa algorithm will be applied to track the #uid interface. The primary reason for this scheme is to prevent the #uid region from #ooding all over the domain, should the #ood-over algorithm be applied in the "rst place when there is no signi"cant structure acting as a containment for the #uid region.
5. Numerical examples In the following numerical examples, we apply the aforementioned method to design internally and externally pressurized structures. It will also be demonstrated that the proposed scheme can be used to simulate the uni-directional loading besides the hydrostatic pressure load. The "lter algorithm, as described in the Petersson and Sigmund paper [17], is applied to prevent the checkerboard and mesh-dependent problems in the optimized topology. 5.1. Externally pressurized structure design In this example, an underwater structure subject to the hydrostatic pressure is optimized. The initial setup of the problem is shown in Fig. 6 where the gray area denotes the "xed #uid region. The "nite element analysis model consists of 1800 quad elements. The structure is "xed near the two end corners at the bottom. The hydrostatic pressure forces come from the top, left and right sides. The thermal stress tensor of the solid region is chosen to be isotropic, i.e., b"!bd , to GH GH simulate the hydrostatic pressure acting in the direction of inward normal of the solid-phase interface. The `lift-upa coe$cient h in the modi"ed SIMP material model is chosen to be 0.5 to prevent intermediate density elements. The solid region shall occupy no more than 30% of the total volume in the design domain. The optimized structure has a near semi-circular topology as shown in Fig. 7(a) and conforms to the shape of underwater structures widely adopted in the engineering practice. Fig. 7(b) illustrates
Fig. 6. The design domain with the "xed #uid region located on top, left and right sides.
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
67
Fig. 7. The underwater structure subject to the hydrostatic pressure. (a): The optimal structural topology. (b): The region of non-#uid phase and the hydrostatic pressure force vectors acting on the structure.
Fig. 8. The civil structure subject to the uni-directional snow weight. (a): The optimal structure with an arch-like topology (b): The region of non-#uid phase and the uni-directional force vectors acting on the structure.
the non-#uid phase region and the corresponding hydrostatic pressure force vectors acting at the interface. It is noted that Hammer and Olho! [12] have also designed a similar structure but with a less circular shape. In the following, we will use the identical problem setting of the previous example (Fig. 6) to design a civil structure subject to the uni-directional load due to the snow's weight by simply modifying the thermal stress tensor b of the solid phase. Instead of choosing the b to be GH GH isotropic, one can set the b to have only one non-zero component, namely, b "0, b "!b, GH b "b "0 to imitate the downward snow weight acting on a house roof in the !x direction. The optimized topology is an arch-like structure as shown in Fig. 8(a). The non-#uid phase region and the downward force vectors are illustrated in Fig. 8(b). It is interesting to note that the structure is similar to the arch computed by Rozvany and Prager [18]. 5.2. Internally pressurized structure design In this example, we will design a cover-like structure subject to the internal pressure acting on the top. The problem setup is shown in Fig. 9(a) with the "xed #uid region in the rectangular shape on the top. The structure is clamped in both endpoints at the top corners. The "nite element analysis model consists of 2400 quad elements. The solid region is to use no more than 30% of the total
Finel=812=Ravi=VVC
68
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
Fig. 9. The internally pressurized container design. (a): The design domain and the initial pressure distribution. (b): The optimized topology.
Fig. 10. The internally pressurized container design with enlarged "xed #uid region. (a): The design domain and the initial pressure distribution. (b): The optimized topology.
volume. Again, the thermal stress tensor of the solid region is chosen to be isotropic to simulate the hydrostatic pressure. The optimized topology is shown in Fig. 9(b). Unlike the cover-like structure obtained in the Hammer and Olho! work [12], the optimized structure bears a distant resemblance to the MBB beam [19] and has a relative small #uid-phase area. The structure is basically composed of two overlapping arch trusses. The larger inverted arch truss spans the two supports and the smaller one is located adjacent to the "xed #uid region to resist the pressure force acting on the top. These two arch trusses are further reinforced by two additional supporting trusses connecting them. To form a cover-like structure, we enlarge the "xed #uid region to form an inverted triangle shape as shown in Fig. 10(a). The optimized topology is a cover-like structure as shown in Fig. 10(b). By extending the "xed #uid-phase further into the design domain, we essentially rule out this occupied area as a candidate region for the structure to reside. It should be noted that the mean compliance of the double-arch structure of Fig. 9(b) is 20.26% smaller than that of the cover-like structure of Fig. 10(b). In other words, the double-arch structure is sti!er when subjected to the hydrostatic pressure.
Finel=812=Ravi=VVC
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
69
Fig. 11. The #uid chamber design. (a): The design domain and the "xed #uid and "xed solid regions. (b): The optimal chamber pro"le with inlet from a channel to a bigger pressure chamber.
5.3. Flow chamber design The last example models the internally pressurized structure with inlet from a channel to a larger chamber as shown in Fig. 11(a). The inlet is designated as the solid region and is not subject to change in the optimization process as shown by the black color, while the gray color denotes the "xed #uid region. The "nite element analysis model consists of 2280 quad elements and the solid region occupies no more than 25% of the total volume. The thermal stress tensor of the solid region is isotropic. The optimized topology is shown in Fig. 11(b). The optimized structure is similar to the result in the Hammer and Olho! work [12]. 6. Conclusion A generalization of the classical topology optimization problem of designing structures subject to design-dependent loads has been presented. A "ctitious thermal force is used to simulate this design-dependent load such as the hydrostatic pressure force. The formulation stays within the context of the classical topology optimization paradigm and uses the "xed "nite element mesh in the optimization process. The major advantage is to circumvent the computational complexity associated with the domain shape variation if the shape of the analysis domain changes. A modi"ed SIMP material model has been proposed to con"ne the pressure force acting on a small range of density contour curves without explicitly constructing the iso-volumetric density curve. Both internally and externally pressurized structures have been designed. It has also been demonstrated that by manipulating the components of the thermal stress tensor, one can simulate a wide variety of interface force loading conditions in addition to the hydrostatic pressure. This work serves as the starting point to design structures involving solid and #uid interaction in which the magnitude of interaction force is a design-dependent quantity itself. Acknowledgements The "nancial support of DARPA for this research under Grant No. DARPA/NIST F000765 is gratefully acknowledged.
Finel=812=Ravi=VVC
70
B.-C. Chen, N. Kikuchi / Finite Elements in Analysis and Design 37 (2001) 57}70
References [1] R.V. Kohn, G. Strang, Optimal design in elasticity and plasticity, Numer. Methods Eng. 22 (1986) 183}188. [2] R.V. Kohn, G. Strang, Optimal design and relaxation of variational problems, Comm. Pure Appl. Math (New York). 39 (1986) 1}25 (Part I), 139}182 (Part II), 353}357 (Part III). [3] M.P. Bends+e, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Eng. 71 (2) (1988) 197}224. [4] K. Suzuki, N. Kikuchi, A homogenization method for shape and topology optimization, Comput. Methods Appl. Mech. Eng. 93 (1991) 291}318. [5] A. Diaz, M.P. Bends+e, Shape optimization of structures for multiple loading situations using a homogenization method, Struct. Optim. 4 (1992) 17}22. [6] M.P. Bends+e, J.M. Guedes, R.B. Haber, P. Pedersen, J.E. Taylor, An analytical model to predict optimal material properties in the context of optimal structural design, Trans. ASME J. Appl. Mech. 61 (4) (1994) 930}937. [7] M.P. Bends+e, A.R. DmH az, R. Lipton, J.E. Taylor, Optimal design of material properties and material distribution for multiple loading conditions, Int. J. Numer. Methods Eng. 38 (7) (1995) 1149}1170. [8] A. Diaz, N. Kikuchi, Solution to shape and topology eigenvalue optimization problems using a homogenization method, Int. J. Numer. Methods Eng. 35 (1992) 1487}1502. [9] Z.-D. Ma, N. Kikuchi, I.H Hagiwara, Structural topology and shape optimization for a frequency response problem, Comput. Mech. 13 (3) (1993) 157}174. [10] M.P. Bends+e, Optimization of Structural Topology, Shape, and Material, Springer, London, UK, 1995. [11] M.P. Bends+e, J. Soko"owski, Shape sensitivity analysis of optimal compliance functionals, Mech. Struct. Mach. 23 (1) (1995) 35}58. [12] V.B. Hammer, N. Olho!. Topology optimization with design dependent loads, WCSMO-3 Bu!alo, 1999, pp. 629}631. [13] M.P. Bends+e, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. 69 (1999) 635}654. [14] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behavior of multiphase materials, J. Mech. Phys. Solids 11 (1963) 127}140. [15] K. Dems, Z. MroH z, Variational approach by means of adjoint systems to structural optimization and sensitivity analysis * I, Int. J. Solids Struct. 19 (8) (1983) 677}692. [16] K.K. Choi, E.J. Haug, V. Komkov, Design Sensitivity Analysis of Structural Systems, Academic Press, London, 1986. [17] J. Petersson, O. Sigmund, Slope constrained topology optimization, Int. J. Numer. Methods Eng. 41 (8) (1998) 1417}1434. [18] G.I.N. Rozvany, W. Prager, New class of structural optimization problems: optimal archgrids, Comput. Methods Appl. Mech. Eng. 19 (1) (1979) 127}150. [19] N. Olho!, M.P. Bends+e, J. Rasmussen, On CAD-integrated structural topology and design optimization, Comput. Methods Appl. Mech. Eng. 89 (1992) 259}279.