INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING

Int. J. Numer. Meth. Engng 2010; 84:733–756

Published online 17 May 2010 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/nme.2922

Topological derivative for multi-scale linear elasticity models

applied to the synthesis of microstructures

S. Amstutz1 , S. M. Giusti2 , A. A. Novotny2 and E. A. de Souza Neto3, ∗, †

1 Laboratoire
d’analyse non-linéaire et géométrie, Faculté des sciences 33,
rue Louis Pasteur 84000 Avignon, France
2 Laboratório Nacional de Computação Científica LNCC/MCT, Coordenação de Matemática Aplicada e

Computacional, Av. Getúlio Vargas 333, 25651-075 Petrópolis RJ, Brazil

3 Civil and Computational Engineering Centre, School of Engineering, Swansea University, Singleton Park,

Swansea SA28PP, U.K.

SUMMARY
This paper proposes an algorithm for the synthesis/optimization of microstructures based on an
exact formula for the topological derivative of the macroscopic elasticity tensor and a level set domain
representation. The macroscopic elasticity tensor is estimated by a standard multi-scale constitutive
theory where the strain and stress tensors are volume averages of their microscopic counterparts over a
representative volume element. The algorithm is of simple computational implementation. In particular,
it does not require artificial algorithmic parameters or strategies. This is in sharp contrast with existing
microstructural optimization procedures and follows as a natural consequence of the use of the topological
derivative concept. This concept provides the correct mathematical framework to treat topology changes
such as those characterizing microstuctural optimization problems. The effectiveness of the proposed
methodology is illustrated in a set of finite element-based numerical examples.Copyright 䉷 2010 John
Wiley & Sons, Ltd.

Received 8 December 2009; Revised 18 March 2010; Accepted 24 March 2010

KEY WORDS: topological derivative; sensitivity analysis; multi-scale modelling; level set domain repre-
sentation; synthesis of microstructures

∗ Correspondence to: E. A. de Souza Neto, Civil and Computational Engineering Centre, School of Engineering,
Swansea University, Singleton Park, Swansea SA28PP, U.K.
† E-mail: [email protected]

Contract/grant sponsor: The Royal Society; contract/grant number: JP080066

Contract/grant sponsor: Brazilian research funding council, CNPq; contract/grant number: 382485/2009-2

Copyright 䉷 2010 John Wiley & Sons, Ltd.

734 S. AMSTUTZ ET AL.

1. INTRODUCTION

The accurate prediction of the constitutive behavior of a continuum body under loading is of
paramount importance in many areas of engineering and science. Until about a decade ago, this
issue has been addressed mainly by means of conventional phenomenological constitutive theories.
More recently, the increasing understanding of the microscopic mechanisms responsible for the
macroscopic response, allied to the industrial demand for more accurate predictive tools, led to the
development and use of the so-called multi-scale constitutive theories. Such theories are currently
a subject of intensive research in applied mathematics and computational mechanics. Their starting
point can be traced back to the pioneering developments reported in [1–6]. Early applications
were concerned with the description of relatively simple micro-scale phenomena often treated
by analytical or semi-analytical methods [7–12]. More recent applications rely often on finite
element-based computational simulations [13, 14] and are frequently applied to more complex
material behavior in areas such as the modelling of human arterial tissue [15], bone [16], the
plastic behavior of porous metals [17] and the microstructural evolution and phase transition in
the solidification of metals [18].
One interesting branching of such developments is the study of the sensitivity of the macroscopic
response to changes in the underlying microstucture. The sensitivity information becomes essential
in the analysis and potential purpose-design and optimization of heterogeneous media. For instance,
sensitivity information obtained by means of a relaxation-based technique has been successfully
used in [19–21] to design microstructural topologies with negative macroscopic Poisson ratio.
Multi-scale models have also been applied with success to the topology optimization of load-bearing
structures in the context of the so-called homogenization approach to topology optimization (see,
for instance, the review paper by Eschenauer and Olhoff [22]) based on the fundamental papers
by Bendsøe and Kikuchi [23] and Żochowski [24]. In such cases, the microscale model acts as a
regularization of the exact problem posed by a material point turning into a hole [25]. The homog-
enization approach has also been applied to microstructural topology optimization problems where
the target is the design of topologies that yield pre-specified or extreme macroscopic response
[26–28]. One of the drawbacks of this methodology, however, is that it often produces designs
with large regions consisting of perforated material. To deal with this problem, a penalization of
intermediate densities is commonly introduced.
In contrast to the homogenization approach, we propose in this paper a microstructural
synthesis/optimization algorithm relying on the mathematical concept of topological derivative
[29–31] combined with a level set domain representation. In this context, a (remarkably simple)
exact formula for the sensitivity of the macroscopic elasticity tensor to the insertion of inclusions
at the micro-scale is used. This analytical formula provides a rigorous first-order approximation to
the variation in the macroscopic elasticity tensor resulting from the insertion of a circular inclusion
of a given phase contrast. It is derived by making use of the notions of topological asymptotic anal-
ysis and topological derivative within the variational formulation of well-established multi-scale
constitutive theory [13, 32], where the macroscopic strain and stress tensors are volume averages
of their microscopic counterparts over a representative volume element (RVE) of material.
The relatively new concept of topological asymptotic expansion allows the exact calculation of
the sensitivity of a given shape functional with respect to infinitesimal domain perturbations such as
the insertion of voids, inclusions, source terms or even cracks. It has proved to be extremely useful
in the treatment of a wide range of problems, including the topology optimization of load-bearing
structures [33–35], inverse analysis [36–38] and image processing [39–41], and is currently a

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 735

subject of great interest in applied mathematics circles [42–45]. In the present context, when
combined with the ideas introduced in [34], the use of the exact topological sensitivity formula
results in a microstructural optimization algorithm of simple computational implementation.
In particular, the algorithm is free from artificial algorithmic parameters and does not require
extra post-processing strategies when new topologies are generated throughout its iterations. This
is only a natural consequence of the use of the notion of topological derivative which provides the
correct mathematical framework for the treatment of the inherently singular changes in topology
that characterize microstructural optimization problems. The effectiveness and robustness of the
proposed algorithm are illustrated in a number of finite element-based numerical examples of
microstructural optimization.
The paper is organized as follows. The multi-scale constitutive framework adopted in the estima-
tion of the macroscopic elasticity tensor is briefly reviewed in Section 2. The concept of topological
asymptotic expansion and topological derivative—upon which the algorithm proposed in this paper
relies—is briefly reviewed in Section 3, where the closed formula for the topological derivative of
the macroscopic elasticity tensor is also presented. The proposed microstructural topology opti-
mization algorithm is presented in Section 4 and the numerical examples are shown in Section 5.
The paper closes in Section 6 with some concluding remarks.

2. MULTI-SCALE MODELLING

The homogenization-based multi-scale constitutive framework presented, among others, in

[13, 14, 32], is adopted here in the estimation of the macroscopic elastic response from the
knowledge of the underlying microstructure. The main idea behind this well-established family
of constitutive theories is the assumption that any point x of the macroscopic continuum (refer to
Figure 1) is associated with a local RVE whose domain  , with boundary * , has characteristic
length L  , much smaller than the characteristic length L of the macro-continuum domain .
The axiomatic structure of this family of theories is described in detail in [46, 47]. Accordingly,
the entire theory can be formally derived on the basis of five basic assumptions: (i) the strain-
averaging relation; (ii) the requirement that the chosen functional set of kinematically admissible
displacement fluctuations of the RVE be a subspace of the minimally constrained space compatible
with the strain-averaging hypothesis; (iii) the mechanical equilibrium of the RVE; (iv) the stress-
averaging relation, and; (v) the Hill–Mandel principle of macro-homogeneity [2, 3], which ensures
the energy consistency between the so-called micro- and macro-scales.

x
y
Lμ
L

Figure 1. Macroscopic continuum with a locally attached microstructure.

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
736 S. AMSTUTZ ET AL.

Assumption (i) uses the concept of homogenization to define the macroscopic strain tensor ε
at a point x of the macroscopic continuum as the volume average of its microscopic counterpart
ε := ∇ s u  over the corresponding RVE:
 
1 1
ε := ∇ u =
s
u  ⊗s n, (1)
V  V *

where V denotes the total volume of the RVE, u  is the microscopic displacement field of the
RVE, n is the outward unit normal to the boundary * , ∇ s is the symmetric gradient operator
and ⊗s denotes the symmetric tensor product between vectors.
For the purposes of the present paper, we shall consider RVEs whose domain consists of a
matrix m , containing inclusions of different materials occupying a domain  (Figure 1). Further,
i

we shall assume the matrix and the inclusions to be modelled as isotropic homogeneous materials.
Hence, the microscopic stress tensor field  () satisfies

 () = C ∇ s , (2)

where C is the fourth-order isotropic elasticity tensor field of the RVE:

E
C := [(1− )I+ (I⊗I)], (3)
1−2

with I and I denoting, respectively, the second- and fourth-order identity tensors and E  and 
Young’s modulus and Poisson ratio fields of the RVE, given by
 m  m
E  if y ∈ m  if y ∈ m 
E  (y) := and  (y) := i . (4)
E  if y ∈ 
i i
 if y ∈  i

If the RVE has more than one inclusion, the parameters E i and i are considered piecewise
constant over i .
Without the loss of generality, u  may be decomposed into a sum

u  (y) = u +u  (y)+
u  (y), (5)

of a constant (rigid) RVE displacement coinciding with the macroscopic displacement u(x), a linear
field u  (y) := εy, and a displacement fluctuation field  u  (y). This, together with (1), (2) and the
Hill–Mandel principle of macro-homogeneity [46, 47] leads to the RVE mechanical equilibrium
problem that consists of finding, for a given macroscopic strain ε, an admissible microscopic
displacement fluctuation field  u  ∈ U , such that
 
 (u  )·∇ s +  (
u  )·∇ s  = 0 ∀ ∈ U , (6)
 

where

 (u  ) = C ∇ s u  ,  (
u  ) = C ∇ s 
u, (7)

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 737

and the (yet to be defined) space U of admissible displacement fluctuations (and virtual kinemati-
  —the minimally constrained
cally admissible displacements) fields of the RVE is a subspace of U
space of kinematically displacements of the RVE compatible with (1):
   
U ⊂ U  := v ∈ [H ( )] :
1 2
v = 0, v ⊗s n = 0 . (8)
 *

Once the problem (6) has been solved, the macroscopic stress tensor  is obtained as the volume
average of the microscopic stress field  (u  ) =  (u  )+ (
u  ) over the RVE, i.e.

1
=  (u  ). (9)
V 

The complete characterization of the multi-scale constitutive model is obtained by defining

the subspace U ⊂ U   of kinematically admissible displacement fluctuations. In general, different
choices produce different macroscopic responses for the same RVE. For example, the choice
U = U  results in a lower bound for the homogenized elastic properties. In the present paper,
the analysis will be focussed on media with periodic microstructure. In this case, the geometry
of the RVE cannot be arbitrary and must represent a cell whose periodic repetition generates the
macroscopic continuum. In addition, the displacement fluctuations must satisfy periodicity on the
boundary of the RVE. Accordingly, we have

U := {   :
u ∈ U u  (y + ) =
u  (y − ) ∀(y + , y − ) ∈ P}, (10)

where P is the set of pairs of points, defined by a one-to-one periodicity correspondence, lying
on opposing sides of the RVE boundary.

2.1. The homogenized elasticity tensor

Crucial to the developments presented in Section 3 is the closed form of the homogenized elasticity
tensor C for the multi-scale model defined in the above. This can be obtained by means of the
methodology suggested in [13], which is based on rewriting problem (6) as a superposition of
linear problems associated with the individual Cartesian components of the macroscopic strain
tensor. The components of the homogenized elasticity tensor C, in the orthonormal basis {ei } of
the Euclidean space, can be written as:

1
(C)ijkl = ( (u kl ))ij . (11)
V 

By virtue of (6), the canonical microscopic displacement field u kl associated with the tensors
εkl = ek ⊗el is the solution of the equilibrium equation

 (u kl )·∇ s  = 0 ∀ ∈ U . (12)


Note that, in view of (5), we have

u kl −u −(ek ⊗el )y =

u kl ∈ U , (13)

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
738 S. AMSTUTZ ET AL.

so that problem (12) is equivalent to finding the field  u kl ∈ U such that
 
 (
u kl )·∇ s +  (ek ⊗s el )·∇ s  = 0 ∀ ∈ U . (14)
 

3. THE TOPOLOGICAL SENSITIVITY OF THE HOMOGENIZED ELASTICITY TENSOR

A closed formula for the sensitivity of the homogenized elasticity tensor of (11) to the nucleation
of a circular inclusion within the RVE is presented in this section. The presented formula is central
to the algorithm for microstructural synthesis/optimization proposed and applied in Sections 4
and 5 and relies on the concepts of topological asymptotic expansion and topological derivative.
These relatively new mathematical concepts have been originally introduced by Sokołowski and
Żochowski [29]. They extend the conventional notion of differentiability and provide the correct
mathematical framework for the treatment of inherently singular topological domain changes such
as those produced by the introduction of an inclusion, as considered in the present context.

3.1. Topological asymptotic expansion. Preliminary concepts

Let  be a shape functional whose value depends on a given domain the topology of which
is parameterized by a non-negative scalar . A nil value of  defines what we refer to as the
original (or unperturbed) domain, so that (0) is the value taken by the shape functional for such a
domain. Any other value >0 defines a domain that differs from the original one by a topological
perturbation of size . Let us assume that  has sufficient regularity so that the following expansion
is possible:
() = (0)+ f ()DT +o( f ()), (15)
where () denotes the value of the functional for the topologically perturbed domain, f () is a
non-negative function such that f () → 0 when  → 0 and o( f ()) contains all terms of higher order
in f (). Expression (15) is referred to as the topological asymptotic expansion of the functional
 and the term DT  of (15) is named the topological derivative of  at the unperturbed domain.

3.2. Topological derivative of the homogenized elasticity tensor

To apply the concepts of topological asymptotic expansion and topological derivative in the present
context of multi-scale linear elasticity models, we begin by defining the so-called original domain
as a given generic RVE of the type described in Section 2. The topological perturbation of the
original domain  consists of the nucleation of a small inclusion of radius  denoted by I .
Formally, the perturbed RVE domain   is obtained by first introducing a circular hole H
of radius  centered at an arbitrary point  y ∈  and then replacing the region of the hole by a
circular inclusion I of different material properties. That is, the perturbed domain is defined
as  = ( \H )∪I (refer to Figure 2).
We shall assume that the microscopic elasticity tensor field of the perturbed RVE is given by

 C in  \H
C := (16)
C in I ,

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 739

ρ
y^
n

Figure 2. Topological perturbation at the microscopic level.

where the parameter 0 defines the contrast between the elasticity tensor of the region of the
original domain  , where the inclusion was inserted and the elasticity tensor of the inclusion
material. We note that this type of perturbation corresponds to a change only in Young’s modulus
of the phases. Now, with this definition at hand, we denote the effective microscopic stress in the
domain  by:

  s  () in  \H
 () = C ∇  = (17)
 () in I .
From the above expression and (11) the components of the macroscopic elasticity tensor corre-
sponding to   can be obtained as

 1  
(C )ijkl = ( (u kl ))ij , (18)
V 

where the canonical microscopic displacement field u kl associated to the topologically perturbed
domain satisfies

 
 (u kl )·∇ s  = 0 ∀ ∈ U (19)


and, analogous to (13), we have

 
u kl −u −(ek ⊗el )y =
u kl ∈ U , (20)

where u kl is the displacement fluctuation field associated with the perturbed domain  and to
the tensor ek ⊗el .
With the above definitions at hand, we may now specialize (15) by identifying the shape
functional  with the homogenized elasticity tensor. The resulting topological asymptotic expansion
of C reads
C = C+ f ()DT C+o( f ()). (21)
This expansion is indeed possible and its derivation for two-dimensional problems is addressed in
the Appendix. The function f in this case is found to be

2
f () = , (22)
V
i.e. it represents the volume fraction of the inserted inclusion. The topological derivative of the
homogenized elasticity tensor at the unperturbed RVE domain, denoted DT C, is a fourth-order

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
740 S. AMSTUTZ ET AL.

tensor field over  that provides a rigorous first-order approximation (in ) to the change in C
resulting from the insertion of a circular inclusion of radius  within the RVE. Its closed form
expression in Cartesian components reads

(DT C)ijkl = H  (u ij )· (u kl ), (23)

where the fields u ij are the solutions to (12) for the unperturbed RVE domain and the fourth-order
tensor H is defined as
 

1 1− 1− (−2)
H := − 4I− (I⊗I) . (24)
E 1+ 1+

with
1+ 3−
= , = . (25)
1− 1+
Remark 1
The above analytical formula for the topological derivative field is particularly suitable for imple-
mentation within finite element-based numerical frameworks. Its computation in this context
(considered in Sections 4 and 5) is extremely simple: Once the finite element approximations to
the vector fields u ij are obtained as solutions of the discretized version of problem (12) for the
original domain, the corresponding discretized version of the fourth-order tensor field DT C can
be promptly calculated according to (23).

Remark 2
Expression (23) allows the exact derivative of any differentiable function of C with respect to the
volume fraction of inclusion to be calculated through the direct application of the conventional
rules of differential calculus. That is, any such function h(C) has exact topological derivative

DT h =
Dh(C), DT C, (26)

with the brackets

·, · denoting the appropriate product between the two terms—the derivative of
h with respect to C and the topological derivative of C. Note, for example, that the properties
of interest such as the homogenized Young’s, shear and bulk moduli as well as the Poisson ratio
are all regular functions of C. This observation together with Remark 1 point strongly to the
suitability of the use of (21) in a finite element-based framework for the synthesis and optimization
of elastic micro-structures based on the minimization/maximization of cost functions defined in
terms of homogenized properties. This is the main aim of the present paper and will be pursued
in Sections 4 and 5.

Remark 3
The tensor H has an explicit dependency on the phase contrast parameter . In the limiting case
with → ∞, corresponding to the insertion of a rigid inclusion, we have

1 −2
H∞ = 4I+ (I⊗I) . (27)
E 

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 741

The limit → 0, on the other hand, corresponds to the insertion of a hole in the RVE. In this case
we have
1
H0 = − [4I−(I⊗I)]. (28)
E
This last expression coincides with the result derived in [48].

4. A TOPOLOGICAL DERIVATIVE-BASED MICRO-STRUCTURE DESIGN ALGORITHM

In this section we propose a microstucture topology design algorithm based on the topological
asymptotic expansion (21) of the homogenized elasticity tensor. The RVE domain here represents
a bi-material microstructure whose constituents properties are defined by a given elasticity tensor
C∗ and phase contrast ∗ so that, as in (16), we have
 ∗
C ∀y ∈ 1
C (y) = (29)
∗ C∗ ∀y ∈ 2 ,

where 1 and 2 denote the domains occupied by materials 1 and 2, respectively. The sought
topology of the RVE is the solution of the general optimization problem stated as

|1 |
Minimize J (1 ) = h(C)+ , (30)
1 ⊂ V

where the cost function J is a scalar-valued functional of the domain 1 , h is a generic (yet to be
defined) function of the homogenized elasticity tensor C and is a fixed Lagrange multiplier that
imposes a constraint on the volume ratio of material 1 (with volume denoted |1 |).
It should be stressed that the design variable in problem (30) is the topology of the domain
1 . Hence, the use of the exact topological sensitivity information provided by the topological
derivative (23) emerges as a natural alternative in the development of a numerical optimization
algorithm to tackle the problem. In this context, by taking into account the comments made in
Remark 2, we see that the topological derivative of the cost function, i.e. the exact derivative of J
with respect to the volume fraction of an inclusion of radius  centred at an arbitrary point y ∈  ,
is given by the simple general formula
DT J (y) =
Dh(C), DT C(y)+ . (31)
At this point, it is worth remarking that the inclusion referred to above is assumed to be made of
material 2 at points y ∈ 1 and of material 1 at points y ∈ 2 so that (31) always measures the
sensitivity of J when materials 1 and 2 are interchanged within the RVE. This allows the algorithm
(described below) to replace, for example, material 1 with material 2 and then replace material
2 back with material 1 and so on at any given region of the RVE. Following this observation,
the computation of (31) is carried out using the expressions (23)–(25) with = ∗ if y ∈ 1 and
= 1/ ∗ if y ∈ 2 .
Having made the above considerations, the topological derivative-based optimization algorithm
devised in [34] stands out as a particularly well-suited choice to solve problem (30). The procedure

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
742 S. AMSTUTZ ET AL.

relies on a level set domain representation [49] and the approximation of the topological optimality
conditions by a fixed point iteration. In particular, the algorithm displays a marked ability to
produce general topological domain changes uncommon to other methodologies based on a level
set representation and has been successfully applied in [34] to topology optimization in the context
of two-dimensional elasticity and flow through porous media. For completeness, the algorithm is
outlined in the following. For further details we refer to [34].
With the adoption of a level set domain representation, the current material 1 domain, 1 , is
characterized by a function  ∈ L 2 ( ) such that

1 = {y ∈  , (y)<0}, (32)

whereas the material 2 domain is defined by

2 = {y ∈  , (y)>0}. (33)

Now, let us consider the topological derivative field DT J of formula (31). According to Amstutz
and Andrä [34], an obvious sufficient condition of local optimality of problem (30) for the class
of perturbations consisting of circular inclusions is

DT J (y)>0 ∀y ∈  . (34)

To devise a level set-based algorithm whose aim is to produce a topology that satisfies (34) it
is convenient to define the function

−DT J (1 )(y) if y ∈ 1
g(y) = (35)
DT J (1 )(y) if y ∈ 2 .

With the above definition and (32), (33) it can be easily established that the sufficient condition
(34) is satisfied if the following equivalence relation between g and the level set function  holds

∃>0 s.t. g = , (36)

or, equivalently,

g, 
 := arccos = 0, (37)
g L 2  L 2

where  is the angle between the vectors g and  in L 2 ( ).

Starting from a given level set function 0 ∈ L 2 ( ) which defines the chosen initial guess for the
optimum topology, the algorithm proposed in [34] produces a sequence (i )i∈N of level set func-
tions (corresponding to a sequence of RVE topologies) which provides successive approximations
to the sufficient condition for optimality (36). The sequence satisfies

0 ∈ L 2 ( ),
(38)
n+1 ∈ co(n , gn ) ∀n ∈ N,

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 743

where co(n , gn ) is the convex hull of {n , gn }. In the actual algorithm, the initial guess 0 is
normalized. With S denoting the unit sphere in L 2 ( ), the algorithm is explicitly given by

0 ∈ S,

1 gn (39)
n+1 = sin((1−n )n )n +sin(n n ) ∀n ∈ N,
sin n gn  L 2

where n ∈ [0, 1] is a step size determined by a line-search in order to decrease the value of the cost
functional J and, by construction of (39)2 , we have that n+1 ∈ S ∀n ∈ N. The iterative process
is stopped when for some iteration the obtained decrease in J is smaller than a given numerical
tolerance. If, at this stage, the optimality condition (36),(37) is not satisfied to the desired degree
of accuracy, i.e. if

n+1 > , (40)

where  is a pre-specified convergence tolerance, then a uniform mesh refinement of the RVE is
carried out and the procedure is continued.

4.1. Finite element implementation

In the computation of C according to (11)–(14), we obtain the fields ũ kl by a finite element
discretization of (14). In particular, the periodicity of the boundary displacement fluctuations is
enforced by a direct approach according to the implementation described in [17]. In the computation
of DT J with (31), a finite element approximation to the topological gradient field DT C is used.
The approximation is obtained by first using the finite element solutions u  kl in (23) and then
smoothing each component field (DT C)ijkl in a standard fashion. That is, the final discretized
version of the tensor field DT C is generated by the finite element shape functions with smoothed
nodal values of (DT C)ijkl . The level set functions n are generated by the same shape functions used
in the finite element approximation of problem (14). In this context, the following approximation is
used for the bi-material micro-structure. The material property associated with 1 (2 ) is assigned
to the nodes with negative (positive) level set function n . The material property within the element
is obtained by a standard interpolation using the element shape functions. Obviously, the resolution
of the RVE domain topology depends directly on the fineness of the adopted mesh.

Remark 4
The present algorithm is of simple implementation. The updated level set function n+1 obtained
at iteration n +1 according to (39) is only a linear combination between the known function n
and the corresponding function gn , which depends on the topological gradient DT J for the known
topology of iteration n. The computation of these quantities is straightforward. Note, in particular,
that no artificial algorithmic parameters or post-processing strategies are required throughout the
iterations. This is in sharp contrast with existing microstructural optimization strategies and follows
as a natural consequence of the use of the concept of topological derivative. This concept provides
the correct mathematical framework for the rigorous treatment of topology changes of the type
present in microstructural optimization problems.

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
744 S. AMSTUTZ ET AL.

5. NUMERICAL EXAMPLES

The application of the above optimization algorithm to the synthesis of microstructures is illustrated
in this section by means of a series of numerical examples.
Following the usual convention we write the homogenized two-dimensional elasticity tensor
C in matrix form:
⎛ ⎞
(C)1111 (C)1122 (C)1112
⎜ ⎟
C = ⎝ (C)1122 (C)2222 (C)2212 ⎠ . (41)
(C)1112 (C)2212 (C)1212

For the case orthotropic symmetry, the effective material properties—Young’s, bulk and shear
moduli as well as the Poisson ratios—are related explicitly to the components of the compliance
tensor C−1 , whose matrix representation reads
⎛ ⎞
1 12
⎛ −1 ⎞ ⎜ − 0 ⎟
(C )1111 (C−1 )1122 0 ⎜ E1 E1 ⎟
⎜ ⎜
⎟ ⎜ 21 ⎟
−1 ⎜ −1 −1 ⎟ = ⎜− 1 ⎟
C = ⎝ (C )1122 (C )2222 0 ⎠ ⎜ E 0 ⎟, (42)
E ⎟
⎜ 2 2 ⎟
0 0 (C−1 )1212 ⎝ 1⎠
0 0
G
where E 1 , E 2 are the effective Young’s moduli along the orthotropy directions e1 and e2 , respec-
tively, G is the effective in-plane shear modulus and 12 and 21 are the effective Poisson ratios
satisfying
21 12
= . (43)
E2 E1
In the examples that follow we shall use the algorithm of Section 4 to synthesize microstructures
with optimized effective properties. To this end, it is convenient to define the following types of
functions h(C):
• The first type of functions is defined as

h(C) = C−1 1 ·2 , (44)

where 1 and 2 are second-order tensors to be defined later. The topological derivative of
the associated cost function J of (30) in this case reads

DT J = −(C−1 (DT C)C−1 )1 ·2 + . (45)

• The second type is defined as

C−1 1 ·2 C−1 2 ·1

h(C) = + . (46)
C−1 1 ·1 C−1 2 ·2

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 745

Figure 3. Microstructural synthesis examples: (a) initial guess for the optimum topology and
(b) initial finite element mesh adopted.

The corresponding topological derivative of J here is

(C−1 (DT C)C−1 )1 ·[(C−1 1 ·1 )2 −(C−1 1 ·2 )1 ]
DT J = −
(C−1 1 ·1 )2

(C−1 (DT C)C−1 )2 ·[(C−1 2 ·2 )1 −(C−1 2 ·1 )2 ]
− + . (47)
(C−1 2 ·2 )2

In all examples the RVE domain is taken as the unit square  = (0, 1)×(0, 1). Young’s modulus
and the Poisson ratio of materials 1 and 2 are, respectively, E 1 = 1, E 2 = 0.01 and 1 = 2 = 0.3.
That is, we consider a phase contrast parameter ∗ = 0.01. The specified convergence tolerance for
the angle  is  = 1◦ . We remark that our numerical experience shows that this choice provides
a particularly stringent criterion producing solutions that satisfy the sufficient condition (36) for
optimality very closely. The adopted initial guess for the optimum topology in all cases is defined
by the level set function
1
0 = [cos2 (
(x − x0 )) cos2 (
(y − y0 ))−0.5], (48)
N
where N is the normalizing constant that ensures that  L 2 = 1. Unless otherwise stated, we choose
(x0 , y0 ) = (0.5, 0.5). The topology of the RVE domain in this case is illustrated in Figure 3(a),
where the dark- and light-coloured areas denote, respectively, materials 1 and 2. Also depicted in
Figure 3 is the initial (regular) mesh discretizing the RVE. It consists of 6400 three-noded linear
triangles and a total of 3281 nodes.

5.1. Basic cases

In this section we consider only functions h(C) of the type (44). In each case the algorithm
converges with the initial mesh to the prescribed accuracy of  = 1◦ . Despite the good accuracy
attained with the initial mesh alone, we perform one uniform mesh refinement step to further
improve the accuracy of the results and the resolution of the optimized RVE topology. The final
meshes have 25 600 elements and 12 961 nodes.

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
746 S. AMSTUTZ ET AL.

Figure 4. Horizontal rigidity maximization: (a) optimum RVE topology and (b) periodic microstructure.

(a) (b)

Figure 5. Horizontal rigidity maximization. Convergence history: (a) cost function and (b) angle .

5.1.1. Horizontal rigidity maximization. This is a trivial case where we choose 1 = 2 = e1 ⊗e1 ,
with e1 denoting the unit vector in the horizontal direction. Then, the corresponding function h(C)
is given by

1
h(C) := (C−1 )1111 = . (49)
E1

Note that the minimization of h(C) in this case corresponds to the maximization of the longitudinal
Young’s modulus E 1 . In addition we choose = 10. The obtained optimized RVE topology is
shown in Figure 4(a). The corresponding periodic microstructure is illustrated in Figure 4(b). The
evolution of the cost function and angle  throughout the iterations of (39)2 is shown in Figure 5.
Here the mesh refinement step has been performed at iteration 9, at which the residual angle  is
sufficiently small. Following the refinement, there is a small increase in  (iteration 10), followed
by a sharp decrease at iteration 11, at which final convergence is attained.

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 747

Figure 6. Bulk modulus maximization: (a) optimized RVE topology and

(b) corresponding periodic microstructure.

5.1.2. Bulk modulus maximization. In this case, we choose 1 = 2 = e1 ⊗e1 +e2 ⊗e2 . The corre-
sponding function h(C) is
1 1
h(C) := (C−1 )1111 +2(C−1 )1122 +(C−1 )2222 = (1−12 )+ (1−21 ). (50)
E1 E2
We choose = 20. The optimum RVE topology obtained is depicted in Figure 6 alongside an
illustration of the corresponding periodic microstructure. The total number of iterations needed to
reach the optimum topology in the present case was 21 (including the interactions after the mesh
refinement step).

5.1.3. Shear modulus maximization. Here we take 1 = 2 = e1 ⊗e2 +e2 ⊗e1 so that h(C) is simply

h(C) = 4(C−1 )1212 . (51)

The minimization of the present function h corresponds to the maximization of the effective shear
modulus G. The Lagrange multiplier is chosen as = 50. The final optimized topology obtained
is shown in Figure 7. A total of 10 interations were required to achieve convergence.

5.2. The Poisson ratio optimization

In this section we present four examples associated with the optimization of the Poisson ratio.
In the first two cases we consider functions h(C) of the form (44) and, in the last two cases,
functions of the type (46). In all four examples, we perform two steps of uniform mesh refinement
during the optimization process. The final meshes contain a total of 102 400 three-noded elements
and 51 521 nodes. The Lagrange multiplier is set to = 0 so that no constraints are imposed on
the volume ratio of the phases.

5.2.1. Minimization of a modified Poisson ratio. In this first case, we set 1 = e1 ⊗e1 and 2 =
−e2 ⊗e2 in (44). The corresponding function h(C) is
12 21
h(C) := −(C−1 )1122 = = . (52)
E1 E2

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
748 S. AMSTUTZ ET AL.

Figure 7. Shear modulus maximization: (a) optimized RVE and (b) corresponding periodic microstructure.

Figure 8. Minimization of a modified Poisson ratio: (a) optimized RVE and

(b) corresponding periodic microstructure.

The final topology, obtained here in 45 iterations, is shown in Figure 8. The corresponding
homogenized elasticity tensor, in matrix representation, is found to be
⎛ ⎞
0.0825 −0.0308 0
⎜ ⎟
C = ⎝ −0.0308 0.0825 0 ⎠, (53)
0 0 0.0105

which yields a negative Poisson ratio,  = −0.3740.

5.2.2. Maximization of a modified Poisson ratio. Here we set 1 = e1 ⊗e1 and 2 = e2 ⊗e2 so that

12 21
h(C) := (C−1 )1122 = − =− . (54)
E1 E2

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 749

Figure 9. Maximization of a modified Poisson ratio: (a) optimized RVE and

(b) corresponding periodic microstructure.

The final topology, obtained in 32 iterations is depicted in Figure 9(a) alongside the resulting
periodic microstructure. The matrix form of the corresponding homogenized elasticity is
⎛ ⎞
0.0469 0.0368 0
⎜ ⎟
C = ⎝ 0.0368 0.0469 0 ⎠, (55)
0 0 0.0098

with a resulting Poisson ratio  = 0.7847.

5.2.3. The Poisson ratio minimization. Here we use 1 = e1 ⊗e1 and 2 = −e2 ⊗e2 in (46).
The resulting function h(C) is

(C−1 )1122 (C−1 )1122

h(C) := − − = −12 −21 . (56)
(C−1 )1111 (C−1 )2222
The obtained optimized topology is shown in Figure 10(a). The matrix representation of the
corresponding C is
⎛ ⎞
0.1939 −0.0669 0
⎜ ⎟
C = ⎝ −0.0669 0.1939 0 ⎠, (57)
0 0 0.0311

which results in the negative Poisson ratio  = −0.3452. In this case, the algorithm does not converge
well initially (the value of J stops decreasing significantly with  ≈ 7◦ ). This is probably due to an
ill-conditioning of the problem. However, by performing another step of uniform mesh refinement
we obtain a better convergence with a final topology having <5◦ and a Poisson ratio  = −0.4118.
The final resulting topology is depicted in Figure 11 together with the convergence history of the
cost function h(C). A total number of 73 iterations were required to achieve convergence with
 = 5 ◦ .

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
750 S. AMSTUTZ ET AL.

Figure 10. Minimization of the Poisson ratio.

(a) (b)

Figure 11. Minimization of the Poisson ratio: (a) final topology and
(b) convergence history of the cost function.

5.2.4. Poisson ratio maximization. In this last example we set 1 = e1 ⊗e1 and 2 = e2 ⊗e2 and
the resulting function h(C) is

(C−1 )1122 (C−1 )1122

h(C) := + = 12 +21 . (58)
(C−1 )1111 (C−1 )2222
The converged optimized topology of the RVE, attained in 35 iterations, is shown in Figure 12(a).
The associated homogenized elasticity tensor is given by
⎛ ⎞
0.1565 0.1363 0
⎜ ⎟
C = ⎝ 0.1363 0.1565 0 ⎠, (59)
0 0 0.1162
which gives a Poisson ratio  = 0.8711.

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 751

Figure 12. Maximization of the Poisson ratio: (a) optimized RVE topology and
(b) corresponding periodic microstructure.

6. CONCLUSION

An algorithm for the topological design of periodic microstructures has been proposed. The algo-
rithm relies crucially on an exact formula for the topological derivative of the homogenized elasticity
tensor and a level set domain representation. The homogenized elasticity tensor is estimated by a
well-established multi-scale constitutive theory in which the macroscopic stress and strain tensors
are volume averages of their microscopic counterparts over a RVE. Its analytical topological deriva-
tive is a fourth-order tensor field over the RVE that provides a rigorous first-order approximation (in
the volume fraction of inclusion) to the change in the homogenized elasticity tensor resulting from
the insertion of a circular inclusion of a given phase contrast parameter. The proposed algorithm
was successfully and efficiently used in several numerical examples of optimum topology design
of bi-material periodic microstructures. We remark that the algorithm is of simple implementation
and, in sharp contrast to existing microstructural optimization methods, does not require the use
of artificial algorithmic parameters or strategies. This is only a natural consequence of the use
of the concept of topological derivative, which provides the correct mathematical framework to
treat problems involving singular changes of topology such as those present in microstructural
optimization problems. The application of the proposed methodology to other types of multi-scale
models is currently in progress and will be reported in forthcoming publications.

APPENDIX A: TOPOLOGICAL DERIVATIVE CALCULATION

For completeness we present here the main steps of the derivation of the topological asymptotic
expansion of the homogenized elasticity tensor for topological perturbations in the form of circular
inclusions. The derivation presented here offers an alternative to that recently reported in [50] and
the final result is an extension of the formula derived in [48] where perturbations in the form of
voids were considered instead.

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
752 S. AMSTUTZ ET AL.

In view of expressions (11) and (18) we have that the difference between the i jkl-components
of the tensors C and C is given by
 
 1   1
(C −C)ijkl =  (u kl )·(ei ⊗e j )−  (u kl )·(ei ⊗e j ), (A1)
V  V 

or, by making use of some simple tensorial relations,

 
1   1
(C −C)ijkl =  (u kl )·∇ s [(ei ⊗e j )y]−  (u kl )·∇ s [(ei ⊗e j )y]. (A2)
V  V 

In view of the additive decomposition (13) and (20) of the displacement fluctuation fields, we have
the identities
 
∇ s [(ei ⊗e j )y] = ∇ s (u ij −
u ij ) = ∇ s (u ij −
u ij ), (A3)
which replaced back into (A2) yield
 
 1   1  
(C −C)ijkl =  (u kl )·∇ (u ij −
s
u ij )−  (u kl )·∇ s (u ij −
u ij ). (A4)
V  V 

u ij and 
By taking into account that both  u ij belong to U and making use of the equilibrium
equations (12) and (19) together with the constitutive relations (2) and (17), expression (A4)
reduces to
 
1   1 
(C −C)ijkl = C ∇ s u kl ·∇ s u ij − C ∇ s u kl ·∇ s u ij , (A5)
V  V 
   
or, equivalently, by using the symmetry relation C ∇ s u kl ·∇ s u ij = C ∇ s u ij ·∇ s u kl , or (C)ijkl =
(C)klij ,

1  
(C −C)ijkl = [ (u ij )− (u ij )]·∇ s u kl . (A6)
V 

Finally, by taking into account the definitions of the original and perturbed microscopic domains
and their corresponding constitutive relations, the above expression can be equivalently written in
terms of an integral over the perturbation I :

 −1 
(C −C)ijkl =  (u ij )·∇ s u kl . (A7)
V I

The first-order asymptotic expansion of the above quantity for an arbitrarily shaped inclusion
can be obtained in a similar way as in [42]. The present problem, however, differs from that of
[42] in the functional spaces and boundary conditions involved. Thus, we briefly present in what
follows the main steps of the derivation. For simplicity we restrict ourselves to circular inclusions.
We start by re-writing (A7) in the equivalent form
  
 −1 s 
(C −C)ijkl =  (u ij )·∇ u kl +
s
 (u ij )·∇ (u kl −u kl ) , (A8)
V I I

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 753

and then replace the first integral with an approximation to obtain

  
 −1 s 
(C −C)ijkl =
  (u ij )(
2 s
y)·∇ u kl (
y)+E1 ()+  (u ij )·∇ (u kl −u kl ) . (A9)
V I

Here and in what follows Ei () denotes remainders whose behavior will be discussed later. Now
 
we consider the solution  (wkl ) of the exterior problem defined by
 
−div( (wkl )) = 0 in I ∪(R2 \H ),
 
[[ (wkl )]]n = −(1− ) (u kl )(
y)n on *I , (A10)
 
 (wkl ) →0 at ∞,

where [[(·)]] denotes the jump of the function (·) across the matrix/inclusion interface *I :

[[(·)]] := (·)|m −(·)|i , (A11)

with subscripts m and i associated with quantities evaluated on the matrix and inclusion boundaries,
respectively. Then, by considering the symmetry of the constitutive tensor in (2) and taking the
 
solution of (A10) as an estimate for the variation  (u kl −u kl ), we obtain the following for the
rightmost integral in (A8):
 
 
 (u ij )·∇ s (u kl −u kl ) = ∇ s u ij · (u kl −u kl )
I I


= ∇ s u ij · (wkl )+E2 ()
I


= ∇ u ij (
s
y)·  (wkl )+E2 ()+E3 (), (A12)
I

where an extra approximation involving ∇ s u ij has been introduced to arrive at the third right-hand
side.

In the present case of a circular inclusion, the nominal stress tensor  (wkl ) has the following
representation in a polar coordinate system (r, ) with origin at the centre 
y of the inclusion:
• for r 
 2 
1 1− 2 1 I II 1−  4
r (r, ) = − (I +II ) − ( − ) 4 −3 cos 2, (A13)
2 1+ r 2 2  
1+ r2 r4
1 1− 2 3 I II 1− 
4
 (r, ) = (I +II ) − ( − ) cos 2, (A14)
2 1+ r 2 2  
1+ r 4
 2 
r 1 I II 1−  4
 (r, ) = − ( − ) 2 2 −3 4 sin 2; (A15)
2 1+ r r

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
754 S. AMSTUTZ ET AL.

• for 0<r <

1 1− 1 1−
r (r, ) = (I +II ) + (I −II )  cos 2, (A16)
2 1+ 2 1+
1 1− 1 1−
 (r, ) = (I +II ) − (I −II )  cos 2, (A17)
2 1+ 2 1+
1 1−
r (r, ) = − (I −II )  sin 2, (A18)
2 1+

where I and II are the eigenvalues of the tensor  (u kl )(
y). In addition, the constants  and 
are given, respectively, by
1+ 3−
= , = . (A19)
1− 1+

Representation (A13)–(A18) for the nominal stress  (wkl ) allows the integral term in (A12) to
be analytically calculated and results in


 (wkl ) = −
2 T (u kl )(
y), (A20)
I

with the scalar value  and the fourth-order tensor T given by:
1− 1 −
=− and T = I+ I⊗I. (A21)
1+ 2 1+ 
By replacing (A20) in (A12) and then substituting the result into (A9), we finally arrive at the
following estimate for the difference between the i jkl-components of the homogenized elasticity
tensors C and C of the perturbed and original RVEs:
 

2 1+ 1 −
(C −C)ijkl = (1− )  (u kl )·∇ s u ij + tr[ (u kl )]tr(∇ s u ij ) (
y)
V 1− 2 1+ 


3
+( −1) Ei (). (A22)
i=1

By following analogous steps to those of [42] it can be shown that


3
Ei () = o(2 ). (A23)
i=1

Then, by comparing (A22), (A23) with the topological asymptotic expansion (21) we identify the
function f according to (22) and the topological derivative of the tensor C as
 
1+ 1 −
(DT C)ijkl = (1− )  (u kl )·∇ u ij +
s s
tr( (u kl ))tr(∇ u ij ) (A24)
1− 2 1+ 
with  and  given by (25). The topological derivative of C can be expressed in compact form
by (23).

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
 MULTI-SCALE LINEAR ELASTICITY MODELS 755

ACKNOWLEDGEMENTS
The collaboration between the Swansea University and LNCC researchers was supported by The Royal
Society (International Joint Project no. JP080066). S.M. Giusti was supported by the Brazilian research
funding council, CNPq (grant no. 382485/2009-2). All support is gratefully acknowledged.

REFERENCES
1. Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behaviour of multiphase materials.
Journal of the Mechanics and Physics of Solids 1963; 11(2):127–140.
2. Hill R. A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids
1965; 13(4):213–222.
3. Mandel J. Plasticité classique et viscoplasticité. CISM Lecture Notes. Springer: Udine, 1971.
4. Bensoussan A, Lions JL, Papanicolau G. Asymptotic Analysis for Periodic Microstructures. North-Holland:
Amsterdam, 1978.
5. Sanchez-Palencia E. Non-Homogeneous Media and Vibration Theory. Lecture Notes in Physics, vol. 127. Springer:
Berlin, 1980.
6. Suquet PM. Elements of homogenization for inelastic solid mechanics. Homogenization Techniques for Composite
Media. Lecture Notes in Physics, vol. 272. Springer: Berlin, 1987.
7. Auriault JL. Effective macroscopic description for heat conduction in periodic composites. International Journal
of Heat and Mass Transfer 1983; 26(6):861–869.
8. Auriault JL, Royer P. Double conductivity media: a comparison between phenomenological and homogenization
approaches. International Journal of Heat and Mass Transfer 1993; 36(10):2613–2621.
9. Gurson AL. Continuum theory of ductile rupture by void nucleation and growth: part I—yield criteria and flow
rule for porous ductile media. Journal of Engineering Materials and Technology, Transactions on ASME 1977;
99(1):2–15.
10. Nemat-Nasser S. Averaging theorems in finite deformation plasticity. Mechanics of Materials 1999; 31(8):493–523.
11. Nemat-Nasser S, Hori M. Micromechanics: Overall Properties of Heterogeneous Materials. North-Holland:
Amsterdam, 1993.
12. Ostoja-Starzewski M, Schulte J. Bounding of effective thermal conductivities of multiscale materials by essential
and natural boundary conditions. Physical Review B 1996; 54(1):278–285.
13. Michel JC, Moulinec H, Suquet P. Effective properties of composite materials with periodic microstructure: a
computational approach. Computer Methods in Applied Mechanics and Engineering 1999; 172(1–4):109–143.
14. Miehe C, Schotte J, Schröder J. Computational micro-macro transitions and overall moduli in the analysis of
polycrystals at large strains. Computational Materials Science 1999; 16(1–4):372–382.
15. Speirs DCD, de Souza Neto EA, Perić D. An approach to the mechanical constitutive modelling of arterial tissue
based on homogenization and optimization. Journal of Biomechanics 2008; 41(12):2673–2680.
16. Oyen ML, Ferguson VL, Bembey AK, Bushby AJ, Boyde A. Composite bounds on the elastic modulus of bone.
Journal of Biomechanics 2008; 41(11):2585–2588.
17. Giusti SM, Blanco PJ, de Souza Neto EA, Feijóo RA. An assessment of the Gurson yield criterion by a
computational multi-scale approach. Engineering with Computers 2009; 26(3):281–301.
18. Celentano DJ, Dardati PM, Godoy LA, Boeri RE. Computational simulation of microstructure evolution during
solidification of ductile cast iron. International Journal of Cast Metals Research 2008; 21(6):416–426.
19. Almgreen RF. An isotropic three-dimensional structure with Poisson’s ratio—1. Journal of Elasticity 1985;
15(4):427–430.
20. Lakes R. Foam structures with negative Poisson’s ratio. Science, AAAS 1987; 235(4792):1038–1040.
21. Lakes R. Negative Poisson’s ratio materials. Science, AAAS 1987; 238(4826):551.
22. Eschenauer HA, Olhoff N. Topology optimization of continuum structures: a review. Applied Mechanics Reviews
2001; 54(4):331–390.
23. Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using an homogenization method.
Computer Methods in Applied Mechanics and Engineering 1988; 71(2):197–224.
24. Żochowski A. Optimal perforation design in 2-dimensional elasticity. Mechanics of Structures and Machines
1988; 16(1):17–33.
25. Belytschko T, Xiao SP, Parimi C. Topology optimization with implicit functions and regularization. International
Journal for Numerical Methods in Engineering 2003; 57:1177–1196.

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme
756 S. AMSTUTZ ET AL.

26. Kikuchi N, Nishiwaki S, Fonseca JSO, Silva ECN. Design optimization method for compliant mechanisms and
material microstructure. Computer Methods in Applied Mechanics and Engineering 1998; 151(3–4):401–417.
27. Sigmund O. Materials with prescribed constitutive parameters: an inverse homogenization problem. International
Journal of Solids and Structures 1994; 31(17):2313–2329.
28. Silva ECN, Fonseca JSO, Kikuchi N. Optimal design of periodic microstructures. Computational Mechanics
1997; 19(5):397–410.
29. Sokołowski J, Żochowski A. On the topological derivative in shape optimization. SIAM Journal on Control and
Optimization 1999; 37(4):1251–1272.
30. Céa J, Garreau S, Guillaume Ph, Masmoudi M. The shape and topological optimizations connection. Computer
Methods in Applied Mechanics and Engineering 2000; 188(4):713–726.
31. Eschenauer HA, Kobelev VV, Schumacher A. Bubble method for topology and shape optimization of structures.
Structural Optimization 1994; 8(1):42–51.
32. Germain P, Nguyen QS, Suquet P. Continuum thermodynamics. Journal of Applied Mechanics, Transactions
ASME 1983; 50(4):1010–1020.
33. Allaire G, de Gournay F, Jouve F, Toader A. Structural optimization using topological and shape sensitivity via
a level set method. Control and Cybernetics 2005; 34(1):59–80.
34. Amstutz S, Andrä H. A new algorithm for topology optimization using a level-set method. Journal of
Computational Physics 2006; 216(2):573–588.
35. Novotny AA, Feijóo RA, Taroco E, Padra C. Topological sensitivity analysis for three-dimensional linear elasticity
problem. Computer Methods in Applied Mechanics and Engineering 2007; 196(41–44):4354–4364.
36. Amstutz S, Horchani I, Masmoudi M. Crack detection by the topological gradient method. Control and Cybernetics
2005; 34(1):81–101.
37. Feijóo GR. A new method in inverse scattering based on the topological derivative. Inverse Problems 2004;
20(6):1819–1840.
38. Hintermüller M, Laurain A. Electrical impedance tomography: from topology to shape. Control and Cybernetics
2008; 37(4):913–933.
39. Aurous D, Masmoudi M, Belaid L. Image restoration and classification by topological asymptotic expansion. In
Variational Formulations in Mechanics: Theory and Applications, Taroco E, de Souza Neto EA, Novotny AA
(eds). CIMNE: Barcelona, Spain, 2007.
40. Hintermüller M, Laurain A. Multiphase image segmentation and modulation recovery based on shape and
topological sensitivity. Journal of Mathematical Imaging and Vision 2009; 35:1–22.
41. Larrabide I, Feijóo RA, Novotny AA, Taroco E. Topological derivative: a tool for image processing. Computers
and Structures 2008; 86(13–14):1386–1403.
42. Amstutz S. Sensitivity analysis with respect to a local perturbation of the material property. Asymptotic Analysis
2006; 49(1–2):87–108.
43. Garreau S, Guillaume Ph, Masmoudi M. The topological asymptotic for pde systems: the elasticity case. SIAM
Journal on Control and Optimization 2001; 39(6):1756–1778.
44. Nazarov SA, Sokołowski J. Asymptotic analysis of shape functionals. Journal de Mathematiques Pures
et Appliquees 2003; 82(2):125–196.
45. Turevsky I, Gopalakrishnan SH, Suresh K. An efficient numerical method for computing the topological sensitivity
of arbitrary-shaped features in plate bending. International Journal for Numerical Methods in Engineering 2009;
79:1683–1702.
46. de Souza Neto EA, Feijóo RA. Variational foundations of multi-scale constitutive models of solid: small and large
strain kinematical formulation. Technical Report No. 16/2006, Laboratório Nacional de Computação Científica
LNCC/MCT, Petrópolis, Brazil, 2006.
47. de Souza Neto EA, Feijóo RA. Axiomatic variational foundations of multi-scale solid constitutive models:
kinematical formulation. In Advanced Computational Materials Modelling. From Classical to Multi-Scale
Techniques, Vaz Jr M, de Souza Neto EA, Muñoz-Rojas P (eds). Wiley: Chichester, 2010.
48. Giusti SM, Novotny AA, de Souza Neto EA, Feijóo RA. Sensitivity of the macroscopic elasticity tensor to
topological microstructural changes. Journal of the Mechanics and Physics of Solids 2009; 57(3):555–570.
49. Osher S, Sethian JA. Front propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi
formulations. Journal of Computational Physics 1988; 78:12–49.
50. Giusti SM, Novotny AA, de Souza Neto EA. Sensitivity of the macroscopic response of elastic microstructures to
the insertion of inclusions. Proceedings of the Royal Society of London A 2009; DOI: 10.1098/rspa.2009.0499.

Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:733–756
DOI: 10.1002/nme