Asymptotics of the stress concentration in high-contrast elastic composites

Because the aim of this paper is to study the asymptotic behavior of $\nabla u$ in the narrow region between two adjacent inclusions, we may without loss of generality restrict our attention to a situation with only two adjacent inclusions. The basic notations used in this paper follow from [10].

We use $x=(x^{\prime},x_{d})$ to denote a point in $\mathbb{R}^{d}$, where $x^{\prime}=(x_{1},\cdots,x_{d-1})$ and $d=2,3$. Let $D$ be a bounded open set in $\mathbb{R}^{d}$ with $C^{2}$ boundary. $D_{1}$ and $D_{2}$ are two disjoint convex open subsets in $D$ with $C^{2,\gamma}$ $(0<\gamma<1)$ boundaries, $\varepsilon$-apart, and far away from $\partial D$. That is,

$$\begin{split}\overline{D}_{1},\overline{D}_{2}\subset D,\quad\varepsilon:=\mbox{dist}(D_{1},D_{2})>0,\quad\mbox{dist}(D_{1}\cup D_{2},\partial D)>\kappa_{0}>0,\end{split}$$

where $\kappa_{0}$ is a constant independent of $\varepsilon$. We also assume that the $C^{2,\gamma}$ norms of $\partial{D}_{1}$, $\partial{D}_{2}$, and $\partial{D}$ are bounded by some positive constant independent of $\varepsilon$. Set

$$\Omega:=D\setminus\overline{D_{1}\cup D_{2}}.$$

Assume that $\Omega$ and $D_{1}\cup D_{2}$ are occupied, respectively, by two different isotropic and homogeneous materials with different Lamé constants $(\lambda,\mu)$ and $(\lambda_{1},\mu_{1})$. Then the elasticity tensors for the background and the inclusion can be written, respectively, as $\mathbb{C}^{0}$ and $\mathbb{C}^{1}$, with

$$C_{ijkl}^{0}=\lambda\delta_{ij}\delta_{kl}+\mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}),$$

and

$$C_{ijkl}^{1}=\lambda_{1}\delta_{ij}\delta_{kl}+\mu_{1}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}),$$

where $i,j,k,l=1,2,\cdots,d$ and $\delta_{ij}$ is the kronecker symbol: $\delta_{ij}=0$ for $i\neq j$, $\delta_{ij}=1$ for $i=j$. Let $u=(u^{1},u^{2},\cdots,u^{d})^{\mathrm{T}}:D\rightarrow\mathbb{R}^{d}$ denote the displacement field. For a given vector-valued function $\varphi=(\varphi^{1},\varphi^{2},\cdots,\varphi^{d})^{\mathrm{T}}$, we consider the following Dirichlet problem for the Lamé system with pieceiwise constant coefficients:

$\displaystyle\begin{cases}\nabla\cdot\left((\chi_{\Omega}\mathbb{C}^{0}+\chi_{D_{1}\cup{D}_{2}}\mathbb{C}^{1})e(u)\right)=0,&\hbox{in}~{}~{}D,\\ u=\varphi,&\hbox{on}~{}~{}\partial{D},\end{cases}$

(2.1)

where $\chi_{\Omega}$ is the characteristic function of $\Omega\subset\mathbb{R}^{d}$, and

$$e(u)=\frac{1}{2}(\nabla u+(\nabla u)^{\mathrm{T}})$$

is the strain tensor.

Assume that the standard ellipticity condition holds for (2.1), that is,

$\displaystyle\mu>0,\quad d\lambda+2\mu>0,\quad\mu_{1}>0,\quad d\lambda_{1}+2\mu_{1}>0.$

For $\varphi\in C^{1,\gamma}(\partial D;\mathbb{R}^{d})$, it is well known that there exists a unique solution $u\in H^{1}(D;\mathbb{R}^{d})$ to the Dirichlet problem (2.1), which is also the minimizer of the energy functional

$$J_{1}[u]:=\frac{1}{2}\int_{\Omega}\left((\chi_{\Omega}\mathbb{C}^{0}+\chi_{D_{1}\cup{D}_{2}}\mathbb{C}^{1})e(u),e(u)\right)dx$$

on

$\displaystyle H^{1}_{\varphi}(D;\mathbb{R}^{d}):=\left\{u\in H^{1}(D;\mathbb{R}^{d})~{}\big{|}~{}u-\varphi\in H^{1}_{0}(D;\mathbb{R}^{d})\right\}.$

As mentioned previously, Li and Nirenberg [41] proved that $\nabla u$ is uniformly bounded with respect to $\varepsilon$. But, in high-contrast composite media, the concentration of $\nabla u$ is a very usual phenomenon when the distance $\varepsilon$ is sufficiently small. In order to investigate the role of $\varepsilon$ in such concentration phenomenon, let us assume that the Lamé constant in $D_{1}\cup D_{2}$ degenerates to infinite and consider this extreme case. To this end, we first introduce the linear space of rigid displacement in $\mathbb{R}^{d}$:

$$\Psi:=\{\psi\in C^{1}(\mathbb{R}^{d};\mathbb{R}^{d})\ |\ \nabla\psi+(\nabla\psi)^{\mathrm{T}}=0\},$$

with a basis $\left\{\psi_{\alpha}~{}|~{}\alpha=1,2,\cdots,\frac{d(d+1)}{2}\right\}$, namely,

$$\left\{~{}e_{i},~{}x_{j}e_{k}-x_{k}e_{j}~{}\big{|}~{}1\leq\,i\leq\,d,~{}1\leq\,j<k\leq\,d~{}\right\},$$

where $e_{1},\cdots,e_{d}$ denote the standard basis of $\mathbb{R}^{d}$. For fixed $\lambda$ and $\mu$, denoting $u_{\lambda_{1},\mu_{1}}$ as the solution of (2.1), then we have [10]

$\displaystyle u_{\lambda_{1},\mu_{1}}\rightarrow u\quad\hbox{in}\ H^{1}(D;\mathbb{R}^{d}),\quad\hbox{as}\ \min\{\mu_{1},d\lambda_{1}+2\mu_{1}\}\rightarrow\infty,$

where $u$ is the unique $H^{1}(D;\mathbb{R}^{d})$ solution of

$\displaystyle\begin{cases}\mathcal{L}_{\lambda,\mu}u:=\nabla\cdot(\mathbb{C}^{0}e(u))=0,\quad&\hbox{in}\ \Omega,\\ u|_{+}=u|_{-},&\hbox{on}\ \partial{D}_{i},i=1,2,\\ e(u)=0,&\hbox{in}~{}~{}D_{i},i=1,2,\\ \int_{\partial{D}_{i}}\frac{\partial u}{\partial\nu}\Big{|}_{+}\cdot\psi_{\alpha}=0,&i=1,2,\alpha=1,2,\cdots,\frac{d(d+1)}{2},\\ u=\varphi,&\hbox{on}\ \partial{D},\end{cases}$

(2.2)

where

$\displaystyle\frac{\partial u}{\partial\nu}\Big{|}_{+}$

$\displaystyle:=(\mathbb{C}^{0}e(u))\vec{n}=\lambda(\nabla\cdot u)\vec{n}+\mu(\nabla u+(\nabla u)^{\mathrm{T}})\vec{n},$

and $\vec{n}$ is the unit outer normal of $D_{i}$, $i=1,2$. Here and throughout this paper the subscript $\pm$ indicates the limit from outside and inside the domain, respectively. The existence, uniqueness and regularity of weak solutions to (2.2) can be referred to the Appendix of [10]. We note that it suffices to consider the problem (2.2) with $\varphi\in C^{0}(\partial D;\mathbb{R}^{d})$ replaced by $\varphi\in C^{1,\gamma}(\partial D;\mathbb{R}^{d})$. Indeed, it follows from the maximum principle [50] that $\|u\|_{L^{\infty}(D)}\leq C\|\varphi\|_{C^{0}(\partial D)}$. Taking a slightly small domain $\widetilde{D}\subset\subset D$, then in view of the interior derivative estimates for Lamé system, we find that $\tilde{\varphi}:=u\big{|}_{\partial\widetilde{D}}$ satisfies

$$\|\tilde{\varphi}\|_{C^{1,\gamma}(\partial\widetilde{D})}\leq C\|u\|_{L^{\infty}(D)}\leq C\|\varphi\|_{C^{0}(\partial D)}.$$

Without loss of generality, we assume that $\|\varphi\|_{C^{0}(\partial D)}=1$ by considering $u/\|\varphi\|_{C^{0}(\partial D)}$ if $\|\varphi\|_{C^{0}(\partial D)}>0$. If $\varphi\big{|}_{\partial D}=0$, then $u\equiv 0$.

We first point out that problem (2.2) has free boundary value feature. Although $e(u)=0$ implies $u$ in $\overline{D}_{i}$ is linear combination of $\psi_{\alpha}$,

$$u=\sum_{\alpha=1}^{d(d+1)/2}C_{i}^{\alpha}\psi_{\alpha}\quad\text{in}~{}~{}\overline{D}_{i},$$

(2.3)

these $d(d+1)$ constants $C_{i}^{\alpha}$ are free, which will be uniquely determined by $u$. We would like to emphasize that this is exactly the biggest difference with the conductivity model [12], where only two free constants need us to handle in any dimension. It is the increase of the number of free contants that makes elastic problem quite difficult to deal with. Therefore, how to determine such many constants is one of main difficulties we need to solve.

Our strategy in spirit follows from [10, 11]. First, by continuity of $u$ across $\partial D_{i}$, we can decompose the solution of (2.2),

$$u(x)=\sum_{i=1}^{2}\sum_{\alpha=1}^{d(d+1)/2}C_{i}^{\alpha}v_{i}^{\alpha}(x)+v_{0}(x),\quad x\in\,\Omega,$$

(2.4)

where $v_{i}^{\alpha},v_{0}\in{C}^{2}(\Omega;\mathbb{R}^{d})$, respectively, satisfying

$$\begin{cases}\mathcal{L}_{\lambda,\mu}v_{i}^{\alpha}=0,&\mathrm{in}~{}\Omega,\\ v_{i}^{\alpha}=\psi_{\alpha},&\mathrm{on}~{}\partial{D}_{i},\\ v_{i}^{\alpha}=0,&\mathrm{on}~{}\partial{D_{j}}\cup\partial{D},~{}j\neq i,\end{cases}\quad i=1,2,~{}\alpha=1,\cdots,d(d+1)/2,$$

(2.5)

and

$$\begin{cases}\mathcal{L}_{\lambda,\mu}v_{0}=0,&\mathrm{in}~{}\Omega,\\ v_{0}=0,&\mathrm{on}~{}\partial{D}_{1}\cup\partial{D_{2}},\\ v_{0}=\varphi,&\mathrm{on}~{}\partial{D}.\end{cases}$$

(2.6)

So

$$\nabla u(x)=\sum_{i=1}^{2}\sum_{\alpha=1}^{d(d+1)/2}C_{i}^{\alpha}\nabla v_{i}^{\alpha}(x)+\nabla v_{0}(x),\quad x\in\,\Omega.$$

(2.7)

To investigate the symptotic behavior of $\nabla u$, we need both the asymptotic formulas of $\nabla v_{i}^{\alpha}$ and the exact value of $C_{i}^{\alpha}$. To solve $C_{i}^{\alpha}$, from the fourth line in (2.2) and the decomposition (2.4), we have the following linear system of these free constants $C_{i}^{\alpha}$,

$$\sum_{i=1}^{2}\sum\limits_{\alpha=1}^{\frac{d(d+1)}{2}}C_{i}^{\alpha}\int_{\partial D_{j}}\frac{\partial v_{i}^{\alpha}}{\partial\nu}\Big{|}_{+}\cdot\psi_{\beta}+\int_{\partial D_{j}}\frac{\partial v_{0}}{\partial\nu}\Big{|}_{+}\cdot\psi_{\beta}=0,$$

(2.8)

where $j=1,\,2$, $\beta=1,\,\cdots,\frac{d(d+1)}{2}$. But these coefficients are all boundary integrals. By integration by parts,

$\displaystyle a_{ij}^{\alpha\beta}:=-\int_{\partial{D}_{j}}\frac{\partial{v}_{i}^{\alpha}}{\partial\nu}\Big{|}_{+}\cdot\psi_{\beta}=\int_{\Omega}\big{(}\mathbb{C}^{0}e(v_{i}^{\alpha}),e(v_{j}^{\beta})\big{)}\ dx.$

(2.9)

Therefore, in order to solve $C_{i}^{\alpha}$ from (2.8), we have to calculate the energy integral on the right hand side of (2.9). This in turn needs to estimate $\nabla v_{i}^{\alpha}$. In fact, even if we can have the asymptotic formulas of $\nabla v_{i}^{\alpha}$, see Theorem 2.2 and Theorem 2.6 below, it is still hard to solve every $C_{i}^{\alpha}$. It seems to be a mission impossible. To avoid this difficulty, we rewrite (2.7) as

$$\nabla{u}=\sum_{\alpha=1}^{d(d+1)/2}\left(C_{1}^{\alpha}-C_{2}^{\alpha}\right)\nabla{v}_{1}^{\alpha}+\nabla u_{b},\quad\mbox{in}~{}\Omega,$$

(2.10)

where

$$u_{b}:=\sum_{\alpha=1}^{d(d+1)/2}C_{2}^{\alpha}(v_{1}^{\alpha}+v_{2}^{\alpha})+v_{0}.$$

(2.11)

This is because the following boundedness estimates for $|\nabla(v_{1}^{\alpha}+v_{2}^{\alpha})|$ and $|\nabla v_{0}|$, together with the boundedness of $C_{2}^{\alpha}$, $\alpha=1,\cdots,d(d+1)/2$, makes $\nabla u_{b}$ is a “good” term, which has no singularity in the narrow region.

Thus, we reduce the establishment of the asymptotics of $\nabla u$ to that of the asymptotics of $\nabla v_{1}^{\alpha}$, $\alpha=1,\cdots,d(d+1)/2$, and to solving $C_{1}^{\alpha}-C_{2}^{\alpha}$, $\alpha=1,\cdots,d(d+1)/2$. These are two main difficulties that we need to solve in this paper. For the former, we separate all singular terms of $\nabla v_{1}^{\alpha}$, up to constant terms, by using a family of improved auxiliary functions, which depend on the parameters of Lamé system and the geometry informations of $\partial D_{1}$ and $\partial D_{2}$. This essentially improves the results in [10, 11], where only estimates of $|\nabla v_{i}^{\alpha}|$ are obtained. For the latter, we need to characterize the coefficients $a_{ij}^{\alpha\beta}$ to solve the big linear system generated by (2.10), and the convergence of $\int_{\partial D_{1}}\frac{\partial u_{b}}{\partial\nu}\Big{|}_{+}\cdot\psi_{\beta}$, $\beta=1,\cdots,d(d+1)/2$; see Proposition 3.1. Finally, we remark that one crucial step in proving the asymptotics of $C_{1}^{\alpha}-C_{2}^{\alpha}$ is to derive the asymptotic expression of $a_{11}^{\alpha\alpha}$, which is of independent interest, see Proposition 2.3 and Remark 2.4 below. To state it precisely, we further fix our domain and notations.

Recalling the assumptions about $D_{1}$ and $D_{2}$, there exist two points $P_{1}\in\partial D_{1}$ and $P_{2}\in\partial D_{2}$, respectively, such that

$$\mbox{dist}(P_{1},P_{2})=\mbox{dist}(\partial D_{1},\partial D_{2})=\varepsilon.$$

By a translation and rotation of coordinates, if necessary, we suppose that

$$P_{1}=(0^{\prime},\frac{\varepsilon}{2})\in\partial D_{1},\quad P_{2}=(0^{\prime},-\frac{\varepsilon}{2})\in\partial D_{2}.$$

Now we further assume that there exits a constant $R$, independent of $\varepsilon$, such that the portions of $\partial D_{1}$ and $\partial D_{2}$ near $P_{1}$ and $P_{2}$, respectively, can be represented by

$\displaystyle x_{d}=\frac{\varepsilon}{2}+h_{1}(x^{\prime})\quad\mbox{and}\quad x_{d}=-\frac{\varepsilon}{2}+h_{2}(x^{\prime}),\quad\mbox{for}~{}|x^{\prime}|<2R.$

Suppose that the convexity of $\partial D_{1}$ and $\partial D_{2}$ is of order $m\geq 2$ ($m\in\mathbb{N}^{+}$) near the origin,

$\displaystyle h_{i}(x^{\prime})=\begin{cases}(-1)^{i+1}\kappa_{i}|x^{\prime}|^{2}+O(|x^{\prime}|^{2+\gamma}),&\quad m=2,\\ (-1)^{i+1}\kappa_{i}|x^{\prime}|^{m}+O(|x^{\prime}|^{m+1}),&\quad m\geq 3,\end{cases}\quad\mbox{for}~{}i=1,2,$

(2.12)

and

$$\|h_{1}\|_{C^{2,\gamma}(B^{\prime}_{2R})}+\|h_{2}\|_{C^{2,\gamma}(B^{\prime}_{2R})}\leq C,$$

(2.13)

where $\kappa_{1}$, $\kappa_{2}$, and $C$ are constants independent of $\varepsilon$. We call these inclusions $m$-convex inclusions. For simplicity, we assume that $\kappa_{1}=\kappa_{2}=\dfrac{\kappa}{2}$. For $0<r\leq\,2R$, set the narrow region between $\partial{D}_{1}$ and $\partial{D}_{2}$ as

$$\Omega_{r}:=\left\{(x^{\prime},x_{d})\in\mathbb{R}^{d}~{}\big{|}~{}-\frac{\varepsilon}{2}+h_{2}(x^{\prime})<x_{d}<\frac{\varepsilon}{2}+h_{1}(x^{\prime}),~{}|x^{\prime}|<r\right\}.$$

By the standard theory for elliptic systems, we have

$$\|\nabla v_{i}^{\alpha}(x)\|_{L^{\infty}(\Omega\setminus\Omega_{R})}\leq\,C,\quad\alpha=1,\cdots,\frac{d(d+1)}{2},~{}i=1,2.$$

(2.14)

Therefore, in the following we only need to deal with the problems in $\Omega_{R}$. To this end, we denote

$$\delta(x^{\prime}):=\varepsilon+h_{1}(x^{\prime})-h_{2}(x^{\prime}),$$

and introduce a scalar auxiliary function $\bar{u}\in C^{2}(\mathbb{R}^{d})$ such that

$$\bar{u}(x)=\frac{x_{d}+\frac{\varepsilon}{2}-h_{2}(x^{\prime})}{\delta(x^{\prime})},\quad\hbox{in}\ \Omega_{2R},$$

(2.15)

$\bar{u}=1$ on $\partial{D}_{1}$, $\bar{u}=0$ on $\partial{D}_{2}\cup\partial{D}$, and $\|\bar{u}\|_{C^{2}(\Omega\setminus\Omega_{R})}\leq\,C$. Next we use the function $\bar{u}$ to generate a family of vector-valued auxiliary functions $u_{1}^{\alpha}$, $\alpha=1,\cdots,d$, which is much finer than before used in [10, 11].

For $d=2$, recalling that

$$\psi_{1}=\begin{pmatrix}1\\ 0\end{pmatrix}\quad\mbox{and}\quad\psi_{2}=\begin{pmatrix}0\\ 1\end{pmatrix},$$

we define $u_{1}^{\alpha}\in C^{2}(\Omega)$ such that $u_{1}^{\alpha}=v_{1}^{\alpha}$ on $\partial\Omega$, and, in $\Omega_{2R}$

$$\begin{split}u_{1}^{1}&:=\bar{u}_{1}^{1}+\tilde{u}_{1}^{1}:=\bar{u}\psi_{1}+\frac{\lambda+\mu}{\lambda+2\mu}f(\bar{u})\,\delta^{\prime}(x_{1})\,\psi_{2},\\ u_{1}^{2}&:=\bar{u}_{1}^{2}+\tilde{u}_{1}^{2}:=\bar{u}\psi_{2}+\frac{\lambda+\mu}{\mu}f(\bar{u})\,\delta^{\prime}(x_{1})\,\psi_{1},\end{split}$$

(2.16)

where $f(\bar{u}):=\dfrac{1}{2}\Big{(}\bar{u}-\frac{1}{2}\Big{)}^{2}-\dfrac{1}{8}$, $f(\bar{u})=0$ on $(\partial D_{1}\cup\partial D_{2})\cap\{|x_{1}|<R\}$, and $\|u_{1}^{\alpha}\|_{C^{2}(\Omega\setminus\Omega_{R})}\leq C$.

For $d=3$, noting that

$$\psi_{1}=\begin{pmatrix}1\\ 0\\ 0\end{pmatrix},\quad\psi_{2}=\begin{pmatrix}0\\ 1\\ 0\end{pmatrix},\quad\mbox{and}\quad\psi_{3}=\begin{pmatrix}0\\ 0\\ 1\end{pmatrix},$$

we define $u_{1}^{\alpha}\in C^{2}(\Omega)$ such that, in $\Omega_{2R}$

$$\begin{split}u_{1}^{\alpha}&:=\bar{u}_{1}^{\alpha}+\tilde{u}_{1}^{\alpha}:=\bar{u}\psi_{\alpha}+\frac{\lambda+\mu}{\lambda+2\mu}f(\bar{u})\,\partial_{x_{\alpha}}\delta\,\psi_{3},\qquad\qquad\alpha=1,2,\\ u_{1}^{3}&:=\bar{u}_{1}^{3}+\tilde{u}_{1}^{3}:=\bar{u}\psi_{3}+\frac{\lambda+\mu}{\mu}f(\bar{u})\,(\partial_{x_{1}}\delta\,\psi_{1}+\partial_{x_{2}}\delta\,\psi_{2}).\end{split}$$

(2.17)

Notice that the corrected terms $\tilde{u}_{1}^{\alpha}$ depend on the Lamé parameters $\lambda$ and $\mu$, which can help us capture all singular terms of $\nabla v_{i}^{\alpha}$, $\alpha=1,\cdots,d$.

We now use these auxiliary functions $u_{1}^{\alpha}$ to obtain the asymptotics of $\nabla v_{1}^{\alpha}$ in the narrow region $\Omega_{R}$, $\alpha=1,\cdots,d$.

Here we would like to emphasize the importance of the introduction of $\tilde{u}_{1}^{\alpha}$. Although $\bar{u}_{1}^{\alpha}$ can be used to obtain the upper bound estimates of $|\nabla v_{1}^{\alpha}|$ by $|\nabla(v_{1}^{\alpha}-\bar{u}_{1}^{\alpha})|\leq\frac{C}{\sqrt{\varepsilon}}$ derived in [27, Corollary 5.2], it as well shows to us it is possible that there remains more other singular terms in $\nabla(v_{1}^{\alpha}-\bar{u}_{1}^{\alpha})$. As it turns out, the appearance of $\nabla\tilde{u}_{1}^{\alpha}$ can find all the singular terms of $\nabla v_{1}^{\alpha}$ and make $|\nabla(v_{1}^{\alpha}-\bar{u}_{1}^{\alpha}-\tilde{u}_{1}^{\alpha})|$ be bounded. The proof is an adaption of the iteration technique with respect to the energy which was first used in [39] and further developed in [10, 11] to obtain the gradient estimates. We leave it to the Appendix.

The asymptotic expression (2.18) is an essential improvement of the estimate $|\nabla(v_{1}^{\alpha}-\bar{u}_{1}^{\alpha})|\leq\frac{C}{\sqrt{\varepsilon}}$. More importantly, it also allows us to obtain the asymptotics of $a_{11}^{\alpha\alpha}$, $\alpha=1,\cdots,d$, defined in (2.9). This kind of formulas is also the key to establish the asymptotics of $C_{1}^{\alpha}-C_{2}^{\alpha}$, $\alpha=1,\cdots,d$. We here give the results for $m=2$. For the general cases $m\geq 3$, see Section 6 below.

For $\alpha=d+1,\cdots,d(d+1)/2$, we use the auxiliary functions as in [10, 11], that is,

$\displaystyle u_{1}^{\alpha}:=\bar{u}\psi_{\alpha}.$

(2.19)

We remark that Theorems 2.2 and 2.6 also hold true for $v_{2}^{\alpha}$ by replacing $\bar{u}$ with $\underline{u}$, where $\underline{u}\in C^{2}(\mathbb{R}^{d})$ is a scalar function satisfying $\underline{u}=1$ on $\partial{D}_{2}$, $\underline{u}=0$ on $\partial{D}_{1}\cup\partial D$, $\underline{u}=1-\bar{u}$ in $\Omega_{2R}$, and $\|\underline{u}\|_{C^{2}(\Omega\setminus\Omega_{R})}\leq\,C$. In this case we denote the auxiliary functions by $u_{2}^{\alpha}$.

Throughout the paper, unless otherwise stated, we use $C$ to denote some positive constant, whose values may vary from line to line, depending only on $d,\kappa_{0},R$, and an upper bound of the $C^{2,\gamma}$ norms of $\partial D_{1},\partial D_{2}$, and $\partial D$, but not on $\varepsilon$. We call a constant having such dependence a universal constant.

As mentioned in Subsection 2.2, $|\nabla u_{b}|$ is uniformly bounded with respect to $\varepsilon$. In order to derive the asymptotics of the gradient of solution with respect to the sufficiently small parameter $\varepsilon>0$, we consider the case when two inclusions touch each other. Let $u^{*}$ be the solution of

$\displaystyle\begin{cases}\mathcal{L}_{\lambda,\mu}u^{*}=0,\quad&\hbox{in}\ \Omega^{*},\\ u^{*}=\sum_{\alpha=1}^{d(d+1)/2}C_{*}^{\alpha}\psi_{\alpha},&\hbox{on}\ \partial D_{1}^{*}\cup\partial D_{2}^{*},\\ \int_{\partial{D}_{1}^{*}}\frac{\partial{u}^{*}}{\partial\nu}\big{|}_{+}\cdot\psi_{\beta}+\int_{\partial{D}_{2}^{*}}\frac{\partial{u}^{*}}{\partial\nu}\big{|}_{+}\cdot\psi_{\beta}=0,&\beta=1,\cdots,d(d+1)/2,\\ u^{*}=\varphi,&\hbox{on}\ \partial{D},\end{cases}$

(3.1)

where

$$D_{1}^{*}:=\{x\in\mathbb{R}^{d}\big{|}x+P_{1}\in D_{1}\},\quad D_{2}^{*}:=D_{2},\quad\mbox{and}\quad\Omega^{*}:=D\setminus\overline{D_{1}^{*}\cup D_{2}^{*}},$$

and the constants $C_{*}^{\alpha}$, $\alpha=1,\cdots,d(d+1)/2$, are uniquely determined by minimizing the energy

$$\int_{\Omega^{*}}\left(\mathbb{C}^{(0)}e(v),e(v)\right)dx$$

in an admissible function space

$$\mathcal{A}^{*}:=\left\{v\in{H}^{1}(D;\mathbb{R}^{d})~{}\big{|}~{}e(v)=0~{}~{}\mbox{in}~{}{D}^{*}_{1}\cup{D}^{*}_{2},~{}\mbox{and}~{}v=\varphi~{}\mbox{on}~{}\partial{D}\right\}.$$

We emphasize that we use the third line of (3.1), which implies that total flux of $u^{*}$ along the boundaries of both two inclusions is zero, to make a distinction with the forth line of (2.2). This kind of limit function for conductivity problem was also used in [24, 26, 29, 36]. We define the functionals of $\varphi$, for $\beta=1,\cdots,d(d+1)/2$,

$\displaystyle b_{1}^{\beta}[\varphi]:=\int_{\partial{D}_{1}}\frac{\partial{u}_{b}}{\partial\nu}\Big{|}_{+}\cdot\psi_{\beta}\quad\mbox{and}\quad b_{1}^{*\beta}[\varphi]:=\int_{\partial{D}_{1}^{*}}\frac{\partial u^{*}}{\partial\nu}\Big{|}_{+}\cdot\psi_{\beta},$

(3.2)

and denote

$$\rho_{d}(\varepsilon)=\begin{cases}\sqrt{\varepsilon},&\quad d=2,\\ |\log\varepsilon|^{-1},&\quad d=3.\end{cases}$$

We will prove that $b_{1}^{*\beta}[\varphi]$ are convergent to $b_{1}^{\beta}[\varphi]$, with rates $\rho_{d}(\varepsilon)$.

Here we need some symmetric assumptions on the domain and the boundary data. First, we assume that

	$\displaystyle({\rm S_{1}}):~{}$	$\displaystyle~{}D_{1}\cup D_{2}~{}\mbox{ and~{}}D~{}\mbox{are~{} symmetric~{} with~{} repect~{} to~{} each~{} axis~{}}x_{i},\mbox{~{}and~{}}$
		$\displaystyle~{}\mbox{ the superplane~{}}\{x_{d}=0\};$

For boundary data $\varphi$, we assume that in dimension two, $\varphi^{1}(x_{1},x_{2})$ is odd with respect to $x_{2}$ and $\varphi^{2}(x_{1},x_{2})$ is even with respect to $x_{2}$; and in dimension three, $\varphi^{i}(x)$ is odd for $i=1,2,3$, that is, for $x\in\partial D$,

	$\displaystyle({\rm S_{2}}):~{}~{}\varphi^{1}(x_{1},x_{2})=-\varphi^{1}(x_{1},-x_{2}),~{}\varphi^{2}(x_{1},x_{2})=\varphi^{2}(x_{1},-x_{2}),\quad\mbox{for}~{}d=2;$
	$\displaystyle({\rm S_{3}}):~{}~{}\varphi^{i}(x)=-\varphi^{i}(-x),\quad\,i=1,2,3,~{}\quad\mbox{for}~{}d=3.$

We shall use $O(1)$ to denote those quantities satisfying $|O(1)|\leq C$, for some constant independent of $\varepsilon$. We assume that for some $\delta_{0}>0$,

$$\delta_{0}\leq\mu,d\lambda+2\mu\leq\frac{1}{\delta_{0}}.$$

Then