Sensitivity of the Second Order Homogenized Elasticity Tensor to Topological Microstructural Changes

The multiscale elasticity model of solids with singular geometrical perturbations of microstructure is considered for the purposes, e.g., of optimum design. The homogenized linear elasticity tensors of first and second orders are considered in the framework of periodic Sobolev spaces. In particular, the sensitivity analysis of second order homogenized elasticity tensor to topological microstructural changes is performed. The derivation of the proposed sensitivities relies on the concept of topological derivative applied within a multiscale constitutive model. The microstructure is topologically perturbed by the nucleation of a small circular inclusion that allows for deriving the sensitivity in its closed form with the help of appropriate adjoint states. The resulting topological derivative is given by a sixth order tensor field over the microstructural domain, which measures how the second order homogenized elasticity tensor changes when a small circular inclusion is introduced at the microscopic level. As a result, the topological derivatives of functionals for multiscale models can be obtained and used in numerical methods of shape and topology optimization of microstructures, including synthesis and optimal design of metamaterials by taking into account the second order mechanical effects. The analysis is performed in two spatial dimensions however the results are valid in three spatial dimensions as well.


Introduction
The study of synthesis and design of materials involving multiscale effects gave rise to a wide interest in Engineering, Mechanics, and Mathematics during the two past decades, and it broadened the application scope, among others structural mechanics, biomechanics, aerospace engineering, wave propagation in solids, and acoustics. The research works on this subject have increased with the emergence of recent experimental and manufacturing techniques, computational methods and tools, and theoretical developments. The various length scales of this type of materials allow the elaboration of multiscale constitutive theories, the so-called theory of homogenization, in order to explain, more accurately than standard phenomenological approaches, their macroscopic response under loading for example. The first Extended author information available on the last page of the article developments has been made for periodic structures [18,30,34,36,45,46,54,61], and since then the framework has not stopped expanding to fit more general and complex models. An introduction to homogenization and related mathematical framework can be found in [23,62].
In this context, the design of the microstructure is a major issue, as well for a mixture of different materials, as for the arrangement of a single one. For example, in [6] and [42] microstructural topologies that produce negative macroscopic Poisson's ratio are obtained with a relaxation-based technique. To find out new microstructures producing exotic behaviours at the macroscopic scale, we can cite the use of classical shape optimization method (see, e.g., [35,60]), based on shape gradient of the desired criterion, with respect to a smooth variation of the boundary. More recently, the combination of shape gradient concept and level-set method (investigated in [53]), has produced several results in structural optimization, the reader can see for instance [2,5,20,52,55,63]. This approach depends deeply on the initial guess for the microstructure, because it does not allow for topology changes. Relaxed formulation based on homogenization theory has been developed in [1,3,4,17], and provide topology variations in certain cases.
Another strategy is based on the concept of topological derivative, which was rigorously introduced in [57]. The idea is to compute a topological asymptotic expansion of the investigated criterion with respect to an infinitesimal topological perturbation of the domain. The reader may find the use of this concept in topology optimization in [5,9,21]. In the framework of homogenized model of elastic materials, the topological derivative of the first order homogenized elasticity tensor has been calculated in [31,32] in the case of void and soft inclusion, respectively, and in [12] in the case of a soft material inclusion, completed with a numerical investigation. More recently the topological derivative of the second-order macroscopic model associated with scalar waves in periodic media has been evaluated in [19], making use of integral equations together with the periodic Green's function.
In the present paper, the elasticity system in periodic media is considered. The microstructure of the underlying material is topologically perturbed by the nucleation of a small circular inclusion endowed with different material properties from the background. In this case, the asymptotic behaviour of the perturbed displacement solutions is computed, allowing the evaluation of the sensitivity of the second-order homogenized elasticity tensor with respect to topological microstructural changes. The resulting topological derivative is given by a symmetric sixth-order tensor field over the RVE (Representative Volume Element) domain that measures how the second-order homogenized elasticity tensor changes when a small circular inclusion is introduced at the microscale level. This information is crucial for the synthesis and optimal design of microsctructures having a macroscopic behaviour depending on the second order derivative of the average displacement. We start by describing in Sect. 2 a homogenization scheme (see [56]), which leads to the formal definition of higher-order homogenized tensors. For this we need the solutions of auxiliary problems posed on the RVE, in the framework of periodic homogenization by asymptotic expansions. Then we undertake in Sect. 3 a perturbation of the RVE by including a new material characterised by a finite contrast, and contained in a small ball. After a substantial introduction to the topological derivative concept through Sect. 3.1, we recall in Sect. 3.2 the formula of the topological asymptotic of the classical first-order homogenized tensor derived in [12,31,32]. We also calculate the topological derivative associated with a simple higher-order homogenized tensor for introducing the method. In Sect. 3.3 we finally derive in details the topological derivative of the second-order homogenized tensor. The paper ends with some concluding remarks in Sect. 4. The proofs of certain lemmas are send back to Appendix B.

Homogenization
This section describes the multi-scale method to calculate the homogenized coefficients of an elasticity problem written for a periodic media, in order to calculate their topological sensitivities with respect to a configurational perturbation in the periodic cell. Let , D ⊂ R 2 be two domains. We denote by D τ the periodic medium consisting of the domain T D, for T > 0, paved with the periodic cell t , for t > 0, where we have set τ = t/T . We define y = x/t and Y = x/T for all x ∈ D τ (see Fig. 1). Let the vector field u τ (x) ∈ R 2 be the displacement, solution of the elasticity system in D τ . We assume that u can be expanded as where the functions u i (Y, ·) are -periodic for i ≥ 0. Using this expansion in the equilibrium and constitutive equations, we obtain a family of auxiliary problems, and we can write u τ thanks to a sequence of tensor fields (H ∇i (y)) i 0 called corrector fields, each one of order i + 2, and thanks to a sequence of macroscopic vector fields (U i (Y )) i 0 assumed to be constant on a cell. This gives where the (i + 2)-order tensor H ∇i operates for all j ≥ 0 on the (i + 1)-order tensor , the above expansion suggests to seek approximations of u τ in the form of truncations with respect to different order of τ of the following form Therefore, we seek the total field u τ as the sum of a macroscopic displacement field U and its i-th derivative weighted by τ i and a corrector field, for 1 ≤ i ≤ k. Calculating the homogenized elastic energy induced by this truncations, we can identify the homogenized tensors. We recall how to obtain formally the auxiliary equations and the corrector fields in the framework of the multi-scale method (see, e.g., [23,28]).

Auxiliary Equations
Let us write the auxiliary problems in their strong formulations. For a load f ∈ L 2 (D τ ), and a Dirichlet data u D ∈ H 1/2 (∂D τ ), the displacement vector field u τ from (2.1) is given by the solution of the following boundary value problem: (2.4) where the second order tensor field σ x (u τ ), called the total stress tensor, is specified throughout the constitutive law (2.5), and it depends linearly on the total strain tensor ε x (u τ ). The right lower index of a differential operator denotes the differentiation variable. We first determine the convention used for classical tensor calculus. Let u and v be two vectors, A and B be two second order tensors, and T be a fourth order tensor, we write is an orthonormal basis of R 2 , using the Einstein summation convention. Then we define where the elasticity tensor C = (C ij kl ) 1≤i,j,k,l≤2 is a fourth order tensor such that for all indices i, j, k, l = 1, 2 we have: C ij kl ∈ L ∞ ( ), is -periodic, C ij kl = C jikl = C klij , and there exists two real numbers 0 < a < b such that |CA| ≤ b|A| for any second order tensor A, and CA · A ≥ a|A| 2 when A is symmetric. Assuming that ∂ x = T −1 (∂ Y + τ −1 ∂ y ), and in view of ansatz (2.1), we can formally write ε x (u τ ) = τ −1 ∇ s y u 0 + ∞ i=0 τ i ε i where ε i := ∇ s Y u i + ∇ s y u i+1 . Let us define in the same way σ i := Cε i . Introducing expansion (2.1) of u τ in the equilibrium equation (2.4), we obtain a sequence of equations relating to the order of τ : Each of these equations is written on a unit cell , with the -periodicity of u i as boundary condition. The corresponding boundary value problems, also named auxiliary problems, can be solved by induction. The first equation Let us rewrite the second equation for all vectors a, b ∈ R 2 . By linearity we can write where the vector fieldũ ij is the solution of the -periodic boundary value problem on the unit cell for the first auxiliary equation: div y C∇ s yũ ij + C(e i ⊗ s e j ) = 0. (2.12) Before solving the third auxiliary problem, let us evaluate the average on the unit cell of equation (2.8c). We assume that f = T −1 F (Y ). This gives us that: where for all tensor fields A and where V = | | denotes the scaled area of the RVE, | | being the Lebesgue measure of . We can show that σ 0 = C h E 0 , where the homogenized tensor C h is constant. Now let us rewrite the equation (2.8c) setting E 1 := ∇ s Y U 1 and K 0 := ∇ Y E 0 , taking into account (2.13), and defining the vector field u ij (y) := (e i ⊗ s e j )y +ũ ij (y). We find div y C∇ s y u 2 + div y (CE 1 ) + [div y (C(ũ ij ⊗ s e k )) + (C∇ s y u ij − C h (e i ⊗ s e j ))e k ]K 0 ij k = 0.
(2.15) Once again by linearity we can write the solution u 2 in the following way where the vector fieldũ ij k is the solution of the -periodic boundary value problem on the unit cell for the second auxiliary equation: div y C∇ s yũ ij k + div y C(ũ ij ⊗ s e k ) + (C∇ s y u ij − C h (e i ⊗ s e j ))e k = 0. (2.17)

First-Order Truncation
We recall that Y = x/T , y = x/t and τ = t/T , with x ∈ R 2 . Motivated by expansion (2.3), we introduce the macroscopic displacement field U(Y ) ∈ R 2 , and the macroscopic deformation is defined as We also introduce the expansion The displacement fieldsũ ij are solutions of the following canonical set of variational problemsũ ij ∈ V : σ y (ũ ij ) · ε y (η) + C(e i ⊗ s e j ) · ε y (η) = 0, ∀η ∈ W, (2.19) where σ y (ũ ij ) = Cε y (ũ ij ) and the spaces W and V are defined as follows

20)
where H 1 per ( ; R 2 ) is the completion of the space of functions in C ∞ (R 2 , R 2 ) which are -periodic, in the norm of H 1 ( ; R 2 ). From these elements we have where u ij := (e i ⊗ s e j )y +ũ ij , (2.23) withũ ij solutions to the set of canonical variational problems (2.19). Then we calculate the average of the elastic energy 1 2 σ x (u τ ) · ε x (u τ ) on the RVE domain , denoted by W h , in order to identify the homogenized elasticity tensors. From now on we set T = 1 for convenience. We find We assume that terms of odd order τ in the expression of the homogenized elastic energy are equal to zero. This property may be proved thanks to the parity of the fields in the case of a centrosymmetric cell (see [56]). Thus, from now on we consider that the unit cell is centrosymmetric.
which defines the two following homogenized elasticity tensors: the fourth order tensor C h = (C h ij kl ) 1≤i,j,k,l≤2 , and the sixth-order tensor D = (D ij kpqr ) 1≤i,j,k,p,q,r≤2 given by and (2.27)

Second-Order Truncation
We introduce the expansion The displacement fieldsũ ij k are solutions of the following canonical set of variational problems Let us calculate the strain tensor of u τ : Same as before, we need to calculate 1 2 σ x (u τ ) · ε x (u τ ) in order to evaluate the average of the elastic energy on the cell and then identify the macroscopic energy law. Since the unit cell is assumed to be centrosymmetric, it turns out that the odd order terms in τ are null, and performing a formal macroscopic integration by parts in order to turn the coupled terms E ij ∂ y k K pqr into K ij k K pqr (see [56]), we calculate where D = (D ij kpqr ) 1≤i,j,k,p,q,r≤2 is the homogenized sixth-order tensor given by By setting η =ũ pqr as test function in (2.29), we obtain the following equality σ y (ũ ij k ) · ε y (ũ pqr ) + C(ũ ij ⊗ s e k ) · ε y (ũ pqr ) = (σ y (u ij ) − C h (e i ⊗ s e j ))e k ·ũ pqr , (2.33) which allows to write (2.32) as . So far, we have defined the homogenized tensors C h , D and D. In the following section, we slightly change the RVE, write the new tensors deriving from the perturbed RVE, and explore their behaviour regarding the size of such a perturbation.

Topological Sensitivity
The topological optimization framework is as follows. The original domain is composed of two phases of isotropic materials, the first one represented by the domain 1 , and the second represented by γ , such that = 1 ∪ γ ∪ γ , where γ = ∂ γ ∩ with ∂ and ∂ γ being Lipschitz continuous. These two phases result in a piecewise constant elasticity tensor denoted by C, which is defined as follows. Let C 0 = 2μI + λI ⊗ I, where the so-called Lamé coefficients μ, λ ∈ R are chosen such that C 0 satisfies the conditions from Sect. 2.1, and with I and I used to denote the second and the fourth order identity tensors, respectively.
We consider the nucleation of a finite number of small disjoint circular inclusions within the unit cell . The inclusions are spatially distributed far from each other and in such a way that the centrosymmetry of the RVE is preserved, which allows for focusing the attention over a single inclusion. From there, is subjected to a perturbation confined in a small circular open set B ρ (ŷ) of radius ρ and center at an arbitrary pointŷ of , such that B ρ (ŷ) ⊂ , and which does not touch the interface γ (see Fig. 2). Then, the region occupied by B ρ (ŷ) is filled by an inclusion with different material property from the background. The material properties of the perturbed domain are characterized by the piecewise constant function γ ρ of the form (3.1) Namely the elasticity tensor is given by γ ρ C in the perturbed domain. Henceforth we leave the lower indices of differential operators behind, because we only deal with y-variable depending fields. Therefore, the topologically perturbed counterparts of problems (2.19) and (2.29) are respectively given bỹ and Due to Korn's inequality, the existence and uniqueness of solutions of problems (2.19), (2.29), (3.2) and (3.3) are ensured on W endowed with the norm · W , which is defined as follows Actually, in view of the properties introduced in Sect. 2.1 and satisfied by the elasticity tensor C, the bilinear form of these problems is symmetric, continuous and coercive. We finally fix for each problem a solution belonging to H 1 per ( ; R 2 ) by choosing the representative which satisfies · = 0 for problems (2.19), (2.29) and (3.2), (3.3), so thatũ ij andũ ij k , as well as u ρ ij andũ ρ ij k belong to V.
Remark 1 According to Lax-Milgram theorem, the existence and uniqueness of solutions of variational problems on (W, · W ) with a symmetric, continuous and coercive bilinear form, are ensured as long as the continuous linear form belongs to the dual space of W, which can be identified with the subspace of the dual space (H 1 per ( ; R 2 )) whose the elements F ∈ (H 1 per ( ; R 2 )) are such that F (c) = 0 for all c ∈ R. All these conditions are satisfied for the problems (2.19), (2.29), (3.2), and (3.3).
As we did in Sect. 2, we can consequently define the topologically perturbed counterparts of the homogenized tensors, denoted as C ρ , D ρ and D ρ . By setting u ρ ij := (e i ⊗ s e j )y +ũ ρ ij , this gives

The Topological Derivative Method
We are interested in the behaviour of the homogenized tensors with respect to the size of the topological perturbation. For this purpose, we will use the concept of topological derivative. It has been rigorously introduced in [57] in the context of heat conduction and elasticity problems. Developments of the theory have been led the past two decades in among others [8,11,13,27,29,40,43,47,49,50,58,59]. Furthermore, the topological derivative was applied in many fields, such as topology optimization [9,10,51], inverse problems [22,33,38,39,44], and image processing [14,16,37]. Let ρ,ŷ be used to represent the topologically perturbed counterpart of . In the particular case of a perforation, for instance, ρ,ŷ = \ B ρ (ŷ). For a given shape functional ρ,ŷ → ψ( ρ,ŷ ) we are looking for the topological asymptotic expansion The sign of the topological derivative indicates whether it is interesting or not regarding the considered criterion ψ to add a small inclusion of material at the pointŷ. To this end, we will calculate the topological derivative of the homogenized tensors in the next section. Before proceeding, let us introduce a truncated domain of the form We fix R and consider small positive parameter ρ → 0 with R > ρ > 0, and such that B R (ŷ) is included in . Note that B R (ŷ) contains the inclusion B ρ (ŷ). The existence of the topological derivatives for the components of homogenized tensors is ensured by the following two lemmas. The proofs of these results are relegated to Appendix B.
Lemma 2 Letũ ij andũ ρ ij be the solutions of original (2.19) and perturbed (3.2) problems, respectively. Then, the following estimates hold true

First-Order Truncation
First, let us consider the expansion of the homogenized tensors C h and D . To this end, we exclusively need the estimates from Lemma 2. The calculations of the topological derivative of C h is well-known, and from [12,31] we have the following result.

Theorem 4
The topological asymptotic expansion of the homogenized elasticity tensor C h is given by which, setting f (ρ) = πρ 2 /V , allows for identifying the topological derivative of any component of C h , namely

17)
where u ij is given by (2.23) and the polarization tensor is defined as

18)
with the parameters α and β given by We are more interested in the topological derivative of D, but the computations performed for D are helpful to understand the method applied afterwards. Because the original and perturbed fieldsũ ij andũ ρ ij are living in the same functional space V, we can perform a direct calculation. Thus the topological asymptotic expansion of the tensor D given by (2.27) is obtained from its definition as follows where the remainder E 1 (ρ) is given by and can be bounded as follows where we have used Lemma 2. Still evoking Lemma 2, we notice that the last integral on the right-hand side of expression (3.20) gives rise to a ρ 2 order term in the asymptotic expansion of D ρ . But we can not use the same arguments to analyse the first integral (second line) on the right-hand side of (3.20). To overcome this difficult, we make use of the classical adjoint method by introducing suitable adjoint sates v kr ij ∈ V for i, j, k, r ∈ {1, 2}, solution of the following set of variational problems: From these elements, we can state the main result of this section, where the details of the calculations can be found in the Appendix A.1.

Theorem 5 The topological asymptotic expansion of tensor D is given by
where P is the polarization tensor defined in (3.18). By setting f (ρ) = πρ 2 /V , the topological derivative of any component of tensor D can be identified, namely

25)
Finally, u ij is given by (2.23),ũ ij are solutions to the set of canonical variational problems (2.19) and the adjoint states v kr ij are solutions to (3.23).

Second-Order Truncation
Now we want to investigate the topological sensitivity of the tensor D involved in the macroscopic elastic energy calculated for the second order truncation (2.31). As we did for D , we perform a direct calculation described in Appendix A.2. By taking into account that ũ pqr = 0 and ũρ pqr = 0, we have where the notation δ(·) ρ = (·) ρ − (·) has been introduced. The remainder E 1 (ρ) is defined as since u ρ ij := (e i ⊗ s e j )y +ũ ρ ij , so that δu ρ ij = δũ ρ ij . From Lemmas 2 and 3, the following estimate for the remainder E 1 (ρ) given by (3.27) holds true (3.28) By taking into account once again Lemmas 2 and 3, we note that the last two integrals in the expansion (3.26) are of order ρ 2 . However, in order to analyse the first fourth integrals in (3.26), the introduction of convenient adjoint states p r ij k ∈ V for i, j, k, r ∈ {1, 2} are required, which are solutions of the following set of variational problems: (3.29) Finally, from these elements we can state the main result of the article, where the details of the calculations can be found in the Appendix A.2.

30)
where P is the polarization tensor defined in (3.18). By setting f (ρ) = πρ 2 /V , we can identify the topological derivative of any component of tensor D, namely

Remark 7
In the particular case of the weak phase material which simulates a void contained in a ball B ρ of radius ρ centered atŷ, with the infinitely small contrast 0 < γ 0 1, the definitions of the macroscopic and cell boundary value problems should be changed. In such a case the loads are applied on the stiff phase only. For this purpose the characteristic function χ of the solid domain is introduced. The topologically perturbed counterpart of this characteristic function is written χ ρ = χ − χ Bρ , χ Bρ being the characteristic function of the ball B ρ . Let ρ 0 > 0, we define for 0 ≤ ρ ≤ ρ 0 ϕ ρ (y) := V χ ρ χ ρ (y) and η ϕ ρ := For the macroscopic problem, we assume that f = T −1 ϕ ρ 0 (y)F (Y ). We select the periodic solutions of cell problems with null averages on the solid phase. Therefore, the first auxiliary problem remains unchanged, and the second auxiliary problem becomes We can calculate the topological derivative of D in the same way as before, using Lemmas 2 and 3 which are still valid. The normalized characteristic function is introduced into the adjoint boundary value problem. The expression of topological derivative of D contains some new terms, (D T D) ij kpqr = −Pσ (u ij ) · (ε(p k pqr ) + (ũ pqr ⊗ s e k )) − Pσ (u pq ) · (ε(p r ij k ) + (ũ ij k ⊗ s e r )) (3.35) Note that the non-uniform contribution of C h in equation (3.33) results in corrections to the topological derivative.

Remark 8
We cover as well the case of three spatial dimensions which is important for applications. In particular, the method of asymptotic analysis performed in R 2 can be extended to R 3 . In three spatial dimensions, the topological expansion of homogenized tensors is obtained by setting f (ρ) = (4/3)πρ 3 /V and replacing the polarization tensor by [7] P = −3βI − (α − β)I ⊗ I, (3.36) with the coefficients α and β redefined as follows , (3.37) where E is the Young modulus and ν the Poisson ratio.
Finally, it is important to note that formula (3.31) can be used to evaluate the topological derivative of any differentiable function of D through the direct application of the conventional calculus rules for composed functions. That is, any such function D → (D) admits the topological derivative of the form (3.38) with the brackets ·, · denoting the appropriate product between the derivative of with respect to D and the topological derivative D T D of D. In order to fix these ideas, let us consider a pair 1 , 2 ∈ R 2 × R 2 × R 2 of third order tensors. Then we obtain the following results, which can be used in numerical methods of synthesis and/or topology design of microstructures analogously to [12]: We consider a function (D) of the form Therefore, according to (3.38), its topological derivative is given by If we set 1 = e i ⊗ e j ⊗ e k and 2 = e l ⊗ e m ⊗ e n , for instance, we get (D) = (D) ij klmn and its topological derivative is given by D T (D) = (D T D) ij klmn . It means that D T (D) actually represents the topological derivative of the components (D) ij klmn of the tensor D.

Remark 10
In [24] a pantographic material is investigated, leading to the first order homogenized tensor C h which is not invertible. In such a case, the addition of the strain-gradient terms K = ∇E in the macroscopic model is needed to predict more accurately the behaviour of the material. It turns out that an useful information arises from the projection of D onto the kernel of C h . This fact points strongly to the suitability of the use of (3.31) in a topology design algorithm for the synthesis and optimization of elastic microstructures based on minimization/maximization of cost functions defined in terms of homogenized properties.

Conclusion
In this paper, the topological derivative of the second order homogenized elasticity tensor with respect to the nucleation of circular inclusions at the microscopic level has been presented in the framework of periodic Sobolev spaces. The sensitivity has been derived in its closed form with the help of appropriate adjoint states. The limiting case associated with the nucleation of a very weak inclusion has also been considered. As expected, the associated topological derivative leads to a sixth order tensor field over the microstructural domain, measuring the sensitivity of the second order homogenized elasticity tensor to topological microstructural changes. Therefore, this information can be used in the context of synthesis and optimal design of metamaterials, for instance, accounting for second order mechanical effects.

A.1 Proof of Theorem 5
First we rewrite below the topological asymptotic expansion (3.20) In order to simplify the second term on the right-hand side of (A.1), we evoke the adjoint method. We start by subtracting (2.19) from (3.2), to obtain where u ρ pq := (e p ⊗ s e q )y +ũ ρ pq . By setting η =ũ ρ pq −ũ pq as test function in the adjoint problem (3.23) for v kr ij , and noting that ũ ρ pq −ũ pq = 0, we obtain After taking η = v kr ij as test function in (A.2), we have From the symmetry of the bilinear forms we conclude that Similarly we have, after replacing the indexes pq by ij in (A.2) and (3.23), that These results lead to Since we assume that the inclusion B ρ is located neither on the interface nor on the boundary, the solutions of elliptic boundary value problems are smooth in B ρ by the elliptic regularity. Finally, from Lemma 2, and by taking into account the Lebesgue differentiation theorem combined with the Eshelby theorem [25,26], we can write the topological derivative with the use of the polarization tensor P (see, e.g., [12,31]), and we deduce Theorem 5.

A.2 Proof of Theorem 6
The topological asymptotic expansion of the tensor D given by (2.32) or alternatively by (2.34) is obtained as follows. We start by rewriting expansion (3.26), namely where the notation δ(·) ρ = (·) ρ − (·) has been introduced. After subtracting (2.29) from (3.3), we obtain ∀η ∈ W σ (δũ ρ ij k ) · ε(η) = (σ (δũ ρ ij ) − (C h ρ − C h )(e i ⊗ s e j ))e k · η − C(δũ ρ ij ⊗ s e k ) · ε(η) (A.9) We set η =ũ pqr in (A.9) and η = δũ ρ ij k in (2.29), the equation forũ pqr , leading to simplification of the terms depending on δũ ρ ij k in (A.8). Reordering the members of the equation in order to gather the terms of the type δũ ρ ij and δũ ρ pq , we find (A.10) Let us set η = δũ ρ ij and η = δũ ρ pq in the adjoint equation (3.29) for p k pqr and p r ij k , respectively. Now, we set in (A.2) η = p r ij k , and η = p k pqr after replacing the indexes ij by pq. After comparing the obtained results, we can conclude that We recall that the inclusion B ρ is located neither on the interface nor on the boundary, so that the solutions of elliptic boundary value problems are smooth in B ρ by the elliptic regularity. Finally, from Lemmas 2 and 3, and by taking into account the Lebesgue differentiation theorem combined with the Eshelby theorem [25,26], we can write the topological derivative with the use of the polarization tensor P (see, e.g., [12,31]), and we deduce Theorem 6.

Appendix B: Proofs of Lemmas 2 and 3
For the convenience of the reader we provide the proofs of auxiliary Lemmas which are used to evaluate of the topological derivatives of homogenized tensors.

B.3 Preliminary Lemmas
We write in this subsection two useful statements for the proof of Lemmas 2 and 3.
We consider an open subset of R 2 , andŷ a fixed arbitrary point of . We denote by B ρ be the small disk of radius ρ, centered atŷ, and take ρ 0 > 0, such that B ρ ⊂ for all 0 < ρ ρ 0 . We have the following results.

Lemma 12
Let η ∈ H 1 ( ; R 2 ). Then we have for all δ > 0 a constant c (δ) > 0 depending on δ and such that for all 0 < ρ ≤ ρ 0 Proof For simplicity we set ρ 0 = 1. Let 0 < ρ ≤ 1 and η ∈ H 1 ( ; R 2 ). We introduce φ ρ : The restriction φ ρ : ∂B 1 → ∂B ρ is also a diffeomorphism, allowing us the following change of variable where c (1) > 0 is the fixed constant given by the Trace theorem applied on H 1 (B 1 ; R 2 ). Once again, a change of variable yields and Then we have The Sobolev Embedding Theorem gives H 1 ( ; R 2 ) → L p ( ; R 2 ) continuously, for all 2 ≤ p < +∞. Then from Hölder inequality with q −1 = 1 − (p/2) −1 , q ∈ (1, +∞] we have where c (p) is the constant depending on p and given by the Sobolev Embedding Theorem. Thus we can carry on the derivation of inequality (B.5) for δ > 0 as small as we want when q goes to 1. We finally conclude taking the square root of the last inequality.
We end this preliminary section presenting convenient notation used to name the norms and seminorms of the needed functional spaces. For an open set ⊂ R 2 , bounded, connected and with Lipschitz continuous boundary, H 1 ( , R 2 ) is endowed with a norm equivalent to the usual one, due to Korn's inequality [48], defined by and we also define in H 1 ( , R 2 ) the seminorm We finally denote the usual norms of L 2 ( , R 2 ), H 1/2 (∂ , R 2 ) and H −1/2 (∂ , R 2 ), respectively by η 0, , η 1/2,∂ , and η −1/2,∂ .

B.4 Lemma 2
Proof of Lemma 2 Let us introduce an ansatz of the form [41] u ρ ij (y) =ũ ij (y) + w ρ ij (y) + z ρ ij (y), (B.10) where w ρ ij is solution to the following exterior problem: where h = −(1 − γ )σ (u ij )(ŷ)n, n being the inward normal vector on ∂B ρ , and · denotes the jump across the interface of the inclusion: The solution w ρ ij of this classical problem is explicitly known [15], and gives rise to these estimates with R defined by (3.9). We want to estimate z ρ ij which compensates the discrepancies introduced by w ρ ij . But defined as in (B.10), we don't have z ρ ij ∈ H 1 per ( ; R 2 ). For having that, let us modify slightly w ρ ij near the boundary ∂ . We define the ring C := {y ∈ | dist(y, ∂ ) < } , (B.14) with ∂C = ∂ ∪ ∂C int , ∂C int = {y ∈ | dist(y, ∂ ) = }, for a fixed > 0 small enough to have B R (ŷ) ⊂ \ C. Then we set where w ρ C,ij is solution to the following boundary value problem Let us introduce the new ansatz wherew ρ ij := w ρ, ij − w ρ, ij , with · defined as in (2.14), so thatw ρ ij ∈ V, and thereforẽ z ρ ij ∈ H 1 per ( ; R 2 ) with z ρ ij = 0. Thus we can calculate the variational problem satisfied bỹ z ρ ij . Let us apply the operator γ ρ σ to (B.17), multiply the obtained expression by ε(η) where η is a test function in H 1 per ( ; R 2 ), and integrate over . Applying the Green formula over C ∩ , ( \ C ∪ B ρ ), and B ρ to the integral term depending onw ρ ij , and in view of (2.19), (3.2), (B.11) and (B.16), we find Because of the regularity of u ij , we have the following estimate In view of the elliptic regularity of the problem (B.16), the continuity of the trace operator, the estimate (B.13), and because we can show that there exists a constant c > 0 such that for all η ∈ H 1 (C, R 2 ) satisfying div σ (η) = 0, we have σ (η)n −1/2,∂ ≤ c |η| 1,C , (B.20) then we have the following estimate ∂ σ (w ρ C,ij )n · η ≤ σ (w ρ C,ij )n −1/2,∂ η 1/2,∂ ≤ c|w ρ C,ij | 1,C η 1, ≤ c w ρ ij 1/2,∂C int η 1, ≤ c w ρ ij 1, R \C η 1, ≤ cρ 2 η 1, . (B.21) In the same way we obtain the estimate O(ρ 2 ) η 1, for the integrals over ∂C int . Noting that the bilinear form of the problem (B.18) is uniformly coercive on W, and becausez ρ ij ∈ V, we can conclude, making use of Poincaré-Wirtinger inequality, that We can finally write the following expansioñ and noting after few calculations that we finally find the estimates (3.10), (3.11) and (3.12).

B.5 Lemma 3
Before proving Lemma 3, let us make some preliminary calculations, for which the results directly ensue from Lemma 2. We denote by a ρ the bilinear form on W of problem (3.3), that is to say We want to estimate, for all η ∈ W, the expression a ρ (ũ ρ ij k , η) − a ρ (ũ ij k , η). According to expressions (3.3) and (2.29), we find Recalling the notation δ(·) ρ = (·) ρ − (·), we start developing the first two terms of the right hand side of this expression.

(B.40)
Now we can introduce the new ansatz: u ρ ij k =ũ ij k + w ρ ij k + z ρ ij k +w ρ ij k − w ρ ij k , (B.41) wherew ρ ij k = w ρ, ij k − w ρ, ij k , so thatw ρ ij k ∈ V. In this way, we effectively havez ρ ij k ∈ H 1 per ( ; R 2 ), with z ρ ij k = 0. As in the proof of Lemma 2, we write the variational problem satisfied byz ρ ij k , applying the bilinear form a ρ toz ρ ij k and a test function η. In view of (B.41) and the notations introduced before, the problem is expressed bỹ z ρ ij k ∈ W : a ρ (z ρ ij k , η) = a ρ (ũ ρ ij k , η) − a ρ (ũ ij k , η) − a ρ (w ρ ij k , η), ∀η ∈ W. Let us apply the Green formula to − γ ρ σ (w ρ ij k ) · ε(η) on domains C ∩ , ( \ C ∪ B ρ ), and B ρ . In view of the definition ofw ρ ij k (B.38), and denoting by n the inward normal to the boundary ∂B ρ , we finally obtain thatz ρ ij k follows the variational problem: In the same way as in Lemma 2, we find that the two first terms on the right-hand side of equation ( Furthermore the solution w ρ ij k of the classical problem (B.38) is explicitly known (see, e.g., [15]), gives rise to the same estimates as those written equations (B.13), and few developments lead as well to