## Abstract

This paper proposes an exact analytical formula for the topological sensitivity of the macroscopic response of elastic microstructures to the insertion of circular inclusions. The macroscopic response is assumed to be predicted by a well-established multi-scale constitutive theory where the macroscopic strain and stress tensors are defined as volume averages of their microscopic counterpart fields over a representative volume element (RVE) of material. The proposed formula—a symmetric fourth-order tensor field over the RVE domain—is a topological derivative which measures how the macroscopic elasticity tensor changes when an infinitesimal circular elastic inclusion is introduced within the RVE. In the limits, when the inclusion/matrix phase contrast ratio tends to zero and infinity, the sensitivities to the insertion of a hole and a rigid inclusion, respectively, are rigorously obtained. The derivation relies on the topological asymptotic analysis of the predicted macroscopic elasticity and is presented in detail. The derived fundamental formula is of interest to many areas of applied and computational mechanics. To illustrate its potential applicability, a simple finite element-based example is presented where the topological derivative information is used to automatically generate a bi-material microstructure to meet pre-specified macroscopic properties.

## 1. Introduction

The ability to predict the macroscopic constitutive response of materials from the knowledge of their microstructure has long been a subject of intensive research in applied and computational mechanics circles. The large body of publications currently available in this field ranges from early fundamental work (e.g. Hashin & Shtrikman 1963; Hill 1965; Mandel 1971; Bensoussan *et al.* 1978; Sanchez-Palencia 1980; Germain *et al.* 1983; Suquet 1987) to more recent numerical simulations mainly based on the finite element method (Michel *et al.* 1999; Miehe *et al.* 1999; Terada *et al.* 2003; Matsui *et al.* 2004; Speirs *et al.* 2008).

Among the various applications of such so-called multi-scale constitutive theories, of particular interest is the study of the sensitivity of the macroscopic response to changes in the microstructure. In this context, the use of sensitivity information has been successfully applied in the topological optimization of microstructures, among others, by Sigmund (1994), Silva *et al.* (1997), Kikuchi *et al.* (1998) and Hyun & Torquato (2001). The approach adopted by these authors relies essentially on the regularization of the actual representative volume element (RVE) topology optimization problem by introducing a fictitious density field of which the elastic material parameters are assumed to be a function. Throughout the optimization iterations, voids are assumed to be located wherever the fictitious density field (the design variable in the regularized optimization problem) falls below a given numerical tolerance. The sensitivity of the macroscopic elastic parameters to topological changes of the RVE (the introduction of voids in this case) is calculated only in an approximate sense. The approximate character of the sensitivity calculation stems from the fact that procedures of this type are based on the conventional concept of derivative, and the introdution of topological changes such as voids is inherently singular and does not fit within the conventional notion of differentiability.

In the present paper, we propose an exact analytical formula for the sensitivity of the macroscopic elasticity tensor to topological microstructural changes. The macroscopic elasticity tensor is assumed to be predicted by a well-established multi-scale constitutive theory based on the volume averaging of the stress and strain tensors over the RVE. Within this constitutive framework, upper and lower bounds for the elastic behaviour can be obtained by assuming, respectively, linear RVE boundary displacements and minimum RVE kinematical constraint compatible with the strain averaging assumption (and resulting in uniform RVE boundary tractions). In addition, the widely used assumption of periodic RVE boundary displacement fluctuations provides an estimate for the response of periodic media. The topological change considered consists of the insertion of a circular isotropically elastic inclusion within the isotropically elastic RVE matrix. Young’s modulus of the inclusion is assumed to be a scalar multiple of Young’s modulus of the matrix. In the limits when the scalar multiple (the phase contrast parameter) tends to zero and infinity, the proposed formula gives the exact topological sensitivity of the macroscopic elasticity tensor to the insertion of a hole and a rigid inclusion, respectively. The derivation of the analytical sensitivity is presented in detail. It relies on the concepts of topological asymptotic expansion and topological derivative (Sokołowski & Żochowski 1999; Céa *et al.* 2000) which provides the correct mathematical framework for the treatment of singularities of the present type. The derived formula is of great potential use in applied and computational mechanics and its final format is remarkably simple. Its potential applicability is illustrated in a finite element-based example where the topological derivative information is used in a very simple algorithm to automatically generate a bi-material microstucture to meet a pre-specified macroscopic behaviour.

This paper is organized as follows. The multi-scale constitutive framework adopted in the estimation of the macroscopic elasticity tensor is briefly reviewed in §2. The main contribution of this paper is presented in §3 with a detailed derivation of the proposed sensitivity formula. The application of the topological derivative to the synthesis of microstructures is shown in §4 and some concluding remarks are made in §5.

## 2. Multi-scale constitutive modelling

The family of constitutive theories upon which we rely for the estimation of the macroscopic elastic properties has been formally presented in a rather general setting by Germain *et al.* (1983) and later exploited, among others, by Michel *et al.* (1999) and Miehe *et al.* (1999) in the computational context. When applied to linearly elastic periodic media, it coincides with the asymptotic expansion-based theory described by Bensoussan *et al.* (1978) and Sanchez-Palencia (1980).

The fundamental idea here is the assumption that any point ** x** of the macroscopic continuum 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,

*Ω*. An axiomatic variational framework for this family of constituive theories is presented in detail by de Souza Neto & Feijóo (2006). Accordingly, the entire theory can be derived from five basic principles: (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 of displacement fluctuations compatible with the strain averaging hypothesis; (iii) the equilibrium of the RVE; (iv) the stress averaging relation; and (v) the Hill–Mandel principle of macro-homogeneity (Hill 1965; Mandel 1971), ensuring the energy consistency between the two scales.

The strain averaging assumption defines the macroscopic strain tensor **E** at a point ** x**∈

*Ω*as the volume average of its microscopic counterpart

**E**

_{μ}over the RVE: 2.1 where

*V*

_{μ}is a total volume of the RVE and

**E**

_{μ}:=∇

^{s}

**u**

_{μ}, with

**u**

_{μ}denoting the microscopic displacement field of the RVE.

For the purposes of the present paper, we consider RVEs whose domain consists of a matrix, , with inclusions of different materials occupying a domain . Futher, we assume the behaviour of the matrix and inclusion to be isotropic linear elastic. The microscopic stress tensor field **T**_{μ} is given by
2.2
where is the fourth-order isotropic elasticity tensor
2.3
with **I** and denoting, respectively, the second- and fourth-order identity tensors and *E* and *ν* Young’s modulus and Poisson’s ratio, given by
2.4
*E*_{i} and *ν*_{i} are assumed constant within each inclusion but may in general vary from inclusion to inclusion.

Without loss of generality, **u**_{μ} may be decomposed as a sum
2.5
of a rigid displacement coinciding with the displacement **u** of ** x**, a field linear in the local RVE coordinate

**(whose origin is assumed without loss of generality to be located at the centroid of the RVE) and a displacement fluct- uation field that, in general, varies with**

*y***. The above renders the split of**

*y***E**

_{μ}, 2.6 as a sum of a homogeneous strain over the RVE equal to the macroscopic strain and a strain fluctuation, , about the average value. Similarly, we have 2.7 where is the stress field associated with the uniform strain and is the stress fluctuation field.

The above considerations together with equations (2.1), (2.2), axioms (ii) and (iii) and the Hill–Mandel principle of macro-homogeneity (Hill 1965) lead to the mechanical equilibrium problem for the RVE which consists of finding, for a given macroscopic strain **E**, a kinematically admissible displacement fluctuation field , such that
2.8
where the (as yet not defined) space of admissible displacement fluctuations (and virtual admissible displacements) of the RVE is a subspace of —the minimally constrained space of kinematically admissible displacement fluctuation fields compatible with equation (2.1)
2.9
with **n** denoting the outward unit normal to the boundary ∂*Ω*_{μ}, ⊗_{s} the symmetric tensor product and [[** v**]] the

*jump*of

**across the matrix/inclusion interface 2.10 With the solution to problem (2.8) at hand, the macroscopic stress is obtained according to stress averaging relation 2.11**

*v*### (a) Classes of multi-scale constitutive models

The characterization of a multi-scale model of the present type is completed with the choice of a suitable space of admissible displacement fluctuations which defines the kinematical constraints assumed over the RVE. The three classical choices are

—

*Linear boundary displacements*. For this class of models, the choice is 2.12—

*Periodic boundary fluctuations*. This assumption is typical of periodic media. The geometry of the RVE here must satisfy certain constraints not needed by the other two classes of models discussed. For polygonal RVE geometries we have that the boundary ∂*Ω*_{μ}is composed of a number of pairs of equally sized subsets {*Γ*^{+}_{i},*Γ*^{−}_{i}} with normals**n**^{+}_{i}=−**n**^{−}_{i}. For each pair {*Γ*^{+}_{i},*Γ*^{−}_{i}} of sides, there is a one-to-one correspondence between points*y*^{+}∈*Γ*^{+}_{i}and*y*^{−}∈*Γ*^{−}_{i}. The periodicity of the structure requires that the displacement fluctuation at any point*y*^{+}coincide with that of the corresponding point*y*^{−}. Hence, the space of displacement fluctuations is defined as 2.13—

*Minimally constrained or uniform boundary traction model*. In this case, 2.14 It can be shown (de Souza Neto & Feijóo 2006) that this assumption results in uniform RVE boundary tractions 2.15

It should be noted that the following holds for the above definitions: 2.16 Then, in the RVE equilibrium problem (2.8), the assumption of linear boundary displacements (the most kinematically constrained of the three) results in the least compliant solution while the assumption of minimal kinematical constraint leads to the most compliant one. In this sense, the choices and provide, respectively, an upper and a lower bound for the estimated macroscopic stiffness of the material.

### (b) The homogenized elasticity tensor

The assumed behaviour of the microscale implies that the macroscopic response is linear elastic. That is, there is a *homogenized elasticity tensor* such that
2.17
A closed form for can be derived, as in Michel *et al.* (1999), based on the representation of equation (2.8) as a superposition of linear problems associated with the Cartesian components of **E**. The resulting formula can be expressed as
2.18
where is the volume average of the microscopic elasticity tensor, and the (generally dependent upon the choice of space ) is given by
2.19
where is the fluctuation stress field associated with the fluctuation displacement field that solves linear variational problem
2.20
for *i*,*j*=1,2 (in the two-dimensional case). In the above, {**e**_{i}} denotes an orthonormal basis for the two-dimensional Euclidean space.

## 3. The topological sensitivity of the homogenized elasticity tensor

In this section, we present the main result of this paper—a closed formula for the sensitivity of the homogenized elasticity tensor (2.18) to the introduction of an infinitesimal circular inclusion centred at an arbitrary point of the RVE domain. The formula is obtained by following analogous steps to those presented by Giusti *et al.* (2009) for the introduction of a circular void, but contains fundamental differences that justify the presentation of the details of its derivation in the following.

We start by providing a brief introduction to the relatively new mathematical concepts of *topological asymptotic expansion*and *topological derivative*. To this end, let *ψ* be a functional whose value depends on a given domain and let *ψ* have sufficient regularity so that the following expansion is possible:
3.1
where *ψ*(0) is the value of the functional for an original (unperturbed) domain and *ψ*(*ε*) denotes the value of the functional for a domain that differs from the original one by a topological perturbation of size *ε*. Note that the original domain is retrieved when *ε*=0. In addition, *f*(*ε*) is a *regularizing function* defined such that with and *o*(*f*(*ε*)) contains all terms of higher order in *f*(*ε*). The right-hand side of equation (3.1) is named the *asymptotic topological expansion* of the functional *ψ* and the term *D*_{T}*ψ* is defined as the *topological derivative* of *ψ* at the unperturbed RVE domain.

The concept of topological derivative was rigorously introduced by Sokołowski & Żochowski (1999). Since then, the notion of topological derivative has proved extremely useful in the treatment of a wide range of problems in mechanics, optimization, inverse analysis and image processing, and has become a subject of intensive research (see, for instance, Céa *et al.* 2000; Garreau *et al.* 2001; Amstutz *et al.* 2005).

### (a) Application to the multi-scale elasticity model

To begin the topological sensitivity analysis in the present context, it is appropriate to define the following functional:
3.2
where **T**^{ε} denotes the macroscopic stress tensor associated with an RVE topologically perturbed by a small inclusion of radius *ε* defined by and **T** is the macroscopic stress tensor associated with the unperturbed domain *Ω*_{μ}. More precisely, the perturbed domain is obtained when a circular hole of radius *ε* is introduced at an arbitrary point and, then, this region is replaced with the circular inclusion made of a different material. Then, the perturbed domain is defined as (refer to figure 1). The asymptotic topological expansion of the functional (3.2) reads
3.3

### (b) Topological derivative calculation

Our purpose here is to derive the closed formula for the topological sensitivity of the macroscopic elasticity tensor (2.18). Then, we start deriving a closed formula for the associated topological derivative *D*_{T}*ψ*, which characterizes the asymptotic expansion (3.3). To this end, we re-define the functional
3.4
where **T**_{με} is the microscopic stress field of the perturbed domain *Ω*_{με}. The microscopic stress field **T**_{με} is given by
3.5
with **E**_{με}=∇^{s}**u**_{με} denoting the microscopic strain field in *Ω*_{με} and constitutive fourth-order tensor given by
3.6
where *γ*∈ℜ^{+} is the contrast parameter defining the ratio between the properties of the original and new material at the location of the perturbation. Clearly, this type of perturbation corresponds to a change only in Young’s modulus of the phases. We remark that the particular shape functional (3.4) subjected to the topological perturbation (3.6) has the sufficient degree of regularity for the existence of the topological derivative. For further details on the existence of the topological derivative for energy shape functionals associated with elastic systems, we refer to Sokołowski & Żochowski (1999) and Nazarov & Sokołowski (2003).

The microscopic displacement field of the perturbed RVE is decomposed as
3.7
where the displacement fluctuation field is the solution of the following variational problem in *Ω*_{με}: Find such that
3.8
where is the space of kinematically admissible displacement fluctuations of the perturbed RVE and is the microscopic stress field, associated with *Ω*_{με}, induced by the macroscopic strain **E**, i.e. .

For the calculation of the topological derivative, we shall adopt the approach presented by Sokołowski & Żochowski (2001) and Novotny *et al.* (2003), whereby the topological derivative is obtained as
3.9

The derivative of the functional with respect to the perturbation parameter *ε* can be seen as the sensitivity of , in the classical sense, to the change in shape produced by a uniform expansion of the inclusion .

### Proposition 3.1.

*Let* *be the functional defined by equation (3.4). Then, the derivative of the functional* *with respect to the small parameter ε is given by*
3.10
*where* *v**is the RVE shape-change velocity field defined in Ω*

_{με}

*and*

**Σ**

_{με}

*is a generalization of the classical Eshelby momentum–energy tensor (*Eshelby 1975

*;*Gurtin 2000

*) of the RVE, given in the present case by*3.11

### Proof.

By making use of Reynolds’ transport theorem (Sokołowski & Zolésio 1992), we obtain the identity
3.12
Next, by using the concept of material derivative of a spatial field, we find that the first term of the above right-hand side integral can be written as
3.13
where the superimposed dot denotes the (total) material derivative with respect to *ε*. Further, note that
3.14
Then, by introducing the above expression into equation (3.13), we obtain
3.15
which, substituted in equation (3.12), gives
3.16
where we have made use of the identity *d**i**v* **v=I**⋅∇**v**. Now, note that by definition of the spaces of virtual displacements, we have This, together with the equilibrium equation (3.8), implies that the first term of equation (3.16) vanishes. Then, a straightforward rearrangement of the above yields equation (3.10). ▪

### Proposition 3.2.

*Let* *be the functional defined by equation (3.4). Then, the derivative of the functional* *with respect to the small parameter ε can be written as*
3.17
*where* *v**is the RVE shape-change velocity field and* **Σ**_{με}*is given by equation (3.11).*

### Proof.

Let us compute the shape derivative of the functional using the following version for the Reynolds’ transport theorem (Sokołowski & Zolésio 1992):
3.18
Next, by using the concept of spatial derivative and equation (2.2), we find that the first term of the above right-hand side integral can be written as
3.19
where the prime denotes the (partial) spatial derivative with respect to *ε*. Further, note that the relation (3.14) gives
3.20
Then, by introducing the above expression into equation (3.19), we obtain
3.21
With the above result, the sensitivity of the functional reads
3.22
Now, note that by definition of the spaces of virtual displacements, we have This, together with the equilibrium equation (3.8), implies that the first term of equation (3.22) vanishes. Then, we obtain
3.23
In view of the tensor relation
3.24
and the divergence theorem, expression (3.23) can be written as
3.25
Finally, as the stress field **T**_{με} is in equilibrium, we have that *d**i**v* **T**_{με}=**0** in *Ω*_{με} so that a straightforward rearrangement of the above yields equation (3.17). ▪

### Corollary 3.3.

*By applying the divergence theorem to the right-hand side of equation (3.10 ), we obtain*
3.26
*As equations (3.17 ) and (3.26 ) remain valid for all velocity fields* , *we have*
3.27
*i.e.* **Σ**_{με} *is a divergence-free field.*

From the above corollary, it follows that the shape derivative of the functional , for the particular case of a uniform expansion of the perturbation—circular inclusion—can be expressed exclusively in terms of an integral over the boundary of the inclusion 3.28

In order to derive an explicit expression for the integrand on the right-hand side of equation (3.28), we consider a curvilinear coordinate system along , characterized by the orthonormal vectors **n** and **t**. Then, we can decompose the stress tensor **T**_{με} and the strain tensor **E**_{με} on the boundary as follows:
3.29

The Neumann boundary condition along , together with equation (3.29)_{1}, gives
3.30

Similar to equation (3.29), the fluctuation displacement field can be decomposed on as 3.31 As a consequence, the continuity condition of along implies 3.32

Alternatively, the above condition can be written in terms of the components of the fluctuation strain tensor in the basis {**n**,**t**} as follows:
3.33
In addition, in view of the additive split (2.6), condition (3.33) establishes the continuity of component *t**t* of the strain tensor **E**_{με}, i.e.
3.34

By taking into account the decompositions (3.29) and (3.31) and the continuity condition (3.30), (3.32) and (3.34), the jump of the Eshelby tensor flux in the normal direction to the boundary of the perturbation can be written as 3.35

Note that, by using the constitutive law (3.5), the jump terms on the right-hand side of the above expression satisfy
3.36
3.37and
3.38
where , , and are the constant and fluctuation part of components and of the stress tensor given by equation (3.29)_{1}.

By introducing the above results into equation (3.35) and taking into account the additive decomposition of the components of the microscopic stress field **T**_{με}, we find that the jump of the Eshelby tensor flux in the normal direction to the boundary has the following representation in terms of the solution inside the perturbation :
3.39

In order to obtain an analytical formula for the boundary integral (3.28), we make use of the classical asymptotic analysis for two-dimensional elasticity problems (see appendix A). Thus, the distribution of the microscopic stress field on boundary can be written as
3.40
with as and the fourth-order tensors and given by
3.41and
3.42
where the constants *α* and *β* are defined as
3.43

With the stress distribution along the boundary shown in equation (3.40) and the result (3.39), we can obtain the topological derivative by evaluating the boundary integral (3.28) analytically. This gives
3.44
Finally, by substituting the above in equation (3.9) and adopting the function *f*(*ε*) as the size (area) of the circular perturbation, we obtain the explicit closed form expression for the topological derivative of *ψ*
3.45
where the fourth-order tensor is defined as
3.46

### (c) The sensitivity of the macroscopic elasticity tensor

Expressions (3.3) and (3.45) promptly lead to the explicit formula for the topological asymptotic expansion of *ψ*
3.47
where *v*(*ε*):=*π**ε*^{2}/*V*_{μ} is the RVE volume fraction of perturbation.

The approach used in §2*b* to obtain a closed form expression for the macroscopic elasticity tensor can be easily extended to derive an analytical formula for its topological sensitivity—the main result of this paper. Accordingly, we write **E**_{μ} as a linear combination of the Cartesian components of **E**
3.48
With the introduction of the above expression in equation (2.2), the microscopic stress tensor **T**_{μ} can be written as
3.49
where **T**_{μij} denotes the microscopic stress field associated with each displacement fluctuation field that solves the variational equation (2.20). Then, by combining equations (3.45) and (3.49), we obtain the following alternative formula for the topological derivative of *ψ*:
3.50
where is the fourth-order symmetric tensor field over *Ω*_{μ} defined by
3.51

Our main result—the asymptotic expansion of the macroscopic elasticity tensor and the corresponding sensitivity formula—is presented in the following. From equations (3.45), (3.47) and (3.50) and by making use of the fact that the macroscopic constitutive response is linear elastic, we obtain
3.52
where is the macroscopic elasticity tensor of the topologically perturbed microstructure, i.e.
3.53
Finally, as equation (3.52) is valid for any **E**, we arrive at explicit expression for the topological expansion of the macroscopic elasticity tensor
3.54

### Remark 3.4.

The topological sensitivity tensor field over *Ω*_{μ}, whose expli- cit expression is given by equation (3.51), provides a first-order accurate measure of how the macroscopic elasticity tensor varies when a topological perturbation (a circular inclusion) is added to the RVE. The value of each Cartesian component at an arbitrary point ** y**∈

*Ω*

_{μ}represents the derivative of the component

*i*

*j*

*k*

*l*of the macroscopic elasticity tensor with respect to the volume fraction

*v*(

*ε*) of a circular inclusion of radius

*ε*inserted at

**. The remarkable simplicity of the closed form sensitivity given by equation (3.51) is to be noted. Once the vector fields have been obtained as solutions of equation (2.20) for the**

*y**original*RVE domain, the sensitivity tensor can be trivially assembled.

### Remark 3.5.

The topological derivative tensor is dependent upon the contrast parameter *γ*. This dependency is made explicit in definition (3.46) of the tensor . Note that the limiting cases and correspond, respectively, to the insertion of a hole (infinitely compliant material) and a rigid inclusion (infinitely stiff material). In such cases, the sensitivity tensor is calculated by means of equation (3.51), respectively, with
3.55
and
3.56
Expression (3.55) coincides with the result derived in Giusti *et al.* (2009) for the sensitivity to the insertion of holes.

## 4. Example of application: microstructure topology synthesis

The explicit formula for the macroscopic elasticity sensitivity derived above is a fundamental result with potential application in different areas of interest. In this section, we show what is probably its most intuitive and straightforward application—the finite element-based automatic synthesis of a microstructure to meet a pre-specified macroscopic response.

Crucial to the proposed application is the definition of a measure of distance between two generic symmetric fourth-order tensors and . Here, we shall adopt simply the square of the Euclidean norm of the difference between the two tensors, i.e. 4.1

Then, for a pre-specified *target*constitutive tensor , we conveniently define the functions *ϕ*(*ε*) and *ϕ*(0), respectively, as
4.2
Next, note that by taking into account the asymptotic expansion (3.54) of the homogenized elasticity tensor and definitions (4.1) and (4.2), we obtain the following asymptotic topological expansion of the function *ϕ*(*ε*):
4.3
where we have identified the topological derivative of the function *ϕ* as
4.4

## Remark 4.1.

Note that the topological derivative of *ϕ* has been derived here without the need for a full topological asymptotic analysis. The derivation has been based on simply introducing equation (3.54) into equation (4.2)_{1} and then collecting the terms of same power of *ε*. This can be done for any function of the with the required degree of regularity.

### (a) A simple microstructural synthesis algorithm

With the above at hand, we are now ready to devise a simple algorithm to synthesize a bi-material microstructure whose macroscopic response is characterized by . Our task here is to produce a sequence of changes to a given initial microstructure topology in order to achieve a final topology for which the value of *ϕ* (a non-negative function) is as close as possible to zero (the zero value corresponds to a perfect match between the target and actual predicted response). A very simple procedure consists in requiring at each iteration of the algorithm that a small inclusion be introduced at points of the RVE where the field *D*_{T}*ϕ* is most negative. That is, we use the sensitivity information provided by the topological derivative field to obtain a feasible descent direction. This approach has been successfully used in the context of topological optimization of elastic structures. For further details, we refer to Giusti *et al.* (2008). The algorithm is summarized below:

—

*Provide*a complete description of the initial guess (initial topology) of the RVE with domain*Ω*_{μ}and the parameter*γ*defining the phase contrast; a target macroscopic constitutive response ; the maximum number of iterations allowed*N*and a numerical convergence tolerance*ϵ*.—

*While**ϕ*>*ϵ*and*j*≤*N*,*do*:*Compute*and*D*^{j}_{Tμ}*ϕ*in the domain*Change*the material properties of the elements for which*D*^{j}_{Tμ}*ϕ*is most negative. The maximum number*n*_{e}of elements whose properties are allowed to change at each iteration is a user-specified fixed parameter.*Set*and

In the above, is the RVE domain at iteration *j*; and and *D*^{j}_{Tμ}*ϕ* are the corresponding quantities evaluated for . The maximum volume fraction of material that may be changed in each iteration is limited by *n*_{e} and the size of the elements, but changes may occur simultaneously at different RVE locations.

### (b) Application: automatic synthesis of a bi-material microstructure

The above algorithm is used here to generate a two-dimensional bi-material microstructure under plane stress conditions. The (normalized) Young modulus of the matrix and inclusion materials is, respectively, and (corresponding to a phase contrast *γ*=100), both having Poisson’s ratio *ν*=0.3. The target two-dimensional elasticity tensor is chosen as
4.5
The starting topology of the RVE consists of a unit square homogeneous matrix with a circular inclusion of diameter 0.50 at the centre of the RVE. In this particular example, the widely used periodicity assumption (2.13) is adopted in the prediction of the macroscopic behaviour. The corresponding macroscopic elasticity is
4.6
To solve the set of linear variational problems required in the computation of and its topological derivative field, we use a standard finite element strategy. A fine structured finite element mesh consisting of 93 080 three-noded (linear) triangles is used to discretize the RVE. The mesh is kept constant and only Young’s modulus of the elements is changed (being assigned the values or ) throughout the iterations according to the criterion set above in the definition of the algorithm. The algorithmic parameter *n*_{e} has been chosen here as 1 per cent of the total number of elements in the mesh. Figure 2 illustrates the evolution of the topology during the automatic synthesis procedure. It shows the initial, an intermediate and the final topology. The dark areas correspond to inclusion material and the light-coloured areas to matrix material. Figure 3 shows the periodic assembly of the final microstructure topology and the corresponding evolution of the distance function *ϕ*. The final topology has macroscopic elasticity tensor
4.7
with corresponding normalized value of the distance function *ϕ*≈0.013. We remark that owing to the simplicity of the derived topological derivative formula, the calculations required by the algorithm are straightforward and of easy computational implementation.

## 5. Conclusions

A fundamental analytical formula for the sensitivity of the two-dimensional macroscopic elasticity tensor to the insertion of circular inclusions has been derived by applying the concept of topological derivative within a well-established multi-scale constitutive modelling framework for linear elasticity. The multi-scale constitutive framework is based on the assumption that the macroscopic strain and stress tensors are defined as volume averages of their microscopic counterparts over an RVE. It allows different predictions of macroscopic behaviour—including an upper and a lower bound for stiffness—to be obtained according to the chosen kinematical constraints imposed upon the RVE. The derived sensitivity—a symmetric fourth-order tensor field over the RVE domain—measures how the estimated macroscopic elasticity tensor changes when a small circular inclusion is introduced at the microscale. This result is fundamental and has potential application in different areas of interest. To illustrate its potential applicability, a very simple algorithm based on the proposed formula and relying on the finite element approximation of the RVE equilibrium problem has been devised to synthesize a microstructure to a pre-specified macroscopic response. The numerical results have shown the successful use of the proposed expression in this context. Further use of the present formula in the context of microstructural optimization is currently under investigation and will be the subject of a forthcoming publication.

## Acknowledgements

The collaboration between the Swansea University and LNCC researchers has been funded by the Royal Society (International Joint Project no. JP080066). S.M.G. has been partly supported by CNPq (Brazilian Council for Scientific and Technological Development—Grant no. 382485/2009-2). All support is gratefully acknowledged.

## Appendix A. Asymptotic analysis

This appendix presents the derivation of the asymptotic formula used in the topological sensitivity analysis developed in §3*b*. We start by considering the following expansion of the stress fluctuation field associated with the solution to problem (3.8) (Sokołowski & Żochowski 1999):
A1
where denotes the solution of the elasticity system (3.8) in the infinite domain , such that the stresses tend to a constant value when . Then, the exterior problem can be written as
A2
where **n** denotes the outward unit normal to the boundary , is the solution of the unperturbed problem (2.8) and is defined in equation (2.7).

In a polar coordinate system (*r*,*θ*) having its origin at the centre of the hole and with the angle *θ* measured with respect to one of the principal directions of , the components of the solution of the partial differential equation (A2) are given by

— Exterior solution (

*r*≥*ε*) A3— Interior solution (0<

*r*<*ε*) A4

where *φ* indicates the angle between principal stress directions associated with the stress fields and . In addition, we denote
A5and
A6
with and representing the principal stresses associated with the displacement fields and of the original (unperturbed) domain *Ω*_{μ}. The constants *α* and *β* used in equations (A3) and (A4) are given by
A7

## Footnotes

- Received September 24, 2009.
- Accepted November 27, 2009.

- © 2010 The Royal Society