Integral constraints in multiple scales problems with a slowly varying microstructure

Abstract

Asymptotic homogenisation is considered for problems with integral constraints imposed on a slowly varying microstructure; an insulator with an array of perfectly dielectric inclusions of slowly varying size serves as a paradigm. Although it is well-known how to handle each of these effects (integral constraints, slowly varying microstructure) independently within multiple scales analysis, additional care is needed when they are combined. Using the flux transport theorem, the multiple scales form of an integral constraint on a slowly varying domain is identified. The proposed form is applied to obtain a homogenised model for the electric potential in a dielectric composite, where the microstructure slowly varies and the integral constraint arises due to a statement of charge conservation. A comparison with multiple scales analysis of the problem with established approaches provides validation that the proposed form results in the correct homogenised model.

1. Introduction

Homogenisation via multiscale asymptotics is one coarse-graining method that can be used to derive the effective properties of composite media [14]. Typically used for periodic microstructure, example applications of the technique include modelling flow in porous media, wave propagation in poroelastic materials, filtration and decontamination processes [4, 7, 14, 15].

The result of the homogenisation process is the reduction of a problem posed on a complicated domain, or with rapidly varying coefficients, to two simpler problems: one ‘cell problem’ describing the microscale variation; and a second ‘homogenised model’ describing the macroscale variation of variables across the whole domain. The technique can be extended to problems with a slowly varying geometry, albeit at the cost of having a cell problem which varies with the macroscale [8]. A mapping depending on the slow spatial scale can be applied to transform a heterogeneous microstructure to an exactly periodic reference configuration [13, 19]. Standard homogenisation can be performed in this reference configuration before inverting the mapping to obtain the homogenised equations featuring spatially dependent coefficients which reflect microstructural variation. A similar approach can be used to treat microstructures with temporal and spatiotemporal variations. Examples of problems using a prescribed mapping include [3, 6, 21] and the mapping can be coupled to the macroscale variables [18]. When the domain is locally periodic and the unit cell has fixed size, transformation to a reference configuration is no longer required as the slow variable features as a parameter in the microscale problem [8]. Examples of this approach are found in [7, 11, 20].

When considering problems with a slowly varying domain, care must be taken in converting Neumann and Robin boundary conditions on microscopic inclusions into multiple scales form. Typically, a level set function is introduced to define the boundary of the inclusion, with the expansion of the normal to the boundary derived by writing the gradient of the level set function in multiple scales form [3, 11, 21]. An alternative approach is to treat the normal as defining the microstructure so that it does not need to be expanded [12, 16, 17]. However, in that case, the position of the boundary itself should be expanded, to ensure that the asymptotic expansions are self-consistent. Many authors do not do this expansion, resulting in homogenised equations which are incorrect. We address this problem in an Appendix, where we attempt to clear up some of the confusion regarding these contrasting approaches.

A second extension of the standard method allows for homogenisation of problems featuring integral constraints [5]. These constraints generally appear as conservation conditions, for example, of charge or momentum, with applications in modelling nematic crystals, radiation in porous media and bubbly liquids [2, 5, 22]. Unlike standard multiple scales, where the macroscale coordinate can be assumed to take a constant value within a given unit cell, it is crucial to account for the small variation in macroscale coordinate along the integration path, since this variation causes a change in flux which affects the parameters in the homogenised model.

In the present work, we aim to combine these two extensions, developing an understanding of how to write integral constraints on a slowly varying domain in multiple scales form. Although this seems like a routine task, we will see that in fact that the answer is not obvious a priori, and the ‘obvious’ approach is incorrect. We use as a paradigm the problem of the electric potential in an insulator interspersed with a periodic array of perfectly dielectric inclusions of slowly varying size. This problem has the advantage that the perfectly dielectric limit can also be taken after a standard homogenisation procedure so that we know what the homogenised model should be. Not all integral constraint problems can be recast in this way.

2. Paradigm problem

We consider the electric potential $\phi$ in a dielectric material, which satisfies Poisson’s equation

(2.1) \begin{equation} \nabla \cdot ( \varepsilon \nabla \phi) = -\rho, \end{equation}

where $\varepsilon$ is the permittivity and $\rho$ is the charge density (which we suppose is given). We consider a material composite comprising an insulator $\Omega _{\mathrm {ex}}$ of constant permittivity $\varepsilon _{\mathrm {ex}}$ with an array of inclusions $\Omega _{\mathrm {in}}$ of constant permittivity $\varepsilon _{\mathrm {in}}$ . At the boundary between the two regions

(2.2) \begin{align} \left [\textbf{n} \cdot ( \varepsilon \nabla \phi) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0,\\[-8pt]\nonumber \end{align}

(2.3) \begin{align} \left [ \phi \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

where $\textbf{n}$ is the (outward-facing) normal to the boundary of the inclusion, and $[\cdot]^{\mathrm {ex}}_{\mathrm {in}}$ represents the jump in the enclosed quantity across the interface. We suppose that the centres of the inclusions are arranged on a regular square lattice ¹ of side $\delta$ , and that the shape of each inclusion varies slowly with (macroscopic) position. ² In figure 1, we illustrate the set-up using circular inclusions of slowly varying radius, but our analysis will be valid for inclusions of any slowly varying shape. We aim to find effective (homogenised) equations in the limit $\delta \rightarrow 0$ .

In the limit $\varepsilon _{\mathrm {in}} \to \infty$ , the inclusions are perfectly dielectric, and the model becomes

(2.4) \begin{align} \nabla \cdot ( \varepsilon _{\mathrm {ex}} \nabla \phi) & = -\rho \qquad \mbox{in }\Omega _{\mathrm {ex}}, \end{align}

(2.5) \begin{align} \nabla \phi & = \textbf{0}\qquad \mbox{in }\Omega _{\mathrm {in}}, \end{align}

with boundary condition

(2.6) \begin{align} \left [ \phi \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}

The potential $\phi$ is constant on each inclusion but may take different values on different inclusions. To close the problem, we need to integrate (2.1) over each inclusion and use (2.2) to give the integral constraint

(2.7) \begin{equation} \int _{{\partial } \Omega _{\mathrm {in}}} \left.\varepsilon _{\mathrm {ex}}\, \textbf{n}\cdot {\nabla } \phi \right |_{\mathrm {ex}}\, {\textrm {d}} S = - \int _{\Omega _{\mathrm {in}}} \rho \, {\textrm {d}} {\textbf{x}}, \end{equation}

where $\Omega _{\mathrm {in}}$ is an individual inclusion and $|_{\mathrm {ex}}$ denotes the evaluation of the integrand on the exterior side of the inclusion boundary.

We will approach the limit $\varepsilon _{\mathrm {in}} \to \infty$ in two different ways. We will first homogenise (2.1)–(2.3) following [3], before taking the limit $\varepsilon _{\mathrm {in}} \to \infty$ in the homogenised model. We will then homogenise (2.4)–(2.7) directly, which will require us to determine how to cast (2.7) in multiple-scales form when the domain $\Omega _{\mathrm {in}}$ is a function of (slow) position.

Figure 1. An example of a 2D composite. Perfectly dielectric inclusions shown in grey lie on a periodic array within an insulator. The inclusions have a radius $a({\textbf{x}})$ which varies slowly across the domain. In this example $a({\textbf{x}}) = 0.3+0.1(x_1+x_2)$ so that the boundaries of the inclusions are $|(x_1,x_2)- \delta (n,m)| = \delta ( 0.3 + 0.1 x_1+0.1x_2)$ , for $n,m \in {\mathbb Z}$ .

2.1 Standard multiple scales

We introduce the fast scale ${\textbf{X}} = {\textbf{x}}/\delta$ , and suppose that $\phi = \phi ({\textbf{x}},{\textbf{X}})$ , treating $\textbf{x}$ and $\textbf{X}$ as independent, with derivatives transforming according to the chain rule

(2.8) \begin{equation} {\nabla } \rightarrow \nabla _{{\textbf{x}}} + \frac {1}{\delta }\nabla _{{\textbf{X}}}. \end{equation}

We remove the indeterminacy that this introduces by imposing that $\phi$ is $\textbf{1}$ -periodic in $\textbf{X}$ . To describe the inclusions we introduce the function $h({\textbf{x}},{\textbf{X}})$ , periodic in $\textbf{X}$ with period 1 in each component, such that the level set $h({\textbf{x}},{\textbf{X}})=0$ defines the boundary of the inclusions. For example, for the slowly varying circular inclusions in figure 1, we may take

(2.9) \begin{equation} h({\textbf{x}},{\textbf{X}}) = |{\textbf{X}}| - a({\textbf{x}}), \end{equation}

for ${\textbf{X}} \in [{-}1/2,1/2]^d$ where $d$ is the dimension, and then extend this function to be periodic in $\textbf{X}$ with period 1 in each component. ³

The normal to the inclusion can then be written in multiple scales form as

(2.10) \begin{equation} \textbf{n} =\frac {{\nabla } h}{|{\nabla } h|} = \frac {\nabla _{{\textbf{X}}} h + \delta \nabla _{{\textbf{x}}} h }{|\nabla _{{\textbf{X}}} h + \delta \nabla _{{\textbf{x}}} h |} = \textbf{n}_0 + \delta \textbf{n}_1 + O(\delta ^2), \end{equation}

where

(2.11) \begin{equation} \textbf{n}_0 = \frac {\nabla _{{\textbf{X}}} h}{|\nabla _{{\textbf{X}}} h|}, \qquad \textbf{n}_1 = \frac {\nabla _{{\textbf{x}}} h}{|\nabla _{{\textbf{X}}} h|} - \frac {(\nabla _{{\textbf{x}}} h \cdot \nabla _{{\textbf{X}}} h) \nabla _{{\textbf{X}}} h}{|\nabla _{{\textbf{X}}} h|^3} \end{equation}

(see Appendix C). Many authors fail to expand the normal in powers of $\delta$ in this way, often leading to incorrect final homogenised equations [4, 16, 17]. Since this question is tangential to our main question (that of how to handle integral constraints), we defer a more thorough discussion to Appendix A.

Within a given unit cell $D$ we denote the region occupied by the inclusion as $D_{\mathrm {in}}({\textbf{x}})$ and that occupied by the insulator as $D_{\mathrm {ex}}({\textbf{x}})$ , that is

\begin{equation*} D_{\mathrm {in}} = \{ {\textbf{X}} \in D\,:\,h({\textbf{x}},{\textbf{X}})\lt 0\}, \qquad D_{\mathrm {ex}} = \{ {\textbf{X}} \in D\,:\, h({\textbf{x}},{\textbf{X}})\gt 0\}. \end{equation*}

Substituting (2.8) and (2.10) into (2.1)–(2.3), expanding

(2.12) \begin{equation} \phi \sim \phi _0({\textbf{x}}, {\textbf{X}}) + \delta \phi _1({\textbf{x}}, {\textbf{X}}) + \cdots, \end{equation}

and equating coefficients of powers of $\delta$ , we find that at leading-order

(2.13) \begin{align} \nabla _{{\textbf{X}}} \cdot ( \varepsilon \nabla _{{\textbf{X}}} \phi _0) &= 0, \end{align}

(2.14) \begin{align} \left [\textbf{n}_0 \cdot ( \varepsilon \nabla _{{\textbf{X}}} \phi _0) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

(2.15) \begin{align} \left [ \phi _0 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with $\phi _0$ $\textbf{1}$ -periodic in $\textbf{X}$ . Thus, $\phi _0$ is constant on the fast scale, so that $\phi _0 = \phi _0({\textbf{x}})$ . At next order, we find

(2.16) \begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} \phi _1 \right) &= 0, \end{align}

(2.17) \begin{align} \left [\textbf{n}_0 \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \right) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

(2.18) \begin{align} \left [ \phi _1 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with $\phi _1$ $\textbf{1}$ -periodic in $\textbf{X}$ , where we have used the fast scale independence of the leading-order potential. The solution is

(2.19) \begin{equation} \phi _1 = \boldsymbol {\Psi }\cdot \nabla _{{\textbf{x}}} \phi _0 + \overline {\phi }_1, \end{equation}

where $\overline {\phi }_1$ is independent of $\textbf{X}$ and $\boldsymbol{\Psi }$ satisfies the cell problem

(2.20) \begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} \boldsymbol{\Psi } \right) &= 0, \end{align}

(2.21) \begin{align} \left [\textbf{n}_0 \cdot \left ( \varepsilon (\nabla _{{\textbf{X}}} \boldsymbol{\Psi } + I)\right) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

(2.22) \begin{align} \left [ \boldsymbol{\Psi } \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with $\boldsymbol{\Psi }$ $\textbf{1}$ -periodic in $\textbf{X}$ , where $I$ is the identity matrix, and uniqueness is achieved by imposing zero mean, for example. Finally, equating coefficients of $\delta ^0$ , we find

(2.23) \begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1)\right) + \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right) &= -\rho, \end{align}

(2.24) \begin{align} \left [\textbf{n}_0 \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1) \right) \right]^{\mathrm {ex}}_{\mathrm {in}} + \left [\textbf{n}_1 \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \right) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

(2.25) \begin{align} \left [ \phi _2 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}

Integrating (2.23) over the unit cell $D$ , applying the divergence theorem to terms involving the fast divergence, and using (2.24) we find

(2.26) \begin{align} & \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \left ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 \right) \cdot \textbf{n}_1 {\textrm {d}} S_{\textbf{X}} - \int _{\partial D_{\mathrm {in}}} \varepsilon _{\mathrm {in}} \left ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 \right) \cdot \textbf{n}_1 {\textrm {d}} S_{\textbf{X}} \nonumber \\& \quad +\int _{D_{\mathrm {ex}}} \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon _{\mathrm {ex}} ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right){\textrm {d}} {\textbf{X}} + \int _{D_{\mathrm {in}} } \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon _{\mathrm {in}} ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right) {\textrm {d}} {\textbf{X}} = -\rho _{\textrm{eff}}, \end{align}

where $\partial D_{\mathrm {in}}$ and $\partial D_{\mathrm {ex}}$ denote the interior and exterior of the inclusion boundary in the unit cell respectively, and the effective charge is given by

(2.27) \begin{equation} \rho _{\textrm{eff}} = \int _{D} \rho \, {\textrm {d}} {\textbf{X}}. \end{equation}

Taking the slow divergence outside the integral using the Reynolds transport theorem, we find

(2.28) \begin{align} & \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \left ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 \right) \cdot \textbf{n}_1 \,{\textrm {d}} S_{\textbf{X}} - \int _{\partial D_{\mathrm {in}}} \varepsilon _{\mathrm {in}} \left ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 \right) \cdot \textbf{n}_1 \,{\textrm {d}} S_{\textbf{X}} \nonumber \\[5pt]& \quad + \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \left ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 \right) \cdot \textbf{V} \cdot \textbf{n}_0 \, {\textrm {d}} S_{\textbf{X}} - \int _{\partial D_{\mathrm {in}}} \varepsilon _{\mathrm {in}} \left ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 \right)\cdot \textbf{V} \cdot \textbf{n}_0 \, {\textrm {d}} S_{\textbf{X}} \nonumber \\[5pt] & \quad +\nabla _{{\textbf{x}}} \cdot \int _{D_{\mathrm {ex}}} \left ( \varepsilon _{\mathrm {ex}} ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right)\, {\textrm {d}} {\textbf{X}} + \nabla _{{\textbf{x}}} \cdot \int _{D_{\mathrm {in}} } \left ( \varepsilon _{\mathrm {in}} ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right) \, {\textrm {d}} {\textbf{X}} = -\rho _{\textrm{eff}}, \end{align}

where the matrix $\textbf{V}$ is the ‘velocity’ of the boundary, i.e. the derivative of position on the boundary with respect to $\textbf{x}$ . Differentiating the equation $h=0$ with respect to $\textbf{x}$ gives

\begin{equation*}\textbf{V}\cdot \nabla _{{\textbf{X}}} h + \nabla _{{\textbf{x}}} h =\textbf{0},\end{equation*}

so that

(2.29) \begin{equation} \textbf{V} \cdot \textbf{n}_0 = -\frac{\nabla _{{\textbf{x}}} h}{|\nabla _{{\textbf{X}}} h|}, \qquad \textbf{V} \cdot \textbf{n}_0 + \textbf{n}_1 = - \frac {(\nabla _{{\textbf{x}}} h \cdot \nabla _{{\textbf{X}}} h) \nabla _{{\textbf{X}}} h}{|\nabla _{{\textbf{X}}} h|^3} = - \frac {(\nabla _{{\textbf{x}}} h \cdot \nabla _{{\textbf{X}}} h)}{|\nabla _{{\textbf{X}}} h|^2}\textbf{n}_0. \end{equation}

Thus, using (2.17), the surface integrals cancel in (2.28), leaving

(2.30) \begin{equation} \nabla _{{\textbf{x}}} \cdot \int _{D_{\mathrm {ex}}} \left ( \varepsilon _{\mathrm {ex}} ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right)\, {\textrm {d}} {\textbf{X}} + \nabla _{{\textbf{x}}} \cdot \int _{D_{\mathrm {in}} } \left ( \varepsilon _{\mathrm {in}} ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right)\, {\textrm {d}} {\textbf{X}} = -\rho _{\textrm{eff}}. \end{equation}

Substituting (2.19), gives, finally, the homogenised problem

(2.31) \begin{equation} \nabla _{{\textbf{x}}} \cdot \left ( {\boldsymbol \varepsilon }_{\textrm{eff}} \nabla _{{\textbf{x}}} \phi _0 \right) = -\rho _{\textrm{eff}}, \end{equation}

where the effective permittivity ${\boldsymbol \varepsilon }_{\textrm{eff}}$ is given by

(2.32) \begin{equation} {\boldsymbol \varepsilon }_{\textrm{eff}} = \int _{D} \varepsilon \left ( I + \nabla _{{\textbf{X}}} \boldsymbol{\Psi } \right) \, {\textrm {d}} {\textbf{X}}. \end{equation}

2.1.1 The limit $\varepsilon _{\mathrm {in}} \to \infty$

As $\varepsilon _{\mathrm {in}} \to \infty$ in the cell problem (2.20)–(2.22) we find

(2.33) \begin{align} \nabla _{{\textbf{X}}}^2 \boldsymbol{\Psi } &= \textbf{0} \quad \text {in} \quad D,\\[-5pt]\nonumber \end{align}

(2.34) \begin{align} \textbf{n}_0 \cdot \left ( \nabla _{{\textbf{X}}} \boldsymbol{\Psi } + I \right) &= \textbf{0} \quad \text {on} \quad \partial D_{\mathrm {in}},\\[-5pt]\nonumber \end{align}

(2.35) \begin{align} \left [ \boldsymbol{\Psi }\right]^{\mathrm {ex}}_{\mathrm {in}} &= \textbf{0}. \end{align}

Thus $\boldsymbol{\Psi } = -{\textbf{X}} +$ constant in $D_{\mathrm {in}}$ , where the constant must be chosen so that $\boldsymbol{\Psi }$ has zero mean. In the effective permittivity (2.32) this gives zero times infinity in the inclusion, so we must manipulate this expression into something more suitable before we take the limit. Switching to index notation, using (2.21), the divergence theorem, and (2.33), we find

(2.36) \begin{eqnarray} {{\varepsilon }_{\textrm{eff}}}_{ij} & =& \int _{D} \varepsilon \delta _{ij} \, {\textrm {d}} {\textbf{X}} + \int _{D} \varepsilon \frac {{\partial }}{\partial X_k} \left ( X_j \frac {{\partial } \Psi _i}{{\partial } X_k} \right)\, {\textrm {d}} {\textbf{X}}\nonumber \\ & = & \int _{D_{\mathrm {in}}} \varepsilon _{\mathrm {in}} \delta _{ij} \, {\textrm {d}} {\textbf{X}} + \int _{D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \delta _{ij} \,{\textrm {d}} {\textbf{X}} + \int _{\partial D_{\mathrm {in}}} \varepsilon _{\mathrm {in}} X_j \frac {{\partial }\Psi _i }{{\partial } X_k } n_k \,{\textrm {d}} S_{\textbf{X}}\nonumber \\&& \mbox{}- \int _{ \partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} X_j \frac {{\partial }\Psi _i }{{\partial } X_k } n_k \,{\textrm {d}} S_{\textbf{X}} + \int _{ \partial D} \varepsilon _{\mathrm {ex}} X_j\frac {{\partial }\Psi _i }{{\partial } X_k } n_k \,{\textrm {d}} S_{\textbf{X}} \nonumber \\ & = & \int _{D_{\mathrm {in}}} \varepsilon _{\mathrm {in}} \delta _{ij} \, {\textrm {d}} {\textbf{X}} + \int _{D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \delta _{ij}\, {\textrm {d}} {\textbf{X}} - \int _{\partial D_{\mathrm {in}}} \varepsilon _{\mathrm {in}} X_j n_i \,{\textrm {d}} S_{\textbf{X}} \nonumber \\&& \mbox{} + \int _{ \partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} X_j n_i \,{\textrm {d}} S_{\textbf{X}} + \int _{ \partial D} \varepsilon _{\mathrm {ex}} X_j \frac {{\partial }\Psi _i }{{\partial } X_k } n_k \,{\textrm {d}} S_{\textbf{X}}\nonumber \\ &=& \varepsilon _{\mathrm {ex}} \left ( \delta _{ij} + \int _{\partial D} X_j \frac {{\partial }\Psi _i }{{\partial } X_k } n_k\,{\textrm {d}} S_{\textbf{X}} \right), \end{eqnarray}

where ${\partial } D$ is the boundary of $D$ , i.e. the boundary of the unit cell on which periodic conditions were imposed. We can now safely take the limit $\varepsilon _{\mathrm {in}} \to \infty$ .

2.2 Multiple scales with integral constraints

We now apply the method of multiple scales directly to the problem (2.4)–(2.7), hoping to retrieve (2.31) with (2.36). Substituting (2.8) into (2.4)–(2.6), expanding as in (2.12), and equating coefficients of powers of $\delta$ we find that at leading-order

(2.37) \begin{align} \nabla _{{\textbf{X}}}\cdot ( \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \phi _0) & = 0 \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}

(2.38) \begin{align} \nabla _{{\textbf{X}}} \phi _0 &=\textbf{0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}

(2.39) \begin{align} \left [ \phi _0 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with $\phi _0$ $\textbf{1}$ -periodic in $\textbf{X}$ . Thus, as before, $\phi _0 = \phi _0({\textbf{x}})$ . At first-order, we find

(2.40) \begin{align} \nabla _{{\textbf{X}}}\cdot ( \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \phi _1) & = 0 \quad \text {in} \quad D_{\mathrm {ex}},\\[-12pt]\nonumber \end{align}

(2.41) \begin{align} \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}},\\[-12pt]\nonumber \end{align}

(2.42) \begin{align} \left [ \phi _1 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with $\phi _1$ $\textbf{1}$ -periodic in $\textbf{X}$ . As in section 2.1, the solution is $\phi _1 = \boldsymbol{\Psi }\cdot \nabla _{{\textbf{x}}} \phi _0 + \overline {\phi }_1$ where $\overline {\phi }_1$ is independent of $\textbf{X}$ and

(2.43) \begin{align} \nabla _{{\textbf{X}}} \cdot ( \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \boldsymbol{\Psi }) &= \textbf{0} \quad \text {in} \quad D_{\mathrm {ex}},\\[-12pt]\nonumber \end{align}

(2.44) \begin{align} \nabla _{{\textbf{X}}} \boldsymbol{\Psi } + I &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}},\\[-12pt]\nonumber \end{align}

(2.45) \begin{align} \left [ \boldsymbol{\Psi } \right]^{\mathrm {ex}}_{\mathrm {in}} &= \textbf{0}, \end{align}

with $\boldsymbol{\Psi }$ $\textbf{1}$ -periodic in $\textbf{X}$ , and we impose

(2.46) \begin{equation} \int _{D} \boldsymbol{\Psi } \, {\textrm {d}} {\textbf{X}} = \textbf{0}. \end{equation}

Equating coefficients of $\delta ^0$ we find

(2.47) \begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1)\right) + \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right) &= - \rho \quad \text {in} \quad D_{\mathrm {ex}},\\[-12pt]\nonumber \end{align}

(2.48) \begin{align} \nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1 &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}},\\[-12pt]\nonumber \end{align}

(2.49) \begin{align} \left [ \phi _2 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}

Integrating (2.47) over the exterior region and applying the divergence theorem to the first term gives

(2.50) \begin{equation} - \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1)\cdot \textbf{n}_0 \, {\textrm {d}} S_{\textbf{X}} + \int _{D_{\mathrm {ex}}} \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right)\, {\textrm {d}} {\textbf{X}} = -\int _{D_{\mathrm {ex}}} \rho \, {\textrm {d}} {\textbf{X}}, \end{equation}

where the integral over the exterior boundary of the unit cell ${\partial } D$ vanishes due to periodicity. To evaluate the surface integral in (2.50), we need to use the integral constraint (2.7).

2.2.1 Dealing with the integral

As discussed in [5], it seems natural to write (2.7) in multiple scales form as

\begin{equation*} \delta ^2\int _{{\partial } \Omega _{\mathrm {in}}} \varepsilon _{\mathrm {ex}}\, \textbf{n}\cdot \left ( \nabla _{{\textbf{x}}} \phi + \frac {1}{\delta }\nabla _{{\textbf{X}}} \phi \right) \, {\textrm {d}} S_{\textbf{X}} = - \delta ^3\int _{\Omega _{\mathrm {in}}} \rho \, {\textrm {d}} {\textbf{X}}, \end{equation*}

but this is incorrect, as it neglects the small variation in the slow variable $\textbf{x}$ around the boundary of the inclusion, which turns out to be crucial. Writing $\textbf{Q} = {\nabla } \phi$ , the approach taken in [5] was to recognise that on the interface ${\textbf{x}} = \hat {{\textbf{x}}} + \delta {\textbf{X}}$ , where $\hat {{\textbf{x}}}$ is the position of the centre of the unit cell and ${\textbf{X}} \in [{-}1/2,1/2]^d$ , expanding

(2.51) \begin{equation} \textbf{Q}({\textbf{x}},{\textbf{X}}) =\textbf{Q}(\hat {{\textbf{x}}}+\delta {\textbf{X}},{\textbf{X}}) = \textbf{Q}(\hat {{\textbf{x}}},{\textbf{X}})+\delta {\textbf{X}} \cdot \nabla _{{\textbf{x}}} \textbf{Q}(\hat {{\textbf{x}}},{\textbf{X}})+\cdots \end{equation}

in the integrand of (2.7). But how should we proceed when the domain and the normal, as well as the integrand, depend on the slow variable $\textbf{x}$ ?

A plausible but incorrect approach

A plausible way to write the integral constraint on a slowly varying domain in multiple scales form seems to be to combine the normal expansion (2.10) with the expansion of the integrand given in (2.51) by writing

(2.52) \begin{equation} \int _{\partial \Omega _{\mathrm {in}}} \textbf{Q} \cdot \textbf{n}\, {\textrm {d}} S \to \delta ^2 \int _{\partial \Omega _{\mathrm {in}}} \left (\textbf{Q}_0 + \delta ( \textbf{Q}_1 + {\textbf{X}} \cdot \nabla _{{\textbf{x}}} \textbf{Q}_0) + \ldots \right) \cdot (\textbf{n}_0 + \delta \textbf{n}_1 + \ldots) \, {\textrm {d}} S_{\textbf{X}}. \end{equation}

This is the approach taken in [2, 9] (on much more complicated problems). ⁴ However, it is not correct, as we now highlight. Writing (2.47) in terms of the flux $\textbf{Q}$ , and integrating over the exterior region, we find

(2.53) \begin{equation} \int _{D_{\mathrm {ex}}} \nabla _{{\textbf{X}}} \cdot \textbf{Q}_1 \, {\textrm {d}} {\textbf{X}} + \int _{D_{\mathrm {ex}}} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 \, {\textrm {d}} {\textbf{X}} = -\int _{D_{\mathrm {ex}}} \rho \, {\textrm {d}} {\textbf{X}}. \end{equation}

Applying the divergence theorem to the first term and using the integral constraint, we find

(2.54) \begin{equation} \int _{\partial D_{\mathrm {ex}}} \textbf{Q}_0 \cdot \textbf{n}_1 \,{\textrm {d}} S_{\textbf{X}} + \int _{\partial D_{\mathrm {ex}}} {\textbf{X}} \cdot \nabla _{{\textbf{x}}} \textbf{Q}_0 \cdot \textbf{n}_0\, {\textrm {d}} S_{\textbf{X}}+ \int _{D_{\mathrm {ex}}} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0\, {\textrm {d}} {\textbf{X}} = -\int _{D_{\mathrm {ex}}} \rho \, {\textrm {d}} {\textbf{X}} - \int _{D_{\mathrm {in}}} \rho \, {\textrm {d}} {\textbf{X}}. \end{equation}

Applying the flux transport theorem to the second term and Reynolds transport theorem to the final term of (2.54), we obtain

(2.55) \begin{equation} \int _{\partial D_{\mathrm {ex}}} \textbf{Q}_0 \cdot \textbf{n}_1 \, {\textrm {d}} S_{\textbf{X}} +\nabla _{{\textbf{x}}} \cdot (\varepsilon _{\textrm {eff}} \nabla _{{\textbf{x}}} \phi _0) = - \rho _{\textrm {eff}}, \end{equation}

after simplification, where we have defined the effective charge as in (2.27). In comparison to (2.31), we have an additional boundary integral, which should not be there.

The correct approach

Figure 2. Schematic of the open surface $\partial \Omega$ changing with slow coordinate.

In identifying the multiple scales form of integral constraints on a periodic domain, the integrand was expanded about a fixed slow position as in (2.51). When the microstructure slowly varies, we must also expand the boundary position about a fixed slow position. Thus, given a slowly varying surface ${\partial }\Omega$ , we look to expand

(2.56) \begin{equation} \int _{\partial \Omega } \textbf{Q} \cdot \textbf{n} \, {\textrm {d}} S = \delta ^2 \int _{\partial \Omega (\hat {{\textbf{x}}} + \delta {\textbf{X}})} \textbf{Q} (\hat {{\textbf{x}}} + \delta {\textbf{X}}, {\textbf{X}})\cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}} \end{equation}

For generality (and for ease of illustration in figure 2), we allow the surface $\partial \Omega$ to be open, though it will usually be the boundary of an inclusion. Expanding the integrand as in [5], we find

(2.57) \begin{equation} \int _{\partial \Omega } \textbf{Q} \cdot \textbf{n} \, {\textrm {d}} S = \delta ^2 \int _{\partial \Omega (\hat {{\textbf{x}}} + \delta {\textbf{X}})} \big {(} \textbf{Q} (\hat {{\textbf{x}}}, {\textbf{X}}) + \delta {\textbf{X}} \cdot \nabla _{{\textbf{x}}} \textbf{Q} (\hat {{\textbf{x}}}, {\textbf{X}})) \big) \cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}} + \cdots. \end{equation}

To project the boundary onto that at $\hat {{\textbf{x}}}$ , we take a similar approach to heuristic derivations of the flux transport theorem, see, for example, [10]. We apply the divergence theorem to the volume $\Omega _{\delta }$ swept out by the surface as we move from $\hat {{\textbf{x}}}$ to $\hat {{\textbf{x}}} + \delta {\textbf{X}}$ (illustrated schematically in figure 2), writing

(2.58) \begin{align} \int _{\partial \Omega (\hat {{\textbf{x}}} + \delta {\textbf{X}})} \textbf{Q} (\hat {{\textbf{x}}}, {\textbf{X}}) \cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}} & = \int _{\Omega _\delta (\hat {{\textbf{x}}})} \nabla _{{\textbf{X}}} \cdot \textbf{Q}(\hat {{\textbf{x}}}, {\textbf{X}}) {\textrm {d}} {\textbf{X}} + \int _{\partial \Omega (\hat {{\textbf{x}}})} \textbf{Q}(\hat {{\textbf{x}}}, {\textbf{X}}) \cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}} \nonumber \\& \quad - \int _{\partial \Omega _\delta (\hat {{\textbf{x}}})} \textbf{Q}(\hat {{\textbf{x}}}, {\textbf{X}})\cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}}, \end{align}

where $\partial \Omega _\delta (\hat {{\textbf{x}}})$ is the surface of $\Omega _\delta$ excluding the ends $\partial \Omega (\hat {{\textbf{x}}})$ and $\partial \Omega (\hat {{\textbf{x}}} + \delta {\textbf{X}})$ . Note that in the integral over ${\partial } \Omega (\hat {{\textbf{x}}})$ , the vector $\textbf{n}$ is the inward normal to $\Omega _\delta (\hat {{\textbf{x}}})$ (see figure 2), which is the reason this term appears with a plus rather than a minus. In the limit $\delta \to 0$ , we can write the volume and surface elements as ${\textrm {d}} {\textbf{X}} = \delta {\textbf{X}} \cdot \textbf{V} \cdot \textbf{n} \, {\textrm {d}} S_X$ and $\textbf{n} {\textrm {d}} S_{\textbf{X}} = -\delta {\textbf{X}} \cdot \textbf{V} \times {\textrm {d}} {\mathbf {s}}$ , respectively, where $\textbf{V}$ is the ‘velocity’ of the boundary as before, and ${\textrm {d}} {\mathbf {s}}$ is the line element of $\Gamma (\hat {{\textbf{x}}})$ . Substituting the form for the volume element of $\Omega _\delta$ and area element of $\partial \Omega _\delta$ into (2.58), we have

(2.59) \begin{align} \int _{\partial \Omega (\hat {{\textbf{x}}} + \delta {\textbf{X}})} \textbf{Q} (\hat {{\textbf{x}}}, {\textbf{X}}) \cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}} & = \int _{\partial \Omega (\hat {{\textbf{x}}})} \delta (\nabla _{{\textbf{X}}} \cdot \textbf{Q}) {\textbf{X}} \cdot \textbf{V} \cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}}\nonumber \\& \quad + \int _{\partial \Omega (\hat {{\textbf{x}}})} \textbf{Q} \cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}} + \int _{\Gamma (\hat {{\textbf{x}}}) } \delta \textbf{Q}\cdot ({\textbf{X}} \cdot \textbf{V}) \times {\textrm {d}} {\mathbf {s}}, \end{align}

Combining (2.59) with (2.57), we obtain the multiple scales form for integral constraints on a slowly varying domain

(2.60) \begin{align} \int _{\partial \Omega } \textbf{Q} \cdot \textbf{n}\, {\textrm {d}} S \rightarrow \delta ^2 \int _{\partial \Omega } \left (\textbf{Q} + \delta {\textbf{X}} \cdot \nabla _{{\textbf{x}}} \textbf{Q} + \delta (\nabla _{{\textbf{X}}} \cdot \textbf{Q}) {\textbf{X}} \cdot \textbf{V} \right)\cdot \textbf{n} \, {\textrm {d}} S_{\textbf{X}} + \delta ^2 \int _{ \Gamma } \delta \textbf{Q}\cdot ({\textbf{X}} \cdot \textbf{V}) \times {\textrm {d}} {\mathbf {s}}. \end{align}

In (2.60), $\textbf{n}$ is the normal to the boundary at fixed $\hat {{\textbf{x}}}$ ; in the example of inclusions with a slowly varying radius, this is given by $\textbf{n}_0$ . Note that, unlike in section 2.1 when approximating (2.2), and perhaps counter-intuitively, we do not need to expand the normal to introduce $\textbf{n}_1$ , or to apply the operator ${\textbf{X}} \cdot \nabla _{{\textbf{x}}}$ to $\textbf{n}_0$ : the perturbation to the normal is already accounted for by the term involving $\textbf{V}$ . An expansion of the normal will only appear in (2.60) when the function $h$ defining the boundary through its level sets is itself a function of $\delta$ (as in Appendix B, for example).

Returning to our problem of dielectric inclusions, the integral constraint (2.7), in multiple scales form, is

(2.61) \begin{equation} \int _{{\partial } D_{\mathrm {ex}}} \left (\textbf{Q} + \delta {\textbf{X}} \cdot \nabla _{{\textbf{x}}} \textbf{Q} + \delta (\nabla _{{\textbf{X}}} \cdot \textbf{Q}) {\textbf{X}} \cdot \textbf{V}\right)\cdot \textbf{n}_0\, {\textrm {d}} S_{\textbf{X}} = -\delta \int _{D_{\mathrm {in}}} \rho \, {\textrm {d}} {\textbf{X}}, \end{equation}

where

\begin{equation*} \textbf{Q} = \frac {1}{\delta } \nabla _{{\textbf{X}}} \phi + \nabla _{{\textbf{x}}} \phi \end{equation*}

and the surface integral is over the exterior surface of the inclusion. Using the expansion (2.12) in (2.61), we find at leading-order that

(2.62) \begin{equation} \int _{{\partial } D_{\mathrm {ex}}}\nabla _{{\textbf{X}}} \phi _0 \cdot \textbf{n}_0\, {\textrm {d}} S_{\textbf{X}} =0, \end{equation}

consistent with $\phi _0 = \phi _0({\textbf{x}})$ . At first-order, we find

(2.63) \begin{equation} \int _{{\partial } D_{\mathrm {ex}}} \textbf{Q}_0 \cdot \textbf{n}_0\, {\textrm {d}} S_{\textbf{X}} =0,\qquad \qquad \textbf{Q}_0 = \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0, \end{equation}

which is consistent with (2.40). Finally, equating coefficients of $\delta$ , and noting that $\nabla _{{\textbf{X}}} \cdot \textbf{Q}_0 = 0$ , we obtain

(2.64) \begin{equation} \int _{{\partial } D_{\mathrm {ex}}} \left (\textbf{Q}_1 + {\textbf{X}} \cdot \nabla _{{\textbf{x}}} \textbf{Q}_0 \right)\cdot \textbf{n}_0\, {\textrm {d}} S_{\textbf{X}} = - \int _{D_{\mathrm {in}}} \rho \, {\textrm {d}} {\textbf{X}},\qquad \qquad \textbf{Q}_1 = \nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1. \end{equation}

Substituting into (2.50) gives

(2.65) \begin{equation} \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} {\textbf{X}} \cdot \nabla _{{\textbf{x}}} \textbf{Q}_0 \cdot \textbf{n}_0 \,{\textrm {d}} S_{\textbf{X}} + \int _{D_{\mathrm {ex}}} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 \, {\textrm {d}} {\textbf{X}} = - \rho _{\textrm{eff}}, \end{equation}

where

(2.66) \begin{equation} \rho _{\textrm{eff}} = \int _{D} \rho \, {\textrm {d}} {\textbf{X}} \end{equation}

is the effective charge as before. Using the transport theorem to take the slow derivatives outside the integral gives

\begin{eqnarray*} \int _{{\partial } D_{\mathrm {in}}} ({\textbf{X}}\cdot \nabla _{{\textbf{x}}}) \textbf{Q}_0\cdot \textbf{n}_0\, {\textrm {d}} S_{\textbf{X}} & =& \int _{{\partial } D_{\mathrm {in}}} X_i \frac {{\partial } Q_{0,j}}{{\partial } x_i}n_{0,j}\, {\textrm {d}} S_{\textbf{X}}\\ & = & \frac {{\partial }}{{\partial } x_i} \int _{{\partial } D_{\mathrm {in}}} X_i Q_{0,j}n_{0,j}\, {\textrm {d}} S_{\textbf{X}} - \int _{{\partial } D_{\mathrm {in}}} Q_{0,i}V_{ij}n_{0,j}\, {\textrm {d}} S_{\textbf{X}}\\ & = & \nabla _{{\textbf{x}}} \cdot \int _{{\partial } D_{\mathrm {in}}} {\textbf{X}} (\textbf{Q}_{0} \cdot \textbf{n}_{0})\, {\textrm {d}} S_{\textbf{X}} - \int _{{\partial } D_{\mathrm {in}}} \textbf{Q}_0 \cdot \textbf{V} \cdot \textbf{n}_{0}\, {\textrm {d}} S_{\textbf{X}}, \end{eqnarray*}

while

\begin{eqnarray*} \int _{D_{\mathrm {ex}}} \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon _{\mathrm {ex}} \textbf{Q}_0\right) \, {\textrm {d}} {\textbf{X}} & = & \nabla _{{\textbf{x}}} \cdot \int _{D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \textbf{Q}_0 \, {\textrm {d}} {\textbf{X}} + \int _{{\partial } D_{\mathrm {in}} } \textbf{Q}_0 \cdot \textbf{V} \cdot \textbf{n}_{0}\, {\textrm {d}} S_{\textbf{X}} \end{eqnarray*}

(since $\textbf{n}_0$ is the outward normal to $D_{\mathrm {in}}$ ). Thus the two surface integrals cancel. Simplifying now as we did to obtain (2.36), we find that (2.65) becomes

(2.67) \begin{equation} \nabla _{{\textbf{x}}} \cdotp \left ( {\boldsymbol \varepsilon }_{\textrm{eff}} \nabla _{{\textbf{x}}} \phi _0 \right) = -\rho _{\textrm{eff}}, \end{equation}

where the effective permittivity tensor is

(2.68) \begin{equation} {\varepsilon _{\textrm{eff}}}_{ij} = \varepsilon _{\mathrm {ex}} \left ( \delta _{ij} + \int _{\partial D} X_i \frac {\partial \Psi _j}{\partial X_k} n_k \, {\textrm {d}} S_{\textbf{X}} \right) \end{equation}

in agreement with (2.36).

3. Paradigm problem: another limit

In section 2, we illustrated how to treat the multiple-scales problem with integral constraints considered in [5] when the domain slowly varies. In that example $\nabla _{{\textbf{X}}}\cdot \textbf{Q}_0=0$ , so that the extra term we took so much care to get correct in (2.61) actually vanishes. In order to provide further validation of (2.61), we construct now a problem in which this term is non-zero by considering the case in which

(3.1) \begin{equation} \int _D \rho \, {\textrm {d}} {\textbf{X}} = 0. \end{equation}

Then $\rho _{\mathrm {eff}}$ in (2.67) is zero, and (given homogeneous boundary conditions on $\phi _0$ ) the solution is simply $\phi _0\equiv 0$ , indicating that the scaling balance between $\rho$ and $\phi$ is not correct. Rescaling $\phi \rightarrow \delta \phi$ leads to

(3.2) \begin{equation} \nabla \cdot ( \varepsilon \nabla \phi) = -\frac {\rho }{\delta }, \end{equation}

(with $\phi, \rho$ of $O(1)$ ) which we can think of as a problem with a large charge density for which the average charge over a unit cell is zero. Boundary conditions (2.2) and (2.3) are unchanged.

In the limit of perfectly dielectric inclusions, where $\varepsilon _{\mathrm {in}} \to \infty$ , we have

(3.3) \begin{align} \nabla \cdot ( \varepsilon _{\mathrm {ex}} \nabla \phi) & = -\frac {\rho }{\delta }\qquad \mbox{in }\Omega _{\mathrm {ex}}, \end{align}

(3.4) \begin{align} \nabla \phi & = 0\qquad \mbox{in }\Omega _{\mathrm {in}}, \end{align}

with continuity (2.6) at the inclusion boundary. The rescaled integral constraint becomes

(3.5) \begin{equation} \int _{{\partial } \Omega _{\mathrm {in}}} \left.\varepsilon _{\mathrm {ex}}\, \textbf{n}\cdot {\nabla } \phi \right |_{\mathrm {ex}}\, {\textrm {d}} S = - \frac {1}{\delta } \int _{\Omega _{\mathrm {in}}} \rho \, {\textrm {d}} {\textbf{x}}. \end{equation}

We perform a similar analysis to section 2: first we take the limit $\varepsilon _{\mathrm {in}} \to \infty$ in the standard multiple scales problem before comparing it with the problem formulated with an integral condition.

3.1 Standard multiple scales

We substitute the multiple scales expansion (2.12) into (3.2) and compare coefficients at each order of $\delta$ . At leading-order, we find that $\phi _0$ is independent of the fast scale. At next order, we have

(3.6) \begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} \phi _1 \right) &= -\rho \quad \text {in} \quad D, \end{align}

(3.7) \begin{align} \left [\textbf{n}_0 \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \right) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

(3.8) \begin{align} \left [ \phi _1 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with $\phi _1$ $\textbf{1}$ -periodic in $\textbf{X}$ . We write

(3.9) \begin{equation} \phi _1 = \boldsymbol{\Psi }\cdot \nabla _{{\textbf{x}}} \phi _0 + \xi + \overline {\phi }_1, \end{equation}

with $\overline {\phi }_1 = \overline {\phi }_1(\textbf{x})$ , where $\boldsymbol \Psi$ takes care of the inhomogeneity due to $\nabla _{{\textbf{x}}} \phi _0$ and $\xi$ takes care of the inhomogeneity due to $\rho$ . We thus obtain two microscale problems: $\boldsymbol{\Psi }$ satisfies (2.20)–(2.22) as before, while $\xi$ satisfies

(3.10) \begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} {\xi } \right) &= -\rho \quad \text {in} \quad D, \end{align}

(3.11) \begin{align} \left [\textbf{n}_0 \cdot \varepsilon \nabla _{{\textbf{X}}} {\xi } \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

(3.12) \begin{align} \left [ {\xi } \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with

(3.13) \begin{equation} \int _{D} \xi \, d{\textbf{X}} = 0. \end{equation}

Note that a solution for $\xi$ exists because of (3.1). Equating coefficients of $\delta ^0$ , we find

(3.14) \begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1)\right) + \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right) &= 0 \quad \text {in} \quad D, \end{align}

(3.15) \begin{align} \left [\textbf{n}_0 \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1) \right) \right]^{\mathrm {ex}}_{\mathrm {in}} + \left [\textbf{n}_1 \cdot \left ( \varepsilon ( \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \right) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

(3.16) \begin{align} \left [ \phi _2 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}

We integrate (3.14) over the unit cell $D$ , applying the divergence theorem to terms involving the fast divergence, the Reynolds transport theorem and use (3.15). Following similar analysis to section 2.1, we obtain the homogenised problem

(3.17) \begin{equation} \nabla _{{\textbf{x}}} \cdot \left ( {\boldsymbol \varepsilon }_{\textrm{eff}} \nabla _{{\textbf{x}}} \phi _0 \right) = -\rho _{\textrm {eff}}, \end{equation}

with effective permittivity ${\boldsymbol \varepsilon }_{\textrm{eff}}$ given by (2.32) and effective charge density,

(3.18) \begin{equation} \rho _{\textrm {eff}} = \nabla _{{\textbf{x}}} \cdotp \int _{D} \varepsilon \nabla _{{\textbf{X}}} \xi {\rm d}\bf X. \end{equation}

3.1.1 Taking the limit $\varepsilon _{{\mathrm {in}}} \to \infty$

In the limit of perfectly dielectric inclusions, the effective permittivity takes the form (2.36). In the limit $\varepsilon _{{\mathrm {in}}} \to \infty$ , the cell problem for $\xi$ becomes

(3.19) \begin{align} \nabla _{{\textbf{X}}} \cdot (\varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \xi) &= - \rho \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}

(3.20) \begin{align} \nabla _{{\textbf{X}}} \xi &= \boldsymbol {0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}

(3.21) \begin{align} \textbf{n}_0 \cdotp \nabla _{{\textbf{X}}} \xi &= 0 \quad \text {on} \quad \partial D_{\mathrm {in}}, \end{align}

(3.22) \begin{align} \left [ \xi \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with $\xi$ $\boldsymbol {1}$ -periodic in $\boldsymbol {X}$ . We switch to index notation to establish the form of the effective charge $\rho _{\textrm {eff}}$ in the limit $\varepsilon _{{\mathrm {in}}} \to \infty$ ,

(3.23) \begin{equation} \begin{split} \rho _{\textrm {eff}} &= \frac {{\partial } }{{\partial } x_i} \int _{D} \varepsilon \frac {{\partial } \xi }{{\partial } X_i} {\textrm {d}} \textbf{X} \\ &= \frac {{\partial } }{{\partial } x_i} \int _{D} \varepsilon \frac {{\partial } }{{\partial } X_j}\left ( X_i \frac {{\partial } \xi }{{\partial } X_j} \right) {\textrm {d}} \textbf{X} - \frac {{\partial }}{{\partial } x_i} \int _{D} \varepsilon X_i \frac {{\partial }^2 \xi }{{\partial } X_j {\partial } X_j} {\textrm {d}} \textbf{X}\\ &= \frac {{\partial } }{{\partial } x_i} \int _{\partial D} \varepsilon _{\mathrm {ex}} X_i \frac {{\partial } \xi }{{\partial } X_j} n_j {\textrm {d}} S - \frac {{\partial }}{{\partial } x_i} \int _{D} \varepsilon X_i \frac {{\partial }^2 \xi }{{\partial } X_j {\partial } X_j} {\textrm {d}} \textbf{X}\\ &= \frac {{\partial } }{{\partial } x_i} \int _{\partial D} \varepsilon _{\mathrm {ex}} X_i \frac {{\partial } \xi }{{\partial } X_j} n_j {\textrm {d}} S + \frac {{\partial }}{{\partial } x_i} \int _{D_{\mathrm {ex}}} X_i \rho {\textrm {d}} \textbf{X}, \end{split} \end{equation}

where we have used (3.11) in going from the second to third line and (3.19)–(3.20) in the third to fourth line.

3.2 Multiple scales with integral constraints

In this section we treat (3.3)–(3.5) directly, writing (3.5) in multiple scales form using (2.60). Substituting (2.12) into (3.3)–(3.5) written in multiple scales form, we find that $\phi _0 = \phi _0(\textbf{x})$ at leading-order. At first-order, we find

(3.24) \begin{align} \nabla _{{\textbf{X}}}\cdot ( \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \phi _1) & = -\rho \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}

(3.25) \begin{align} \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}

(3.26) \begin{align} \left [ \phi _1 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

(3.27) \begin{align} \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} & = - \int _{D_{\mathrm {in}}} \rho {\textrm {d}} \textbf{X}, \end{align}

with $\phi _1$ $\textbf{1}$ -periodic in $\textbf{X}$ . As in section 3.1, we write $\phi _1 = \boldsymbol{\Psi }\cdot \nabla _{{\textbf{x}}} \phi _0 + \xi + \overline {\phi }_1$ where $\overline {\phi }_1 = \overline {\phi }_1(\textbf{x})$ and $\boldsymbol{\Psi }$ is the solution to (2.43)–(2.45). The second cell function $\xi$ satisfies

(3.28) \begin{align} \nabla _{{\textbf{X}}} \cdot ( \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \xi) &= -\rho \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}

(3.29) \begin{align} \nabla _{{\textbf{X}}} \xi &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}

(3.30) \begin{align} \left [ \xi \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with

(3.31) \begin{equation} \int _{D} \xi {\textrm {d}}\textbf{X } = 0. \end{equation}

Equating coefficients at next order, we have

(3.32) \begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1)\right) + \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right) &= 0 \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}

(3.33) \begin{align} \nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1 &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}

(3.34) \begin{align} \left [ \phi _2 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}

with

(3.35) \begin{align} \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1) \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} + \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} {\textbf{X}}\cdotp \nabla _{{\textbf{x}}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} \nonumber \\ + \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \cdotp (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \textbf{X } \cdotp \textbf{V} \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} = 0. \end{align}

Integrating (3.32) over the exterior region, applying the divergence theorem to the first term and substituting (3.35), we have

\begin{equation*} \begin{split} \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} {\textbf{X}}\cdotp \nabla _{{\textbf{x}}} (\nabla _{{\textbf{X}}} \phi _1 &+ \nabla _{{\textbf{x}}} \phi _0) \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} + \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \cdotp (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \textbf{X } \cdotp \textbf{V}\cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} \\ &+ \int _{D_{\mathrm {ex}}} \nabla _{{\textbf{x}}} \cdot \left ( \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right)\, {\textrm {d}} {\textbf{X}} = 0. \end{split} \end{equation*}

Applying the transport theorem to the first and final integrals gives

\begin{equation*} \begin{split} \nabla _{{\textbf{x}}} \cdotp \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} {\textbf{X}} (\nabla _{{\textbf{X}}} \phi _1 &+ \nabla _{{\textbf{x}}} \phi _0) \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} - \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \cdotp ( {\textbf{X}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)) {\cdotp \bf V} \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} \\ + \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \cdotp (\nabla _{{\textbf{X}}} \phi _1 &+ \nabla _{{\textbf{x}}} \phi _0) \textbf{X } \cdotp \textbf{V}\cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} + \nabla _{{\textbf{x}}} \cdot \int _{D_{\mathrm {ex}}} \left ( \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right)\, {\textrm {d}} {\textbf{X}} \\ &+ \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) {\cdotp \bf V} \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} = 0. \end{split} \end{equation*}

Expanding the divergence in the second integral, we find some of the boundary terms cancel, leaving

(3.36) \begin{equation} \nabla _{{\textbf{x}}} \cdotp \int _{\partial D_{\mathrm {ex}}} \varepsilon _{\mathrm {ex}} {\textbf{X}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0) \cdotp \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} + \nabla _{{\textbf{x}}} \cdot \int _{D_{\mathrm {ex}}} \left ( \varepsilon _{\mathrm {ex}} (\nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0)\right)\, {\textrm {d}} {\textbf{X}} = 0. \end{equation}

We substitute (3.9), using the divergence theorem to take the first integral into the exterior region and use (3.28)–(3.30) to obtain

(3.37) \begin{equation} \nabla _{{\textbf{x}}} \cdot (\varepsilon _{\textrm {eff}} \nabla _{{\textbf{x}}} \phi _0) = - \rho _{\textrm {eff}} \end{equation}

where

(3.38) \begin{equation} \varepsilon _{\mathrm {eff} ij} = \varepsilon _{{\mathrm {ex}}} \left ( \delta _{ij} + \int _{\partial D} X_i \frac {{\partial } \Psi _j}{{\partial } X_k } \textbf{n}_{0k} {\textrm {d}} S_{\textbf{X}} \right). \end{equation}

The effective charge is given by

(3.39) \begin{equation} \rho _{\textrm {eff}} = \nabla _{{\textbf{x}}}\cdotp \left ( \int _{D_{\mathrm {ex}}}\rho \textbf{X} {\textrm {d}} \textbf{X} + \int _{\partial D} \varepsilon _{\mathrm {ex}} \textbf{X} \nabla _{{\textbf{X}}} \xi \cdot \textbf{n}_0 {\textrm {d}} S_{\textbf{X}} \right). \end{equation}

Thus, we have recovered (3.17) in the limit of perfectly dielectric inclusions, providing further evidence that the divergence term present in (2.61) is correct.

4. Discussion

We have outlined how to combine integral constraints with a slowly varying microstructure, when performing asymptotic homogenisation using the method of multiple scales. Our main result is equation (2.60), which shows how to write an integral constraint in multiple scales form when the (fast) domain of the integral is a function of the slow scale. Essentially, the rest of the manuscript is a justification of this equation, showing that it leads to the correct homogenised model for an example in which that model can be identified using a more standard approach. Some problems involving integral constraints, especially those in which different physics holds in the inclusions, do not arise as a limit of a more standard problem, and such an approach is not available. These problems can be handled using equation (2.60).