Composite Material Design of Two-Dimensional Structures Using The Homogenization Design Method

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2001; 50:20312051

Composite material design of two-dimensional
structures using the homogenization design method
D. Fujii
;
, B. C. Chen and N. Kikuchi
Department of Mechanical Engineering and Applied Mechanics; University of Michigan; Ann Arbor;
Michigan 48109-2125; U.S.A.
SUMMARY
Composite materials of two-dimensional structures are designed using the homogenization design method.
The composite material is made of two or three dierent material phases. Designing the composite material
consists of nding a distribution of material phases that minimizes the mean compliance of the macrostructure
subject to volume fraction constraints of the constituent phases, within a unit cell of periodic microstructures.
At the start of the computational solution, the material distribution of the microstructure is represented as
a pure mixture of the constituent phases. As the iteration procedure unfolds, the component phases sepa-
rate themselves out to form distinctive interfaces. The eective material properties of the articially mixed
materials are dened by the interpolation of the constituents. The optimization problem is solved using the
sequential linear programming method. Both the macrostructure and the microstructures are analysed using
the nite element method in each iteration step. Several examples of optimal topology design of composite
material are presented to demonstrate the validity of the present numerical algorithm. Copyright ? 2001 John
Wiley & Sons, Ltd.
KEY WORDS: composite material; homogenization method; topology optimization; microstructures; numerical
algorithm
1. INTRODUCTION
The homogenization design method (HDM) has been proposed by Bendse and Kikuchi [1] and
Suzuki and Kikuchi [2] in order to obtain the optimum shape and topology of the continuum
structure. In this method, the boundary shape optimization problem is converted into the material
optimization problem in the extended design domain based on the approach of Murat and Tartar
[3]. It is assumed that the extended design domain is composed of porous materials, and the
optimum topology and shape of the macrostructure are obtained by optimizing the size and the
direction of the material holes in the microstructure (see Figure 1).
Correspondence to: D. Fujii, Department of Environmental and Ocean Engineering, The University of Tokyo. 7-3-1
Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
E-mail: dfujii@nasl.t.u-tokyo.ac.jp
Contract=grant Sponsor. Steel Material Club Corporation Aggregate, Japan
Received 4 February 1999
Copyright
?
2001 John Wiley & Sons, Ltd. Revised 22 September 1999
2032 D. FUJII, B. C. CHEN AND N. KIKUCHI
Figure 1. A two-dimensional structure with porous microstructures.
In this work, there is no need to limit the material microstructure to a particular porous pattern,
and it can be of any arbitrary shape and topology. In addition, it is also possible to design the
material microstructure which realizes the optimum components of the elasticity tensor of the
macrostructure, if the optimum components of that are already known. Several researchers have
addressed the problem of designing the microstructure of the composite material which realizes
the optimum elasticity constants of the macrostructure based on this idea. Sigmund et al. [47]
and Fonseca [8] formulated this problem as an optimization problem by searching for the lightest
microstructure with the prescribed thermoelastic properties specied as equality constraints. Swan
and Arora [8] and Swan and Kosaka [10] designed inelastic composites to enhance stiness and
strength. Silva et al. [11; 12] and Gibiansky and Torquato [13; 14] designed the microstructure
of piezocomposites with improved performance characteristics. Haslinger and Dvorak [15] treated
this problem as an identication problem.
However, the loading and the boundary conditions of the macrostructure change the optimum
characteristics of the macrostructure actually. In some cases, it is dicult to determine the optimum
characteristics of the macrostructure. The simultaneous analysis of the macrostructure and the
microstructure based on the original idea of Bendse and Kikuchi [1] is necessary for such cases.
Therefore, in this paper, it is assumed that the design domain is composed of the composite
material with two or three material phases instead of the porous material, and the material layout
of unit cell of the microstructure is chosen as the design object. In other words, in the present
method, the topology of the unit cell is free to change to any pattern, while the microstructure
pattern is parameterized as a square region with a sizeable rectangular hole in the original idea of
Bendse and Kikuchi [1]. As an initial attempt, the simple problem in which the macrostructure
is uniform is handled in this paper because the number of the design variables becomes enormous
if the microstructures assume non-uniform distribution within the macroscopic design domain.
The present design procedure consists of nding a distribution of material phases that minimizes
the mean compliance of the macrostructure subject to volume fraction constraints of the constituent
phases, within a unit cell of periodic microstructures. Concretely, the material of each element
of the unit cell with nite element discretization is assigned from two or three material phases.
However, in this discretion optimization, the number of the combination becomes enormous, as the
number of division of the nite element increases, and it is hopeless to solve the optimization prob-
lem using random search methods such as, genetic algorithms or simulated annealing methods [7].
Copyright ? 2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 50:20312051
COMPOSITE MATERIAL DESIGN OF TWO-DIMENSIONAL STRUCTURES 2033
Therefore, in this paper, the mixing rule used in the analysis of macrostructure (e.g. Bendse [16],
Rozvany et al. [17], Mlejnek and Schirrmacher [18], Yang and Chuang [19]), is applied to the
microstructure design as one can nd in literature [415].
Although the mixing rules for two materials mixture involving linear elastic solid and void phases
have been proposed, the mixing rule for three materials mixture is not found except for Sigmund
and Torquato [7]. The three-material-phase design is more advantageous than the two-material-
phase design when the temperature characteristics, permeability, lightening, etc. are considered.
The mixing rule for two materials and void mixture proposed by Sigmund and Torquato [7] is
composed of the classical Voigt mixing rule for two solid materials and the power law for void
and solid ones. However, if we use the Voigt mixing rule, the mixture of the two material phases
dose not separate in the nal stage of the iteration procedure for the optimization, as Swan and
Kosaka [20] has reported. Therefore, the new mixing rule of the three materials mixture based on
the power laws is proposed in this paper.
It is well known that the checkerboarding material distribution is easily obtained if the underlying
mixing rule based on density approach is employed [20; 21]. In order to avoid the checkerboarding
solutions, it is necessary to use the high-order nite element methods [7; 22; 23], or to use ltering
methods [20; 2427]. In this paper, the low-order nite element method for the analysis of the
microstructure is used in order to save calculation time, and a novel ltering method, based on
the concept of gravity, is proposed.
The optimization problem is solved using sequential linear programming. The method for
relaxing the volume fraction constraints of the constituent phases in initial stage of the itera-
tions is used in order to avoid the convergence to the local minima. Some examples of optimum
topology design of composite materials for two-dimensional structures are presented to demonstrate
the validity of the proposed ltering method and the present numerical algorithm.
The paper is organized in the following way. In Section 2 we describe the formulation of the
optimization problem and its discretization. The topology optimization procedure is described in
Section 3. In Section 4, the checkerboarding problem is discussed, and a new ltering method
is presented. The performance of the proposed ltering method is demonstrated with an example.
In Section 5, the eectiveness of the present algorithm is shown by some examples which have
dierent boundary and loading conditions.
2. FORMULATION OF THE OPTIMIZATION PROBLEM
2.1. Optimization problem in the homogenization design method
In the homogenization design method (HDM) proposed by Bendse and Kikuchi [1], it is assumed
that the extended domain is composed of the porous material, and the size and the direction of
the holes of the microstructure which minimize the mean compliance of the macrostructure are
obtained. This optimization problem can be expanded in the more general form as follows:
max
{D
H
ij
}
_
min
v
0
K
(v
0
)
_
; K={v
0
| v
0
=0 on
D
} (1)
(v
0
) =
1
2
_
D
U(v
0
)
T
D
H
U(v
0
) dD
_
T
v
T
0
t d (2)
Figure 2. A two-dimensional structure with composite material microstructures.
where is the total potential energy, v
0
the displacement vector of the macrostructure, U the strain
tensor, D
H
the homogenized elasticity matrix, D
H
ij
the components of D
H
, and t the boundary
traction vector (see Figure 1). Here, it is assumed that the body force can be neglected.
The optimization problem in Equation (1) searches the optimal components D
H
ij
which maximize
the minimum potential energy of the macrostructure. In the problem in which loading condition is
given, Equation (1) becomes the minimization problem of the compliance. On the other hand, for
the problem in which enforcement displacement is given, it becomes the maximization problem of
the compliance. In the problem in search of the topology of the macrostructure, D
H
becomes the
function of the size and direction of the hole in the porous microstructure, namely D
H
(a; b; ).
In this paper, the microstructure of the extended design domain D is composed of the uniform
periodic arrangement of the unit cell as shown in Figure 2, and the unit cell is made by the layout
of 2 or 3 kinds of material phases. In this problem, D
H
is decided by the layout of material phases
in the unit cell. Therefore, the layout of material phases in the unit cell is the design object of
this problem.
2.2. Homogenized minimum potential energy
The homogenized elasticity matrix D
H
is obtained by the homogenization method [28]. In this
method, two kinds of co-ordinate systems, the global macro-co-ordinates x and the local micro-co-
ordinates y, are used, as shown in Figure 2. The relationship y =x= is established between the two
co-ordinate systems, where is the representative length of the microstructure, and is assumed
to be small enough compared with the dimension of the macrostructure. And the microstructure
is supposed to have the periodicity. In this case, the displacement eld can be approximated by
an asymptotic expansion
v
(x) =v(x; y) v
0
(x) + v
1
(x; y) (3)
where the superscript indicates the existence of the microstructure, v
0
is the average displacement
of the unit cell, v
1
is the rst-order variation from the average displacement.
The strain and stress tensor also can be approximated by
U(v
) =999
x
v
999
x
v
0
+ 999
y
v
1
=U
0
+ U
1
(4)
where U
0
is the strain due to the average displacement, U
1
is the strain due to the rst deviation
of the displacement from the average displacement of the unit cell, and 999
x
; 999
y
are the dierential
operators which are expressed in the two dimensional problem as follows:
999
x
=
_
_
@=@x
1
0
0 @=@x
2
@=@x
2
@=@x
1
_
_
; 999
y
=
_
_
@=@y
1
0
0 @=@y
2
@=@y
2
@=@y
1
_
_
(5)
The stress is dened by
A(v
) =DU(v
) D(999
x
v
0
+ 999
y
v
1
) (6)
where D is the elasticity matrix.
Because v
1
is the rst-order variation from the average displacement, it is possible that v
1
is
considered to be proportional for the strain U
0
:
v
1
= E(y)U
0
(x) = E(y)999
x
v
0
(x) (7)
where E is the 2 3 matrix in the two-dimensional problem, and it indicates a kind of proportion-
ality constant. E is called the characteristic displacement function of the microstructure, because it
is the displacement for the unit strain of the macrostructure.
Substituting Equation (7) into Equations (4) and (6), the total potential energy is expressed as
(v
) =
1
2
_
D
(999
x
v
0
)
T
(I 999
y
E)
T
D(I 999
y
E)(999
x
v
0
) dD
T
v
T
0
t d (8)
The following equation
lim
0
_
D
(x) dD
=
_
D
_
1
|Y|
_
Y
(x; y) dY
_
dD (9)
is considered, and by 0
(v
) =
1
2
_
D
(999
x
v
0
)
T
1
|Y|
_
Y
(I 999
y
E)
T
D(I 999
y
E) dY(999
x
v
0
) dD
_
T
v
T
0
t d (10)
is obtained where I is the identity matrix, and |Y| the area of unit cell. We dene
D
Y
=
1
|Y|
_
Y
(I 999
y
E)
T
D(I 999
y
E) dY (11)
then Equation (10) is expressed as follows:
(v
) =
1
2
_
D
(999
x
v
0
)
T
D
Y
(999
x
v
0
) dD
_
T
v
T
0
t d (12)
From Equation (12), the minimum potential energy is obtained as follows:
min
v
[(v
)] min
v
0
;E
[(v
0
)] =
H
(u
0
) = min
v
0
_
1
2
_
D
(999
x
v
0
)
T
D
H
(999
x
v
0
) dD
_
T
v
T
0
t d
_
(13)
where
H
(u
0
) is the homogenized minimum potential energy, and D
H
is the homogenized elasticity
matrix dened as follows:
D
H
= min
E
Periodic
D
Y
= min
E
Periodic
_
1
|Y|
_
Y
(I 999
y
E)
T
D(I 999
y
E) dY
_
(14)
Here, (v
0
) in Equation (13) is equal to Equation (2).
2.3. Discretization for the macrostructure
In order to solve the optimization problem dened by Equation (1), the minimum potential energy
dened by Equation (13) is discretized using the nite element method.
The displacement eld in an element is expressed using the nodal displacements of the element
as follows:
v
0
=N
0
d
L
e
(15)
where N
0
is the shape function, and d
L
e
is the nodal displacement vector with respect to the local
co-ordinate system. From Equation (15),
U
0
=999
x
v
0
=B
0
d
L
e
(16)
is obtained, where B
0
is the straindisplacement matrix. And the potential energy for an element
is as follows:
e
=
1
2
d
L
T
e
_
D
e
B
T
0
D
H
B
0
dD
e
d
L
e
d
L
T
e
_
Te
N
T
0
t d
e
(17)
The displacement vector d
L
e
is transformed to d
e
with respect to the global co-ordinate system
by
d
L
e
=T
e
d
e
(18)
where T
e
is the coordinate transformation matrix. Substitute Equation (18) into Equation (17),
then
e
=
1
2
d
T
e
T
T
e
_
D
e
B
T
0
D
H
B
0
dD
e
T
e
d
e
d
e
T
T
e
_
Te
N
T
0
t d
e
(19)
When the total potential energy is minimum,
e
=0 is consistent, and
K
0e
d
e
=f
e
(20)
where
K
0e
= T
T
e
_
D
e
B
T
0
D
H
B
0
dD
e
T
e
(21)
f
e
= T
T
e
_
Te
N
T
0
t d
e
(22)
Figure 3. Discretization of a unit cell.
The following global stiness equation is obtained by assembling the element stiness equation
dened by Equation (20) for all elements.
K
0
d =f (23)
In the problem in which the loading condition is prescribed, Equation (23) is solved considering
the geometrical boundary condition, and using this solution, the minimum potential energy is
obtained as follows:
H
(u
0
) =
1
2
d
T
K
0
d d
T
f =
1
2
d
T
K
0
d (24)
Therefore, Equation (1) can be written as follows:
min
{D
H
ij
}
[d
T
K
0
d] (25)
2.4. Discretization for the microstructure
The homogenized elasticity matrix D
H
dened by Equation (14) is also discretized by the nite
element method. The unit cell is divided by the nite elements as shown in Figure 3. The charac-
teristic displacement function E and its strain 999
y
E in an element is approximated using the nodal
displacement function X
e
of the element as follows
E = N
1
X
e
(26)
999
y
E = B
1
X
e
(27)
where N
1
and B
1
are the shape function and the straindisplacement matrix for the microstructure,
respectively. Using Equations (26) and (27), D
Y
dened by Equation (11) is expressed as follows:
D
Y
=
1
|Y
e
|
_
Y
e
(I B
1
X
e
)
T
D(I B
1
X
e
) dY
e
(28)
The necessary condition for the minimizing of Equation (28) is
X
e
__
Y
e
B
1
DdY
e

_
Y
e
B
1
DB
1
dY
e
X
e
_
=0 (29)
From Equation (29), the following equation is obtained:
K
1e
X
e
=q
e
(30)
where
K
1e
=
_
Y
e
B
1
DB
1
dY
e
(31)
q
e
=
_
Y
e
B
1
DdY
e
(32)
The following equation is obtained by assembling Equation (30) for all elements in the unit
cell.
K
1
X=q (33)
Equation (33) is solved considering the periodicity condition of the unit cell, and using this
solution, D
H
dened by Equation (14) is calculated by the following equation:
D
H
=

elements
_
Y
e
(I B
1
X
e
)
T
D(I B
1
X
e
) dY
e
(34)
3. TOPOLOGY OPTIMIZATION PROCEDURE
The optimization problem dened by Equations (1) and (2) in the discretized form is stated in
Equations (25) and (34). The design object {D
H
ij
} in Equation (25) is calculated by Equation (34).
In this paper, the unit cell is composed of two- or three-phase material structures (see Figure 4).
Therefore, the optimization problem consists of assigning the material phase to each element in
the unit cell from two or three material phases. In other words, it is to nd the topology of each
material phase in the unit cell. Here, it is assumed that the material properties in each phase are
linear elastic and isotropic.
3.1. Constitutive mixing rules
In this paper, the density approach [1719] proposed for the topology optimization analysis of the
macrostructure is adopted for the analysis of the microstructure. In the density approach, the solid
material and void are mixed articially using the power law. This mixing rule also can be applied
to two kinds of solid materials.
It is known that if one uses the Voigt mixing rule (the parallel spring model, see Figure 5),
the mixed material is not decomposed in the nal solution, that is to say, intermediate density
(gray scale) is remained in the topology solution. On the other hand, if we use the Reuss mixing
rule (the serial spring model, see Figure 5), the mixed material is decomposed clearly. However,
in the case of the mixture of solid material and void, the variation of the stiness of the mixed
material is extremely steep (it will be shown in Figure 6(b) later). The mixing rule with power
law corresponds to the intermediate one between Voigt and Reuss mixing rules. As an alternative,
the hybrid rule of the Voigt and Reuss mixing rules was proposed by Swan and Kosaka [20].
Though the hybrid rule is more eective for the elastoplastic problem, the power law is adopted
in this paper because this rule is simpler than the hybrid rule for the three material phase problem.
Figure 4. Design domain for topology optimization of composite material microstructures. Each square
represents one nite element which consists of either phase A; B or C.
Figure 5. Voigt mixing rule and Reuss mixing rule, where E
A
and E
B
are the stiness of A and B materials.
So far, we have addressed the mixing rule for two-phase material. However, in the design
of composite material, three-phase material design is more advantageous when the temperature
characteristics, permeabilities, lightening, ecology, etc. are considered. Sigmund and Torquato [7]
addressed the design of materials using three kinds of materials, viz. two solid materials and void.
They used the mixing rule which combines the power law and the Voigt mixing. However, as it
has been mentioned above, it is easy to have overlapping between the constituent phases in the
nal solution when the Voigt mixing rule is used. Therefore, in this paper, instead of using the
material rule based on the Voigt mixing, a new mixing rule, which used the combines power laws
for three-phase material design is proposed.
Figure 6. Relationship between E and a for the variation of s in the case of: (a) E
A
=72:4 GPa; E
B
=3:0 GPa
and b =1; (b) E
A
=3:0 Gpa; E
B
=100 Pa and b =1.
In considering the microstructures made of constituent phases AC, the Youngs moduli for
each constituent phases are denoted by E
A
; E
B
; E
C
. In this case, the Youngs modulus E of the
mixed material is dened using the power law as follows:
E =r{pE
A
+ (1 p)E
B
} + (1 r){qE
C
+ (1 q)E
B
} (35)
where p; q and r are dened by
p = a
s
if E
B
6E
A
(36)
p = 1 (1 a)
s
if E
B
E
A
(37)
q = a
s
if E
B
6E
C
(38)
q = 1 (1 a)
s
if E
B
E
C
(39)
r = b
s
if {qE
C
+ (1 q)E
B
} 6{pE
A
+ (1 p)E
B
} (40)
r = 1 (1 b)
s
if {qE
C
+ (1 q)E
B
}{pE
A
+ (1 p)E
B
} (41)
and 0 6a 61 and 0 6b 61. In addition, the power s in Equations (36)(41) is the factor used
to penalize the intermediate properties. For a special case when b =1, the problem degenerates
into a two-phase design composed of phases A and B. As one can observe from above equations,
if a =1; b =1 then E =E
A
; if a =0 then E =E
B
; and if a =1; b =0 then E =E
C
.
Figures 6(a) and 6(b) show the relationships between E and a (b =1) for dierent sets of s.
Figure 6(a) shows the relationships for the case when the two constituent phases are of comparable
rigidity, E
A
=72:4 GPa; E
B
=3:0 GPa. Figure 6(b) shows the relationships for the case when phase
B is almost negligible when compared with phase A, E
A
=3:0 GPa; E
B
=100 Pa. In these gures,
the Voigt and Reuss mixing rules are also shown respectively,
E =aE
A
+ (1 a)E
B
(42)
E =
E
A
E
B
(1 a)E
A
+ aE
B
(43)
where the Voigt rule corresponds to the case of s =1.
As it can be seen from these gures, the Reuss mixing rule is eective for penalizing the inter-
mediate value of a when the Youngs moduli are comparable in magnitude. It is not appropriate,
however, when the magnitudes of constituent phases are of dierent scales. In contrast, the power
function dened in Equations (36)(41) is eective for both cases.
On the other hand, the Voigt mixing rule is applied for the Poissons ratio as follows
=b{a
A
+ (1 a)
B
} + (1 b){a
C
+ (1 a)
B
} (44)
where
A
;
B
;
C
are the Poissons ratios of constituent phases A, B and C, respectively.
3.2. Formulation of the topology optimization problem
By using the aforementioned mixing rule, it is possible to analyse simultaneously, the topology
of the macrostructure and the topology of the microstructure in theory, because this optimization
problem is the extension of the method proposed by Bendse and Kikuchi [1]. However, it is
impractical, since the number of the design variables becomes enormous. Therefore, in this paper,
it is assumed that the structure consists of a single uniform composite material as shown in
Figure 4. In this case, the number of design variables equals the number of nite elements of the
unit cell for two-material-phase design, and the twice of that for three-material-phase design.
The design variables which prescribe the material phase of the i th element are denoted by
a
i
; b
i
, where i =1; 2; : : : ; N, and N denotes the number of the nite elements of the unit cell.
The constraint conditions of this optimization problem are the volume fractions of the constituent
materials:
f
A
=
N
i=1
a
i
b
i
N
; f
B
=
N
i=1
(1 a
i
)
N
; f
C
=1 f
A
f
B
(45)
where f
A
; f
B
; f
C
are the volume fractions of constituent phases AC. Here, it is assumed that the
area of all nite elements of the unit cell is equal. In this case, the optimization problem dened
by Equation (25) is rewritten as follows:
min
Q L
[C =d
T
K
0
d]; Q ={a
1
; a
2
; : : : ; a
N
; b
1
; b
2
; : : : ; b
N
}
L = {Q|f
min
A
6f
A
6f
max
A
; f
min
B
6f
B
6f
max
B
; 06a
i
61; 06b
i
61; i =1; : : : ; N}
(46)
where C denotes the mean compliance, f
max
; f
min
are the lower and upper limits of the volume
fractions.
3.3. Sequential linear programming method
In this paper, the sequential linear programming method is used to solve the optimization problem
dened by Equation (46). Figure 7 shows the ow chart of the optimization procedure. To begin
with, the initial values of design variables are given. Next, the characteristic displacement function
X is calculated from Equation (33), and D
H
is calculated by Equation (34). The mean compliance
C is obtained by solving Equation (23).
Next, it is necessary to calculate the sensitivities of the compliance and the volume fractions.
The sensitivity of the mean compliance is calculated by the following equation:
@C
@
i
=d
T
@K
0
@
i
d =

elements
d
L
T
e
_
D
e
B
T
0
@D
H
@
i
B
0
dD
e
d
L
e
(47)
where
i
is the ith component of Q, other variables are already dened in Section 2, and the
sensitivity of D
H
is calculated from Equation (34) as follows:
@D
H
@
i
=
_
Y
e
(I B
1
X
e
)
T
@D
@
i
(I B
1
X
e
) dY
e
(48)
The simplex method is used for solving the linear programming method, and the optimum solution
of the increment value of Q, viz. Q, is obtained. Such calculation is repeated until all components
of Q are less than
tol
.
The moving limit is imposed in Q as follows:
6
i
6 (i =1; : : : ; 2N) (49)
In this paper, the initial value of the moving limit is set to 0.5, and we change the value by
=1:05
k
(k denotes the times of iteration) for k650, and =1:05
50
=1:1
k50
for k 50. And, the
tolerance
tol
is set to 0.005. In this case, the maximum number of iterations for optimization
is 68.
The constraints on the volume fraction of the constituent phases are applied in a step-wise
fashion. It is dicult to nd the global optimal solution if the volume fractions are set to the
pre-dened bounds at the start of the optimization iteration. Therefore, for the relaxation we set
these limit values as
f
min
A
=

f
A
_
1 {1 (k 1)=n
c
}
f
max
A
=

f
A
_
1 + {1 (k 1)=n
c
}
; (k6n
c
+ 1)
f
min
B
=

f
B
_
1 {1 (k 1)=n
c
}
f
max
B
=

f
B
_
1 + {1 (k 1)=n
c
}
; (k6n
c
+ 1)
f
min
A
= f
max
A
=

f
A
; (k n
c
+ 1)
f
min
B
= f
max
B
=

f
B
; (k n
c
+ 1)
= min
_
1=

f
A
1; 1=

f
B
1; 2
(50)
where

f
A
;

f
B
;

f
c
are the volume fractions specied, n
c
denotes an appropriate number and 16n
c
6n
(n is the total number of iterations). In the numerical examples in this paper, we set n
c
to 30.
4. CHECKERBOARDING PROBLEM
4.1. Checkerboarding solutions
In this paper, the four-node nite element proposed by Wilson and Tayler [29; 30] is used for
the analysis of two-dimensional macrostructure, and the four-node isoparametric nite element is
used for the analysis of microstructure. The checkerboarding solutions are obtained in the present
approach, because the density approach with low order nite element is used for the analysis of
the microstructure.
A numerical example is solved by the present approach in order to demonstrate this phenomenon.
Figure 8 shows the example of a plate model with shear loads. The boundary condition, the loading
condition and the nite element division of this model are shown in this gure. The thickness t
Figure 7. Flowchart of the optimization
procedure.
Figure 8. A plate model subjected to
shear forces.
Table I. The properties of material phases.
E (GPa)
Grey material Cast epoxy resin 3.0 0.25
Black material E-glass bre 72.4 0.15
White material Void 10
7
0.25
of the plate is 0.01 m. This problem is solved under the plane stress condition. In this analysis,
the power s in Equations (36)(41) is set to 6 for the two-phase design and 5 for the three-phase
design. The symmetrical condition with respect to the two axes y
1
; y
2
and the diagonal line of the
unit cell is imposed. The volume fractions of the constituent phases are set to

f
A
=0:5 and

f
B
=0:5
for two-phase materials, and

f
A
=0:4;

f
B
=0:4 and

f
C
=0:2 (void) for three-phase materials. We
chose the E-glass bre and resin epoxy used by Swan and Kosaka [10] as the two-component
material phases. In the three-phase design problem, a very compliant material is used to simulate
the void phase. The material properties each material phase are shown in Table I.
Figure 9 shows the solutions of the topology optimization of composite material of the plate for
the two-phase design and the three-phase design, respectively. These gures show the distribution
of the material phases in a unit cell, where the mesh division is 40 40. As shown in these
gures, a lot of checkerboard patterns have appeared in these solutions.
4.2. Perimeter control approach
For the purpose of preventing such checkerboard solutions, several methods have been proposed
by Haber et al. [24], Beckers [25], Duysinx [26], Swan and Kosaka [20], Petersson and Sigmund
[27], and it was summarized by Sigmund and Petersson [21]. According to Sigmund and Petersson
[21], the perimeter control approach [2426] and the mesh-independent ltering approach [27], are
eective to prevent such numerical instabilities in the topology optimization.
Figure 9. Optimal distribution of material phases in a unit cell for two- and three-phase designs:
(a) two-phase design; (b) three-phase design.
Figure 10. Values of the perimeter control function for some typical patterns of two-phase materials.
In the perimeter control approach, the checkerboarding solution is prevented by restricting the
length of the perimeter which appears by the dierence of the material density between adjacent
elements. If the discretized perimeter control approach proposed by Beckers [25] is applied to the
present approach, the perimeter control function P is dened by
P =
N
i=1
p
i
; p
i
=
m
i
j=1
l
ij
_
A
i

A
j
B
i

B
j
C
i

C
j
_
(51)
where
A
i
=a
i
b
i
;
B
i
=1 a
i
;
C
i
=1
A
i

B
i
, and m
i
denotes the number of elements which
have a common side with the ith element (m
i
=4 for an interior elements, m
i
=3 or 2 for edge
or corner elements). Figure 10 shows the values of the perimeter control function p
i
for some
typical patterns in the case of two-phase materials, where m
i
=4; l
ij
=1.
However, there are two problems in the perimeter control method. The rst is that the setting of
the constraint value of the perimeter length is dicult. If the constraint value of the perimeter length
is too small, the convergent solution may not be obtained. If it is too large, the checkerboarding
solution may occur. Many trials and errors are required in order to nd the appropriate constraint
Figure 11. Values of the gravity function for some typical patterns of two-phase materials.
value. The second problem is that the mixed condition of the material phases is allowed in this
method. It can be seen from gure that the pattern of 5 grey-scale elements yields the identical
value of p
i
as the pattern of 5 black or white elements. That is to say, the penalty is not imposed
on the grey-scale pattern. Therefore, if we use this method, the articially mixed material phases
are often not decomposed clearly in the solution. In the topology analysis of the macrostructure, the
solution of the mixed phases may be admitted since they can be interpreted as a composite structure
in the microscale. But in the analysis of the microstructure, the solution of the mixing phases is
not favourable since one runs the risk of interpreting the microstructure of microstructure in the
regions of the mixed constituent phases.
4.3. Proposed lter
In this paper, a novel ltering method for addressing the shortcomings of the perimeter control
method is proposed. In this method, the following gravity control function G is used instead of
the perimeter control function dened by Equation (51):
G =
N
i=1
g
i
; g
i
=
m
i
j=1
_
A
i

A
j
+
B
i

B
j
+
C
i

C
j
_
(52)
where
A
i
=a
i
b
i
;
B
i
=1a
i
;
C
i
=1
A
i

B
i
, and m
i
denotes the number of elements which have
a common side with the ith element. This function is based on the concept of gravity. Because the
gravity force of the particles next to each other is dened by
i
j
=r
2
ij
, where
i
;
j
are the masses
of the particles, and r
ij
denotes the distance of the particles. When the size of the nite elements
in the unit cell are all equal, the distances between the ith element and the elements which have a
common side with the ith element are all equal. Therefore, in Equation (52), r
2
ij
has been removed.
Figure 11 shows the values of the gravity control function g
i
for some typical patterns in the case
of two-phase materials. It can be seen from Figure 11 that the penalty is imposed on the grey
scale in this case, and the penalty is also imposed on the checkerboard pattern as well as the
perimeter control method.
In this method, if the magnitude of the gravity control function increases, both the checkerboard
and the grey-scale pattern decrease. Therefore, we add the gravity control function with a weighting
ratio to the objective function as follows:
f=C + w
g

C(G=G
ini
)
2
(53)
Figure 12. Optimal distribution of material phases in a unit cell for two- and three-phase designs:
(a) two-phase design; (b) three-phase design.
where C and

C are the mean compliance (C =

C in each step), but

C is just scaling constant
and independent of the design variables (i.e. @

C=@
i
=0). Another scaling constant, G
ini
is the
value of the gravity function calculated at the initial iteration. The coecient w
g
[0; 1] represents
the weighting ratio between the mean compliance and the gravity function. The power coecient
( =2) is set in order to emphasize the eect of the gravity control function. The sensitivity of this
multi-objective function is calculated by
@f
@
i
=
@C
@
i
2w
g

C
_
G=G
2
ini
_
@G
@
i
(54)
The numerical example shown in Figure 8 is solved to demonstrate the eectiveness of the
proposed ltering method. Figure 12 shows the optimum topology solutions of the microstructure.
The weighting ratios are set to 0.5 for two-phase design and 0.3 for three-phase design. It can be
seen that clear topologies are obtained by using the proposed ltering method.
5. DESIGN EXAMPLES
In this section, some numerical examples are presented to demonstrate the applicability of the
presented algorithm to the design of composite materials within the context of two-dimensional
structures. We chose the E-glass bre and resin epoxy used by Swan and Kosaka [5] as the two-
component material phases. In the three-phase design problem, a very compliant material is used
to simulate the void phase. The properties of these material phases are shown in Table I. In the
following analysis, the power s in Equations (36)(41) is set to 6 for the two-phase design and
5 for the three-phase design. The volume fractions of the constituent phases are set to

f
A
=0:5
and

f
B
=0:5 for two-phase materials, and

f
A
=0:4;

f
B
=0:4 and

f
C
=0:2 (void) for three-phase
materials. The mesh division of the unit cell is 40 40. These problems are solved under the
plane stress condition.
Figure 13 shows the shear plate models which have dierent boundary conditions. The horizontal
displacement of the right-hand edge is constrained in Case 1, and it is free in Case 2. The nite
Figure 13. Plate models subjected to shear forces with dierent boundary conditions.
Figure 14. Optimal distribution of material phases in four unit cells for two-phase design.
element division of the macrostructure is shown in these gures. The symmetrical condition with
respect to the two axes y
1
; y
2
and the diagonal line of the unit cell is imposed.
Figure 14 shows the solutions of the topology optimization of composite material for two-
material-phase design. Figure 15 shows the solutions for three-material-phase design. These gures
show the distribution of the material phase in four unit cells. The values of the mean compliance
and the weighting ratio in Equation (53) are shown in these gures. Tables II and III show the
components of the elasticity matrix in these solutions. It can be seen from these gures that the
topology of Case 1 is clearly dierent from that of Case 2. Tables II and III show that the stiness
for the normal stresses of Case 2 is greater than that of Case 1. It also shows that the topology
for two-phase design is similar to that of three-phase design, and the value of the shear stiness
D
H
33
is greater than 9.5 GPa which was obtained by Swan and Kosaka [5].
Figure 15. Optimal distribution of material phases in four unit cells for three-phase design.
Table II. Components of the elasticity matrix D
H
for two-phase design.
D
H
11
(GPa) D
H
12
(GPa) D
H
33
(GPa)
Case 1 18.2 9.72 10.9
Case 2 21.9 5.25 7.29
Table III. Components of the elasticity matrix D
H
for three-phase design.
D
H
11
(GPa) D
H
12
(GPa) D
H
33
(GPa)
Case 1 12.0 7.21 7.62
Case 2 14.9 4.07 4.95
Figure 16 shows the model of another design example. In this case, the loading condition is
changed, and the boundary conditions are the same as the previous example. The symmetrical
condition with respect to the two axes y
1
; y
2
is also imposed, but the symmetrical condition with
respect to the diagonal line is removed in this case. Figure 17 shows the solutions of the topology
optimization of composite material for two-material-phase design. Figure 18 shows the solutions
for three-material-phase design. Tables IV and V show the components of the elasticity matrix in
these solutions. It can be seen from these gures that the topology of Case 1 is also dierent from
Case 2. Table IV and V show that the stiness in the y
2
direction is greater than that in the y
1
direction in Case 1, but the stiness in the y
1
direction is greater than that in the y
2
direction in
Case 2. It is also seen that the shear stiness of Case 1 is greater than that of Case 2.
6. CONCLUSIONS
In this paper, the microstructural design algorithm of two-dimensional structures with composite
materials was presented. Using this algorithm, it is possible to design both two- and three-phase
microstructures. Using the articial mixture assumption, the intermediate Youngs modulus and
Figure 16. Plate models subjected to normal forces with dierent boundary conditions.
Figure 17. Optimal distribution of material phases in four unit cells for two-phase design.
Poissons ratio were made through the interpolation of two- or three-phase materials. The power-
law function was used for the interpolation of the Youngs modulus, and the linear function was
used for the interpolation of the Poissons ratio. The mean compliance of a macrostructure was
chosen as the objective function, and the volume fractions of each material phase were given as
the constraint conditions. To solve the optimization problem, the sequential linear programming
method with move limit was used. A method for the relaxation of a fractions constraint condition
was used to prevent local minima. The ltering method based on the concept of gravity was used
to prevent checkerboards. In this method, the gravity control function is added to the objective
function with weight.
It was shown by the example of the shear plate that the proposed ltering method is eective for
preventing the checkerboards. The eectiveness of the present algorithm for the design of composite
Figure 18. Optimal distribution of material phases in four unit cells for three-phase design.
Table IV. Components of the elasticity matrix D
H
for two-phase design.
D
H
11
(GPa) D
H
12
(GPa) D
H
22
(GPa) D
H
33
(GPa)
Case 1 17.9 9.01 19.9 10.1
Case 2 25.5 5.92 18.5 8.15
Table V. Components of the elasticity matrix D
H
for three-phase design.
D
H
11
(GPa) D
H
12
(GPa) D
H
22
(GPa) D
H
33
(GPa)
Case 1 12.5 4.78 14.3 5.72
Case 2 16.4 3.87 11.9 4.82
material was demonstrated by some numerical examples which have dierent boundary or loading
conditions. That is to say, the distinctive topology of the microstructures of the composite material
can be obtained using the present algorithm. The eciency of the calculation of the present
algorithm was also veried because the number of iterations required for the optimization procedure
was less than 68.
Finally, we note that the algorithm presented here can be directly extended to three-dimensional
problems.
ACKNOWLEDGEMENTS
This research was partially supported by Steel Material Club Corporation Aggregate in Japan. Its support is
gratefully acknowledged.
REFERENCES
1. Bendse MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Computer
Methods in Applied Mechanics and Engineering 1988; 71:197224.
2. Suzuki K, Kikuchi N. A homogenization method for shape and topology optimization. Computer Methods in Applied
Mechanics and Engineering 1991; 93:291318.
3. Murat F, Tartar L. Calcul des variations et homogeneisation, in: les methodes de lhomogeneisation: theorie et
applications en physique. Coll. de la Dir. des Etudes et Recherches de Electricite de France. Eyrolles, Paris, 1985;
319370.
4. Sigmund O. Material with prescribed constitutive parameters: an inverse homogenization problem. International Journal
of Solids and Structures 1994; 31(17):23132329.
5. Sigmund O. Tailoring materials with prescribed elastic properties. Mechanics of Materials 1995; 20(4):351368.
6. Sigmund O. Design of material structures using topology optimization. Report No. 502, Danish Center for Applied
Mathematics and Mechanics, Technical University of Denmark, 1995; 4350.
7. Sigmund O, Torquato S. Design of materials with extreme thermal expansion using a three-phase topology optimization
method. Journal of the Mechanics and Physics of Solids 1997; 45(6):10371067.
8. Fonseca JSO. Design of microstructures of periodic composite materials. PhD Thesis, The University of Michigan,
1997.
9. Swan CC, Arora JS. Topology design of material layout in structured composites of high stiness and strength.
Structural Optimization 1997; 13:4559.
10. Swan CC, Kosaka I. Homogenization-based analysis and design of composites. Computers and Structures 1997;
64:603621.
11. Silva, ECN, Kikuchi N. Optimal design of piezoelectric microstructures. Computer Mechanics 1997; 19(5):398410.
12. Silva ECN, Fonseca JSO, Kikuchi N. Optimal design of periodic piezocomposites. Computer Methods in Applied
Mechanics and Engineering 1998; 159:4977.
13. Gibiansky LV, Torquato S. On the use of homogenization theory to design optimal piezocomposites for hydrophone
applications. Journal of the Mechanics and Physics of Solids 1997; 45(5):689708.
14. Gibiansky LV, Torquato S. Optimal design of 13 composite piezoelectrics. Structural Optimization 1997; 13:2328.
15. Haslinger J, Dvorak J. Optimum composite material design. RAIRO Modelisation Mathematiqueet Analyse Numerique
1995; 29(6):657686.
16. Bendse MP. Optimal shape design as a material distribution problem. Structural Optimization 1989; 1:193202.
17. Rozvany GIN, Zhou M, Birker T, Sigmund O. Topology optimization using iterative continuum type optimality criteria
(COC) methods for discretized systemes. In Topology Design of Structures. Bendse MP, Soares CA (eds), Springer:
Berlin, 1992; 273286.
18. Mlejnek HP, Schirrmacher R. An engineers approach to optimal material distribution and shape nding. Computer
Methods in Applied Mechanics and Engineering 1993; 106:126.
19. Yang RJ, Chuang CH. Optimal topology design using linear programming. Computers and Structures 1994; 52:
265275.
20. Swan, CC, Kosaka I. Voigt-Reuss topology optimization for structures with linear elastic material behaviors.
International Journal for Numerical Methods in Engineering 1997; 40:30333057.
21. Sigmund O, Petersson J. Numerical instabilities in topology optimization: A survey on procedures dealing with
checkerboards, mesh-dependencies and local minima. Structural Optimization 1998; 16:6875.
22. Jog CS, Habar RB. Stability of nite element models for distributed-parameter optimization and topology design.
Computer Methods in Applied Mechanics and Engineering 1996; 130:203226.
23. Diaz A, Sigmund O. Checkerboard patterns in layout optimization. Structural Optimization 1995; 10:4045.
24. Haber RB, Bendse MP, Jog C. A new approach to variable-topology shape design using a constraint on the perimeter.
Structural Optimization 1996; 11:112.
25. Beckers M. Optimisation de structures en variables disct3tes. Th3se de Doctorat, Universite de Liege, Ing3nieur civil
3lectricien-m3canicien(a3rospatiale), 1997; 108135.
26. Duysinx P. Layout optimization: a mathematical programming approach. DACAMM Report, No.540, Technical
University of Denmark, 1997.
27. Petersson J, Sigmund O. Slope constrained topology optimization. International Journal for Numerical Methods in
Engineering 1998; 41:14171434.
28. Lions JL, Some Methods in the Mathematical Analysis of Systems and Their Control. Science Press: Beijing, China,
Gordon and Breach Science Publishers: New York, 1981.
29. Wilson EL, Taylor RL, Doherty WP, Ghaboussi J. Incompatible displacement models. In Numerical and Computer
Methods in Structural Mechanics, Fenves ST et al. (eds), Academic Press: New York, 1973; 4357.
30. Taylor RL, Beresford PJ, Wilson EL. A non-conforming element for stress analysis. International Journal for
Numerical Methods in Engineering 1976; 10:12111220.

Composite Material Design of Two-Dimensional Structures Using The Homogenization Design Method

Uploaded by

Copyright:

Available Formats

Composite Material Design of Two-Dimensional Structures Using The Homogenization Design Method

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Composite Material Design of Two-Dimensional Structures Using The Homogenization Design Method

Uploaded by

Copyright:

Available Formats

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING

Int. J. Numer. Meth. Engng 2001; 50:20312051

You might also like