Optimal structural design considering flexibility

Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

www.elsevier.com/locate/cma

Optimal structural design considering ¯exibility Shinji Nishiwaki *, Seungjae Min 1, Jeonghoon Yoo 2, Noboru Kikuchi Department of Mechanical Engineering and Applied Mechanics, The University of Michigan, Ann Arbor, MI 48109-2125, USA Received 5 May 1998; received in revised form 26 April 2000

Abstract A topology optimization method based on the homogenization method (the homogenization design method) has been successfully applied to a variety of optimization problems such as the sti€ness maximization problem and the eigen-frequency maximization problem. In this study, a methodology to obtain the optimal structure design considering ¯exibility is developed as a new extension of the homogenization design method. First, ¯exibility is formulated using the mutual energy concept. Second, a new multi-objective function is proposed to obtain optimal solutions incorporating ¯exibility and sti€ness. Next, the topology optimization procedure is constructed using the homogenization method and sequential linear programming (SLP). Finally, some examples are presented to con®rm that the methodology presented here can provide ¯exible structures satisfying the problem speci®cations. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction We shall discuss the design optimization method which determines a topology considering structural ¯exibility. In general structural design, the sti€est structure has been considered optimal. The objective function for the static design problem is usually selected as the maximization of the sti€ness or the minimization of the mean compliance at the portion where we apply forces. The reason for selecting these physical values is that the sti€est structure can provide high stability to maintain the shape of the structure. However, the structure can provide higher performance or additional function if the ¯exibility is implemented in the appropriate portions in the structures. This is true even though the structural design requires sti€ness in some portions. A typical example which can obtain higher performance using ¯exibility is an automotive body structure. In general, four design criteria must be considered in automotive body design: drivability, durability, crashworthiness, and comfort (see details in [1]). Fukushima et al. [2] and Yang and Chahande [3] applied the topology optimization method to the design of an automotive body structure by setting the objective function as the minimization of the mean compliance in the multi-loading case. This is because the ®rst three design criteria can be improved to certain values using this objective function. However, these criteria can be further improved by increasing the ¯exibility of the body structure if the direction of the deformation can be speci®ed. That is, drivability will improve if the deformation mode of the automotive structure posed by the ¯exibility can automatically adjust the alignment of the automobile. Durability will also improve

*

Corresponding author. Present address: Structural Dynamics Laboratory, Toyota Central R&D, Labs., Inc., Nagakute, Aichi-Gun, Aichi Prefecture 480-1192, Japan. Tel.: +81-561-63-4115; fax: +81-561-63-6119. E-mail address: [email protected] (S. Nishiwaki). 1 Present address: School of Mechanical Engineering, Hanyang University, Seoul, 133-791, Korea. 2 Present address: Dept. of Mechanical Engineering, Yonsei University, Seoul, 120-749, Korea. 0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 0 ) 0 0 3 2 9 - 7

4458

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

further by relaxing the crucial stress concentration due to the speci®ed ¯exibility, and crashworthiness will also improve by specifying the deformation mode which determines the route to eciently absorb crash energy in the automotive body. Thus, by implementing the appropriate ¯exibility which can specify the direction of deformation, we can obtain an automotive body with higher performance. On the other hand, we can obtain an additional function such as the mechanical function if ¯exibility which can specify the direction of deformation is implemented in the appropriate portions in the structures. A typical example of a structure with mechanical function is a compliant mechanism (strictly speaking, the distributed compliant mechanism explained later). It is designed to be ¯exible in order to achieve a speci®ed motion as a mechanism. That is, the ¯exibility implemented in the structure provides the intentional direction of the deformation for the motion as a mechanism. As noted by Howell and Midha [4], ``The advantages of compliant mechanisms include simpli®cation of manufacturing and assembly, as well as reduction in cost, weight wear, backlash, noise, and need for lubrication''. Fully compliant mechanisms, which are composed only of ¯exible structures, are classi®ed into two categories, as Ananthasuresh et al. [5] described: lumped compliant mechanisms and distributed compliant mechanisms. In lumped compliant mechanisms, the ¯exibility is provided in localized areas, that is, at ¯exural pivots. On the other hand, in distributed compliant mechanisms, the ¯exibility is distributed through the mechanism. The idea of incorporating ¯exibility into a mechanism is not new. The earliest e€ort for the analysis and synthesis of compliant mechanisms was made by Burns and Crossley [6]. Later, several synthesis and analysis methods were developed by many researchers. Shoup and McLarnan [7] investigated the several types of compliant mechanisms having lower pairs. Sevak and McLarnan [8] analyzed compliant mechanisms using the nonlinear ®nite method, and used the optimization technique for synthesis. Midha and his associates combined and extended their preliminary ideas, and developed the kinematic synthesis approach. In this method, the basic design con®guration obtained by the knowledge from traditional rigid body kinematics is converted to a partially compliant mechanism or a lumped compliant mechanism. Her and Midha [9] presented a methodology for obtaining all possible compliant mechanisms from a given rigid body kinematic chain. Murphy et al. [10] extended this methodology to a topological synthesis approach using graph theory. Howell and Midha [4,11] developed an analysis and synthesis method for lumped compliant mechanisms using a pseudo-rigid-body model. This kinematic synthesis approach, however, is limited to the design of lumped compliant mechanims. Another major approach is the continuum synthesis approach which uses the topology optimization technique for mechanical structures. Sigmund [12] ®rst proposed a design approach based on the truss topology optimization method. Frecker et al. [13] also developed a design approach using the truss topology optimization method and a new multi-objective function. Sigmund [14] extended the design approach using continuum topology optimization based on the density approach and an objective function called the mechanical advantage. Larsen et al. [15] also developed a similar design approach using the continuum topology optimization and a di€erent objective function that minimizes the error in obtaining prescribed values of the geometrical and mechanical advantages. On the other hand, Ananthasuresh [16] applied the topology optimization technique based on the homogenization method to the continuum synthesis approach. However, his results seem to be rather the maximum sti€ness design than the mechanism design because of an inappropriate formulation of the multi-objective function. Several methods have been proposed for the topology optimization. Among them, the density approach and the homogenization design method are most commonly used. In the density approach, each value of elasticity tensor is assumed to be a function of the penalized material density with an exponent parameter (see details in [17]). Since this formulation is easy for the numerical implementation, this density approach has been used for applications such as the microstructure design [18] and the minimum compliance problem [19]. Recently, Swan and Kosaka [20,21] developed a new method based on hybrid combinations of classical Reuss (compliant) and Voigt (sti€) mixing rules. This method is attractive because this method does not need a penalization parameter which cannot be theoretically determined, and a function of the penalized material density is theoretically derived. On the other hand, the homogenization design method has been developed based on the mathematical theory. This method was ®rst proposed by Bendsùe and Kikuchi [22]. Prior to this, Cheng and Olho€ [23]

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4459

reported important nonsmooth characteristics about the optimal design of solid elastic plates. Lurie et al. [24,25] answered this problem with the G-convergence theory. Kohn and Strang [26] introduced the general relaxation concept to deal with an ill-posed variational problem in the optimal design. They discussed the relation between relaxation and homogenization in the optimal design. Murat and Tartar [27] introduced the characteristic function to deal with the topology optimization problem, and they also pointed out the necessity of the homogenization method for the relaxation. Rozvany et al. [28] investigated the implication of this relaxation method by designing perforated elastic plates. Bendsùe and Kikuchi [22] constructed a methodology for topology optimization using the concept of the ®xed design domain and the special microstructure including an empty rectangle. This method was applied to many problems such as the minimum compliance problem [29], the eigen-frequency problem [30±32], the frequency response problem [31,32], and the buckling problem [33,34]. In this paper, a methodology to obtain the optimal structure design considering ¯exibility is developed using the homogenization design method. In Section 2, the ¯exibility and sti€ness are formulated using the mutual energy concept. Their sensitivities with respect to a design variable are also derived. In Section 3, the homogenization design method is brie¯y discussed. That is, the concept of the ®xed design domain is introduced, and the procedure of the homogenization method for relaxation of the design domain is explained. In Section 4, three types of design optimization problems are formulated: the unconstrained single ¯exibility case, the constrained single ¯exibility case, and the multi-¯exibility case. In Section 5, for each case, the multi-objective functions are constructed in order to obtain appropriate optimal solutions in the Pareto optima. In Section 6, the optimization algorithm is constructed using these multi-objective functions, the homogenization method, and sequential linear programming (SLP). In Section 7, some twodimensional and three-dimensional design examples in three types of design optimization problems are presented in order to investigate the con®gurations of the optimal solutions. These examples con®rm that the methodology presented here can provide ¯exible structures satisfying the problem speci®cations. 2. Formulation of ¯exibility and sti€ness The ¯exibility and sti€ness of a structure are formulated using the mutual energy concept [35,36]. Consider a linear elastic body occupying a three-dimensional domain, X, as shown in Fig. 1. Suppose that the elastic body is ®xed at boundary Cd . Now, consider the two cases, Case (a) and Case (b). In each case, the elastic body is separately subjected to a di€erent boundary traction at a di€erent boundary. That is, in Case (a), the elastic body is subjected to boundary traction t1 at boundary Ct1 , and in Cast (b), it is subjected to boundary traction t2 at boundary Ct2 . Body forces applied to the elastic body are ignored for simplicity in the formulation. The displacement ®elds are u1 ˆ fu11 ; u12 ; u13 gt in Case (a), and u2 ˆ fu21 ; u22 ; u23 gt in Case (b). Here, we introduce the linear form (the mean compliance) de®ned by Z Z Lk …ul † ˆ tk
ul dC ˆ tik uli dC; ul 2 V ; …2:1† Ctk

Ctk

Fig. 1. An elastic body subjected to two tractions.

4460

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

where k and l are either 1 or 2, and V is the admissible linear space such that V ˆ fv ˆ vi ei : vi 2 H 1 …X† with v ˆ 0 on Cd ; i ˆ 1; . . . ; 3g: The bilinear form is also introduced as Z Z t a…u; v† ˆ e…v† Ee…u† dX ˆ Eijkl ekl …u†eij …v† dX; X

X

with linearized strains,   1 oui ouj eij …u† ˆ ‡ ; 2 oxj oxi

…2:2†

…2:3†

…2:4†

where Eijkl is the elasticity tensor of the linear elastic body. The displacements satisfy, respectively, the following principles of virtual work: a…u1 ; v1 † ˆ L1 …v1 † 2

2

2

2

a…u ; v † ˆ L …v †

for u1 2 V

8 1

…2:5†

2

8 2

…2:6†

for u 2 V

v 2V; v 2V;

Furthermore, in the sense of the mutual energy concept, the displacements also satisfy, respectively, the following principles of virtual work: a…u1 ; v2 † ˆ L1 …v2 † 2

1

2

1

a…u ; v † ˆ L …v †

for u1 2 V

8 2

…2:7†

2

8 1

…2:8†

for u 2 V

v 2V; v 2V:

Therefore, these four equations of the principle of virtual work are generalized in the following equation: a…uk ; v1 † ˆ Lk …v1 †

for uk 2 V

8 l

v 2V;

…2:9†

where k and l are either 1 or 2. We assume that the boundary traction t1 in Case (a) is an applied force, and the boundary traction t2 in Case (b) is a unit dummy load. Then, the mutual mean compliance, L2 …u1 †, described as Z 2 1 t2
u1 dC ˆ a…u2 ; u1 †; …2:10† L …u † ˆ Ct2

can be interpreted as the measure of deformation at boundary Ct2 when we apply boundary traction t1 at boundary Ct1 . That is, the mutual mean compliance, L2 …u1 †, is interpreted as how ¯exible or sti€ boundary Ct2 is when boundary traction t1 is applied at boundary Ct1 . Therefore, if L2 …u1 † is increased or maximized in an optimization problem, and is suciently large, then sucient ¯exibility in the direction of t2 with respect to t1 can be obtained. Note that the mutual mean compliance, L2 …u1 †, must be at least positive in order to obtain the deformation in the speci®ed direction, although a larger mutual mean compliance provides a higher ¯exibility performance. Conversely, if the absolute value of L2 …u1 † is decreased or minimized in an optimization problem and suciently small, then sucient mutual sti€ness in the direction of t2 with respect to t1 can be obtained. Note that this mutual sti€ness behaves like a constraint between two boundaries rather than like the ordinary sti€ness because this sti€ness is not obtained at Ct2 if boundary traction t1 is not applied at boundary Ct1 even though a boundary traction is applied at Ct2 . If only Case (a) is considered, the ordinary mean compliance, L1 …u1 †, described as Z t1
u1 dC ˆ a…u1 ; u1 †; …2:11† L1 …u1 † ˆ Ct1

can be interpreted as the measure of deformation at boundary Ct1 when we apply the boundary traction t1 at boundary Ct1 . That is, this mean compliance is interpreted as the measure of the stifness at boundary Ct1 Therefore, if the mean compliance, L1 …u1 †, is decreased or minimized in an optimization problem, and is suciently small, sucient sti€ness can be obtained at boundary Ct1 .

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4461

Next, the sensitivities of the mutual mean compliance, L2 …u1 †, and the mean compliance, L1 …u1 †, with respect to a design variable A are formulated. These derivatives are used in the optimization procedure which is composed of sequential linear programming (SLP). The sensitivity of Lk …u1 † is derived, where k and l are either 1 or 2. Consider the mutual total potential energy de®ned by 1 F …vk ; vl † ˆ a…vk ; vl † 2

1 k l L …v † 2

1 1 k L …v †: 2

…2:12†

In the case where k is equal to l, minimizing mutual total potential energy with respect to vk and vl corresponds to the principle of virtual work de®ned by Eqs. (2.5) and (2.6). In the case where k is not equal to l, minimizing mutual total potential energy with respect to vk and vl corresponds to the principle of virtual work de®ned by Eqs. (2.7) and (2.8). Since the following relation is obtained at equilibrium using Eq. (2.9): Lk …u1 † ˆ a…uk ; ul † ˆ a…ul ; uk † ˆ Ll …uk †;

…2:13†

the mutual total potential energy satis®es the following equation at equilibrium: 1 F …uk ; ul † ˆ a…uk ; ul † 2

1 k l L …u † 2

1 l k L …u † ˆ 2

1 k l L …u † ˆ 2

1 l k L …u †: 2

…2:14†

Since the displacement ®elds uk and ul are a function of design variable A, F is a function of uk ; ul and A at equilibrium. We assume that the shape of the domain, X, does not change even though the value of the design variable A changes. That is, we assume that the domain X, is not expressed by a function of the design variable A. Taking the ®rst variation of F with respect to uk ; ul and A at equilibrium, dF …uk ; ul ; A† ˆ DF …uk †…duk † ‡ DF …ul †…dul † ‡ DF …A†…dA†;

…2:15†

where DF …u†…v† is the G-di€erential such that o F …u ‡ nv† n ! 0 on

DF …u†…v† ˆ lim yields

      l k 1 1 k oul 1 k ou l l l ou k L dF …u ; u ; A† ˆ a u ; dA ‡ du dA ‡ du ‡ a u ; dA ‡ du 2 2 2 oA oA oA   Z 1 l ouk 1 t oE L e…uk † dX dA: dA ‡ duk ‡ e…ul † 2 2 X oA oA k

l

Since vl is arbitrary in Eq. (2.9), by setting vl ˆ …oul =oA†dA ‡ dul , we have   l   l ou k ou l k l dA ‡ du dA ‡ du ˆ 0: a u; L oA oA

…2:16†

…2:17†

Similarly, by switching the superscripts k and l, and by setting vk ˆ …ouk =oA†dA ‡ duk , we also have   k   k ou l ou k l k dA ‡ du dA ‡ du ˆ 0: a u; L …2:18† oA oA Substituting Eqs. (2.17) and (2.18) into Eq. (2.16), ®nally, we have Z 1 oE k k l e…u † dX dA: dF …u ; u ; A† ˆ e…ul †t 2 X oA

…2:19†

4462

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Using Eq. (2.14), the sensitivity of the mutual mean compliance, Lk …ul †, is obtained by Z oLk …ul † t oE ˆ e…uk † dX: e…ul † oA oA X

…2:20†

By setting k ˆ 2 and l ˆ 1, the sensitivity of the mutual mean compliance, L2 …u1 †, is obtained as follows: Z oL2 …u1 † oE 2 ˆ e…u † dX: e…u1 †t …2:21† oA oA X Similarly, by setting k ˆ l ˆ 1, the sensitivity of the ordinary mean compliance, L1 …u1 †, is obtained as follows: Z oL1 …u1 † t oE ˆ e…u1 † dX: e…u1 † …2:22† oA oA X The formulation of the sensitivity in Eq. (2.22) is identical to the formulation which Suzuki and Kikuchi [29] introduced when constructing the minimal mean compliance design problem. 3. Homogenization design method 3.1. Concept of the ®xed design domain The concept of the ®xed design domain [22,37] is brie¯y discussed in this section. Consider the design problem determining the boundary of the design domain Xd by minimizing or maximizing objective functions. The key idea of this method is the introduction of a ®xed and extended design domain D including the orginal design domain Xd , a priori, and the utilization of the following characteristic function introduced by Murat and Tartar [27] to describe any of the boundaries of the original design domain Xd ,  1 if x 2 Xd ; vX …x† ˆ …3:1† 0 if x 2 D nXd ; and the equilibrium equation is rewritten as follows: Z Z Z t e…v† vX Ee…u† dD ˆ vX b
v dD ‡ t
v dC for u 2 V D

D

Ct

8

v2V;

…3:2†

where E is the elasticity tensor, u the displacement ®eld, b the body force, and t is the boundary traction applied at boundary Ct . It is clear that by extending the domain Xd into D, we construct the equilibrium equation of a structure distributing a new elasticity tensor, vX E, in the extended design domain D. Note that both the original design domain Xd and the extended design domain D have the same boundary conditions in the same portions. Since the characteristic function can have only the discretized value, either 0 or 1 in an in®nitely small area, this function is very discontinuous everywhere. That is, elasticity tensor, vX E, is at most de®ned in L1 …D†. In other words, the optimal con®guration may have solid material distribution in an in®nitely small area while it may have voids around this area such as in porous composite media. This property makes numerical treatment of the structural optimization problem dicult. To overcome this problem, the homogenization method is applied. In this method, the structural optimization problem is interpreted as a problem of distributing the porous media or the composite material which has in®nite microscale voids as shown in Fig. 2. To represent these porous media, a microstructure composed of solid material and void is used. The scheme, which replaces the wildly discontinuous tensor, vX E, with the smooth function in the global sense, is equivalently treated as the homogenization procedure for the design domain. That is, the relaxed physical properties such as the elasticity tensor are obtained as the homogenized properties. Finally, the extended design domain D is recon®gured using the homogenized properties which have numerical smoothness.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4463

Fig. 2. The homogenization design method.

3.2. Homogenization method Consider the microstructures shown in Fig. 3. In the two-dimensional case, the microstructure is formed inside an empty rectangle in a unit cell, where a; b, and h are regarded as the design variables. In order to develop a complete void, both a and b must be 1, whereas for solid material a and b must be 0. The variable h represents the rotation of the unit cell. In the three-dimensional case, the microstructure is formed inside an empty rectangular box in a unit cell, where a; b; c; /; w, and h are regarded as the design variables. In order to develop a complete void, a; b and c must all be 1, whereas for solid material a; b, and c must be 0. The variables /, w, and h represent the three-dimensional rotations of the unit cell. Note that the optimal solutions obtained by the rectangle type microstructure are sub-optimal, and the convergence with respect to grid re®nement is still not proved. We may use the rank-n type microstructure (e.g. see [17,38]) to obtain the true optimal solutions and to guarantee the convergency. However, the rankn type microstructure provides many gray-scale porous regions implying a composite status in the optimal con®gurations. Several methods such as the perimeter control method [39] are proposed to eliminate grayscale regions from the optimal con®gurations by paying penalty to describe 0±1 sub-optimal solutions. However, these methods do not work well for all design problems [40]. On the other hand, the rectangle type microstructure provides a clear con®guration which does not have many gray regions in the engineering sense. This is because the rectangle type microstructure has a penalty characteristic between the elasticity tensor and the density of the microstructure which can adjust an isotropic response during the optimization process. Furthermore, we also numerically con®rmed that the optimal solutions are largely independent of the grid re®nement as shown in the ®gures of examples described in this paper.

Fig. 3. Microstructures for the relaxation of the design domain. (a) Two-dimensional case, (b) three-dimensional case.

4464

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Here, the procedure of the homogenization method is brie¯y explained. The elasticity tensor, vX E, is described as E e . Suppose that the elasticity tensor, E e , has the Y-periodic characteristic in order to occupy the extended design domain D with the in®nite microstructures shown in Fig. 3. Using parameter e which represents the periodicity and is assumed to be very small, the local coordinate y in the microstructure is de®ned by y ˆ x=e:

…3:3†

Using this local coordinate, the elasticity tensor, E e , is described as E e …x† ˆ E…x; y†:

…3:4†

Because of the Y-periodic unit cell of the microstructure, the elasticity tensor, E e , has the following periodic characteristic: E…x; y† ˆ E…x; y ‡ Y†:

…3:5†

First, we calculate the homogenized elasticity tensor, E H , in the case where the angle h is set to 0 in the two-dimensional microstructure, or where the angles /; w, and h are set to 0 in the three-dimensional microstructure. To obtain this homogenized elasticity tensor, the characteristic deformations, v…x; y†, are calculated using the following equation: Z Z t t ey …v† E…x; y†ey …v…x; y†† dY ˆ ey …v† E…x; y† dY for 8 v 2 Vy ; …3:6† Y

Y

where t



ey …v† ˆ

ov1 ov2 ov3 1 oy1 oy2 oy3 2



ov2 ov3 ‡ oy3 oy2

     1 ov3 ov1 1 ov1 ov2 ‡ ‡ ; 2 oy1 oy3 2 oy2 oy1

and Vy is the admissible space de®ned in the unit cell Y such that Vy ˆ fv ˆ vi ei : vi 2 H 1 …Y † j v is Y-periodic in the unit cell Yg:

…3:7†

After obtaining the characteristic deformation, v…x; y†, the homogenized elasticity tensor, E H , is computed by EH ˆ

1 jYj

Z Y

E…x; y†…I

ey …v†† dY;

…3:8†

where jYj stands for the area of the unit cell. Next, when the unit cell is rotated by angel h in the two-dimensional case, the homogenized elasticity tensor, E G , is computed by E G ˆ R…h†t E H R…h†:

…3:9†

When the unit cell is rotated by angles /; w, and h in the three-dimensional case, the homogenized elasticity tensor, E G , is computed by t

t

t

E G ˆ R…h† R…w† R…/† E H R…/†R…w†R…h†;

…3:10†

where R is the rotation matrix. Thus, the homogenized elasticity tensor, E G , is determined by the microscopic design variables a; b; and h in the two-dimensional case, and a; b; c; /; w, and h in the three-dimenisonal case.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4465

4. Formulation of optimization problems In this section, the three types of the design optimization problems of ¯exible structures are formulated: the unconstrained single ¯exibility case, the constrained single ¯exibility case, and the multi-¯exibility case. In the unconstrained single ¯exibility case, a ¯exible structure is designed to have only a single ¯exibility without any relative deformation constraints. In the constrained single ¯exibility case, a ¯exible structure is designed to have relative deformation constraints at the portions where the ¯exible structure is supposed to deform, or where we apply the traction. In the multi- ¯exibility case, a ¯exible structure is designed to have more than a single ¯exibility. 4.1. Unconstrained single ¯exibility case First, the design of ¯exible structures in the unconstrained single ¯exibility case is discussed. Suppose that the original design domain of a ¯exible structure Xd is ®xed at boundary Cd , and is subjected to the applied traction t1 at boundary Ct1 as shown in Fig. 4(a). We also consider the extended design domain D including Xd . Now, we intend to design a ¯exible structure which starts to deform in the speci®ed direction t2 at boundary Ct2 in order to work as a compliant mechanism or obtain higher performance as a structure when we apply traction t1 at boundary Ct1 . To implement this function in the ¯exible structure, one must take into account the kinematic (motion) requirement and the structural requirement simultaneously. For the kinematic requirement, the ¯exible structure must have sucient ¯exibility, which provides sucient deformation along a direction speci®ed by a dummy load t2 at Ct2 when traction t1 is applied at Ct1 as shown in Fig. 4(a). It is obtained by maximizing the mutual mean compliance posed by traction t2 at boundary Ct2 , and the displacement ®eld u1 due to traction t1 . Indeed, this kinematic requirement is identical to our problem speci®cations mentioned above. Note that if we use the ¯exible structure as a compliant mechanism, the ¯exible structure may deform considerably. However, we can ignore the e€ect of the large deformation in the topology design phase. That is, we can design compliant mechanisms based on the small linear deformation assumption because our goal in this design phase is that the compliant mechanism qualitatively deforms in the desired direction of motion. This assumption is appropriate only if we take into account the qualitative characteristics of the mechanism function. Conversely, if we consider the quantitative performance of the deformation in the mechanism, we must include large deformation analysis in compliant mechanism design. Since our goal is to ®nd the optimum design that qualitatively provides the initial motion of the mechanism in the speci®ed direction, the topology of the compliant mechanism can be designed based on the small deformation assumption. However, we also note that the design con®guration must be subsequently examined to determine whether this con®guration generates the speci®ed large deformation which we originally intended in the topology design phase. In the following section, the design con®guration will be examined using the large displacement ®nite element analysis.

Fig. 4. Speci®cations for a ¯exible structure design.

4466

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Next, we discuss the structural requirement. If we only take into account the kinematic function for the design of the ¯exible structure, ideally, no structure in the design domain may be the optimal con®guration. That is, even though our main goal is implementation of sucient ¯exibility in the structure, the ¯exible structure must have sucient sti€ness in the appropriate portions in order to form structure. We de®ne this sucient sti€ness as the structural requirement. For this structural requirement, we impose two di€erent sti€nesses at two di€erent boundaries. One is the sucient sti€ness at boundary Ct1 as shown in Fig. 4(b) when we apply traction t1 at boundary Ct1 . This sucient sti€ness works to maintain the shape of the ¯exible structure when the ¯exible sturcture deforms due to traction t1 . For example, in a compliant mechanism design, this sti€ness is required to maintain the structural con®guration in a part where a force from an actuator or a hand is applied, while a ¯exible structure is working as a mechanism due to the kinematic requirement. It is obtained by minimizing the mean compliance at boundary Ct1 posed by traction t1 while boundary Ct2 is ®xed since the ¯exible structure is supposed to be imposed by the reaction force of the object. The other sucient sti€ness is one at boundary Ct2 as shown in Fig. 4(c). This sucient sti€ness works to maintain the shape of the ¯exible structure against the reaction force by the object. In a compliant mechanism design, this sti€ness works after a ¯exible structure deforms to contact a workpiece. That is, this sti€ness is required to keep operating on a workpiece without a local deformation at a part where the reaction force by the workpiece is imposed, while we are applying a force. It is obtained by minimizing the mean compliance at boundary Ct2 posed by traction t2 while boundary Ct1 is ®xed because the direction of the reaction force is opposite to that of the dummy load t2 and ¯exible structure must be imposed by the applied traction t1 . To be satis®ed with the kinematic requirement as shown in Fig. 4(a), and two structural requirements in Fig. 4(b) and (c), we formulate the following multi-optimization problem: Z 2 1 maximize L …u † ˆ t2
u1 dC a;b;c;/;w; and h Ct2 Z t3
u3 dC and …4:1† minimize L3 …u3 † ˆ a;b;c;/;w; and h 1 Ct Z t4
u4 dC minimize L4 …u4 † ˆ a;b;c;/;w; and h

Ct2

subject to t3 ˆ t1 4

t ˆ

…4:2† 2

t;

…4:3†

1

1

1

1

1

for u 2 V

…a† 8 1

…a†

;

…4:4†

2

2

2

2

2

for u 2 V

…a† 8 2

…a†

;

…4:5†

3

3

3

3

3

for u 2 V

…b† 8 3

…b†

;

…4:6†

4

4

4

4

4

…c† 8 1

…c†

;

…4:7†

a…u ; v † ˆ L …v † a…u ; v † ˆ L …v † a…u ; v † ˆ L …v † a…u ; v † ˆ L …v †

for u 2 V

v 2V v 2V v 2V

v 2V

0 6 a 6 1; 0 6 b 6 1;

…4:8† …4:9†

06c61

…4:10†

…three-dimensional case†; Z Z g…a; b; c† ˆ q dX Xs ˆ q dX Xd

D

Xs 6 0;

…4:11†

where q ˆ 1 ab in the two-dimensional case, and q ˆ 1 abc in the three-dimensional case. Xs is the total volume constraint of the solid material forming the porous structure, and V …a† ˆ fv ˆ vi ei : vi 2 H 1 …D† with v ˆ 0 on Cd g; V

…b†

V

…c†

…4:12†

1

…4:13†

1

…4:14†

ˆ fv ˆ vi ei : vi 2 H …D† with v ˆ 0 on Cd and Ct2 g; ˆ fv ˆ vi ei : vi 2 H …D† with v ˆ 0 on Cd and Ct1 g:

Note that boundaries Ct3 and Ct4 are set identical to boundaries Ct1 and Ct2 , respectively.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4467

4.2. Constrained single ¯exibility case Next, we shall discuss the case where the ¯exible structure is required to have a relative deformation constraint in one of two portions: the portion where we need to have the deformation for the kinematic requirement, corresponding to boundary Ct2 in Fig. 5 or the portion where we apply the traction, corresponding to boundary Ct1 in Fig. 6. Since we only specify the initial direction of the translational deformation of the ¯exible structure, if only the mutual mean compliance, L2 …u1 †, is maximized to obtain sucient ¯exibility for the kinematic requirement, then the ¯exible structure may provide an additional translational deformation which may cause the rotational motion to maximize the ¯exibility. That is, even though the initial direction of deformation can be speci®ed, the next direction of deformation after the ¯exible structure deforms may not be determined, a priori. Therefore, if we intend to specify only the desired translational deformation, which is sometimes required for operations such as the compliant mechanisms, the relative energy constraint based on the mutual energy, the mutual sti€ness, must be imposed on the ¯exible structure. Suppose that the ¯exible structure must deform along a direction speci®ed by dummy load t2 at boundary Ct2 continuously by the translational motion where traction t1 is applied at boundary Ct1 . Here, the original design domain of a ¯exible structure Xd is ®xed at boundary Cd . Let us consider the threedimensional case shown in Fig. 5(b). The problem formulation in the two-dimensional case shown in Fig. 5(a) is obtained by simplifying the formulation in the three-dimensional case. We introduce the two

Fig. 5. Constrained single ¯exibility case (constraint is added at the portion where we need to have the deformation for the kinematic requirement). (a) Two-dimensional case, (b) three-dimensional case.

Fig. 6. Constrained single ¯exibility case (constraint is added at the portion where we apply traction). (a) Two-dimensional case, (b) three-dimensional case.

4468

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

dummy loads, t5 and t6 , whose directions are perpendicular to the direction of t2 in order to specify the deformed motion as shown in Fig. 5(b). Using these dummy loads, The translational deformation in the direction of t2 can be obtained by solving the problem described equivalently as follows. The ¯exible structure must deform in the direction of dummy load t2 , and must be ®xed in the direction of dummy loads t5 and t6 as a kinematic function at boundary Ct2 where traction t1 is applied at boundary Ct1 . These constraints in the directions of t5 and t6 can be achieved by minimizing the absolute values of the additional two mutual mean compliances. That is, the following optimization problem is considered in addition to the optimization problems formulated in Eqs. (4.1)±(4.11): Z Z minimize jL5 …u1 †j ˆ t5
u1 dC ˆ t5
u1 dC ˆ a…u5 ; u1 † and a;b;c;/;w; and h 5 Ct C2 …4:15† Z Z t 6 1 6 1 6 1 6 1 t
u dC ˆ t
u dC ˆ a…u ; u † minimize jL …u †j ˆ a;b;c;/;w; and h 6 2 Ct

Ct

subject to …t5
t2 † ˆ 0; 6

2

5

6

…4:16†

…t
t † ˆ 0; …t
t † ˆ 0; a…u5 ; v5 † ˆ L5 …v5 † ˆ a…u6 ; v6 † ˆ L6 …v6 † ˆ

…4:17† Z Ct5

Z

Ct6

t5
v5 dC ˆ t6
v6 dC ˆ

Z Ct2

Z

Ct2

…4:18† t5
v5 dC

for u5 2 V …a† 8 v5 2 V …a† ;

…4:19†

t6
v6 dC

for u6 2 V …a† 8 v6 2 V …a† ;

…4:20†

where the operator (
) denotes the dot product of two vectors, and boundaries Ct5 and Ct6 are set identical to boundary Ct2 . Note that the ordinary ®xed condition applied in FEM must not be utilized since that ¯exible structure itself must deform only with the speci®ed motion even though we do not apply any ®xed conditions of boundary Ct2 . If we use the ordinary ®xed condition in the design phase, we must implement the constraint in practice to obtain the e€ect of the ®xed condition speci®ed in the design phase. Also note that the problem formulation in the two-dimensional case is obtained by setting t6 to 0. Similarly, if we intend to specify the deformation only in the desired translational direction at the portion where the external force is applied, corresponding to boundary Ct1 in Fig. 6, the relative energy constraint, the mutual sti€ness, must also be imposed on the ¯exible structure at boundary Ct1 . A typical example which requires the constraint at boundary Ct1 is the design of the automotive suspension system. In this design, one of the critical design criteria is controlling the motion or the deformation at the portion where the external forces are applied, corresponding to boundary Ct1 . Note that we took into account the sti€ness at boundary Ct1 as a structural requirement. However, since this sti€ness is for forming or maintaining the ¯exible structure where boundary Ct2 is ®xed, we do not specify the direction of the deformation due to this sti€ness at boundary Ct1 . Suppose that the ¯exible structure must deform in the direction of traction t1 at boundary Ct1 continuously by the translational motion. Here, the original design domain of a ¯exible structure Xd is ®xed at boundary Cd . Let us consider the three-dimensional case shown in Fig. 6(b). The problem formulation in the two-dimensional case shown in Fig. 6(a) is obtained by simplifying the formulation in the threedimensional case. We also introduce the two dummy loads, t7 and t8 , whose directions are perpendicular to the direction of t1 as shown in Fig. 6(b). The translational deformation in the direction of t1 can be obtained by solving the problem described equivalently as follows. The ¯exible structure must deform in the direction of dummy load t2 at boundary Ct2 , and must be ®xed in the direction of dummy loads t7 and t8 at boundary Ct1 , where traction t1 is applied at boundary Ct1 . These constraints in directions t7 and t8 can be also achieved by minimizing the absolute values of the additional two mutual mean compliances. That is, the following optimization problem is considered in addition to the optimization problems formulated in Eqs. (4.1)±(4.11):

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Z Z 7 1 1 1 minimize jL …u †j ˆ t
u dC ˆ t
u dC ˆ a…u7 ; u1 † and C7 a;b;c;/;w; and h Ct7 Z Z t t8
u1 dC ˆ t8
u1 dC ˆ ja…u8 ; u1 †j minimize jL8 …u1 †j ˆ a;b;c;/;w; and h 8 1 7

4469

1

Ct

…4:21†

Ct

subject to …t7
t1 † ˆ 0; 8

1

7

8

…4:22†

…t
t † ˆ 0; …t
t † ˆ 0; a…u7 ; v7 † ˆ L7 …v7 † ˆ a…u8 ; v8 † ˆ L8 …v8 † ˆ

…4:23† Z Ct7

Z

Ct8

t7
v7 dC ˆ t8
v8 dC ˆ

Z Ct1

Z

Ct1

…4:24† t7
v7 dC

for u7 2 V …a† 8 v7 2 V …a† ;

…4:25†

t8
v8 dC

for u8 2 V …a† 8 v8 2 V …a† ;

…4:26†

where boundaries Ct7 and Ct8 are set identical to boundary Ct1 . Note that the problem formulation in the two-dimensional case is obtained by setting t8 to 0. 4.3. Multi-¯exibility case Some ¯exible structures require more than one ¯exibility to increase various performance criteria or to provide two or more additional functions. For example, in this case of automotive body structure design, in general, two or more performance criteria, such as drivability, durability, and crashworthiness, must be considered. These criteria can be increased if the automotive body structure has ¯exibilities in the appropriate directions at the appropriate portions. However, the ¯exibility mode for one of the performance criteria is usually di€erent from those of any other performance criteria. Therefore, we must take into account the multi-¯exibility case in order to increase all the performance criteria. Moreover, in the case of compliant mechanisms, if we implement two or more ¯exibilities, the compliant mechanisms can have more than one function as a kinematic function. We shall consider the case where the ¯exible structure is required to have n different ¯exibility modes. Suppose that the original design domain of a ¯exible structure Xd is ®xed at boundary Cd , and the ¯exible structure starts to deform in the speci®ed direction, i t2 , at boundary i Ct2 when we apply traction i t1 at boundary i Ct1 …i ˆ 1; . . . ; n† as shown in Fig. 7. We also consider the extended design domain D including Xd . For each ¯exibility case, we must take into account both the kinematic (motion) requirement (Fig. 7, Case(a)) and the structural requirement (Fig. 7, Case(b) and Case(c)), which are identical in the

Fig. 7. Multi-¯exibility case.

4470

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

unconstrained single ¯exibility case. The kinematic requirement is satis®ed by maximizing the mutual mean compliance posed by traction i t2 at boundary i Ct2 and displacement ®eld i u1 due to traction i t1 . The structural requirement is satis®ed by minimizing both the mean compliance at boundary i Ct1 posed by traction i t1 while boundary i Ct2 is ®xed and the mean compliance at boundary i Ct2 posed by traction i t2 while boundary i Ct1 is ®xed. Here, we introduce the extended linear form for the multi-¯exibility case de®ned by Z m k m l m k m l L … u†ˆ t
u dC; m ul 2 V : …4:27† mC tk

Following the same optimization formulation of each requirement as was done in the unconstrained single ¯exibility case, we formulate the multi-objective optimization problem as follows: Z i 2 i 1 maximize i L2 …i u1 † ˆ t
u dC; a;b;c;/;w; and h i Ct 2 Z i 3 t
u3 dC and …4:28† minimize i L3 …i u3 † ˆ a;b;c;/;w; and h i Ct 1 Z i 4 i 4 t
u dC for i ˆ 1; . . . ; n minimize i L4 …u4 † ˆ a;b;c;/;w; and h

i Ct 2

subject to i 3

t ˆ i t1

i 4

t ˆ

…4:29†

i 2

t

…4:30†

i 1 i 1

i 1 i 1

i 1

i

…a† 8i 1

i

…a†

;

…4:31†

i 2 i 2

i 2 i 2

i 2

i

…a† 8i 2

i

…a†

;

…4:32†

i 3 i 3

i 3 i 3

i 3

i

…b† 8i 3

i

…b†

;

…4:33†

i 4 i 4

i 4 i 4

i 4

i

…c† 8i 4

i

…c†

;

…4:34†

a… u ; v † ˆ L … v † for u 2 V a… u ; v † ˆ L … v † for u 2 V a… u ; v † ˆ L … v † for u 2 V a… u ; v † ˆ L … v † for u 2 V

v 2 V v 2 V v 2 V

v 2 V

for i ˆ 1; . . . ; n, and 0 6 a 6 1; 0 6 b 6 1;

…4:35† …4:36†

06c61

…4:37†

…three-dimensional case†; Z Z g…a; b; c† ˆ q dX Xs ˆ q dX Xd

D

Xs 6 0;

…4:38†

where i

V …a† ˆ fv ˆ vi ei : vi 2 H 1 …D† with v ˆ 0 on Cd g;

i

V …b† ˆ fv ˆ vi ei : vi 2 H 1 …D† with v ˆ 0 on Cd and i Ct2 g;

i

V …c† ˆ fv ˆ vi ei : vi 2 H 1 …D† with v ˆ 0 on Cd and i Ct1 g:

and

Note that boundaries i Ct3 and i Ct4 are set identical to boundaries i Ct1 and i Ct2 , respectively. As we discussed in Section 2, the necessary condition for the mutual mean compliance is that each mutual mean compliance must be at least positive in order to obtain deformation in the speci®ed direction, although a larger mutual mean complicance provides a higher ¯exibility performance. Therefore, if we have trade-o€ relations among the n maximization problems of the mutual mean compliances, i L2 …i u1 † for i ˆ 1; . . . ; n, all the mutual mean compliances must be positive in the optimization procedure. This operation can be accomplished by maximizing the smallest of the mutual mean compliances, or placing the largest weight on the smallest one in the multi-objective functions.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4471

5. Formulation of multi-objective functions 5.1. Pareto optimality of multi-objective functions In general, one optimal solution is obtained in the single objective optimization problem. However, in the multi-objective probelm, there exists a series of optimal solutions called a set of Pareto optima, or noninferior solutions. Several methods have been developed in order to obtain the Pareto optima. Typical methods include the weighting method, the e-constraint method, and goal programming (see details in [41,42]). Among these the weighting method has been most widely used because of the simplicity of its formulation. In this method, the multi-objective functions are de®ned by the summation of objective functions by the weighting coecients. As Koski [43] explained, by changing the coecients, we have Pareto optima if each objective function has convexity characteristics. In the case of ¯exible structure design, as discussed in Section 4, we have to take into account both the maximization of the mutual mean compliance and the minimization of the mean compliances. The problem in this optimization is that all Pareto optima do not have physical meaning. That is, since the mutual mean compliance has a possiblity of being in®nite, some optimal con®gurations on the Pareto optima provide the full void in the structure numerically when the mutual mean compliance is large. Even though the mutual mean compliance has a ®nite value, we may have optimal solutions on the Pareto optima which do not have continuity in the structure since higher ¯exibility may disconnect structures as shown in Fig. 8. Therefore, although we can obtain a series of optimal solutions along the Pareto optima using the weighting method with di€erent weighting coecients, some of them do not produce practical structures. Furthermore, our goal of design for the ¯exible structure is to ®nd an appropriate optimal solution which has physical meaning on the Pareto optima, not to search all possible Pareto optima. In order to ®nd the appropriate optimal solution, we constructed new multi-objective functions as described in the next subsection. In this subsection, we discuss the Pareto optimality of a solution obtained by a new multi-objective function. Suppose that we have objective functions f1 …x†; . . . ; fp …x† which must be simultaneously minit mized with inequality constraints g1 …x† 6 0; . . . ; gm …x† 6 0 for the design variables x ˆ ‰x1 ; . . . ; xn Š . This multi-objective optimization problem is described as follows [42,43]: 2

3 f1 …x† 6 . 7 7 minimize f …x† ˆ 6 4 .. 5 x fp …x† subject to

Fig. 8. Disconnected structure due to higher ¯exibility.

…5:1†

4472

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

2

3 g1 …x† 6 7 g…x† ˆ 4 ... 5 6 0;

…5:2†

gm …x†

We de®ne the Lagrangian by t

L1 ˆ lt f …x† ‡ kt g…x†

t

for l ˆ ‰l1 ; . . . ; lp Š and k ˆ ‰k1 ; . . . ; km Š ;

…5:3†

where l and k are the Lagrange multipliers. Using this Lagrangian, the necessary condition for the Pareto optimality is described as follows: ^ is a noninferior solution, there must exist l ^ and ^k such that If x oL1 t og…x† t of …x† ^ ^ ˆl ‡k ˆ 0 for i ˆ 1; . . . ; n; …5:4† oxi xˆ^x oxi xˆ^x oxi xˆ^x oL1 ˆ gj …^ x† 6 0 okj kˆ^k

for j ˆ 1; . . . ; m;

oL1 ^ kj ˆ k^j gj …^ x† ˆ 0 okj kˆ^k

…5:5†

for j ˆ 1; . . . ; m;

…5:6†

^ P 0; l

…5:7†

^ k P 0:

…5:8†

^ is a proper noninferior solution, then we have the following equality instead of Eq. (5.8): If x ^ k > 0:

…5:9†

Next, we consider the optimization problem de®ned as follows: minimize f R …f …x†† ˆ f R …f 1 …x†; . . . ; f p …x††

…5:10†

x

subject to 2

3 g1 …x† 6 7 g…x† ˆ 4 ... 5 6 0;

…5:11†

gm …x†

where fR stands for a multi-objective function composed of the functions, f1 …x†; . . . ; fp …x†. We de®ne the Lagrangian by L2 ˆ fR …f …x†† ‡ kt g…x†

t

for k ˆ ‰k1 ; . . . ; km Š :

…5:12†

Using this Lagrangian, the necessary condition for the optimality (KKT-condition) is described as follows: ^ is an optimal solution, there must exit ^ If x k such that p X oL2 ofR …f …x†† ofj …x† t og…x† ^ ˆ ‡k ˆ 0 for i ˆ 1; . . . ; n; …5:13† ofj …x† oxi xˆ^x oxi xˆ^x oxi xˆ^x jˆ1 oL2 ˆ gj …^ x† 6 0 okj kˆ^k

for j ˆ 1; . . . ; m;

…5:14†

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

oL2 k^j ˆ k^j gj …^ x† ˆ 0 okj kˆ^k

for j ˆ 1; . . . ; m;

4473

…5:15†

^ k P 0:

…5:16†

Comparing Eqs. (5.13)±(5.16) with Eqs. (5.4)±(5.8), if ofR …f …x††=ofj …x† for j ˆ 1; . . . ; p is not negative, Eqs. (5.13)±(5.16) are identical to Eqs. (5.4)±(5.8) because we can ®nd lj corresponding to ofR …f …x††=ofj …x† for j ˆ 1; . . . ; p. This means that if an optimal solution is obtained by the multi-objective fucntion, and ofR …f …x††=ofj …x† is not negative for j ˆ 1; . . . ; p, then this optimal solution is one of the Pareto optima. Note that even though Eqs. (5.4)±(5.8) are the necessary conditions, if we consider only the noninferior solutions, the identi®cation of Eqs. (5.13)±(5.16) to Eqs. (5.4)±(5.8) is sucient to guarantee that this optimal solution is one of the Pareto optima. This is due to the fact that the noninferior solutions which are satis®ed with Eqs. (5.4)±(5.8) are satis®ed with the sucient condition (the second-order condition). Therefore, if the multi-objective function has the condition ofR …f …x††=ofj …x† > 0 for 8 x and for j ˆ 1; . . . ; p and if we obtain an optimal solution using this multi-objective function, then this optimal solution is one of the Pareto optima. Typical functions which are satis®ed with this condition are the multiplication of two or more objective functions, the logarithm function, and the exponential function. Using a multi-objective function composed of these functions, we obtain one of the Pareto optima. 5.2. Formulation of new multi-objective functions and their characteristics In this subsection, we consider the multi-objective function to obtain the appropriate optimal solution which has physical meaning. 5.2.1. Unconstrained single ¯exibility case First, we formulate the multi-objective function in the case of unconstrained single ¯exibility case. As we explained in Section 4.1, we must simultaneously take into account the three objective functions: the mutual mean compliance, L2 …u1 †, must be maximized and the two mean compliances, L3 …u3 † and L4 …u4 †, must be minimized as indicated in Eq. (4.1). To ®nd the incorporating optimal solution (a noninferior solution), we formulate the following multi-objective function: maximize fUD ˆ

a;b;c;/;w; and h

ws

L3 …u3 †

L2 …u1 † ; ‡ …1 ws †L4 …u4 †

…5:17†

or more generally maximize fUL ˆ W log…L2 …u1 ††

a;b;c;/;w; and h

…1

W † log…ws L3 …u3 † ‡ …1

ws †L4 …u4 ††;

…5:18†

where W is the weighting coef®cient determining the weight between the mutual mean compliance L2 …u1 † and two mean compliances L3 …u3 † and L4 …u4 †. ws is the weighting coef®cient determining the weight between mean compliance L3 …u3 † and mean compliance L4 …u4 †. Note that when we use the multi-objective function formulated in Eq. (5.18), we must ®nd an appropriate initial solution where L2 …u1 † is positive because the logarithm function in this multi-objective function can only deal with the positive value. If a simple initial solution such as the solution composed of the uniform values of design variables is used, L2 …u1 † may be negative. To avoid this problem, we can use, for example, the optimal solution obtained by Eq. (5.17) as the initial solution of the optimization problem using the multi-objective function Eq. (5.18). Since these multiobjective functions are con®gured with the division operator or the logarithm function, the optimal solution obtained by Eq. (5.17) or Eq. (5.18) is one of the Pareto optima from our discussion in Section 5.1. Note that the division operator in the case of this multi-objective optimization problem is equivalent to the multiplication operator explained in Section 5.1, because the mutual mean compliance is maximized and two mean compiances are minimized while all functions are minimized in Section 5.1.

4474

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Taking the variations of fUD in Eq. (5.17) and fUL in Eq. (5.18) yields dfUD

L2 …u1 † ˆ …ws L3 …u3 † ‡ …1 ws †L4 …u4 ††

dfUL ˆ W

dL2 …u1 † L2 …u1 †

…1

W†



dL2 …u1 † L2 …u1 †

ws dL3 …u3 † ‡ …1 ws L3 …u3 † ‡ …1

ws dL3 …u3 † ‡ …1 ws L3 …u3 † ‡ …1

 ws †dL4 …u4 † ; ws †L4 …u4 †

ws †dL4 …u4 † : ws †L4 …u4 †

…5:19†

…5:20†

Eqs. (5.19) and (5.20) imply that the equivalent weighting coecient of the small perturbation of the mutual mean compliance, L2 …u1 †, is proportional to the inverse of L2 …u1 †. Therefore, as L2 …u1 † becomes larger, these objective functions provide L2 …u1 † with the smaller value as the wieghting coecient. That is, if L2 …u1 † tends to in®nity, its equivalent coecient automatically tends to zero and prevents this optimization problem from being ill-conditioned. Thus, we can have an appropriate optimal solution which has physical meaning on the Pareto optima using either Eq. (5.17) or Eq. (5.18). 5.2.2. Constrained single ¯exibility case We formulate the multi-objective function in the case of the constrained single ¯exbility case here. As we explained in Section 4.2, we consider two portions as implementation of a relative deformation constraint; the portion where we need to have the deformation for the kinematic requirement, corresponding to boundary Ct2 in Fig. 5 or the portion where we apply the traction, corresponding to boundary Ct1 in Fig. 6. Let us consider the three-dimensional case shown in Fig. 5(b). The multi-objective function in the twodimensional case can be obtained by setting t6 to 0 in the three-dimensional case. In the three-dimensional case, we must simultaneously take into account the ®ve objective functions: the mutual mean compliance, L2 …u1 †, must be maximized; the two mean compliances, L3 …u3 † and L4 …u4 †, must be minimized; and the absolute value of two mutual mean compliances, L5 …u1 † and L6 …u1 †, must be minimized. To ®nd the incorporating optimal solution (a noninferior solution), we formulate the following multiobjective function: maximize fCD ˆ  a;b;c;/;w; and h 2 ws L3 …u3 † ‡ …1

L2 …u1 † 2

2

ws †L4 …u4 † ‡ wM1 L5 …u1 † ‡ wM2 L6 …u1 †

2

12 ;

…5:21†

or more generally maximize fCL ˆ W log…L2 …u1 ††

a;b;c;/;w; and h

1 …1 2

 2 W † log ws L3 …u3 † ‡ …1

 2 2 2 ws †L4 …u4 † ‡ wM1 L5 …u1 † ‡ wM2 L6 …u1 † ;

…5:22†

where wM1 and wM2 are the weighting coecients for the relative deformation constraint. Next, we consider the three-dimensional case shown in Fig. 6(b). The multi-objective function the twodimensional case can be obtained by setting t8 to 0 in the three-dimensional case. In the three-dimensional case, we also must simultaneously take into account the four objective functions: the mutual mean compliance, L2 …u1 †, must be maximized; the two mean compliances, L3 …u3 † and L4 …u4 †, must be minimized; and absolute value of two mutual mean compliances, L7 …u1 † and L8 …u1 †, must be minimized. Therefore, we formulate a multi-objective function as follows: maximize fCD ˆ  a;b;c;/;w; and h ws L3 …u3 †2 ‡ …1

L2 …u1 † ws †L4 …u4 †2 ‡ wM1 L7 …u1 †2 ‡ wM2 L8 …u1 †2

12 ;

…5:23†

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4475

or more generally maximize fCL

a;b;c;/;w; and h

ˆ W log…L2 …u1 ††

 2 W † log ws L3 …u3 † ‡ …1

1 …1 2

 2 2 2 ws †L4 …u4 † ‡ wM1 L7 …u1 † ‡ wM2 L8 …u1 † ; …5:24†

Note that, since these multi-objective functions are also con®gured with the division operator or the logarithm function, the optimal solution obtained by these functions is one of the Pareto optima. 5.2.3. Multi-¯exibility case As we discussed in Section 4.3, if we need to implement n ¯exibilities in the ¯exible structure for the kinematic (motion) requirement, the necessary condition for the mutual mean compliance i L2 …i u1 † for i ˆ 1; . . . ; n corresponding to the ith ¯exibility) is that each mutual mean compliance must be at least positive in order to obtain deformation in the speci®ed direction, although a larger mutual mean compliance provides a higher ¯exibility performance. Furthermore, we also need to implement 2n stiffnesses to the ¯exible structure for the structural requirement. This requirement is satis®ed by minimizing the mean compliances …i L3 …i u3 † and i L4 …i u4 † for i ˆ 1; . . . ; n correspondign to the two ith stiffnesses). To ®nd the incorporating optimal solution which is satis®ed with all the requirements, we formulate the following multi-objective function using an extension of the Kreisselmeier±Steinhauser function (KS-functions) [44]: , n n X X
1 nL i L2 …i u1 † 1i 3 i 3 i 4 i 4 maximize fMD ˆ log e w2i L … u † ‡ w2i …5:25† s s L …u † ; L a;b;c;/;w; and h n iˆ1 iˆ1 or more generally n X 1 e L log n iˆ1

maximize fML ˆ W log

a;b;c;/;w; and h

…1

W † log

n X iˆ1

! nL i L2 …i u1 †

1i 3 i 3 w2i L …u † s

‡

i 4 i 4 w2i s L …u †




! ;

…5:26†

where nL isPthe positive and constant value, and wis is the weighting coecient for the mean compliances 2n such that iˆ1 wis ˆ 1: To verify the necessary condition, we consider the characteristics of the numerator of Eq. (5.25), fF ˆ

n X 1 e L log n iˆ1

nL i L2 …i u1 †

:

…5:27†

Taking the variation of Eq. (5.27) yields Pn dfF ˆ

Li 2 i 1 …u †

di L2 …i u1 †e n L Pn nL i L2 …i u1 † iˆ1 e

iˆ1

:

…5:28†

Eq. (5.28) implies that the equivalent coecient of the small perturbation of the mutual mean Pn weighting Li 2 i 1 Li 2 i 1 compliance, i L2 …i u1 †; is e n L … u † = iˆ1 e n L … u † , and this coecient becomes larger as the mutual compliance, i L2 …i u1 †, becomes smaller. That is, in the optimization procedure, the function fF provides the smaller mutual compliance, i L2 …i u1 †, for i ˆ 1; . . . ; n, with the larger equivalent coecient. Therefore, in the ®rst stage of the optimization procedure, this function can make all mutual mean compliances positive. Moreover, since the function fF is composed of the logarithm function and the exponential functions, Eqs. (5.25) and (5.26) are con®gured with the division operator or the logarithm function, along with the function fF , and the following inequality is satis®ed:

4476

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504 Li 2 i 1

ofF e n L …u † ˆ > 0; P n nL i L2 …i ui † oL2 …i u1 † iˆ1 e

…5:29†

the optimal solution obtained by Eq. (5.25) or Eq. (5.26) is one of the Pareto optima from our discussion in Section 5.1 In the denominator of Eq. (5.25), we have 2n weighting coecients, wis for i ˆ 1; . . . ; 2n, in order to determine the weights for the 2n mean compliances, i L3 …i u3 † and L4 …i iu4 †, for i ˆ 1; . . . ; n. In some designs of ¯exible structures which provide additional functions suchs as compliant mechanisms, we must set each weighting coecient wis because we need larger sti€ness at some portions in the ¯exible structures. However, in other ¯exible structures design such as automotive body design, since the ¯exible structure must have sti€ness that is almost equivalent for the structural requirement, 2n mean compliances …i L3 …i u3 †, and …i L4 …i u4 † for i ˆ 1; . . . ; n† must be set as equivalently as possible. In order to achieve this requirement, we formulate the following objective functions using the KS-function: , n n   X X 1 1 Si 3 i 3 Si 4 i 4 nL i L2 …i u1 † maximize fMD ˆ log e log en L … u † ‡ e n L … u † ; …5:30† L S a;b;c;/;w;and h n n iˆ1 iˆ1 or more generally maximize fML ˆ W log

a;b;c;/;w; and h

…1

n X 1 e L log n iˆ1

W † log

! nL i L2 …i u1 †

n  X 1 Si 3 i 3 log en L … u † ‡ e S n iˆ1

nS i L4 …i u4 †

! 

;

…5:31†

where nS is the positive and constant value. Note that the con®guration of the denominator in Eq. (5.30) is the same as the con®guration of the KS-function. Since the KS-function is also composed of the logarithm function and the exponential functions, and Eqs. (5.30) and (5.31) are con®gured with the division operator or the logarithm function, the optimal solution obtained by Eq. (5.30) or Eq. (5.31) is one of the Pareto optima from our discussion in Section 5.1.

6. Numerical implementation In this section, we shall discuss the numerical implementation in order to obtain the optimal design of ¯exible structures. Fig. 9 shows a ¯owchart of the optimization procedure. There are six steps in the perinteration loop of the optimization. 6.1. Calculation of the homogenized coecients for the elasticity tensor First, the homogenized coecients for elasticity tensor E H are calculated using the ®nite element method. The numerical procedure to obtain homogenized elasticity tensor E H has been discussed by Guedes and Kikuchi (see details in [45]). Using the ®nite element method, the numerical values of homogenized elasticity E H are calculated at six discretized points between 0 and 1 of microscopic design variables a and b in the two-dimensional case, or a; b, and c in the three-dimensional case. Therefore, the numerical values of E H are obtained at 36 points in the two-dimensional case and at 216 points in the threedimensional case. Next, E H is interpolated using Bezier curves. The numerical values obtained by the ®nite element method are used as control points to construct the Bezier curves. Furthermore, the sensitivities of E H , with respect to the microscopic variables, are also interpolated using Bezier curves. An advantage of using Bezier curves is that we have smooth curves which are not necessarily oscillating like the higher-order Lagrange polynomials. These interpolated sensitivities are used in order to calculate the sensitivities of the objective functions in each optimization problem.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4477

Fig. 9. Flowchart of optimization procedure.

6.2. Discretization of design domain using FEM As we discussed in Section 3.1, we deal with the extended design domain D in order to obtain the topological optimal structure. To solve the optimization problem numerically, this extended design domain D is discretized by the ®nite elements. In each ®nite element, we approximate that the con®guration of the microstructure is the same inside the element. That is, in each element, we have three design variables, ai ; bi ; and hi , in the two-dimensional case, and six design varibales, ai ; bi ; ci ; /i ; wi , and hi , in the threedimensional case for i ˆ 1; . . . ; n, where n is the number of the elements. Therefore, in the whole extended design domain D we have a total of 3n design variables in the two-dimensional case, and 6n design variables in the three-dimensional case. For the design variables, ai ; bi , and ci , we set upper bounds, aupp ; bupp ; and cupp , respectively, which are speci®ed as suf®ciently large positive numbers, but less than one in order to avoid singularity in the structural analysis. Therefore, in the actual optimization process, the design variables, ai ; bi , and ci , for i ˆ 1; . . . ; n, are bounded as follows: 0 6 ai 6 aupp < 1; 0 6 bi 6 bupp < 1;

…6:1† …6:2†

0 6 ci 6 cupp < 1 …in the three-dimensional case†:

…6:3†

In the two-dimensional analysis, four-node isoparametric full integration elements with bilinear shape functions are used. In the process of constructing the element sti€ness matrix, the full integration sti€ness matrix is obtained in the following equation [46]: RI CORR K FI ; e ˆ Ke ‡ Ke

…6:4†

RI where K FI e is the full integration element sti€ness matrix, K e the rank-de®cient element sti€ness matrix by CORR reduced integration, and K e is the corrective element sti€ness matrix for the hourglass control. This corrective element sti€ness matrix is calculated without any arti®cial parameters to control hourglass modes. The details of this method were discussed by Koh and Kikuchi (see [46]).

4478

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

In the three-dimensional analysis, a new element which uses the selective integration method and a bubble node is used. The full integration (FI) method has been widely used because stability is always achieved and approximate solution always converge to the exact one, as the number of elements in®nitely increase. However, the FI method involves several diculties such as shear locking for the bending dominated problems and volumetric locking problem for nearly incompressible materials. This is because the FI method overestimates an element sti€ness matrix [47]. Moreover, the FI method requires many computational operations to construct an element matrix. The reduced integration (RI) method is useful in order to avoid these problems, but su€ers instabilities such as the so-called hourglass modes which are due to the rank-de®ciency of an element matrix. The selective reduced integration (SRI) method has been developed to improve the accuracy of the approximate solutions in the problems involved with the FI method. Malkus and Hughes [48] developed the SRI scheme and examined the numerical results obtained by SRI. Hughes [49] re®ned this scheme and proposed the general SRI procedure for anisotropic incompressible material. Liu et al. [50] developed a uni®ed approach with hourglass and volumetric locking control. Koh and Kikuchi [46] developed a method to control the hourglass modes with analytical correction terms (K CORR in Eq. (6.4)) which are obtained by the SRI scheme. e In the design of ¯exible structure, the shear locking problem is critical since the ¯exibility provides the structures with the large bending or torsion, especially in the three-dimensional case. Therefore, the shear stress due to the bending or torsion must be correctly estimated during the optimization procedure. Thus, we developed a new element using SRI and a bubble node. Suppose a hexagonal element as shown in Fig. 10. This element is composed of 8 nodes and 1 bubble node. The bubble node is located at the center of the element. This bubble node is required to improve the accuracy of the solutions in the case of bending dominated problems and incompressible materials. First, we decompose the isotropic elasticity tensor E (homogenized elasticity tensor E G in the optimization procedure) into the following ®ve terms: E ˆ E N ‡ E V ‡ E yz ‡ E zx ‡ E xy ;

…6:5†

where E N is the normal-deviatoric part of the elasticity tensor, E V the volumetric part of the elasticity tensor, and E yz ; E zx , and E xy are the shear parts of the elasticity tensor. These are composed of the following matrix: 2

1 m0 =…1 m0 † 6 m =…1 m † 1 0 6 0 6 6 m0 =…1 m0 † m0 =…1 m0 † E…1 m0 † 6 EN ˆ …1 ‡ m0 †…1 2m0 † 6 0 0 6 6 4 0 0 0 0

m0 =…1 m0 † m0 =…1 m0 † 1 0 0 0

Fig. 10. A hexagonal element.

0 0 0 0 0 0

0 0 0 0 0 0

3 0 07 7 7 07 7; 07 7 7 05 0

…6:6†

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

2

d2 d1 d2 0 0 0

d1 6 d2 6 6 d2 E V ˆ E6 60 6 40 0 2

0 6 60 6 60 E yz ˆ E6 60 6 60 4 0 2

0 6 60 6 60 E zx ˆ E6 60 6 60 4 0

d2 d2 d1 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

…6:7†

0 0 0 0 0 0

3 0 7 07 7 07 7; 07 7 07 5 0

0 0 0 0 0 0

0 0 0 0 1=2…1 ‡ m† 0

3 0 7 07 7 07 7; 07 7 07 5 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 1=2…1 ‡ m† 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

2 60 60 6 60 E xy ˆ E6 60 6 60 4 0

3 0 07 7 07 7; 07 7 05 0

4479

…6:8†

…6:9†

3 0 7 7 0 7 7 0 7; 7 0 7 7 0 5 1=2…1 ‡ m†

…6:10†

where m0 is the threshold to decompose the normal part of the elasticity tensor into normal-deviatoric part E N and volumetric part E V , E the Young's modulus, and m is Poisson's ratio. d1 and d2 are expressed as follows: d1 ˆ

1 m …1 ‡ m†…1 2m†

1 m0 …1 ‡ m0 †…1 2m0 †;

…6:11†

d2 ˆ

m …1 ‡ m†…1

m0 …1 ‡ m0 †…1

…6:12†

2m†

2m0 †

;

Next, we construct the element matrix, K e , con®gured with 9 nodes (8 nodes and 1 bubble node) using SRI are as follows: Z Ke ˆ ˆ

1

Z

1

Z

1

1

1

1

iˆ1

jˆ1

kˆ1

2 X 2 X 2 X

‡

2 X iˆ1

Bt EBJ dn dg df B t E N BJ dn dg dfj…n;g;f†ˆ…ai ;aj ;ak † ‡ 8B t E V BJ j…n;g;f†ˆ…0;0;0†

4B t E yz BJ j…n;g;f†ˆ…ai ;0;0† ‡

2 X jˆ1

4Bt E zx BJ j…n;g;f†ˆ…0;aj ;0† ‡

2 X kˆ1

4B t E xy BJ j…n;g;f†ˆ…0;0;ak †;

…6:13†

4480

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

where matrix B is the strain±displacement matrix obtained by the shape function, J the and the p pJacobian,  notations ai ; aj and ak denote 1 3 if the subscripts i; j, or k are equal to 1, and 1 3 if the subscripts i; j; or k are equal to 2. After obtaining this element matrix, we condense each part of the 9-node element matrix into the corresponding parts of the 8-node element matrix using the Guyan reduction in order to eliminate the 1 bubble node (see [51] for details of the element con®gurations). Using these elements, we solve the FE problems in the optimization procedure. After obtaining the strains and stresses in each case of the optimization problem, we compute the mutual mean compliance, the mean compliances, the objective function, and their sensitivities with respect to the microscopic variables ai ; bi , and ci , in each element. 6.3. Optimization method The optimality criteria method has been used in the topology optimization (e.g. [22,29,31,32,52]) based on the ®xed grid method. This is because the optimality criteria method can deal with numerous design variables. However, Ma et al. [32] pointed out that the standard optimality criteria method cannot be applied to the dynamic problem. To overcome this problem Ma et al. introduced a convex approximation, which is the basis of sequential convex programming such as CONLIN [53] and MMA [54]. Thus, in the optimality criteria method, the heuristics determining the update rule must be adequately adjusted for di€erent objective functions or di€erent optimization problems. Recently, the mathematical programming method has been used to solve the topology optimization problem (e.g. [19,55,56]) based on the ®xed grid method. This is because the mathematical programming method is ¯exible enough to set any optimization problem. That is, it can theoretically deal with any objective function. Moreover, it is easy to implement by combining the commerical packages, especially for SLP. In this study, we use SLP as the optimizer. This is because SLP can deal with a variety of objective functions, it can handle numerous design variables (more than 10,000), and it requires only sensitivities of the objective functions and constraints, although the fast convergence cannot be expected. SLP uses a trust region scheme because, as Haftka and G urdal [57] noted, this linearized approximation is available in a neighborhood of the design variables. Therefore, SLP usually sets the move limits and approximates the feasible design domain to the limited region. There are many ways to set the move limits. Haftka and G urdal [57] suggested that the move limits be set in the range of 10±30% of the design variables. However, they pointed out that this choice is reasonable only if a design variable is not exceedingly small. Thomas et al. [58] discussed a way to set the move limits in the case of shape optimization. They also pointed out that if the move limits are based on a percentage of the current design variable value, the design variable, will never be able to cross zero. Even if we only consider the positive value of the design variable, we have a similar problem when the design variable is exceedingly small as Hafka and G urdal [57] pointed out. That is, in this case, the move limits will lock the design variable. To overcome this diculty, Thomas et al. set move limit DA of design variable A in the following way: DA ˆ Max…fA; DAmin †;

…6:14†

where f is a constant move ratio, and DAmin is the minimum move limit. On the other hand, Chen [59] and Kirsch [60] developed the methods to calculate e€ective move limits. They also veri®ed these methods using the shape optimization problem. In this research, the move limits are simply set using Eq. (6.14) in order to avoid the locking of the move limits in the case of exceedingly small design variables. f is set to 0.1 in the two-dimensional case, and 0.15 in the three-dimensional case. DAmin is set to 1% of the maximum value of the design variable, i.e., 0.01. When the objective function reaches the proximity of the optimal value, f is decreased to half of the original value. In the linearized optimization problem, a package of the simplex method, DSPLP from the SLATEC library [61], is used.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4481

6.4. Filtering for the checkerboard Checkerboard patterns often appear when the ®xed grid method is used in the topology optimization. They are de®ned as a periodic pattern of high and low values of the density in the checkerboard shown in Fig. 11(b). Jog et al. [62] explained that if four-node isoparametric elements in the two-dimensional case, or eightnode isoparametric elements in the three-dimensional case are used, checkerboard patterns appear because of the insuciency of their degrees of freedom. Jog et al. and Bendsùe [17] claimed that this insuciency is explained in a way similar to the violation of the so-called Babuska±Brezzi or LBB condition. Dõaz and Sigmund [63] explained that if four-node isoparametric elements with bilinear shape functions are used, checkerboard patterns appear because the pattern of checkerboard shown in Fig. 11(b) is arti®cially sti€er in terms of the strain energy than the pattern in which each element has a uniform density of 1/2 as shown in Fig. 11(a). Kikuchi et al. [64] discussed the checkerboard mechanism using the eigenvalue analysis of an element. They mentioned that in the shearing mode and volumetric mode, the pattern of checkerboard has much higher eigenvalues than the pattern of the uniform density. Therefore, in the case of the topology optimization problem which maximizes the global sti€ness, if elements are imposed by loads along with either the shearing mode or the volumetric mode, the checkerboard patterns will appear in these elements. Furthermore, they suggested that there are two types of methods to eliminate the checkerboard patterns. The ®rst is introduction of an appropriate ®ltering scheme, and the second is use of appropriate ®nite elements which possess either equal or smaller eigenvalue in the checkerboard patterns than in the uniform density patterns. In this study, the modi®ed ®ltering scheme based on Bendsùe's method [17] is developed to eliminate the checkerboard patterns. The algorithm of this ®ltering scheme is described as follows: 1. Construct one group consisting of the four contiguous ®nite elements in the two-dimensional cases, or eight contiguous ®nite elements in the three-dimensional cases, shown in Fig. 12. The shape of all contiguous ®nite elements is assumed to be the same. The microscopic variables a and b of the ith element for i ˆ 1; . . . ; 4 in the two-dimensional case are described as ai and bi , respectively. The density qi of the ith element is calculated by qi ˆ 1

ai bi

for i ˆ 1; . . . ; 4:

…6:15†

Similarly, in the three-dimensional case, the microscopic variables a; b; and c of the ith element for i ˆ 1; . . . ; 8 are described as ai ; bi ; and ci ; respectively. The density qi of the ith element is calculated by qi ˆ 1

ai bi ci

for i ˆ 1; . . . ; 8:

…6:16†

2. Search for the maximum and minimum values of qi among the elements. The maximum value is set to qmax , and the minimum value is set to qmin .

Fig. 11. Uniform density pattern and checkerboard pattern. (a) Uniform density pattern, (b) checkerboard pattern.

4482

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 12. Contiguous ®nite elements for checkerboard elimination. (a) Two-dimensional case, (b) three-dimensional case.

3. In the two-dimensional case, calculate the ®ltered density qi for i ˆ 1; . . . ; 4 using the following equations: 1 q1 ˆ …3q1 ‡ q2 ‡ q3 q4 †; 4 1 q2 ˆ …q1 ‡ 3q2 q3 ‡ q4 †; 4 1 q3 ˆ …q1 q2 ‡ 3q3 ‡ q4 †; 4 1 q4 ˆ … q1 ‡ q2 ‡ q3 ‡ 3q4 †: 4

…6:17† …6:18† …6:19† …6:20†

Similarly, in the three-dimensional case, calculate the ®ltered density qi for i ˆ 1; . . . ; 8 using the following equations: 1 q1 ˆ …7q1 ‡ q2 ‡ q3 q4 ‡ q5 q6 q7 ‡ q8 †; 8 1 q2 ˆ …q1 ‡ 7q2 q3 ‡ q4 q5 ‡ q6 ‡ q7 q8 †; 8 1 q3 ˆ …q1 q2 ‡ 7q3 ‡ q4 q5 ‡ q6 ‡ q7 q8 †; 8 1 q4 ˆ … q1 ‡ q2 ‡ q3 ‡ 7q4 ‡ q5 q6 q7 ‡ q8 †; 8 1 q5 ˆ …q1 q2 q3 ‡ q4 ‡ 7q5 ‡ q6 ‡ q7 q8 †; 8 1 q6 ˆ … q1 ‡ q2 ‡ q3 q4 ‡ q5 ‡ 7q6 q7 ‡ q8 †; 8 1 q7 ˆ … q1 ‡ q2 ‡ q3 q4 ‡ q5 q6 ‡ 7q7 ‡ q8 †; 8 1 q8 ˆ …q1 q2 q3 ‡ q4 q5 ‡ q6 ‡ q7 ‡ 7q8 †: 8

…6:21† …6:22† …6:23† …6:24† …6:25† …6:26† …6:27† …6:28†

4. Search for the maximum and minimum values of qi among the elements. The maximum value is set to qmax , and the minimum value is set to qmin . 5. If the following inequalities are satis®ed in each group consisting of the four contiguous ®nite elements in the two-dimensional case, or the eight contiguous ®nite elements in the three-dimensional case, replace the density of each element with qi : Otherwise, the density is not ®ltered.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

qmax < qmax

and

qmin > qmin :

4483

…6:29†

Note that only if these equalities are satis®ed, will this group have a checkerboard pattern. 6. In the elements in which the density is ®ltered, calculate the equivalent microscopic variables as follows: ai ˆ gi ai ; bi ˆ gi bi ; ci ˆ gi ci …three-dimensional case†;

…6:30† …6:31† …6:32†

1=2 1=3 where gi ˆ ……1 q†=ai bi † in the two-dimensional case, and gi ˆ ……1 q†=ai bi ci † in the three-dimensional case. If severe checkerboard patterns are recognized in the optimal con®guration, this ®ltering scheme is used to eliminate them.

6.5. Update schemes of angles The angles /i ; wi ; and hi must be updated to maximize each multi-objective function as well. However, these angles are practically updated to minimize the mean compliances, that is, to maximize the sti€ness in the ¯exible structures for simplicity of computation. This is because if the mean compliances are minimized, any multi-objective function discussed in this research tends to be increased to almost the maximum value. To minimize the mean compliances, the direction of these angles are updated to the principal direction of the stress (see details in Pedersen [65]) using the multi-loading criteria which are proposed by Suzuki and Kikuchi [29] in the minimum mean compliance (the maximum global sti€ness) problem. That is, the angles are updated using the following equations. In the unconstrained single ¯exibility case and the constrained single ¯exibility case, minimize Max…L3 …u3 †; L4 …u4 ††; /;w; and h

…6:33†

In the multi-¯exibility case, minimize Max …i L3 …i u3 †; i L4 …i u4 ††; /;w; and h iˆ1;...;n

…6:34†

where n is the number of ¯exibilities implemented in the ¯exible structures. 6.6. Penalization of density The optimal con®gurations obtained by the homogenization design method may contain gray-scale portions where the value of density q is greater than 0 and less than 1, and the gray-scale portions mean the existence of the microscale structure even though the rectangle type microstructure is used. Since these portions are dicult to interpret as a real structure to be manufactured, or are dicult to provide candidate con®gurations for the shape optimization phase, the existence of gray-scale portions are not appropriate, although the con®guration having the gray-scale distributions are optima. Several methods have been introduced to control the number of macroscale holes in the solutions or to penalize the optimal con®gurations (see [39,52,64]). The principal idea is to push the density distribution of the microstructure toward either 0 or 1 by solving the sub-optimization problem using a penalty function. In this research, we obtain the penalized value A~ of a design variable A using the following logistic function fp …A†: er…A‡s† er…A A~ ˆ fp …A† ˆ r…A‡s† e ‡ er…A

s† s†

;

…6:35†

where A stands for ai or bi in the two-dimensional case, and ai ; bi ; or ci in the three-dimensional case for i ˆ 1; . . . ; n; here, n is the number of the elements. Using the logistic function, a design variable is penalized to either 1 or 0 during the optimization process. Parameters r and s are not related to the design variable. Parameter r stands for the slope of the logistic function. This must be suf®ciently large for the penalization

4484

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

of a design variable. Parameter s is a shift parameter which determines whether a design variable rises to 1 or falls to 0. Using the optimal con®gurations as the initial con®gurations, we solve the following optimization problem to eliminate gray-scale portions: maximize f

…6:36†

s

subject to Equilibrium equations; Z q dX Xs DXs 6 0;

…6:37† …6:38†

Xd

where f is a multi-objective function discussed in Section 5, and DXs is the additional value to relax the total volume constraint Xs . This additional volume is required since the penalization of a design variable to one may violate the original total volume constraint Xs : The penalization scheme using the logistic function is achieved during the optimization in Eqs. (6.36)±(6.38). 7. Numerical examples In this section, several numerical examples are presented to examine the con®gurations. In all the examples, the properties of the isotropic material correspond to Young's modulus ˆ 100 and Poisson's ratio ˆ 0.3. As we explained in Section 6.2, in the two-dimensional analysis, four-node isoparametric full integration elements with bilinear shape functions are used. In the three-dimensional analysis, a new element which uses the selective integration method and a bubble node is used. In all cases, the applied force is assumed to be a unit load. 7.1. E€ect of mesh size The ®rst example is used to examine the uniqueness of the optimal topology con®guration when ®nite element meshes are uniformly re®ned, while other conditions are ®xed. Fig. 13 shows the two-dimensional design domain speci®ed as a 60  30 rectangle with a ®xed support boundary on the left-hand side. In this two-dimensional design domain, unconstrained single ¯exibility is implemented. Using the multiobjective function by Eq. (5.17) in the unconstrained single ¯exibility case, the deformation at point P2 in the direction of dummy load F 2 is to be maximized when the external force, F 1 , is applied at point P1 , while

Fig. 13. Design domain for a ¯exible structure.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4485

the sti€nesses at both points P1 and P2 are to be maximized. That is, the mutual mean compliance de®ned by the two loads, F 1 and F 2 , is to be maximized for the kinematic requirement, while the two mean compliances at both points P1 and P2 are to be minimized for the structural requirement. Loads F 1 and F 2 correspond to tractions t1 and t2 in Section 4.1, respectively. Three ®nite element models are utilized for this veri®cation. One is a coarse mesh using 60  30 ®nite elements, another is an intermediate mesh using 80  40 ®nite elements, and the third is a ®ne mesh using 120  60 ®nite elements. The total volume constraint of the material Xs is considered to be 360, which is 20% of the volume of the whole design domain. The initial values of microscopic design variables ai and bi are set to 0.90, and the value of hi is set to 0.0 for i ˆ 1; . . . ; the number of elements (in all elements). The weighting coecient, ws ; in Eq. (5.17) is set to 0.4 in all mesh sizes. In all cases, we use the ®ltering scheme to eliminate the checkerboard patterns. Fig. 14 illustrates the convergence history of the optimization process in the case of 60  30 ®nite elements. Notice that the mutual mean compliance, L2 …u1 †; increases and the mean compliances, L3 …u3 †; and L4 …u4 †; decrease simultaneously as the multi-objective function de®ned by Eq. (5.17) increases. Moreover, at the optimal point, the volume constraint de®ned by Eq. (4.11) becomes active.

Fig. 14. Convergence history (60  30 ®nite elements).

4486

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 15. Optimal con®gurations in the unconstrained single ¯exibility case (using the multi-objective function de®ned by Eq. (5.17), Xs ˆ 360 (20%), ws ˆ 0:4). (a) 60  30 elements, (b) 80  40 elements, (c) 120  60 elements.

Fig. 15 shows the results of the optimal con®gurations with respect to the various meshes. It is clear that the optimal topology con®gurations using the three di€erent mesh sizes are similar, although slight differences in shape are observed in some portions. This means that the optimal topology con®guration is not a€ected by the size of mesh. 7.2. E€ect of the weighting coecient W In this section, we examine optimal con®gurations obtained by the other multi-objective function de®ned by Eq. (5.18) in the single ¯exibility case. That is we investigate the e€ect of the weighting coecient W on

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4487

the optimal con®gurations using the same design model shown in Fig. 13. Using the multi-objective functions formulated by Eq. (5.18) in the unconstrained single ¯exibility case, the mutual mean compliance de®ned by the two loads, F 1 and F 2 , is to be maximized for the kinematic requirement, while the two mean compliances at both points P1 and P2 are to be minimized for the structural requirement. The total volume constraint of the material Xs is considered to be 360, which is 20% of the volume of the whole design domain. The design domain is discretized using 1800 ®nite elements. The initial values of microscopic design variables ai and bi are set to 0.90 and the value of hi is set to 0.0 for i ˆ 1; . . . ; the number of elements (in all elements). In all cases, the weighting coef®cient, ws ; in Eq. (5.18) is set to 0.4, and the ®ltering scheme is used to eliminate the checkerboard patterns. By using these initial conditions, we examine the relation between the optimal con®gurations and the weighting coef®cient W. Table 1 shows the values of the mutual mean compliance and the two mean compliances at the optimal point with respect to the weighting coecient W. It is clear that as W inreases, the mutual mean compliance and the two mean compliances increase. That is, as W increases, the ¯exibility increases, while the two stiffnesses decrease. Therefore, by changing W, we can adjust their values using Eq. (5.18). Fig. 16 shows the optimal con®gurations when the weighting coecient W is set to 0.1, 0.3, and 0.7. Note that the optimal con®guration in the case of W ˆ 0:5 is identical to the optimal con®guration obtained by Eq. (5.17), which is shown in Fig. 15(a). In the case of W ˆ 0:1, the optimal con®guration is similar to the optimal con®guration obtained by the minimum compliance problem. As W increases from 0.1 to 0.3, the portions connecting with the ®xed left-hand side go down so that the ¯exible structure obtains higher ¯exibility. Comparing Fig. 15(a) …W ˆ 0:5† and Fig. 15(b) …W ˆ 0:3†, the optimal con®guration in the case of W ˆ 0:5 is composed of the smooth curves. On the other hand, in the case of Fig. 15(c) (W ˆ 0:7), clear optimal con®guration is not obtained. This is because higher ¯exibility disconnects structure. Thus, in this design model, weighting coef®cient W must be less than or equal to 0.5 in order to obtain an appropriate ¯exible structure. Indeed, the multi-objective function formulated by Eq. (5.17) provides a critical optimal con®guration which has structural continuity in the ¯exibile structure since Eq. (5.17) is made identical to Eq. (5.18) by setting W to 0.5. 7.3. Design of compliant gripper and e€ect of total volume constraint Xs We examine here the relation between the optimal con®guration and the total volume constraint Xs by designing a compliant gripper. Fig. 17 shows a half view of the two-dimensional design domain where boundary condiions and speci®cations are as indicated. As shown in Fig. 17, the left-hand side boundary of the design domain is ®xed to support the gripper, and the symmetry boundary condition is posed at the bottom boundary. The function of the gripper is to (1) deform along the direction of dummy force F 2 in order to grasp a workpiece at point P2 when the external force, F 1 , is applied at point P1 for the kinematic requirement, and (2) hold the workpiece while the external force is continuously applied for the structural requirement. Using the multi-objective functions formulated by Eq. (5.17) in the unconstrained single ¯exibility case, the mutual mean compliance de®ned by the two loads, F 1 and F 2 , is to be maximized, while the mean compliances at both points P1 and P2 are to be minimized in order to satisfy the required function. Loads F 1 and F 2 correspond to tractions t1 and t2 in Section 4.1, respectively. In this problem, 1600 ®nite elements Table 1 Values of mutual mean compliance and two mean compliances at the optimal point with respect to weighting coecient W Weighting coecient, W

Mutual mean compliance, L2 …u1 †

Mean compliance, L3 …u3 †

Mean compliance, L4 …u4 †

0.1 0.3 0.5 0.7 0.9

47.93 143.42 150.20 191.27 555.20

0.305 0.367 0.453 1.212 50.16

0.102 0.113 0.156 0.439 16.970

4488

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 16. Optimal con®gurations in the unconstrained single ¯exibility case (using the multi-objective function de®ned by Eq. (5.18), Xs ˆ 360 (20%), ws ˆ 0:4). (a) W ˆ 0:1, (b) W ˆ 0:3, (c) W ˆ 0:7.

are used. The initial values of microscopic design variables ai and bi are set to 0.90 and the value of hi is set to 0.0 for i ˆ 1; . . . ; the number of elements (in all elements). In all cases, the weighting coecient, ws , in Eq. (5.17) is set to 0.4, and the ®ltering scheme is used to eliminate the checkerboard patterns. Three values of the total volume constraint Xs , 320, 400, and 480, which correspond to 20%, 25%, and 30% of the volume of the whole design domain, are examined. Fig. 18 shows the optimal topology con®gurations of the compliant gripper in the case of Xs ˆ 320; 400; and 480. Comparing these ®gures, the topology con®gurations are almost the same, although some differences in the optimal shape are observed. However, the material density increases near point P2 as the

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4489

Fig. 17. Design domain for a compliant gripper.

total volume constraint, Xs ; increases because the two mean compliances, L3 …u3 † and L4 …u4 †, decrease. This implies that it is necessary to have sucient material in order to obtain the sti€ness required to resist the external force at point P1 , and to hold a work-piece at point P2 for the structural requirement. On the other hand, the topology con®guration on the left-hand side, which connects to the ®xed boundary portions, does not change so much with respect to the total volume constraint. This implies that this con®guration provides sucient ¯exibility in order to satisfy the kinematic requirement. In summary, it is likely that the optimal topology con®guration is a€ected by changing the total volume constraint. To eliminate the intermediate gray-scale portion in Fig. 18(a), we apply the penalization method which is formulated as the optimization problem shown in Eqs. (6.36)±(6.38). In this case, parameter r in the logistic function is set to 50.0, and an additional volume constraint DXs is set to 160.0. The multi-objective function in Eq. (5.17) is maximized when s is equal to 0.42. Fig. 19 shows the penalized optimal con®guration. It is clear that using the penalization method presented here, gray-scale portions are eliminated. As we discussed in Section 4.1, we design the compliant mechanism based on the small linear deformation even though the compliant mechanism may deform considerably. Therefore, the design con®guration must be examined to determine whether this con®guration generates the speci®ed large deformation which we originally intended in the topology design phase. Next, we examine the performance of the optimal con®guration of the compliant gripper using nonlinear large displacement FE analysis. Fig. 20 shows the ®nite element model constructed based on the optimal con®guration shown in Fig. 18(a) (Xs ˆ 320 (20%)). This model is discretized to 1557 elements with 1853 nodes. The properties of the isotropic material corresponds to Young's modulus ˆ 100 and Poisson's ratio ˆ 0.3. Nonlinear commercial code ABAQUS with the large displacement option is used for the nonlinear analysis. Fig. 21 shows the deformation, and Fig. 22 shows the von-Mises stress distribution of the compliant gripper where an external force of 0.1 unit is applied at point P1 . In Fig. 22, the darker color indicates higher stress distribution. It is clear that the compliant gripper deforms according to the problem speci®cations from Fig. 21. That is, the portion around point P2 goes down in order to grasp a workpiece when the external force is applied at point P1 : The stress concentration is observed in two portions where the ¯exible structure has higher ¯exibility (see Fig. 22). 7.4. Design of compliant gripper in the constrained ¯exibility case and applicability of the ®ltering scheme In this section, we examine the design of the ¯exible structure in the case where the ¯exible structure is required to have a relative deformation constraint in the portion where we need to have the deformation for the kinematic requirement, corresponding to boundary Ct2 . That is, we examine the speci®cation of the deformation using a relative deformation constraint de®ned by Eq. (4.15) in the ¯exible structure by designing the compliant gripper which deforms continuously in the speci®ed translational direction. Fig. 23 shows a half view of the design domain. As shown in Fig. 23, we consider an additional dummy load, F 5 ,

4490

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 18. Optimal con®gurations of the complaint gripper (using the multi-objective function de®ned by Eq. (5.17), in the unconstrained single ¯exibility case, ws ˆ 0:4). (a) Xs ˆ 320, (b) Xs ˆ 400, (c) Xs ˆ 480.

which is perpendicular to dummy load F 2 specifying the direction of deformation at point P2 in addition to the same boundary conditions and speci®cations as the design domain shown in Fig. 17. Using the multi-objective function formulated by Eq. (5.22) in the constrained single ¯exibility case, the mutual mean compliance de®ned by the two loads, F 1 and F 2 , is to be maximized while the mean compliances at both points P1 and P2 ; and the absolute value of the mutual mean compliance de®ned by the two loads, F 1 and F 5 , are to be minimized in order to satisfy the required function, where loads F 1 and F 2 corresponds to tractions t1 and t2 in Section 4.1, respectively, and load F 5 corresponds to traction t5 in

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4491

Fig. 19. The penalized optimal con®guration of the complaint gripper (using the multi-objective function de®ned by Eq. (5.17) in the unconstrained single ¯exibility case, Xs ˆ 320 (20%), ws ˆ 0:4; r ˆ 50:0; s ˆ 0:42; DXs ˆ 160:0).

Fig. 20. Finite element model of the compliant gripper for nonlinear large displacement analysis.

Fig. 21. Deformed shape of the compliant gripper (large displacement analysis).

Section 4.2. The initial values of microscopic design variables ai and bi are set to 0.90 and the value of hi is set to 0.0 for i ˆ 1; . . . ; the number of elements (in all elements). The total volume constraint Xs is set to 320, which corresponds to 20% of the volume of the whole design domain. The weighting coecients, ws ; wM1 , and wM2 in Eq. (5.22) are set to 0.4, 1.0, and 0.0, respectively. First, the design domain is

4492

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 22. Von-Mises stress distribution of the compliant gripper (large displacement analysis).

Fig. 23. Design domain of a compliant gripper which deforms continuously in the speci®ed translational direction.

descretized using 1600 ®nite elements, and the ®ltering scheme is not used. This is because the optimal con®guration shown later has a complex topology to satisfy the problem speci®cation, the ®ltering scheme makes the optimal con®guration ambiguous using only 1600 ®nite elements, and the severe checkerboard patterns are not observed. Table 2 shows the values of the two mutual mean compliances and the two mean compliances at the optimal point with respect to weighting coecient W. Their values at the optimal point in the unconstrained single ¯exibility case shown Fig. 18(a) are also indicated for comparison. It is clear that in the constrained single ¯exibility case, the absolute value of the mutual mean compliance, L5 …u1 †, is much smaller than in the unconstrained single ¯exibility case although the mutual mean compliance, L2 …u1 †, is not Table 2 Values of mutual mean compliances and mean compliances at the optimal point with respect to weighting coecient W Weighting coecient W (a) Constrained single ¯exibility case 0.2 0.3 0.4 (b) Unconstrained single ¯exibility case

Mutual mean compliance 2

1

Mean compliance 5

1

L …u †

L …u †

L3 …u3 †

L4 …u4 †

0.815 4.702 6.709

)0.212 )0.019 )0.016

0.246 0.321 0.384

0.340 0.349 0.496

15.788

)12.436

0.230

0.277

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4493

so large (at most less than half of that in the unconstrained ¯exibility case). This means that the deformation at point P2 is well constrained with L5 …u1 †. Furthermore, as W increases, the mutual mean compliance, L2 …u1 † and the mean compliances, L3 …u3 † and L4 …u4 †, increase. That is, as W increases, the ¯exibility increases, while the two stiffnesses decrease. In fact, we could not obtain a clear optimal con®guration in the case of W ˆ 0:4 due to the lack of stiffness. Fig. 24 shows the optimal con®guration in the case of W ˆ 0:3. It is clear that the optimal con®guration in the constrained single ¯exibility case is quite di€erent from the optimal con®guration in the unconstrained ¯exibility case shown in Fig. 18(a). This optimal con®guration can be interpreted as a practical image structure shown in Fig. 25. Next, using the ®ne mesh, we con®rm the applicability of the ®ltering scheme. The design domain is descretized using the 6400 ®nite elements, and the ®ltering scheme is used to obtain the optimal con®guration. Fig. 26 shows the optimal con®guration using the ®ne mesh. By using ®ne mesh, we obtain the clear complex topology using the ®ltering scheme. Thus, the ®ne mesh eliminates the ambiguity of the complex optimal con®guration such as the large gray-scale portions. 7.5. Design of suspension-like ¯exible structure with a relative deformation constraint In this section, we consider the design of the ¯exible structure in the case where the ¯exible structure is required to have a relative deformation constraint in the portion where we apply the traction corresponding

Fig. 24. Optimal con®guration of the compliant gripper (using the multi-objective function de®ned by Eq. (5.22) in the constrained single ¯exibility case, Xs ˆ 320 (20%), W ˆ 0:3; ws ˆ 0:4; wM1 ˆ 1:0; wM2 ˆ 0:0).

Fig. 25. A practical image of optimal con®guration of the compliant gripper (using the multi-objective function de®ned by Eq. (5.22) in the constrained single ¯exibility case, Xs ˆ 320 (20%), W ˆ 0:3; ws ˆ 0:4; wM1 ˆ 1:0; wM2 ˆ 0:0).

4494

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 26. Optimal con®guration of the compliant gripper (using the multi-objective function de®ned by Eq. (5.22) in the constrained single ¯exibility case and ®ne mesh, Xs ˆ 320 (20%), W ˆ 0:3; ws ˆ 0:4; wM1 ˆ 1:0; wM2 ˆ 0:0).

to boundary Ct1 . That is, we examine the speci®cation of the deformation using a relative deformation constraint de®ned by Eq. (4.21) in the ¯exible structure by designing the suspension-like structure which deforms continuously in the speci®ed translational direction. Fig. 27 shows the two-dimensional design domain speci®ed as a 60  20 rectangle with a ®xed support boundary on the left-hand side. Now we suppose that point P1 is the portion where the entire input force is applied in the suspension-like ¯exible structure, and P2 is the portion that supports the absorber system. The necessary condition for the suspension-like ¯exible structure is that the portion around P1 must deform in the vertical direction to adjust the automotive body alignment. Using the multi-objective function formulated by Eq. (5.24) in the constrained single ¯exibility case, the mutual mean compliance de®ned by the two loads, F 1 and F 2 , is to be maximized while the mean compliances at both points P1 and P2 , and the absolute value of the mutual mean compliance de®ned by the two loads, F 1 and F 7 , are to be minimized in order to satisfy the required function, where loads F 1 and F 2 correspond to traction t1 and t2 in Section 4.1, respectively, and load F 7 corresponds to traction t7 in Section 4.2. The design domain is discretized using 4800 (120  40) ®nite elements. The initial values of microscopic design variables ai and bi are set to 0.90 and value of hi is set to 0.0 for i ˆ 1; . . . ; the number of elements (in all elements). The total volume constraint Xs is set to 240 which corresponds to 20% of the volume of the whole design domain. The weighting coecients, W ; ws ; wM1 , and wM2 in Eq. (5.24) are also set to 0.2, 0.4, 1.0, and 0.0, respectively. The ®ltering scheme is also used in this problem. Fig. 28 shows the optimal con®guration obtained by the multi-objective function equation (5.24). In order to compare it with the unconstrained single ¯exibility case, the optimal con®guration obtained by the multi-objective function equation (5.17) by setting ws to 0.4 is also described. Comparing the two ®gures,

Fig. 27. Design domain for a suspension-like ¯exible structure.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4495

Fig. 28. Optimal con®guration of the suspension-like ¯exible structure (Xs ˆ 240 (20%)). (a) The constrained single ¯exibility case, (b) the unconstrained ¯exibility case.

these topologies are quite di€erent. Thus, by implementing a relative deformation constraint, we obtain a di€erent topology for a ¯exible structure. 7.6. Design of ¯exible structure implementing two ¯exibilities Here we examine multi-¯exibility implementation in a ¯exible structure using a simple model. Fig. 29 shows the design domain where boundary conditions and speci®cations are indicated. As shown in Fig. 29, top and bottom boundaries of the design domain are ®xed, and two ¯exibilities are implemented in this design domain using the multi-objective function formulated by Eq. (5.30) in the multi¯exibility case. That is, while the sti€nesses at points 1 P1 ; 1 P2 ; 2 P1 ; and 2 P2 are to be maximized, two deformations are to be maximized: the deformation at point 1 P2 in the direction of dummy load 1 F 2 is to be maximized when the external force, 1 F 1 , is applied at point 1 P1 , and the deformation at point 2 P2 in the direction of dummy load 2 F 2 is to be maximized when the external force, 2 F 1 ; is applied at point 2 P1 . Using the multi-objective function formulated by Eq. (5.30) in the multi-¯exibility case, the mutual mean compliances de®ned by two sets of two loads, 1 F 1 and 1 F 2 ; 2 F 1 and 2 F 2 , are to be maximized for the kinematic requirement, while the four mean compliances at points 1 P1 ; 1 P2 ; 2 P1 ; and 2 P2 are to be minimized for the structural requirement, where loads 1 F 1 ; 1 F 2 ; 2 F 1 ; and 2 F 2 correspond to tractions 1 t1 ; 1 t2 ; 2 t1 ; and 2 t2 in Section 4.3, respectively. The design domain is discretized using 1600 ®nite elements. The initial values of microscopic design variables ai and bi are set to 0.90 and the value of hi is set to 0.0 for i ˆ 1; . . . ; the number of elements (in all elements). The total volume constraint Xs is set to 480 which corresponds to 30% of the volume of the whole design domain. The ®ltering scheme is also used in this problem. First, we examine the relation between the optimal solutions and parameter nL ; while nS is set to 5.0. Table 3 shows the values of the two mutual mean compliances and the four mean compliances at the optimal point with respect to parameter nL :

4496

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 29. Design domain for a ¯exible structure in the multi-¯exibility case.

Table 3 Values of mutual mean compliances and mean compliances at the optimal point with respect to parameter nL …nS ˆ 5:0† Parameter nL (a) Initial solution

Mutual mean compliance

Mean compliance

1 2 1 1

2 2 2 1

1 3 1 3

1 4 1 4

2 3 2 3

2 4 2 4

1.034

)1.034

1.000

1.000

0.561

0.561

3.431 3.066 3.213

0.419 0.267 0.264

0.326 0.269 0.257

0.442 0.323 0.330

0.478 0.417 0.411

L…u†

(b) Optimal solutions 1.0 4.027 3.0 3.054 5.0 3.078

L…u†

L…u†

L…u†

L…u†

L…u†

It is clear that the values of the two mutual compliances become positive at any choice of nL even though all the initial values are not positive. That is, an extension of the KS-function provides the appropriate optimal solution in the Pareto optima. Comparing cases of nL ˆ 3:0 and 5.0 with 1.0, the di€erence between the two mutual mean compliances is smaller. Thus, by setting nL to an appropriate larger value, we obtain almost the same mutual mean compliances. Fig. 30 shows the optimal con®gurations in the cases of nL ˆ 1:0 and nL ˆ 5:0, where nS is set to 5.0. It is clear that the topologies of the optimal con®gurations are almost the same even though a small di€erence in shape is observed. We also examine the performance of the optimal topology con®guration of the ¯exible structure using linear FE analysis. Fig. 31 shows the ®nite element model constructed based on the optimal con®guration for the case nL ˆ nS ˆ 5:0 shown in Fig. 30(b). This model is descretized to 2384 elements with 2791 nodes. Fig. 32 shows deformed shape of the ¯exible structure. It is also clear that the ¯exible structure deforms according to the problem speci®cations. That is, the ¯exible structure deforms with the speci®ed two ¯exibilities. 7.7. E€ect of relative deformation constraint in the three-dimensional case The e€ect of a relative deformation constraint on the three-dimensional design is con®rmed in this section. That is, by designing a three-dimensional simple ¯exible structure, we examine the speci®cation of the deformation using a relative deformation constraint at either the portion where we need to have the

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4497

Fig. 30. Optimal con®gurations in the multi-¯exibility case (using the multi-objective function de®ned by Eq. (5.30), Xs ˆ 480 (30%), nS ˆ 5:0). (a) nL ˆ 1:0, (b) nL ˆ 5:0.

deformation for the kinematic requirement, corresponding to boundary Ct2 in Fig. 5 or the portion where we apply the traction, corresponding to boundary Ct1 in Fig. 6. Fig. 33 shows the design domain speci®ed as a 24  12  8 rectangular box with a ®xed support boundary on the left-hand side. The design domain is discretized using 2304 (24  12  8) ®nite elements. The initial values of microscopic design variables ai ; bi ; and ci are set to 0.90, and values of /i ; wi ; and hi are set to 0.0 for i ˆ 1; . . . ; the number of elements (in all elements). The total volume constraint Xs is set to 6192.2, which corresponds to 30% of the volume of the whole design domain. The ®ltering scheme is not used, since we do not have severe checkerboard patterns. First, we examine the unconstrained single ¯exibility case. Using the multi-objective function formulated by Eq. (5.17) in the unconstrained single ¯exibility case, the deformation at point P2 in the direction of dummy load F 2 is to be maximized when the external force, F 1 ; is applied at point P1 ; while the sti€nesses at both points P1 and P2 are to maximized. That is, the mutual mean compliance de®ned by the two loads, F 1 and F 2 ; is to be maximized for the kinematic requirement, while the two mean compliances at both points P1 and P2 are to be minimized for the structural requirements. Loads F 1 and F 2 corresponds to tractions t1 and t2 in Section 4.1, respectively. The weighting coecient, ws , in Eq. (5.17) is set to 0.4. Fig. 34 shows the

4498

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 31. Finite element model of the multi-¯exibility structure for linear analysis.

optimal con®guration in the unconstrained single ¯exibility case. In this ®gure and the ®gures shown later, the elements which have a density value of greater than or equal to 0.5 are indicated. Next, we examine the constrained single ¯exibility case in which relative deformation constraints are speci®ed at point P2 . That is, the ¯exible structure is designed to have only the speci®ed translational deformation in the portion where we need to have the deformation for the kinematic requirement, corresponding to boundary Ct2 . Using the multi-objective function formulated by Eq. (5.22) in the constrained single ¯exibility case, the mutual mean compliance de®ned by the two loads, F 1 and F 2 , is to be maximized while the mean compliances at both points P1 and P2 , and the absolute values of the two mutual mean compliances de®ned by two sets of the two loads, F 1 and F 5 ; F 1 and F 6 ; are to be minimized in order to satisfy the required function. Dummy loads F 5 and F 6 correspond to tractions t5 and t6 in Section 4.2, respectively, and F 2 ; F 5 ; and F 6 ; and perpendicular. The weighting coecients, W ; ws ; wM1; and wM2 in Eq. (5.22) are set to 0.4, 0.4, 1.0, and 1.0, respectively. Fig. 35 shows the optimal con®guration in the constrained single ¯exibility case. Comparing this ®gure with Fig. 34, the topology of the optimal con®gurations are quite di€erent. That is, by specifying the translational displacement using a relative deformation constraint, we obtain a quite di€erent topology of the design. Finally, we examine the constrained single ¯exibility case in which relative deformation constraints are speci®ed at point P1 : That is, the ¯exible structure is designed to have only the speci®ed translational deformation in the portion where we apply the force, corresponding to boundary Ct1 . Using the multi-objective function formulated by Eq. (5.24) in the constrained single ¯exibility case, the mutual mean compliance de®ned by the two loads, F 1 and F 2 ; is to be maximized while the mean compliances at both points P1 and P2 , and the absolute values of the two mutual mean compliances de®ned by two sets of the two loads, F 1 and F 7 ; F 1 and F 8 ; are to be minimized in order to satisfy the required function. Dummy loads F 7 and F 8 correspond to tractions t7 and t8 in Section 4.2, respectively, and F 1 ; F 7 ; and F 8 are perpendicular. The weighting coecients, W ; Ws ; wM1 , and wM2 in Eq. (5.24) are also set to 0.4, 0.4 , 1.0, and 1.0, respectively. Fig. 36 shows the optimal con®guration in this constrained single ¯exibility case. It is clear that the topology of the optimal con®guration shown in Fig. 36 is also quite di€erent from that shown in Fig. 34. 7.8. Design of three-dimensional compliant gripper We examine here the design of a three-dimensional compliant gripper. Fig. 37 shows a quarter view of the three-dimensional design domain where boundary conditions and speci®cations are as indicated. The left-hand side boundary is ®xed to support the compliant gripper, and the symmetry boundary condition is posed at both the bottom boundary and the back boundary.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4499

Fig. 32. Deformed shape of the multi-¯exibility structure. (a) An external force of 0.1 unit applied at point 1 P1 , (b) an external force of 0.1 unit applied at point 2 P1 .

The design domain is discretized using 4608 ®nite elements. The initial values of microscopic design variables, ai ; bi ; and ci are set to 0.95, and values of /i ; wi ; and hi are set to 0.0 for i ˆ 1; . . . ; the number of elements (in all elements). The total volume constraint Xs is set to 256, which corresponds to 20% of the volume of the whole design domain. The ®ltering scheme is not used in order to avoid the ambiguity of the optimal con®gurations. First, the optimal con®guration in the unconstrained single ¯exibility case is obtained using the multiobjective functions formulated by Eq. (5.17). In order to satisfy the function of the compliant gripper discussed in Section 7.3, the mutual mean compliance de®ned by two loads F 1 at point P1 and F 2 at point P2 is to be maximized, while the mean compliances at both points P1 and P2 are to be minimized. Loads F 1 and

4500

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 33. Design domain for a simple three-dimensional ¯exible structure.

Fig. 34. Optimal con®guration of the three-dimensional ¯exible structure (using the multi-objective function in Eq. (5.17) in the unconstrained single ¯exibility case, Xs ˆ 619:2 (30%), ws ˆ 0:4).

Fig. 35. Optimal con®guration of the three-dimensional ¯exible structure (using the multi-objective function in Eq. (5.22) in the constrained single ¯exibility case, Xs ˆ 619:2 (30%), W ˆ 0:4; ws ˆ 0:4; wM1 ˆ 1:0; wM2 ˆ 1:0).

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4501

Fig. 36. Optimal con®guration of the three-dimensional ¯exible structure (using the multi-objective function in Eq. (5.24) in the constrained single ¯exibility case, Xs ˆ 619:2 (30%), W ˆ 0:4; ws ˆ 0:4; wM1 ˆ 1:0; wM2 ˆ 1:0).

Fig. 37. Design domain for a three-dimensional compliant gripper.

Fig. 38. Optimal con®guration of the three-dimensional compliant gripper (using the multi-objective function in Eq. (5.17) in the unconstrained single ¯exibility case, Xs ˆ 256 (20%), ws ˆ 0:4).

4502

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

Fig. 39. Optimal con®guration of the three-dimensional compliant gripper (using the multi-objective function in Eq. (5.22) in the constrained single ¯exibility case, Xs ˆ 256 (20%), W ˆ 0:3; ws ˆ 0:4; wM1 ˆ 1:0; wM2 ˆ 0:0).

F 2 correspond to tractions t1 and t2 in Section 4.1, respectively. The weighting coecient, ws , in Eq. (5.17) is set to 0.4. Fig. 38 shows the optimal con®guration in the unconstrained single ¯exibility case. Next, the optimal con®guration in the constrained single ¯exibility case is obtained. we specify the desired translational deformation with a relative deformation constraint de®ned by the external force F 1 and dummy load F 5 . Note that, due to the symmetry boundary condition, we only need dummy load F 5 to specify the translational deformation. Using the multi-objective function formulated by Eq. (5.22) (but, L6 …u1 † is not considered), the mutual mean compliance de®ned by the two loads, F 1 and F 2 , is to be maximized while the mean compliances at both points P1 and P2 , and the absolute value of the mutual mean compliance de®ned by the two loads, F 1 and F 5 , are to be minimized. Load F 5 corresponds to traction t5 in Section 4.2. The weighting coecients, W ; ws ; wM1 ; and wM2 , in Eq. (5.22) are set to 0.3, 0.4, 1.0, and 0.0, respectively. Fig. 39 shows the optimal con®guration in this constrained single ¯exibility case. Comparing these three-dimensional results with the two-dimensional results in Sections 7.3 and 7.4, the topologies of the optimal con®guration have a similar tendency. That is, speci®cation of the translational deformation provides the more complex topology. 8. Conclusion In this research, a methodology to obtain in the optimal structure design considering ¯exibility has been developed using the homogenization design method. First, ¯exibility was formulated using the mutual energy concept, and this sensitivity with respect to a design variable was also formulated. Next, three types of design optimization problems were categorized as follows: the unconstrained single ¯exibility case, the constrained single ¯exibility case, and the multi-¯exibility case. In each case, the optimization problem was formulated, and an appropriate multi-objective function was introduced in order to obtain an appropriate optimal solution which has physical meaning in the Pareto optima. Using these multi-objective functions, the homogenization method and for the relaxation of the ®xed design domain, and SLP, the optimization algorithm was constructed. Finally, some design examples in three types of design optimization problems were presented to examine the uniqueness and characteristics of the optimal con®guration, the relation between the optimal con®guration and the parameters in the multi-objective functions, the e€ect of the total volume constraints on the optimal con®guration, and the applicability to the compliant mechanism design. Thus, we con®rmed that a structure design considering ¯exibility can be accomplished by the methodology presented here. Acknowledgements The authors are supported partially by grant DMI 9622261 from the National Science Foundation and TOYOTA Central R&D Labs., Inc. They are grateful for this support.

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

4503

References [1] M.M. Kamal, J.A. Wolf, in: Modern Automotive Structural Analysis, Van Nostrand Reinhold, New York, 1981. [2] J. Fukushima, K. Suzuki, N. Kikuchi, Application to car bodies: generalized layout design of three dimensional shells, in: G.I.N. Rozvany (Ed.), in: Optimization of Large Structural Systems I, Kluwer Academic Publishers, Dordrecht, 1993, pp. 177±191. [3] R.J. Yang, A.I. Chahande, Automotive applications of topology optimization, Struct. Optim. 9 (1995) 245±249. [4] L.L. Howell, A. Midha, A loop-closure theory for the analysis and synthesis of compliant mechanisms, J. Mech. Des., Trans. ASME 118 (1996) 121±125. [5] G.K. Ananthasuresh, S. Kota, N. Kikuchi, Strategies for systematic synthesis of compliant MEMS, in: Proceedings of the 1994 ASME Winter Annual Meeting, Chicago, IL, DSC-55-2, 1994, pp. 677±686. [6] R.H. Burns, F.R.E. Crossley, Structural permutations of ¯exible link mechanisms, ASME Paper No. 66-MECH-5, 1966. [7] T.E. Shoup, C.W. McLarnan, A survey of ¯exible link mechanisms having lower pairs, J. Mech. 6 (1971) 97±105. [8] N.M. Sevak, C.W. McLarnan, Optimal synthesis of ¯exible link mechanisms with large static de¯ections, ASME Paper No. 74-DET-83, 1974. [9] I. Her, A. Midha, A compliance number concept for compliant mechanisms, and type synthesis, J. Mech., Transmissions, Automat. Des., Trans. ASME 109 (1987) 348±355. [10] M.D. Murphy, A. Midha, L.L. Howell, The topological synthesis of compliant mechanisms, in: Proceedings of the Third National Conference on Applied Mechanisms and Robotics, vol. 2, Cincinnati, Ohio, 1994, Paper No. 99. [11] L.L. Howell, A. Midha, A method for the design of compliant mechanisms with small-length ¯exural pivots, J. Mech. Des., Trans. ASME 116 (1994) 280±289. [12] O. Sigmund, Some inverse problem in topology design of materials and mechanisms, in: Proceedings of the IUTAM Symposium on Optimization of Mechanical Systems, Stuttgart, Germany, 1995, pp. 26±31. [13] M.I. Frecker, N. Kikuchi, S. Kota, Optimal synthesis of compliant mechanisms to satisfy kinematics and structural requirements ± preliminary results, in: Proceedings of the ASME Design Engineering Technical Conferences and Computers in Engineering Conference, Irvibe, California, 1994, 96-DETC/DAC-1417. [14] O. Sigmund, On the design of compliant mechanisms using topology optimization, Dan. Center Appl. Math. Mech. 535 (1996) 1±28. [15] U.D. Larsen, O. Sigmund, S. Bouswstra, Design and fabrication of compliant micromechanisms and structures with negative Poisson's ratio, in: Proceedings of the IEEE Ninth Annual International Workshop on Micro Electro Mechanical Systems, An Investigation of Micro Structures, Sensors, Actuators, Machines and Systems, San Diego, California, 1996, pp. 365±371. [16] G.K. Ananthasuresh, S. Kota, Y. Gianchandani, Systematic synthesis of microcompliant mechanisms ± preliminary results, in: Proceedings of the Third National Conference on Applied Mechanisms and Robotics, Cincinnati, Ohio, 1993, Paper 82. [17] M.P. Bendsùe, in: Optimization of Structural Topology, Shape, and Material, Springer, Berlin, 1995. [18] O. Sigmund, Materials with prescribed constitutive parameters: an inverse homogenization problem, Int. J. Sloids Struct. 31 (1994) 2313±2329. [19] R.J. Yang, C.H. Chuang, Optimal topology design using linear programming, Comput. Struct. 52 (1994) 265±275. [20] C.C. Swan, I. Kosaka, Voigt±Reuss topology optimization for structures with linear elastic material behaviours, Int. J. Numer. Methods Engrg. 40 (1997) 3033±3057. [21] C.C. Swan, I. Kosaka, Voigt±Reuss topology optimization for structures with nonlinear elastic material behaviours, Int. J. Numer. Methods Engrg. 40 (1997) 3785±3814. [22] M.P. Bendsùe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg. 71 (1988) 197±224. [23] K. Cheng, N. Olho€, An investigation concerning optimal design of solid elastic plates, Int. J. Solids Struct. 17 (1981) 305±323. [24] K.A. Lurie, A.V. Cherkaev, A.V. Fedorov, Regularization of optimal design problems for bars and plates, J. Optim. Theory Appl. 37 (1982) 499±521 (Part 1), 523±543 (Part 2). [25] K.A. Lurie, A.V. Cherkaev, A.V. Fedorov, On the existence of solutions to some problems of optimal design for bars and plates, J. Optim. Theory Appl. 42 (1984) 247±281. [26] R.V. Kohn, G. Strang, Optimal design and Relaxation of Variational Problems, Commun. Pure Appl. Math. 39 (1986) 1±25 (Part 1), 139±182 (Part 2), 353±377 (Part 3). [27] F. Murat, L. Tartar, Optimality conditions and homogenization, in: A. Marino, L. Modica, S. Spagnolo, M. Degiovannni (Eds.), in: Nonlinear Variational Problems, Pitman Publishing Program, Boston, 1985, pp. 1±8. [28] G.I.N. Rozvany, T.G. Ong, W.T. Seto, R. Sandler, N. Olho€, M.P. Bendsùe, Least-weight design of perforated elastic plates, Int. J. Solids Struct. 23 (1987) 521±536 (Part 1), 537±550 (Part 2). [29] K. Suzuki, N. Kikuchi, A homogenization method for shape and topology optimization, Comput. Methods Appl. Mech. Engrg. 93 (1991) 291±318. [30] A. Dõaz, N. Kikuchi, Solutions to shape and topology eigenvalue optimization problems using a homogenization method, Int. J. Numer. Methods Engrg. 35 (1992) 1487±1502. [31] Z.-D. Ma, N. Kikuchi, I. Hagiwara, Structural topology and shape optimization for a frequency response problem, Comput. Mech. 13 (1993) 157±174. [32] Z.-D. Ma, N. Kikuchi, H.-C. Cheng, Topological design for vibrating structures, Comput. Methods Appl. Mech. Engrg. 121 (1995) 259±280. [33] M.M. Neves, H. Rodrigues, J.M. Guedes, Generalized topology design of structures with a buckling load criterion, Struct. Optim. 10 (1995) 71±78.

4504

S. Nishiwaki et al. / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4457±4504

[34] S. Min, N. Kikuchi, Optimal reinforcement design of structures under the buckling loads using the homogenization design method, Struct. Engrg. Mech. 5-5 (1997) 565±576. [35] R.T. Shield, W. Prager, Optimal structural design for given de¯ection, J. Appl. Math. Phys., ZAMP 21 (1970) 513±523. [36] N. Huang, On principle of stationary mutual complementary energy and its application to optimal structural design, J. Appl. Math. Phys., ZAMP 22 (1971) 609±620. [37] J.T. Oden, N. Kikuchi, Theory of variational ineqalities with application to problems of ¯ow through porous media, Int. J. Engrg. Sci. 18 (1980) 1173±1284. [38] M.P. Bendsùe, Optimal shape design as a material distribution problem, Struct. Optim. 1 (1989) 193±202. [39] R.B. Haber, C.S. Jog, M.P. Bendsùe, A new approach to variable-topology shape design using a constraint on perimeter, Struct. Optim. 11 (1996) 1±12. [40] D. Fujii, N. Kikuchi, Improvement of numerical instabilities in topology optimization using the SLP method, Struct. Optimization 19 (2000) 113±121. [41] J. Koski, Multicriteria optimization in structural design: state of the art, in: Proceedings of the ASME Advance in Design Automation, Albuquerque, New Mexico, 1993, DE-vol. 65-1, 1993, pp. 621±629. [42] W. Stadler, Fundamentals of multicriteria optimization, in: W. Stadler (Ed.), in: Multicriteria Optimization in Engineering and in the Sciences, Plenum Press, New York, 1988, pp. 1±25. [43] J. Koski, Defectiveness of weighting method in multicriteria optimization of structures, Commun. Appl. Numer. Methods 1 (1985) 333±337. [44] H.P. Mlejnek, R. Schirrmacher, An engineer's approach to optimal material distribution and shape ®nding, Comput. Methods Appl. Mech. Engrg. 106 (1993) 1±26. [45] J.M. Guedes, N. Kikuchi, Preprocessing and postprocessing for materials based on the homogenization method with adaptive ®nite element methods, Comput. Methods Appl. Mech. Engrg. 83 (1990) 143±198. [46] B.C. Koh, N. Kikuchi, New improved hourglass control for bilinear and trilinear elements in anisotropic linear elasticity, Comput. Methods Appl. Mech. Engrg. 65 (1987) 1±46. [47] Y.Y. Zhu, S. Cescotto, Uni®ed and mixes formulation of the 8-node hexahedral elements by assumed strain method, Comput. Methods Appl. Mech. Engrg. 129 (1996) 177±209. [48] D.S. Malkus, T.J.R. Hughes, Mixed ®nite element methods ± reduced and selective integration technique: a uni®cation of concepts, Comput. Methods Appl. Mech. Engrg. 15 (1978) 63±81. [49] T.J.R. Hughes, Generalization of selective integration procedures to anisotropic and nonlinear media, Int. J. Numer. Methods Engrg. 15 (1980) 1413±1418. [50] W.K. Liu, J.S.-J. Ong, R.A. Uras, Finite element stabilization matrics ± a uni®cation approach, Comput. Methods Appl. Mech. Engrg. 53 (1985) 13±46. [51] J. Yoo, Structural optimization in magnetic ®elds using the homogenization design method, Ph.D. dissertation at the University of Michigan, 1999, pp. 12±46. [52] M. Zhou, G.I.N. Rozvany, The COC algorithm, Part II: Topological geometrical and generalized shape optimization, Comput. Methods Appl. Mech. Engrg. 89 (1991) 309±336. [53] C. Fleury, V. Braibant, Structural optimization: a new dual method using mixed variables, Int. J. Numer. Methods Engrg. 23 (1986) 409±428. [54] K. Svanberg, The method of moving asymptotes ± a new method for structural optimization, Int. J. Numer. Methods Engrg. 24 (1987) 359±373. [55] L.H. Tenek, I. Hagiwara, Static and vibrational shape and topology optimization using homogenization and mathematical programming, Comput. Methods Appl. Mech. Engrg. 109 (1993) 143±154. [56] P. Duysinx, Layout optimization: a mathematical programming approach, Dan. Center Appl. Math. Mech. 540 (1997) 1±30. [57] R.T. Haftka, Z. G urdal, Elements of Structural Optimization, third revised and expanded ed., Kluwer Academic Publishers, Dordrecht, 1992. [58] H.L. Thomas, G.N. Vanderplaats, Y.-K. Shyy, A study of move limit adjustment strategies in the approximation concepts approach to structural synthesis, in: Proceedings of tne Fourth AIAA/USAF/NASA/OAI Symposium on Multidisciplinary Analysis and Optimization, Cleveland, Ohio, 1992, pp. 507±512. [59] T.-Y. Chen, Calculation of the move limits for the sequential linear programming method, Int. J. Numer. Methods Engrg. 36 (1993) 2661±2679. [60] U. Kirsch, E€ective move limits for approximate structural optimization, J. Struct. Engrg. 123 (1997) 210±217. [61] R. Hanson, K. Hiebert, A sparse linear programming, Technical Report SAND81-0297, Sandia National Laboratories, 1981. [62] C.S. Jog, R.B. Haber, M.P. Bendsùe, Topology design with optimized self-adaptive materials, Int. J. Numer. Methods Engrg. 37 (1994) 1323±1350. [63] A. Diaz, O. Sigmund, Checkerboard patterns in layout optimization, Struct. Optim. 10 (1995) 40±45. [64] N. Kikuchi, S. Hollister, J. Yoo, A concept of image-based integrated CAE for production engineering, in: Proceedings of the International Symposium on Optimization and Innovative Design in JSME, Tokyo, Japan, 1997, pp. 75±90. [65] P. Pedersen, On optimal orientation of orthotropic materials, Struct. Optim. 1 (1989) 101±106.