Optimal Heat Distribution Using Asymptotic Analysis Techniques

Abstract

In this chapter, we consider the optimization problem of a heat distribution on a bounded domain Ω containing a heat source at an unknown location ω⊂Ω. More precisely, we are interested in the best location of ω allowing a suitable thermal environment. For this propose, we consider the minimization of the maximum temperature and its L2 mean oscillations. We extend the notion of topological derivative to the case of local coated perturbation and we perform the asymptotic expansion of the considered shape functionals. In order to reconstruct the location of ω, we propose a one-shot algorithm based on the topological derivative. Finally, we present some numerical experiments in two dimensional case, showing the efficiency of the proposed method.

1. Introduction

The concept of topological derivative is a powerful tool for solving shape optimization problems constrained by partial differential equations. The method has a great potential of applications in the field of non-destructive control. In this chapter, the topological derivative is applied in the context of optimization of a heat distribution. More precisely, we consider the problem of locating circular coated inclusions in order to get an appropriate layout with minimized maximum temperature distribution. This problem can be encountered in the design of current carrying multicables and in some devices in hybrid and electric cars. The mathematical problem is similar to the mixture of materials with different conduction properties extensively studied in the case of two materials; see for instance [1, 2, 3, 4, 5, 6].

The topological derivative measure the sensitivity of a given shape functional with respect to the insertion of a small hole inside the domain. More precisely, we consider a domain $\mathrm{\Omega }\subset {\mathbb{R}}^2$ and a cost functional $\mathcal{J}\left(\mathrm{\Omega }\right)=j\left(\mathrm{\Omega },u_{\mathrm{\Omega }}\right)$, were $u_{\mathrm{\Omega }}$ is the state variable, i.e. a solution of a given partial differential equation in $\mathrm{\Omega }$. For $\epsilon >0$, let ${\mathrm{\Omega }}_{\epsilon }=\mathrm{\Omega }\backslash \overline{\left(x_0+\mathit{\epsilon D}\right)}$ be the domain obtained by removing a small part $\overline{\left(x_0+\mathit{\epsilon D}\right)}$ from $\mathrm{\Omega }$, at a location $x_0\in \mathrm{\Omega }$, and $D$ is a fixed bounded subset of ${\mathbb{R}}^2$ with $\left(0,0\right)\in D$.

Then, the shape functional $\mathcal{J}\left({\mathrm{\Omega }}_{\epsilon }\right)$ associated with the topologically perturbed domain, admits the following topological asymptotic expansion

$$\mathcal{J}\left({\mathrm{\Omega }}_{\epsilon }\right)=\mathcal{J}\left(\mathrm{\Omega }\right)+f\left(\epsilon \right)g\left(x_0\right)+o\left(f\left(\epsilon \right)\right),$$E1

where $\mathcal{J}(\mathrm{\Omega }$ is the shape functional associated to the unperturbed domain $\mathrm{\Omega }$ and $f\left(\epsilon \right)$ is a positive function such that $f\left(\epsilon \right)\to 0$ as $\epsilon \to 0$. The function $x_0\to g\left(x_0\right)$ is called topological gradient of $\mathcal{J}$ at $x_0$. Note that the topological derivative is defined through the limit passage $\epsilon \to 0$. However, according to (1), it can be used as a descent direction in an optimization process similar to any gradient-based method. This concept has been applied for geometrical inverse problems [7], in linear isotropic elasticity [8], and in the context of solving some design problems in steady-state heat conduction [9]. Another situation addressed in [10, 11, 12], consists in studying the influence of the insertion of a small inhomogeneity which is nonempty, but whose constitutive parameters are different from those of the background medium. This approaches was successfully investigate in the case of inhomogeneities conductor materials. For more details, we refer the reader to the recent works on shape reconstruction and stability analysis for some imaging problems with the topological derivative [13, 14, 15].

In this chapter, we extend the asymptotic analysis with respect to the insertion of local small inhomogeneity to the case of the insertion of local small coated inclusion. An asymptotic expansion of a given shape functional is derived with the help of a relevant adjoint method. The computed topological derivative allows us to solve numerically the minimization problem.

The chapter is organized as follows. In Section 2 we present the model problem and the shape optimization formulation. In Section 3 and 4 we perform the asymptotic expansions of the shape functional. The numerical results are presented in Section 5. The paper ends with some concluding remarks in Section 6.

2. The model problem

Let $\mathrm{\Omega }$ be a bounded domain in ${\mathbb{R}}^2$ with Lipschitz boundary $\mathrm{\partial \Omega }$ and let $\omega$ be an open subset of $\mathrm{\Omega }$ composed of two different subdomains ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ where the subset ${\mathrm{\Omega }}_2$ is surrounded by the subset ${\mathrm{\Omega }}_1$. We denote by ${\mathrm{\Gamma }}_2≔\partial {\mathrm{\Omega }}_2$ and ${\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2≔\partial {\mathrm{\Omega }}_1$ as depicted in Figure 1 and we set ${\mathrm{\Omega }}_0≔\mathrm{\Omega }\backslash \overline{{\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2}$.

Figure 1.

The domain $\mathrm{\Omega }={\mathrm{\Omega }}_0\cup {\mathrm{\Omega }}_1\cup {\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2\cup {\mathrm{\Omega }}_2$.

Throughout the chapter, we consider piecewise constant thermal conductivity

$$\sigma ={\sigma }_0{\chi }_{{\mathrm{\Omega }}_0}+{\sigma }_1{\chi }_{{\mathrm{\Omega }}_1}+{\sigma }_2{\chi }_{{\mathrm{\Omega }}_2},$$E2

where ${\sigma }_0,{\sigma }_1,{\sigma }_2\in {\mathbb{R}}_{+}^{\ast }$ and ${\chi }_E$ denotes the indicator function of the set $E$. We assume further that there exists two constants $c_0$, $c_1$ such that

$$0<c_0\le {\sigma }_0,{\sigma }_1,{\sigma }_2\le c_1.$$

For a given source term $f\in L^2\left(\mathrm{\Omega }\right)$ and the Dirichlet data $g\in H^{1/2}\left(\mathrm{\partial \Omega }\right)$, the temperature $u_{\omega }$ satisfies the following problem

$$\left\{\begin{array}{ll}-div\left(\sigma \nabla u_{\omega }\right)=f& \text{in}\mathrm{\Omega }\\ 〚u_{\omega }〛=0& \mathrm{on}{\mathrm{\Gamma }}_i,i=1,2,\\ 〚\sigma {\partial }_n{u}_{\omega }〛=0& \mathrm{on}{\mathrm{\Gamma }}_i,i=1,2,\\ u_{\omega }=g& \mathrm{on}\mathrm{\partial \Omega }.\end{array}\right.$$E3

For simplicity we take $g=0$ by choosing a lifting function $G\in H^2\left(\mathrm{\Omega }\right)$, $G=g$ in $\mathrm{\partial \Omega }$ and modifying the left hand side that we still denote $f$. Then, the weak solution to problem (3) is defined by:

$$\text{Find}{u}_{\omega }\in H_0^1\left(\mathrm{\Omega }\right)\text{such\hspace{0.17em}that}a\left(u_{\omega },v\right)=l\left(v\right),\forall v\in H_0^1\left(\mathrm{\Omega }\right),$$E4

where

$$a\left(u_{\omega },v\right)={\int }_{\mathrm{\Omega }}\sigma \nabla u_{\omega }\cdot \nabla v\mathit{dx},\text{and}l\left(v\right)={\int }_{\mathrm{\Omega }}\mathit{fv}\mathit{dx}.$$

The existence and uniqueness of the weak solution $u_{\omega }$ follows from the Lax-Milgram Lemma.

We consider the following shape minimization problem

$$\mathit{\text{Determine\hspace{0.17em}the\hspace{0.17em}position\hspace{0.17em}of}}\omega \subset \mathrm{\Omega }\mathit{\text{to\hspace{0.17em}obtain\hspace{0.17em}an\hspace{0.17em}appropriate\hspace{0.17em}layout\hspace{0.17em}temperature}}.$$E5

To deal with numerical computation of problem (5), we consider two shape functionals. The first shape functional corresponds to the maximum temperature

$$M\left(\omega ,u_{\omega }\right)=\parallel u_{\omega }{\parallel }_{L^{\infty }\left(\mathrm{\Omega }\right)}.$$

Since the functional $M$ is not differentiable, we can not use the topological derivative framework to perform the sensitivity analysis. Thus we use the shape functional $J_p$, for large $p\ge 2$ instead of the functional $M$:

$$J_p\left(\omega ,u_{\omega }\right)=\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|u_{\omega }\right|}^p\mathit{dx}.$$

The second shape functional may appear as a particular case of $J_p$ but has its own physical and mathematical interests, corresponds to the minimization of the $L^2$ mean oscillations of the temperature:

$$K\left(\omega ,u_{\omega }\right)≔{\int }_{\mathrm{\Omega }}{\left(u_{\omega }-\frac{1}{\mid \mathrm{\Omega }\mid }{\int }_{\mathrm{\Omega }}{u}_{\omega }\mathit{dx}\right)}^2\mathit{dx}.$$

Then, the optimization problems read:

$$\left\{\begin{array}{l}minimize{J}_p\left(\omega ,u_{\omega }\right)≔\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|u_{\omega }\right|}^p\mathit{dx}\\ \text{subject\hspace{0.17em}to}\omega \in \mathcal{O}\text{and}{u}_{\omega }\text{solves\hspace{0.17em}problem}\left(3\right),\end{array}\right.$$E6

and

$$\left\{\begin{array}{l}minimizeK\left(\omega ,u_{\omega }\right)≔{\int }_{\mathrm{\Omega }}{\left(u_{\omega }-\frac{1}{\mid \mathrm{\Omega }\mid }{\int }_{\mathrm{\Omega }}{u}_{\omega }\mathit{dx}\right)}^2\mathit{dx}\\ \text{subject\hspace{0.17em}to}\omega \in \mathcal{O}\text{and}{u}_{\omega }\text{solves\hspace{0.17em}problem}\left(3\right).\end{array}\right.$$E7

Where $\mathcal{O}$ is the admissible set:

$$\mathcal{O}=\left\{\omega \subset \mathrm{\Omega }:P\left(\omega ,\mathrm{\Omega }\right)<\infty \right\},$$

and $P\left(\omega ,\mathrm{\Omega }\right)$ is the relative perimeter:

$$P\left(\omega ,\mathrm{\Omega }\right)≔sup\left\{{\int }_{\omega }\nabla \cdot \phi \mathit{dx}:\varphi \in C_c^1\left(\mathrm{\Omega },{\mathbb{R}}^2\right),\parallel \phi {\parallel }_{\infty }\le 1\right\}.$$

3. Topological derivatives

Now, we assume that the domain ${\mathrm{\Omega }}_2≔x_0+\mathit{\alpha \epsilon B}$ and ${\mathrm{\Omega }}_1$ is such that ${\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2=x_0+\mathit{\epsilon B}$, where $B$ is the unit ball in ${\mathbb{R}}^2$, $\epsilon >0$ and $0<\alpha <1$. We rewrite ${\mathrm{\Omega }}_1^{\epsilon }$ and ${\mathrm{\Omega }}_2^{\epsilon }$ instead of ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$. This allows to perform an asymptotic expansion of the shape functional $J_p\left({\omega }_{\epsilon }\right)$ where ${\omega }_{\epsilon }≔{\mathrm{\Omega }}_2^{\epsilon }\cup {\mathrm{\Omega }}_1^{\epsilon }$. We also introduce ${\mathrm{\Gamma }}_2^{\epsilon }≔\partial {\mathrm{\Omega }}_2^{\epsilon }$ and ${\mathrm{\Gamma }}_1^{\epsilon }$ the outer boundary of ${\mathrm{\Omega }}_1^{\epsilon }$.

In the perturbed domain, the state $u_{\epsilon }$ is solution to the following problem:

$$\left\{\begin{array}{lll}-div\left({\sigma }_{\epsilon }\nabla u_{\epsilon }\right)& =f_{\epsilon }\text{in}\mathrm{\Omega }& \\ u_{\epsilon }& =0\mathrm{on}\mathrm{\partial \Omega }\end{array}\right.$$E8

where

$${\sigma }_{\epsilon }={\sigma }_2{\chi }_{{\mathrm{\Omega }}_2^{\epsilon }}+{\sigma }_1{\chi }_{{\mathrm{\Omega }}_1^{\epsilon }}+{\sigma }_0{\chi }_{\mathrm{\Omega }\backslash \overline{{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }}},\text{and}{f}_{\epsilon }=f_2{\chi }_{{\mathrm{\Omega }}_2^{\epsilon }}+f_1{\chi }_{{\mathrm{\Omega }}_1^{\epsilon }}+f_0{\chi }_{\mathrm{\Omega }\backslash \overline{{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }}}.$$

The functions $f_i$ are in $L^2$. The variational formulation associated with the problem (8) is defined by:

$$\left\{\begin{array}{l}\text{find}{u}_{\epsilon }\in H_0^1\left(\mathrm{\Omega }\right),\text{such\hspace{0.17em}that}\\ a_{\epsilon }\left(u_{\epsilon },v\right)=l_{\epsilon }\left(v\right),\forall v\in H_0^1\left(\mathrm{\Omega }\right),\end{array}\right.$$E9

where

$$a_{\epsilon }\left(u_{\epsilon },v\right)={\int }_{\mathrm{\Omega }}{\sigma }_{\epsilon }\nabla u_{\epsilon }\cdot \nabla v\mathit{dx},\text{and}{l}_{\epsilon }\left(v\right)={\int }_{\mathrm{\Omega }}{f}_{\epsilon }v\mathit{dx}.$$

We denote by $u_0$ the background solution of the following problem:

$$-div\left({\sigma }_0\nabla u_0\right)=f_0,\text{in}\mathrm{\Omega }{u}_0=0\mathrm{on}\mathrm{\partial \Omega }.$$E10

The following proposition describes in an abstract framework the adjoint method we will use to derive the asymptotic expansion of a given shape functional. For more details the reader is referred to [16] and the references therein.

Proposition 1 Let $\mathcal{H}$ be a Hilbert space. For all parameter $\epsilon \in [0,{\epsilon }_0[,{\epsilon }_0>0$, consider a function $u_{\epsilon }\in \mathcal{H}$ solving a variational problem of the form

$$a_{\epsilon }\left(u_{\epsilon },v\right)=l_{\epsilon }\left(v\right)\forall v\in \mathcal{H},$$E11

where $a_{\epsilon }$ is a bilinear form and $l_{\epsilon }$ is a linear form on $\mathcal{H}$. For all $\epsilon \in [0,{\epsilon }_0[$, consider a functional $J_{\epsilon }:\mathcal{H}\to \mathbb{R}$ that is Fréchet differentiable at $u_0$. Assume that the following hypotheses are satisfied.

${\mathbf{H}}_{\mathbf{1}}$ There exist a scalar function $f\left(\epsilon \right)\ge 0$ and two numbers $\mathit{\delta a}$, $\mathit{\delta l}$ such that

$$\left(a_{\epsilon }-a_0\right)\left(u_0,v_{\epsilon }\right)=f\left(\epsilon \right)\mathit{\delta a}+o\left(f\left(\epsilon \right)\right),$$E12

$$\left(l_{\epsilon }-l_0\right)\left(v_{\epsilon }\right)=f\left(\epsilon \right)\mathit{\delta l}+o\left(f\left(\epsilon \right)\right),$$E13

$$\underset{\epsilon \to 0}{lim}f\left(\epsilon \right)=0,$$E14

where $v_{\epsilon }\in \mathcal{H}$ is an adjoint state satisfying

$$a_{\epsilon }\left(\phi ,v_{\epsilon }\right)=-{\mathit{DJ}}_{\epsilon }\left(u_0\right)\phi ,\forall \phi \in \mathcal{H}.$$E15

${\mathbf{H}}_{\mathbf{2}}$ There exist two numbers $\delta J_1$ and $\delta J_2$ such that

$$J_{\epsilon }\left(u_{\epsilon }\right)=J_{\epsilon }\left(u_0\right)+{\mathit{DJ}}_{\epsilon }\left(u_0\right)\left(u_{\epsilon }-u_0\right)+f\left(\epsilon \right)\delta J_1+o\left(f\left(\epsilon \right)\right),$$E16

$$J_{\epsilon }\left(u_0\right)=J_0\left(u_0\right)+f\left(\epsilon \right)\delta J_2+o\left(f\left(\epsilon \right)\right).$$E17

Then the first variation of the cost function with respect to $\epsilon$ is given by

$$J_{\epsilon }\left(u_{\epsilon }\right)-J_0\left(u_0\right)=f\left(\epsilon \right)\left(\mathit{\delta a}-\mathit{\delta l}+\delta J_1+\delta J_2\right)+o\left(f\left(\epsilon \right)\right).$$E18

3.1 Application to the model problem

In this subsection, we will give explicitly the variations $\mathit{\delta a},\mathit{\delta l},\delta J_1,\delta J_2$ and we perform the asymptotic expansion of the shape functional $J_p$. Analogously, we can derive in the same manner the asymptotic expansion of the shape functional $K$.

3.1.1 Variation of the bilinear form

In this subsection, we look on the asymptotic analysis of the variation

$$\left(a_{\epsilon }-a_0\right)\left(u_0,v_{\epsilon }\right)={\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_1-{\sigma }_0\right)\nabla u_0\cdot \nabla v_{\epsilon }\mathit{dx}+{\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left({\sigma }_2-{\sigma }_0\right)\nabla u_0\cdot \nabla v_{\epsilon }\mathit{dx}.$$E19

Let us first look at the behavior of the adjoint state $v_{\epsilon }$ solution of the following boundary value problem

$$\left\{\begin{array}{ll}div\left({\sigma }_{\epsilon }\nabla v_{\epsilon }\right)=u_{\omega }{\left|u_{\omega }\right|}^{p-2}& \text{in}\mathrm{\Omega },\\ v_{\epsilon }=0& \mathrm{on}\mathrm{\partial \Omega }.\end{array}\right.$$E20

Since $u_{\omega }$ is Hölder continuous, $u_{\omega }{\left|u_{\omega }\right|}^{p-2}$ is at least in $L^2\left(\mathrm{\Omega }\right)$. Therefore problem (20) has a unique solution $v_{\epsilon }\in H_0^1\left(\mathrm{\Omega }\right)$.

We split in (19) $v_{\epsilon }$ into $v_0+\left(v_{\epsilon }-v_0\right)$ and introducing the” small” terms (which will be checked later)

$${\mathcal{E}}_1\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_1-{\sigma }_0\right)\left(\nabla u_0\cdot \nabla v_0-\nabla u_0\left(x_0\right)\cdot \nabla v_0\left(x_0\right)\right)\mathit{dx},$$E21

$${\mathcal{E}}_2\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left({\sigma }_2-{\sigma }_0\right)\left(\nabla u_0\cdot \nabla v_0-\nabla u_0\left(x_0\right)\cdot \nabla v_0\left(x_0\right)\right)\mathit{dx},$$E22

we obtain

$$\begin{array}{l}\left(a_{\epsilon }-a_0\right)\left(u_0,v_{\epsilon }\right)={\mathit{\pi \epsilon }}^2\left(\left(1-{\alpha }^2\right)\left({\sigma }_1-{\sigma }_0\right)+{\alpha }^2\left({\sigma }_2-{\sigma }_0\right)\right)\nabla u_0\left(x_0\right)\cdot \nabla v_0\left(x_0\right)\\ +{\mathcal{F}}_1\left(\epsilon \right)+{\mathcal{F}}_2\left(\epsilon \right)+{\mathcal{E}}_1\left(\epsilon \right)+{\mathcal{E}}_2\left(\epsilon \right),\end{array}$$E23

where

$${\mathcal{F}}_1\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_1-{\sigma }_0\right)\nabla u_0\cdot \nabla \left(v_{\epsilon }-v_0\right)\mathit{dx},$$E24

$${\mathcal{F}}_2\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left({\sigma }_2-{\sigma }_0\right)\nabla u_0\cdot \nabla \left(v_{\epsilon }-v_0\right)\mathit{dx}.$$E25

We will now study the asymptotic of ${\mathcal{F}}_1$ and ${\mathcal{F}}_2$. Introducing the variation ${\tilde{v}}_{\epsilon }≔v_{\epsilon }-v_0$, we obtain from (20) that ${\tilde{v}}_{\epsilon }$ solves

$$\left\{\begin{array}{ll}\mathrm{\Delta }{\tilde{v}}_{\epsilon }=0& \text{in}{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }\cup \left(\mathrm{\Omega }\backslash \overline{{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }}\right),\\ 〚\sigma {\partial }_{\nu }{\tilde{v}}_{\epsilon ,}〛=-\left({\sigma }_1-{\sigma }_0\right)\nabla v_0\cdot \nu & \mathrm{on}{\mathrm{\Gamma }}_1^{\epsilon },\\ 〚\sigma {\partial }_{\nu }{\tilde{v}}_{\epsilon }〛=-\left({\sigma }_2-{\sigma }_1\right)\nabla v_0\cdot \nu & \mathrm{on}{\mathrm{\Gamma }}_2^{\epsilon },\\ {\tilde{v}}_{\epsilon }=0& \mathrm{on}\mathrm{\partial \Omega }.\end{array}\right.$$E26

We set $V≔\nabla v_0\left(x_0\right)$ and we approximate ${\tilde{v}}_{\epsilon }$ by the solution $h_{\epsilon }^V$ of the auxiliary problem

$$\left\{\begin{array}{ll}\mathrm{\Delta }{h}_{\epsilon }^V=0& \text{in}{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }\cup \left({\mathbb{R}}^2\backslash \overline{{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }}\right),\\ 〚\sigma {\partial }_{\nu }{h}_{\epsilon }^V〛=\left({\sigma }_0-{\sigma }_1\right)V\cdot \nu & \mathrm{on}{\mathrm{\Gamma }}_1^{\epsilon },\\ 〚\sigma {\partial }_{\nu }{h}_{\epsilon }^V〛=\left({\sigma }_1-{\sigma }_2\right)V\cdot \nu & \mathrm{on}{\mathrm{\Gamma }}_2^{\epsilon },\\ h_{\epsilon }^V\to 0\mathrm{at}& \infty .\end{array}\right.$$E27

By shifting the coordinate system, we can assume for simplicity that $x_0=0$. For our case, we can compute explicitly the function $h_{\epsilon }^V$ using polar coordinates:

$$h_{\epsilon }^V\left(x\right)=\left\{\begin{array}{ll}\left(\beta +\gamma \right)V\cdot x& \text{if}x\in {\mathrm{\Omega }}_2^{\epsilon },\\ \mathit{\beta V}\cdot x+\gamma {\left(\mathit{\alpha \epsilon }\right)}^2\frac{V\cdot x}{{\left|x\right|}^2}& \text{if}x\in {\mathrm{\Omega }}_1^{\epsilon },\\ \left({\mathit{\beta \epsilon }}^2+\gamma {\left(\mathit{\alpha \epsilon }\right)}^2\right)\frac{V\cdot x}{{\left|x\right|}^2}& \text{if}x\in {\mathbb{R}}^2\backslash \overline{{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }},\end{array}\right.$$E28

where

$$\beta ≔\frac{\left({\sigma }_1+{\sigma }_2\right)\left({\sigma }_0-{\sigma }_1\right)-{\alpha }^2\left({\sigma }_1-{\sigma }_2\right)\left({\sigma }_0-{\sigma }_1\right)}{\left({\sigma }_1+{\sigma }_2\right)\left({\sigma }_1+{\sigma }_0\right)+{\alpha }^2\left({\sigma }_1-{\sigma }_2\right)\left({\sigma }_0-{\sigma }_1\right)},$$

and

$$\gamma ≔\frac{2{\sigma }_0\left({\sigma }_1-{\sigma }_2\right)}{\left({\sigma }_1+{\sigma }_2\right)\left({\sigma }_1+{\sigma }_0\right)+{\alpha }^2\left({\sigma }_1-{\sigma }_2\right)\left({\sigma }_0-{\sigma }_1\right)}.$$

Its gradient is given by

$$\nabla h_{\epsilon }^V=\left\{\begin{array}{ll}\left(\beta +\gamma \right)V& \text{if}x\in {\mathrm{\Omega }}_2^{\epsilon },\\ \mathit{\beta V}+\gamma {\left(\mathit{\alpha \epsilon }\right)}^2\left[\frac{V}{{\left|x\right|}^2}-2V\cdot x\frac{x}{{\left|x\right|}^4}\right]& \text{if}x\in {\mathrm{\Omega }}_1^{\epsilon },\\ \left({\mathit{\beta \epsilon }}^2+\gamma {\left(\mathit{\alpha \epsilon }\right)}^2\right)\left[\frac{V}{{\left|x\right|}^2}-2V\cdot x\frac{x}{{\left|x\right|}^4}\right]& \text{if}x\in {\mathbb{R}}^2\backslash \overline{{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }}.\end{array}\right.$$E29

Denoting

$${\mathcal{E}}_3\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_1-{\sigma }_0\right)\nabla u_0\cdot \nabla \left({\tilde{v}}_{\epsilon }-h_{\epsilon }^V\right)\mathit{dx},$$E30

$${\mathcal{E}}_4\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_1-{\sigma }_0\right)\left(\nabla u_0-\nabla u_0\left(x_0\right)\right)\cdot \nabla h_{\epsilon }^V\mathit{dx},$$E31

$${\mathcal{E}}_5\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left({\sigma }_2-{\sigma }_0\right)\nabla u_0\cdot \nabla \left({\tilde{v}}_{\epsilon }-h_{\epsilon ,}^V\right)\mathit{dx},$$E32

and

$${\mathcal{E}}_6\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left({\sigma }_2-{\sigma }_0\right)\left(\nabla u_0-\nabla u_0\left(x_0\right)\right)\cdot \nabla h_{\epsilon }^V\mathit{dx}.$$E33

Then we obtain

$$\begin{array}{c}{\mathcal{F}}_1\left(\epsilon \right)={\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_1-{\sigma }_0\right)\nabla u_0\left(x_0\right)\cdot \nabla h_{\epsilon }^V\mathit{dx}+{\mathcal{E}}_3\left(\epsilon \right)+{\mathcal{E}}_4\left(\epsilon \right)\\ =\left({\sigma }_1-{\sigma }_0\right)\nabla u_0\left(x_0\right)\cdot {\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\nabla h_{\epsilon }^V\mathit{dx}+{\mathcal{E}}_3\left(\epsilon \right)+{\mathcal{E}}_4\left(\epsilon \right).\end{array}$$

Using polar coordinates and integrating by parts, yields

$${\mathcal{F}}_1\left(\epsilon \right)=\pi \left(1-{\alpha }^2\right){\epsilon }^2\beta \left({\sigma }_1-{\sigma }_0\right)\nabla u_0\left(x_0\right)\cdot V+{\mathcal{E}}_3\left(\epsilon \right)+{\mathcal{E}}_4\left(\epsilon \right),$$

$$\begin{array}{cc}{\mathcal{F}}_2\left(\epsilon \right)& ={\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left({\sigma }_2-{\sigma }_0\right)\nabla u_0\left(x_0\right)\cdot \nabla h_{\epsilon }^V\mathit{dx}+{\mathcal{E}}_5\left(\epsilon \right)+{\mathcal{E}}_6\left(\epsilon \right)\\ & ={\mathit{\pi \alpha }}^2{\epsilon }^2\left(\beta +\gamma \right)\left({\sigma }_2-{\sigma }_0\right)\nabla u_0\left(x_0\right)\cdot V+{\mathcal{E}}_5\left(\epsilon \right)+{\mathcal{E}}_6\left(\epsilon \right).\end{array}$$

After rearrangement, we get

$$\left(a_{\epsilon }-a_0\right)\left(u_0,v_{\epsilon }\right)={\mathit{\pi \epsilon }}^2\mathrm{\Lambda }\nabla u_0\left(x_0\right)\cdot \nabla v_0\left(x_0\right)+\sum _{i=1}^6{\mathcal{E}}_i,$$

where

$$\begin{array}{l}\mathrm{\Lambda }≔\left(1-{\alpha }^2\right)\left(1+\beta \right)\left({\sigma }_1-{\sigma }_0\right)+{\alpha }^2\left({\sigma }_2-{\sigma }_0\right)\left(1+\beta +\gamma \right)\\ =\frac{2\left(1-{\alpha }^2\right){\sigma }_0\left({\sigma }_1-{\sigma }_0\right)\left({\sigma }_1+{\sigma }_2\right)+4{\alpha }^2{\sigma }_0{\sigma }_1\left({\sigma }_2-{\sigma }_0\right)}{\left({\sigma }_1+{\sigma }_2\right)\left({\sigma }_1+{\sigma }_0\right)-{\alpha }^2\left({\sigma }_1-{\sigma }_2\right)\left({\sigma }_1-{\sigma }_0\right)}.\end{array}$$E34

3.1.2 Variation of the linear form

Let us now turn to the asymptotic analysis of the variation

$$\left(l_{\epsilon }-l_0\right)\left(v_{\epsilon }\right)={\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left(f_1-f_0\right)v_{\epsilon }\mathit{dx}+{\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left(f_2-f_0\right)v_{\epsilon ,}\mathit{dx}.$$E35

We can rewrite (35) as

$$\left(l_{\epsilon }-l_0\right)\left(v_{\epsilon }\right)={\mathit{\pi \epsilon }}^2\left[\left(1-{\alpha }^2\right)\left(f_1-f_0\right)+{\alpha }^2\left(f_2-f_0\right)\right]v_0\left(x_0\right)+{\mathcal{E}}_7\left(\epsilon \right)+{\mathcal{E}}_8\left(\epsilon \right),$$

where

$${\mathcal{E}}_7\left(\epsilon \right)={\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left(f_1-f_0\right){\tilde{v}}_{\epsilon }\mathit{dx},{\mathcal{E}}_8\left(\epsilon \right)={\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left(f_2-f_0\right){\tilde{v}}_{\epsilon }\mathit{dx}.$$

Again, it will be proved that ${\mathcal{E}}_7$ and ${\mathcal{E}}_8$ are small terms. Consequently we set

$$\mathit{\delta l}=\pi \left[\left(1-{\alpha }^2\right)\left(f_1-f_0\right)+{\alpha }^2\left(f_2-f_0\right)\right]v_0\left(x_0\right).$$

3.1.3 Variation of the cost function

Expression of $J_{p,\epsilon }\left(u_{\epsilon }\right)-J_{p,\epsilon }\left(u_0\right)$. For simplicity of the calculus, we assume that $p$ is even, then we have

$$\begin{array}{cc}{J}_{p,\epsilon }\left(u_{\epsilon }\right)-J_{p,\epsilon }\left(u_0\right)& =\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|u_{\epsilon }\right|}^p\mathit{dx}-\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|u_0\right|}^p\mathit{dx}\\ & =\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|\left(u_{\epsilon }-u_0\right)+u_0\right|}^p-\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|u_0\right|}^p\mathit{dx}\\ & =\frac{1}{p}\sum _{k=0}^p\left(\begin{array}{c}p\\ k\end{array}\right){\int }_{\mathrm{\Omega }}{u}_0^{p-k}{\left(u_{\epsilon }-u_0\right)}^k\mathit{dx}-\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|u_0\right|}^p\mathit{dx}\\ & =\frac{1}{p}\sum _{k=2}^p\left(\begin{array}{c}p\\ k\end{array}\right){\int }_{\mathrm{\Omega }}{u}_0^{p-k}{\left(u_{\epsilon }-u_0\right)}^k\mathit{dx}+{\int }_{\mathrm{\Omega }}{u}_0^{p-1}\left(u_{\epsilon }-u_0\right)\mathit{dx}.\end{array}$$

Therefore

$$J_{p,\epsilon }\left(u_{\epsilon }\right)-J_{p,\epsilon }\left(u_0\right)-{\mathit{DJ}}_{p,\epsilon }\left(u_0\right)\left(u_{\epsilon }-u_0\right)={\mathcal{E}}_9\left(\epsilon \right),$$

where

$${\mathcal{E}}_9\left(\epsilon \right)≔\frac{1}{p}\sum _{k=2}^p\left(\begin{array}{c}p\\ k\end{array}\right){\int }_{\mathrm{\Omega }}{u}_0^{p-k}{\left(u_{\epsilon }-u_0\right)}^k\mathit{dx}.$$

We will prove in the next section that ${\mathcal{E}}_9\left(\epsilon \right)=o\left({\epsilon }^2\right)$, and thus $\delta J_1=0$.

Expression of $J_{p,\epsilon }\left(u_0\right)-J_{p,0}\left(u_0\right)$. We have

$$J_{p,\epsilon }\left(u_0\right)-J_{p,0}\left(u_0\right)=\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|u_0\right|}^p\mathit{dx}-\frac{1}{p}{\int }_{\mathrm{\Omega }}{\left|u_0\right|}^p\mathit{dx}=0.$$

Consequently $\delta J_2=0$. Now, we are ready to state the main result of this paper.

Theorem 1.1 The topological asymptotic expansion of the functional $J$ with respect to the insertion of small coated inclusion ${\omega }_{\epsilon }$ is given by

$$J_{p,\epsilon }\left(u_{\epsilon }\right)-J_{p,0}\left(u_0\right)={\epsilon }^{2}G\left(x_0\right)+o\left({\epsilon }^2\right),$$

where

$$G\left(x_0\right)=\pi \mathrm{\Lambda }\nabla u_0\left(x_0\right)\cdot \nabla v_0\left(x_0\right)+\pi \left[\left(1-{\alpha }^2\right)\left(f_1-f_0\right)+{\alpha }^2\left(f_2-f_0\right)\right]v_0\left(x_0\right),$$E36

and

$$\mathrm{\Lambda }≔\frac{2\left(1-{\alpha }^2\right){\sigma }_0\left({\sigma }_1-{\sigma }_0\right)\left({\sigma }_1+{\sigma }_2\right)+4{\alpha }^2{\sigma }_0{\sigma }_1\left({\sigma }_2-{\sigma }_0\right)}{\left({\sigma }_1+{\sigma }_2\right)\left({\sigma }_1+{\sigma }_0\right)-{\alpha }^2\left({\sigma }_1-{\sigma }_2\right)\left({\sigma }_1-{\sigma }_0\right)}.$$

Remark 1 When ${\sigma }_1={\sigma }_2$ and $\alpha =0$, the topological derivative defined in (36) becomes

$$G\left(x_0\right)=2{\mathit{\pi \rho \sigma }}_0\nabla u_0\left(x_0\right)\cdot \nabla v_0\left(x_0\right)+\pi \left(f_1-f_0\right)v_0\left(x_0\right),\rho =\frac{{\sigma }_1-{\sigma }_0}{{\sigma }_1+{\sigma }_0}.$$E37

Expression (37) is known in the literature when the inclusion $\omega$ is an homogenous disk; see for instance ([17], Thm 4.3).

4. Estimates of the remainders

In this section the estimation for the remainders on the topological asymptotic expansion are presented. The results are derived by using simple arguments from functional analysis.

4.1 Preliminary lemmas

Lemma 1

i. For any vector $V\in {\mathbb{R}}^2$, $x_0\in \mathrm{\Omega }$ and positive radius $R$, we have

$$\parallel h_{\epsilon }^V{\parallel }_{L^2\left(\mathrm{\Omega }\right)}=O\left({\epsilon }^{3/2}\right),$$

$$\parallel \nabla h_{\epsilon }^V{\parallel }_{L^p\left(\mathrm{\Omega }\right)}=O\left({\epsilon }^{2/p}\right)\forall p>1,$$

$$\parallel h_{\epsilon }^V{\parallel }_{L^p\left(\mathrm{\Omega }\backslash B\left(x_0,R\right)\right)}+\parallel \nabla h_{\epsilon }^V{\parallel }_{L^p\left(\mathrm{\Omega }\backslash B\left(x_0,R\right)\right)}=O\left({\epsilon }^2\right)\forall p\ge 1.$$

ii. Given a function $\psi :\mathrm{\Omega }\to {\mathbb{R}}^2$ which is $\theta$ Hölder continuous ($0<\theta <1$) in a neighborhood of $x_0$ and consider the solution $w_{\epsilon }$ of the system:

$$\left\{\begin{array}{ll}-div\left({\sigma }_{\epsilon }\nabla w_{\epsilon }\right)=0& \text{in}{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }\cup \left(\mathrm{\Omega }\backslash \overline{{\mathrm{\Omega }}_1^{\epsilon }\cup {\mathrm{\Omega }}_2^{\epsilon }}\right),\\ 〚\sigma {\partial }_{\nu }{w}_{\epsilon ,}〛=\left({\sigma }_0-{\sigma }_1\right)\psi \cdot \nu & \mathrm{on}{\mathrm{\Gamma }}_1^{\epsilon },\\ 〚\sigma {\partial }_{\nu }{w}_{\epsilon }〛=\left({\sigma }_1-{\sigma }_2\right)\psi \cdot \nu & \mathrm{on}{\mathrm{\Gamma }}_2^{\epsilon },\\ w_{\epsilon }=0& \mathrm{on}\mathrm{\partial \Omega }.\end{array}\right.$$E38

Then, we have

$$\parallel w_{\epsilon }-h_{\epsilon }^{\psi \left(x_0\right)}{\parallel }_{H^1\left(\mathrm{\Omega }\right)}=o\left(\epsilon \right).$$E39

Proof. The estimates $i)$ in Lemma 1 follow directly from (28) and (29). Now, we prove the second part. Let $\phi \in H_0^1\left(\mathrm{\Omega }\right)$ be an arbitrary test function, then from (38), we have

$${\int }_{\mathrm{\Omega }}{\sigma }_{\epsilon }\nabla w_{\epsilon }\cdot \nabla \phi \mathit{dx}=\left({\sigma }_0-{\sigma }_1\right){\int }_{{\mathrm{\Gamma }}_1^{\epsilon }}\psi \cdot \mathit{\nu \phi }\mathit{ds}+\left({\sigma }_1-{\sigma }_2\right){\int }_{{\mathrm{\Gamma }}_2^{\epsilon }}\psi \cdot \mathit{\nu \phi }\mathit{ds}.$$

Using Green’s formula together with (27), we obtain

$${\int }_{\mathrm{\Omega }}{\sigma }_{\epsilon }\nabla h_{\epsilon }^{\psi \left(x_0\right)}\cdot \nabla \phi \mathit{dx}=\left({\sigma }_0-{\sigma }_1\right){\int }_{{\mathrm{\Gamma }}_1^{\epsilon }}\psi \left(x_0\right)\cdot \mathit{\nu \phi }\mathit{ds}+\left({\sigma }_1-{\sigma }_2\right){\int }_{{\mathrm{\Gamma }}_2^{\epsilon }}\psi \left(x_0\right)\cdot \mathit{\nu \phi }\mathit{ds}.$$

Denote ${\mathrm{\Theta }}_{\epsilon }≔w_{\epsilon }-h_{\epsilon }^{\psi \left(x_0\right)}$. It follows that

$${\int }_{\mathrm{\Omega }}{\sigma }_{\epsilon }\nabla {\mathrm{\Theta }}_{\epsilon }\cdot \nabla \phi \mathit{dx}=\left({\sigma }_0-{\sigma }_1\right){\int }_{{\mathrm{\Gamma }}_1^{\epsilon }}\left(\psi -\psi \left(x_0\right)\right)\cdot \mathit{\nu \phi }\mathit{ds}+\left({\sigma }_1-{\sigma }_2\right){\int }_{{\mathrm{\Gamma }}_2^{\epsilon }}\left(\psi -\psi \left(x_0\right)\right)\cdot \mathit{\nu \phi }\mathit{ds}.$$E40

Using the change of variable, the $\theta$-Hölder continuity of $\psi$ in the vicinity of $x_0$ and the trace theorem, we get for $\epsilon$ small enough

$$\begin{array}{ll}\left|{\int }_{{\mathrm{\Gamma }}_2^{\epsilon }}\left(\psi -\psi \left(x_0\right)\right)\cdot \mathit{\nu \phi }\mathit{ds}\right|& =\epsilon \left|{\int }_{{\mathrm{\Gamma }}_2}\left(\psi \left(\mathit{\epsilon x}\right)-\psi \left(x_0\right)\right)\cdot \mathit{\nu \phi }\left(\mathit{\epsilon x}\right)\mathit{ds}\right|\\ & \le c{\epsilon }^{1+\theta }\parallel \phi \left(\mathit{\epsilon x}\right){\parallel }_{H^{1/2}\left({\mathrm{\Gamma }}_2\right)}\\ & \le c{\epsilon }^{1+\theta }\parallel \phi \left(\mathit{\epsilon x}\right){\parallel }_{H^1\left({\mathrm{\Omega }}_2\right)}\\ & \le c{\epsilon }^{\theta }\parallel \phi {\parallel }_{L^2\left({\mathrm{\Omega }}_2^{\epsilon }\right)}+c{\epsilon }^{\theta +1}\parallel \nabla \phi {\parallel }_{L^2\left({\mathrm{\Omega }}_2^{\epsilon }\right)}.\end{array}$$

From the Hölder inequality and the Sobolev imbedding theorem, we obtain

$$\parallel \phi {\parallel }_{L^2\left({\mathrm{\Omega }}_2^{\epsilon }\right)}\le c{\epsilon }^{\frac{1}{p}}\parallel \phi {\parallel }_{L^{\frac{2p}{p-1}}\left({\mathrm{\Omega }}_2^{\epsilon }\right)}\le c{\epsilon }^{1/p}\parallel \phi {\parallel }_{H^1\left(\mathrm{\Omega }\right)},\text{for}\mathrm{any}p>1.$$

Therefore

$$\left|{\int }_{{\mathrm{\Gamma }}_2^{\epsilon }}\left(\psi -\psi \left(x_0\right)\right)\cdot \mathit{\nu \phi }\mathit{ds}\right|\le c\left({\epsilon }^{\theta +1/p}+{\epsilon }^{\theta +1}\right)\parallel \phi {\parallel }_{H^1\left(\mathrm{\Omega }\right)}.$$

Analogously, we can prove that

$$\left|{\int }_{{\mathrm{\Gamma }}_1^{\epsilon }}\left(\psi -\psi \left(x_0\right)\right)\cdot \mathit{\nu \phi }\mathit{ds}\right|\le c^{\prime }\left({\epsilon }^{\theta +1/p}+{\epsilon }^{\theta +1}\right)\parallel \phi {\parallel }_{H^1\left(\mathrm{\Omega }\right)}.$$

From (40) and the first part of Lemma 1, we deduce that

$$\left|{\int }_{\mathrm{\Omega }}{\sigma }_{\epsilon }\nabla {\mathrm{\Theta }}_{\epsilon }\cdot \nabla \phi \mathit{dx}\right|\le {c^{″}}^{}\left({\epsilon }^{\theta +1/p}+{\epsilon }^{\theta +1}+{\epsilon }^2\right)\parallel \phi {\parallel }_{H^1\left(\mathrm{\Omega }\right)}.$$E41

Choosing $\phi =\mathrm{\Theta }$ and $p\in \left(1,\frac{1}{1-\theta }\right)$ in (41), yield

$$\parallel \mathrm{\Theta }{\parallel }_{H^1\left(\mathrm{\Omega }\right)}=o\left(\epsilon \right),$$

and the proof is completed.

Lemma 2 We have the following estimates

$$\parallel u_{\epsilon }-u_0{\parallel }_{H^1\left(\mathrm{\Omega }\right)}=O\left(\epsilon \right),$$E42

$$\parallel v_{\epsilon }-v_0{\parallel }_{H^1\left(\mathrm{\Omega }\right)}=O\left(\epsilon \right),$$E43

$$\parallel u_{\epsilon }-u_0{\parallel }_{L^2\left(\mathrm{\Omega }\right)}=o\left(\epsilon \right),$$E44

$$\parallel v_{\epsilon }-v_0{\parallel }_{L^2\left(\mathrm{\Omega }\right)}=o\left(\epsilon \right).$$E45

Proof. From the Poincaré inequality, we deduce that

$${\int }_{\mathrm{\Omega }}{\left|u_{\epsilon }-u_0\right|}^2\mathit{dx}\le C\left({\int }_{\mathrm{\Omega }}{\left|\nabla \left(u_{\epsilon }-u_0\right)\right|}^2\mathit{dx}\right),$$E46

for some constant $C$ independent of $\epsilon$. Then, it suffices to show that

$${\int }_{\mathrm{\Omega }}{\left|\nabla \left(u_{\epsilon }-u_0\right)\right|}^2\mathit{dx}\le C{\epsilon }^2.$$

From (9), we obtain immediately that

$$a_{\epsilon }\left(u_{\epsilon }-u_0,v\right)=-\left(a_{\epsilon }-a_0\right)\left(u_0,v\right),\forall v\in H_0^1\left(\mathrm{\Omega }\right).$$E47

According to (19), we get

$$\left(a_{\epsilon }-a_0\right)\left(u_0,v\right)={\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_1-{\sigma }_0\right)\nabla u_0\cdot \nabla v\mathit{dx}+{\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left({\sigma }_2-{\sigma }_0\right)\nabla u_0\cdot \nabla v\mathit{dx}.$$

Using the fact that $\nabla u_0$ is uniformly bounded on ${\mathrm{\Omega }}_1^{\epsilon }$ and ${\mathrm{\Omega }}_2^{\epsilon }$, we obtain

$$\begin{array}{l}\mid \left(a_{\epsilon }-a_0\right)\left(u_0,v\right)\mid \le \mid {\sigma }_1-{\sigma }_0\mid \underset{{\mathrm{\Omega }}_1^{\epsilon }}{sup}{\left|\nabla u_0\Vert {\mathrm{\Omega }}_1^{\epsilon }\right|}^{1/2}\parallel \nabla v{\parallel }_{L^2\left(\mathrm{\Omega }\right)}\\ +\mid {\sigma }_2-{\sigma }_0\mid \underset{{\mathrm{\Omega }}_2^{\epsilon }}{sup}{\left|\nabla u_0\Vert {\mathrm{\Omega }}_2^{\epsilon }\right|}^{1/2}\parallel \nabla v{\parallel }_{L^2\left(\mathrm{\Omega }\right)}\\ \le \mathit{C\epsilon }\parallel \nabla v{\parallel }_{L^2\left(\mathrm{\Omega }\right)},\end{array}$$

and from (47), we obtain

$${\int }_{\mathrm{\Omega }}{\left|\nabla \left(u_{\epsilon }-u_0\right)\right|}^2\mathit{dx}\le C{\epsilon }^2.$$

This proves the asymptotic formula (42). Analogously we derive the estimate $\left(43\right)$. The proof of (44) and (45) follows straightforwardly from [4, Lemma 9.3].

4.2 Asymptotic behavior of the remainders

In this subsection, we shall prove that ${\mathcal{E}}_i\left(\epsilon \right)=o\left({\epsilon }^2\right)$ for $i=1\dots 9$. We have

$${\mathcal{E}}_2\left(\epsilon \right)={\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left(\nabla u_0\cdot \nabla v_0-\nabla u_0\left(x_0\right)\cdot \nabla v_0\left(x_0\right)\right)\mathit{dx}.$$

Using the regularity of $u_0$ and $v_0$ near $x_0$ and Taylor-Lagrange expansion, we straightforwardly obtain ${\mathcal{E}}_2\left(\epsilon \right)\le c{\epsilon }^3,$ and thus ${\mathcal{E}}_2\left(\epsilon \right)=o\left({\epsilon }^2\right)$. Similarly, we can prove that ${\mathcal{E}}_1\left(\epsilon \right)=o\left({\epsilon }^2\right).$ Let’s now prove that ${\mathcal{E}}_4\left(\epsilon \right)=o\left({\epsilon }^2\right)$ and ${\mathcal{E}}_6\left(\epsilon \right)=o\left({\epsilon }^2\right)$. We have

$${\mathcal{E}}_4\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_1-{\sigma }_0\right)\left(\nabla u_0-\nabla u_0\left(x_0\right)\right)\cdot \nabla h_{\epsilon }^V\mathit{dx},$$

and

$${\mathcal{E}}_6\left(\epsilon \right)≔{\int }_{{\mathrm{\Omega }}_2^{\epsilon }}\left({\sigma }_2-{\sigma }_0\right)\left(\nabla u_0-\nabla u_0\left(x_0\right)\right)\cdot \nabla h_{\epsilon }^V\mathit{dx}.$$

Using Cauchy-Schwarz inequality, yields

$$\begin{array}{l}\mid {\mathcal{E}}_6\left(\epsilon \right)\mid \le \mid {\sigma }_2-{\sigma }_0\mid {\left({\int }_{{\mathrm{\Omega }}_2^{\epsilon }}{\left|\left(\nabla u_0-\nabla u_0\left(x_0\right)\right)\right|}^2\mathit{dx}\right)}^{1/2}{\left({\int }_{{\mathrm{\Omega }}_2^{\epsilon }}{\left|\nabla h_{\epsilon }^V\right|}^2\mathit{dx}\right)}^{1/2}\\ \le \mid \left(\beta +\gamma \right)\left({\sigma }_2-{\sigma }_0\right)\mid \parallel V\parallel \alpha \sqrt{\pi }\epsilon {\left({\int }_{{\mathrm{\Omega }}_2^{\epsilon }}{\left|\left(\nabla u_0-\nabla u_0\left(x_0\right)\right)\right|}^2\mathit{dx}\right)}^{1/2}.\end{array}$$

From the regularity of $u_0$ near $x_0$ and Taylor-Lagrange expansion, we obtain the bound $\mid {\mathcal{E}}_6\left(\epsilon \right)\mid \le c{\epsilon }^{5/2}.$ Analogously, we can show that ${\mathcal{E}}_4\left(\epsilon \right)=o\left({\epsilon }^2\right).$ Let’s now focus on ${\mathcal{E}}_3$ and ${\mathcal{E}}_5\left(\epsilon \right)$. Using Hölder inequality, we obtain

$$\begin{array}{cc}\left|{\mathcal{E}}_3\left(\epsilon \right)\right|& =\left|{\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\left({\sigma }_0-{\sigma }_1\right)\nabla u_0\cdot \nabla \left({\tilde{v}}_{\epsilon }-h_{\epsilon }^V\right)\mathit{dx}\right|\\ & \le \mid {\sigma }_0-{\sigma }_1\mid \underset{x\in \mathrm{\Omega }}{sup}{\left|\nabla u_0\left(x\right)\Vert {\mathrm{\Omega }}_1^{\epsilon }\right|}^{1/2}{‖\nabla \left({\tilde{v}}_{\epsilon }-h_{\epsilon }^V\right)‖}_{L^2\left({\mathrm{\Omega }}_1^{\epsilon }\right)}\\ & \le \mathit{C\epsilon }{‖\nabla \left({\tilde{v}}_{\epsilon }-h_{\epsilon }^V\right)‖}_{L^2\left({\mathrm{\Omega }}_1^{\epsilon }\right)}.\end{array}$$

From Lemma 1, we deduce that

$$\parallel \nabla \left({\tilde{v}}_{\epsilon }-h_{\epsilon }^V\right){\parallel }_{L^2\left({\mathrm{\Omega }}_1^{\epsilon }\right)}\le \parallel \nabla \left({\tilde{v}}_{\epsilon }-h_{\epsilon }^V\right){\parallel }_{L^2\left(\mathrm{\Omega }\right)}=o\left(\epsilon \right).$$

Therefore, we conclude that ${\mathcal{E}}_3\left(\epsilon \right)=o\left({\epsilon }^2\right)$. By the same techniques, we prove that

$${\mathcal{E}}_5\left(\epsilon \right)=o\left({\epsilon }^2\right).$$

Let’s now check that ${\mathcal{E}}_7\left(\epsilon \right)=o\left({\epsilon }^2\right)$ and ${\mathcal{E}}_8\left(\epsilon \right)=o\left({\epsilon }^2\right)$. We have

$$\mid {\mathcal{E}}_7\left(\epsilon \right)\mid \le c{\int }_{{\mathrm{\Omega }}_1^{\epsilon }}\mid v_{\epsilon }-v_0\mid \mathit{dx}.$$

From the Hölder inequality, we obtain

$$\mid {\mathcal{E}}_7\left(\epsilon \right)\mid \le c{\epsilon }^{2/r}\parallel v_{\epsilon }-v_0{\parallel }_{L^s\left({\mathrm{\Omega }}_1^{\epsilon }\right)}\le c{\epsilon }^{2/r}\parallel v_{\epsilon }-v_0{\parallel }_{L^s\left(\mathrm{\Omega }\right)},$$

for all $r,s\in \left[1,+\infty \right]$ satisfying $\frac{1}{r}+\frac{1}{s}=1$. Due to Lemma 2, we conclude that ${\mathcal{E}}_7\left(\epsilon \right)=o\left({\epsilon }^2\right)$. Similarly, we obtain

$${\mathcal{E}}_8\left(\epsilon \right)=o\left({\epsilon }^2\right).$$

Let’s now focus on

$${\mathcal{E}}_9\left(\epsilon \right)≔\frac{1}{p}\sum _{k=2}^p\left(\begin{array}{c}p\\ k\end{array}\right){\int }_{\mathrm{\Omega }}{u}_0^{p-k}{\left(u_{\epsilon }-u_0\right)}^k\mathit{dx}.$$

The Hölder inequality combined with the estimates in Lemma 2, yield

$${\mathcal{E}}_9\left(\epsilon \right)=o\left({\epsilon }^2\right).$$

5. Numerical experiments

For the numerical computation of the state and the adjoint, we use the finite element method. The computational domain $\mathrm{\Omega }$ is the unit disk centered at the origin. The term source and the conductivities are given by

$$f_0=0,f_1=0,f_2=20\text{and}{\sigma }_0=1,{\sigma }_1=0.5,{\sigma }_2=30.$$

All the numerical computations are done under Matlab R2018a. To solve the optimization problem, we apply a fast non-iterative algorithm, based on the following steps:

Non-iterative algorithm.

1. Compute the solution of the problem (8) and the solution of the adjoint problem (20) in the unperturbed domain.

2. Compute the topological gradient $G$ in (36).

3. Find the point $x_0$ where the topological derivative is the most negative.

4. Locate the inclusion $\omega$ at the point $x_0$.

5.1 Example 1

In this example, the Dirichlet boundary condition is given by $g=sin\left(\theta \right),\theta \in \left[0,2\pi \right]$. We present some numerical experiments using the functional $J_p$ for different values of $p$ ($p=\mathrm{2,10,50}$) and the functional $K$. The coated inclusion $\omega$ to be located in order to minimize the objective function, is composed of two concentric disks ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ with radius $r_1=0.2$ and $r_2=0.1$. We compute the position $x_0$ of $\omega$ using the proposed algorithm.

In Figures 2–5, the $x_1$-coordinate of the center of the inclusion is fixed to zero and the $x_2$-coordinate $\le 0$. We observe that the shape functional is decreasing with respect to the variation of the second coordinate $x_2$. Figures 6–8 show the image of topological gradient and the image of temperature distribution after the minimization process. Figure 9 shows the image of the topological gradient corresponding to the shape functional $J_{50}$. For $p\ge 50$, we observe that the shape functional tends to zero and the image of the topological gradient is most negative on the boundary $\mathrm{\partial \Omega }$. A suitable boundary condition that takes into account the effect of radiation and convection could be relevant in this case to solve properly the minimization problem.

Figure 2.

Values of the objective function $J_2$ for variation of the $x_2$-coordinate of the center of the inclusion.

Figure 3.

Values of the objective function $J_{10}$ for variation of the $x_2$-coordinate of the center of the inclusion.

Figure 4.

Values of the objective function $J_{50}$ for variation of the $x_2$-coordinate of the center of the inclusion.

Figure 5.

Values of the objective function $K$ for variation of the $x_2$-coordinate of the center of the inclusion.

Figure 6.

On the left the topological derivative of the functional $J_2$ and on the right the temperature distribution relative to the position $x_0=\left(-0.0115,-0.6915\right)$ of the coated inclusion given by Algorithm.

Figure 7.

On the left the topological derivative of the functional $J_{10}$ and on the right the temperature distribution. $x_0=\left(-0.0035,-0.9109\right)$ is the position given by Algorithm. In order to locate the inclusion far from the boundary $\mathrm{\partial \Omega }$, we have taken an approximation of $x_0$, that is $x_a=\left(0,-0.75\right)$ for the optimization process.

Figure 8.

On the left the topological derivative of the functional $K$ and on the right the temperature distribution with respect the position $x_0=\left(-0.0046,-0.5850\right)$ given by Algorithm.

Figure 9.

The topological derivative of the functional $J_{50}$.

5.2 Example 2

In this example, the Dirichlet boundary condition is given by $g=sin\left(2\theta \right),\theta \in \left[0,2\pi \right]$. We show some numerical experiments in the case of two coated inclusions using the proposed algorithm corresponding to the shape functional $K$.

Figure 10 depicts the image of the topological gradient of the shape functional $K$ and the image of temperature distribution after the optimization process.

Figure 10.

On the left the topological derivative of the functional $K$ and on the right the temperature distribution when the positions of the coated inclusions are given by Algorith.

6. Conclusion

In this chapter, we have considered an optimization problem for a heat distribution in order to get a suitable thermal environment. We have performed the asymptotic expansion of the proposed shape functions with respect to the insertion of small coated inclusions. We have used a one shot algorithm based on the topological derivative to solve numerically the optimization problem. Numerical results are presented showing the efficiency of method. The proposed method can be extend to solve practical problems, like the thermal optimization of electrical multicables.