1. Introduction
Homogenisation via multiscale asymptotics is one coarse-graining method that can be used to derive the effective properties of composite media [Reference Hornung14]. 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 [Reference Burridge and Keller4, Reference Dalwadi, Griffiths and Bruna7, Reference Hornung14, Reference Luckins, Breward, Griffiths and Wilmott15].
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 [Reference Dalwadi8]. A mapping depending on the slow spatial scale can be applied to transform a heterogeneous microstructure to an exactly periodic reference configuration [Reference Hornung13, Reference Peter19]. 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 [Reference Bruna and Chapman3, Reference Collis, Brown, Hubbard and O’Dea6, Reference Richardson and Chapman21] and the mapping can be coupled to the macroscale variables [Reference Peter and Böhm18]. 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 [Reference Dalwadi8]. Examples of this approach are found in [Reference Dalwadi, Griffiths and Bruna7, Reference Fatima, Arab, Zemskov and Muntean11, Reference Ray, Elbinger and Knabner20].
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 [Reference Bruna and Chapman3, Reference Fatima, Arab, Zemskov and Muntean11, Reference Richardson and Chapman21]. An alternative approach is to treat the normal as defining the microstructure so that it does not need to be expanded [Reference Holmes12, Reference Penta, Ambrosi and Shipley16, Reference Penta, Ambrosi and Quarteroni17]. 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 [Reference Chapman and McBurnie5]. 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 [Reference Bennett, D’Alessandro and Daly2, Reference Chapman and McBurnie5, Reference Rooney, Please and Howison22]. 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
\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
\begin{align} \left [\textbf{n} \cdot ( \varepsilon \nabla \phi) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0,\\[-8pt]\nonumber \end{align}
\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 latticeFootnote
1
of side
$\delta$
, and that the shape of each inclusion varies slowly with (macroscopic) position.Footnote
2
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
\begin{align} \nabla \cdot ( \varepsilon _{\mathrm {ex}} \nabla \phi) & = -\rho \qquad \mbox{in }\Omega _{\mathrm {ex}}, \end{align}
\begin{align} \nabla \phi & = \textbf{0}\qquad \mbox{in }\Omega _{\mathrm {in}}, \end{align}
with boundary condition
\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
\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 [Reference Bruna and Chapman3], 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.
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
\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
\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.Footnote
3
The normal to the inclusion can then be written in multiple scales form as
\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
\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 [Reference Burridge and Keller4, Reference Penta, Ambrosi and Shipley16, Reference Penta, Ambrosi and Quarteroni17]. 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
\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
\begin{align} \nabla _{{\textbf{X}}} \cdot ( \varepsilon \nabla _{{\textbf{X}}} \phi _0) &= 0, \end{align}
\begin{align} \left [\textbf{n}_0 \cdot ( \varepsilon \nabla _{{\textbf{X}}} \phi _0) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\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
\begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} \phi _1 \right) &= 0, \end{align}
\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}
\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
\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
\begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} \boldsymbol{\Psi } \right) &= 0, \end{align}
\begin{align} \left [\textbf{n}_0 \cdot \left ( \varepsilon (\nabla _{{\textbf{X}}} \boldsymbol{\Psi } + I)\right) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\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
\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}
\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}
\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
\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
\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
\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
\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
\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
\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
\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
\begin{align} \nabla _{{\textbf{X}}}^2 \boldsymbol{\Psi } &= \textbf{0} \quad \text {in} \quad D,\\[-5pt]\nonumber \end{align}
\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}
\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
\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
\begin{align} \nabla _{{\textbf{X}}}\cdot ( \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \phi _0) & = 0 \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}
\begin{align} \nabla _{{\textbf{X}}} \phi _0 &=\textbf{0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}
\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
\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}
\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}
\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
\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}
\begin{align} \nabla _{{\textbf{X}}} \boldsymbol{\Psi } + I &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}},\\[-12pt]\nonumber \end{align}
\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
\begin{equation} \int _{D} \boldsymbol{\Psi } \, {\textrm {d}} {\textbf{X}} = \textbf{0}. \end{equation}
Equating coefficients of
$\delta ^0$
we find
\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}
\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}
\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
\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 [Reference Chapman and McBurnie5], 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 [Reference Chapman and McBurnie5] 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
\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
\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 [Reference Bennett, D’Alessandro and Daly2, Reference Daly and Roose9] (on much more complicated problems).Footnote
4
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
\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
\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
\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
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
\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 [Reference Chapman and McBurnie5], we find
\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, [Reference Davis, Snider and Davis10]. 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
\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
\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
\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
\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
\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
\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
\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
\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
\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
\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
\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 [Reference Chapman and McBurnie5] 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
\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
\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
\begin{align} \nabla \cdot ( \varepsilon _{\mathrm {ex}} \nabla \phi) & = -\frac {\rho }{\delta }\qquad \mbox{in }\Omega _{\mathrm {ex}}, \end{align}
\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
\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
\begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} \phi _1 \right) &= -\rho \quad \text {in} \quad D, \end{align}
\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}
\begin{align} \left [ \phi _1 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
with
$\phi _1$
$\textbf{1}$
-periodic in
$\textbf{X}$
. We write
\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
\begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} {\xi } \right) &= -\rho \quad \text {in} \quad D, \end{align}
\begin{align} \left [\textbf{n}_0 \cdot \varepsilon \nabla _{{\textbf{X}}} {\xi } \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\begin{align} \left [ {\xi } \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
with
\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
\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}
\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}
\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
\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,
\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
\begin{align} \nabla _{{\textbf{X}}} \cdot (\varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \xi) &= - \rho \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}
\begin{align} \nabla _{{\textbf{X}}} \xi &= \boldsymbol {0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}
\begin{align} \textbf{n}_0 \cdotp \nabla _{{\textbf{X}}} \xi &= 0 \quad \text {on} \quad \partial D_{\mathrm {in}}, \end{align}
\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$
,
\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
\begin{align} \nabla _{{\textbf{X}}}\cdot ( \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \phi _1) & = -\rho \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}
\begin{align} \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0 &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}
\begin{align} \left [ \phi _1 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\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
\begin{align} \nabla _{{\textbf{X}}} \cdot ( \varepsilon _{\mathrm {ex}} \nabla _{{\textbf{X}}} \xi) &= -\rho \quad \text {in} \quad D_{\mathrm {ex}}, \end{align}
\begin{align} \nabla _{{\textbf{X}}} \xi &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}
\begin{align} \left [ \xi \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
with
\begin{equation} \int _{D} \xi {\textrm {d}}\textbf{X } = 0. \end{equation}
Equating coefficients at next order, we have
\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}
\begin{align} \nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1 &= \textbf{0} \quad \text {in} \quad D_{\mathrm {in}}, \end{align}
\begin{align} \left [ \phi _2 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
with
\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
\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
\begin{equation} \nabla _{{\textbf{x}}} \cdot (\varepsilon _{\textrm {eff}} \nabla _{{\textbf{x}}} \phi _0) = - \rho _{\textrm {eff}} \end{equation}
where
\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
\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).
Acknowledgements
The authors would like to acknowledge one of the anonymous reviewers, who encouraged them to address some of the confusion regarding the expansion of the normal, which led to the development of Appendix A.
Financial support
A.K. thanks the BBSRC for support under grant BB/M011224/1.
Competing interest
None.
Appendix A. Perturbations of the normal
There are two main ways of describing the microstructure which are prevalent in the literature. Either the position vector of interface is given, so that the interface is described by an equation of the form
\begin{equation} {\textbf{X}} = \textbf{r}({\textbf{x}},s) \end{equation}
say, where
$s$
parameterises the interface, or the interface is given implicitly as the zero-level set of a given function,
\begin{equation} h({\textbf{x}},{\textbf{X}})=0, \end{equation}
say Reference Dalwadi[8, Reference Richardson and Chapman21]. In fact often the former is written as [Reference Holmes12]
\begin{equation*} {\textbf{X}} = \textbf{r}({\textbf{x}},{\textbf{X}}),\end{equation*}
but this is an implicit equation for the interface and is really just the latter with
\begin{equation*} h({\textbf{x}},{\textbf{X}}) = {\textbf{X}} - \textbf{r}({\textbf{x}},{\textbf{X}}).\end{equation*}
Since the functions
$\textbf{r}$
or
$h$
are given and independent of
$\delta$
, many authors suppose that the normal to the microstructure interface
$\textbf{n}$
is also independent of
$\delta$
Reference Burridge and Keller[4, Reference Holmes12, Reference Penta, Ambrosi and Shipley16, Reference Penta, Ambrosi and Quarteroni17]. However, this is a mistake. This is due to the fact that (using the definition
${\textbf{x}} = \delta {\textbf{X}}$
), the microstructure depends on
$\delta$
even though
$h$
and
$\textbf{r}$
do not, and is illustrated in our calculation in Appendix C, and the example which follows. This is a subtle point, which perhaps explains why this is such a common mistake.
Whether the microstructure is described by a parameterised boundary (A.1) or using a level set function (A.2) is irrelevant – if handled correctly the two approaches will give the same answer. If the microstructure is exactly periodic, so that the geometry of the unit cell does not vary with the macroscale, then the normal
$\textbf{n}$
is independent of the macroscale, and can be taken to be independent of
$\delta$
. But when the microstructure varies with macroscopic position the normal will generically depend on
$\delta$
and it is crucial to expand it in powers of
$\delta$
. The level-set approach was developed to facilitate this expansion.
Both approaches above treat the microstructure (i.e. the shape of the inclusion) to be given, and then
$\textbf{n}$
to be derived from this shape. Anytime that a quantity is derived rather than given as part of the problem, it should be expanded in powers of
$\delta$
. An alternative point of view is to not describe the microstructure, but treat the normal
$\textbf{n}$
as the given function, rather than deriving it from
$h$
or
$\textbf{r}$
. In this case
$\textbf{n}$
can indeed be taken to be independent of
$\delta$
. But if
$\textbf{n}$
is to be given, then the shape of the microstructure has to be derived. If
$\textbf{n}$
is assumed independent of
$\delta$
, then the position of the interface
$\textbf{r}$
(or the level set function
$h$
) must be expanded in powers of
$\delta$
, with higher-order terms chosen to be consistent with
$\textbf{n}$
being independent of
$\delta$
. Then all boundary conditions should be expanded onto the leading-order interface, as we do in Appendix B. This is a painful process, and nowhere is it done in any literature we could find, and certainly not in Reference Holmes[12, Reference Penta, Ambrosi and Shipley16, Reference Penta, Ambrosi and Quarteroni17]. These references all assume that both
$\textbf{n}$
and the shape of the microstructure can be chosen to be independent of
$\delta$
, which is a mistake. In section B.3 we show that fixing
$\textbf{n}$
and expanding
$h$
in the present problem leads to the same homogenised equations as those found by fixing
$h$
and expanding
$\textbf{n}$
, albeit with considerably more work.
To illustrate the discussion above with an example, let us consider in two dimensions a microstructure comprising an array of circles with macroscopically varying radius.
A.1 Example: slowly varying circular inclusions
Using (A.1), for circular inclusions of slowly varying radius
$a({\textbf{x}})$
we have
\begin{equation} \textbf{r}(s) = a({\textbf{x}})(\cos s, \sin s), \end{equation}
where
$s$
parameterises the surface. It is crucial to remember here that
$\textbf{x}$
depends on
$s$
also. Then the tangent vector is proportional to
\begin{eqnarray*} \textbf{r}'(s) &=& a({\textbf{x}})(-\sin s, \cos s) + {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a(\cos s, \sin s)\\ & = & (-a \sin s + {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a \cos s,a \cos s + {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a \sin s). \end{eqnarray*}
Thus the normal is
\begin{equation*} \textbf{n} = \frac {(a \cos s + {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a \sin s,a \sin s - {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a \cos s)}{|(a \cos s + {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a \sin s,a \sin s - {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a \cos s)|}. \end{equation*}
Now since
${\textbf{x}} = \delta {\textbf{X}}$
,
${\textbf{x}}'(s) = \delta a(-\sin s,\cos s) + O(\delta ^2)$
. Then
\begin{equation*} |(a \cos s + {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a \sin s,a \sin s - {\textbf{x}}'(s)\cdot {\nabla }_{{\textbf{x}}} a \cos s)| = a + O(\delta ^2).\end{equation*}
Thus
\begin{equation} \textbf{n} = (\cos s, \sin s) + \delta ( (-\sin s,\cos s)\cdot {\nabla }_{{\textbf{x}}} a) ( \sin s,- \cos s). \end{equation}
The correction to the normal is minus the projection of
${\nabla }_{{\textbf{x}}} a$
onto the tangent
$(-\sin s,\cos s)$
.
Now let us do the same thing via the level set method (A.2). We have
\begin{equation*}h = |{\textbf{X}}|-a({\textbf{x}}).\end{equation*}
Then
\begin{eqnarray*} \textbf{n} & = & \frac {{\nabla }_{{\textbf{X}}} h }{|{\nabla }_{{\textbf{X}}} h|} + \delta \left (\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} \right) +O(\delta ^2). \end{eqnarray*}
Now
\begin{equation*} {\nabla }_{{\textbf{X}}} h = \frac {{\textbf{X}}}{|{\textbf{X}}|} = \frac {{\textbf{X}}}{a}, \qquad {\nabla }_{{\textbf{x}}} h = - {\nabla }_{{\textbf{x}}} a, \end{equation*}
giving
\begin{eqnarray*} \textbf{n} & = & \frac {{\textbf{X}}}{a} + \delta \left (- {\nabla }_{{\textbf{x}}} a + \left ( {\nabla }_{{\textbf{x}}} a\cdot \frac {{\textbf{X}}}{a}\right)\frac {{\textbf{X}}}{a} \right) +O(\delta ^2). \end{eqnarray*}
The correction is
$-{\nabla }_{{\textbf{x}}} a$
minus its projection onto the normal, which is the same as the projection onto the tangent. Evaluating explicitly on
${\textbf{X}} = a(\cos s,\sin s)$
gives
\begin{eqnarray*} \textbf{n} & = & (\cos s, \sin s) + \delta ( (-\sin s,\cos s)\cdot {\nabla }_{{\textbf{x}}} a) ( \sin s,- \cos s) \end{eqnarray*}
in agreement with (A.4).
If, in contrast, we assume that the normal is given and independent of
$\delta$
, then, using (A.1) and expanding
$\textbf{r}$
as
\begin{equation*} \textbf{r}(s) = a({\textbf{x}})(\cos s, \sin s) + \delta \textbf{r}_1(s) + \cdots,\end{equation*}
we have
\begin{eqnarray*} \textbf{r}'(s) &=& a({\textbf{x}})(-\sin s, \cos s) + \delta a \left ((-\sin s,\cos s)\cdot {\nabla }_{{\textbf{x}}} a\right)(\cos s, \sin s) + \delta \textbf{r}_1'(s) +O(\delta ^2). \end{eqnarray*}
For the normal to have no
$O(\delta)$
correction, we require
\begin{equation*} \textbf{r}_1'(s) + a \left ((-\sin s,\cos s)\cdot {\nabla }_{{\textbf{x}}} a\right)(\cos s, \sin s)\end{equation*}
to be parallel to
$(-\sin s, \cos s)$
, i.e.
\begin{equation*} (\cos s, \sin s) \cdot \textbf{r}_1'(s) = -a (-\sin s,\cos s)\cdot {\nabla }_{{\textbf{x}}} a. \end{equation*}
In fact, this correction to the interface is exactly what is needed if the inclusions are exactly circles, rather than the approximate circles described by (A.3). However, this perturbation to the position of the boundary must be taken into account when asymptotically expanding variables in powers of
$\delta$
to develop a homogenised model, as we do in Appendix B. This is quite a complicated process; the approach of fixing the microstructure and expanding the normal is far easier.
A.2 Implication of this mistake
In many published works, the failure to expand either the position of the interface or the normal to it in powers of
$\delta$
when the microstructure varies with macroscopic position leads to erroneous homogenised equations containing spurious terms.
Let us illustrate by re-examining the classical work of Reference Burridge and Keller[4]. It is not unambiguously clear whether Burridge and Keller intend for their analysis to be valid when the microstructure varies macroscopically. They certainly intend the microscopic properties to vary macroscopically, but whether they intend geometry to be included as such a property could be open to debate. They do, however, take the characteristic function
$\chi _s$
of the solid domain to be a function of both the slow and fast scales, which leads us to believe that they were considering slowly varying domains; certainly, this is how the manuscript is usually interpreted, with many authors applying their results (erroneously) to macroscopically varying microstructures.
Let us see how to correct the work so that it does apply to slowly varying microstructures. The normal appears only in equation (6e ) of [Reference Burridge and Keller4], which expresses continuity of stress across the fluid/solid interface:
\begin{equation} \textbf{n} \cdot \sigma - \textbf{n} \cdot \tau = 0, \qquad \mbox{on }{\partial } D_f = {\partial } D_s, \end{equation}
At leading order, this gives
\begin{equation} \textbf{n}_0 \cdot \sigma _0 - \textbf{n}_0 \cdot \tau _0 = 0, \qquad \mbox{on }{\partial } D_f = {\partial } D_s, \end{equation}
but at first order, the equation should be
\begin{equation} \textbf{n}_0 \cdot \sigma _1 - \textbf{n}_0 \cdot \tau _1+\textbf{n}_1 \cdot \sigma _0 - \textbf{n}_1 \cdot \tau _0 = 0, \qquad \mbox{on }{\partial } D_f = {\partial } D_s; \end{equation}
Burridge and Keller Reference Burridge and Keller[4] miss out on the terms proportional to
$\textbf{n}_1$
. For a microstructure periodic on the fast scale, with a unit volume unit cell (so that
$V_f=\int _{D_f}\, {\textrm {d}} \textbf{y}$
is still the volume fraction of fluid), their equation (20) is
\begin{equation} \overline {{\nabla }_{{\textbf{x}}} \cdot \sigma _0} + V_f \omega ^2 \rho _f {\textbf{u}}_0 + \omega ^2 \rho _f \bar {\textbf{w}} = - \int _{{\partial } D_f} \sigma _1 \cdot \textbf{n}_0\, {\textrm {d}} \textbf{y} \end{equation}
where we have not substituted
$\overline {{\nabla }_{{\textbf{x}}} \cdot \sigma _0} = V_f {\nabla }_{{\textbf{x}}} \cdot \sigma _0$
as in [Reference Burridge and Keller4]. Their equation (21) is
\begin{equation} \overline {{\nabla }_{{\textbf{x}}} \cdot \tau _0} + \omega ^2 \bar {\rho }_s {\textbf{u}}_0 = \int _{{\partial } D_s} \tau _1 \cdot \textbf{n}_0\, {\textrm {d}} \textbf{y}. \end{equation}
Adding these equations [Reference Burridge and Keller4] claims that the surface integrals cancel, but in fact, this is not the case. Correcting their equation (27) gives, instead
\begin{equation} - \omega ^2 \bar {\rho }_s {\textbf{u}}_0 - \omega ^2 \rho _f \bar {\textbf{w}} = \overline {{\nabla }_{{\textbf{x}}} \cdot \tau _0} +\overline {{\nabla }_{{\textbf{x}}} \cdot \sigma _0} - \int _{D_f}\sigma _0 \cdot \textbf{n}_1 \, {\textrm {d}} \textbf{y}+\int _{{\partial } D_s} \tau _0 \cdot \textbf{n}_1 \, {\textrm {d}} \textbf{y}. \end{equation}
Now
\begin{equation*}{\nabla }_{{\textbf{x}}} \cdot \overline{\tau _0} =\overline {{\nabla }_{{\textbf{x}}} \cdot \tau _0} - \int _{{\partial } D_s} \tau _0 \cdot \textbf{V} \cdot \textbf{n}_0\, {\textrm {d}} S_y, \qquad {\nabla }_{{\textbf{x}}} \cdot \overline{\sigma _0} =\overline {{\nabla }_{{\textbf{x}}} \cdot \sigma _0} +\int _{{\partial } D_f} \sigma _0 \cdot \textbf{V} \cdot \textbf{n}_0\, {\textrm {d}} S_y. \end{equation*}
Thus (27) becomes
\begin{align} - \omega ^2 \bar {\rho }_s {\textbf{u}}_0 - \omega ^2 \rho _f \bar {\textbf{w}} = &{\nabla }_{{\textbf{x}}} \cdot \overline{\tau _0} + \int _{{\partial } D_s} \tau _0 \cdot \textbf{V} \cdot \textbf{n}_0\, {\textrm {d}} S_y +{\nabla }_{{\textbf{x}}} \cdot \overline{\sigma _0}\nonumber \\& -\int _{{\partial } D_f} \sigma _0 \cdot \textbf{V} \cdot \textbf{n}_0\, {\textrm {d}} S_y -\int _{D_f}\sigma _0 \cdot \textbf{n}_1 \, {\textrm {d}} \textbf{y} + \int _{{\partial } D_s} \tau _0 \cdot \textbf{n}_1 \, {\textrm {d}} \textbf{y}. \end{align}
As we have seen,
\begin{equation*} \textbf{V}\cdot \textbf{n}_0 + \textbf{n}_1 = -\frac {{\nabla }_{{\textbf{x}}} h \cdot {\nabla }_{\textbf{y}} h}{|{\nabla }_{\textbf{y}} h|^2}\textbf{n}_0.\end{equation*}
Since
$\sigma _0 \cdot \textbf{n}_0= \tau _0 \cdot \textbf{n}_0$
the remaining surface integrals cancel, leaving
\begin{equation} - \omega ^2 \bar {\rho }_s {\textbf{u}}_0 - \omega ^2 \rho _f \bar {\textbf{w}} = {\nabla }_{{\textbf{x}}} \cdot \overline{\tau _0} +{\nabla }_{{\textbf{x}}} \cdot \overline{\sigma _0}= {\nabla }_{{\textbf{x}}} \cdot \overline{\tau _0} -{\nabla }_{{\textbf{x}}} (V_f p_0). \end{equation}
Note that (Reference Burridge and Keller4) do not ever attempt to write
$\overline {{\nabla }_{{\textbf{x}}} \cdot \tau _0}$
in terms of
${\nabla }_{{\textbf{x}}} \cdot \overline{\tau _0}$
except for the macroscopically uniform case, so that their final system of equations (35)–(36) is not closed. If they did so, using Reynold’s transport theorem, then there would be additional surface integral terms in their homogenised equation which should not be present, as written explicitly in equations (5.47)–(5.48) of (Reference Penta, Ambrosi and Shipley16), for example. Note also that correctly writing
$\overline {{\nabla }_{{\textbf{x}}} \cdot \tau _0}$
in terms of
${\nabla }_{{\textbf{x}}} \cdot \overline {\tau _0}$
brings
$V_f$
inside the gradient in the last term. This might look strange to those familiar with two-phase flow (it appears that a uniform pressure in the fluid can generate stress gradients in the solid), but remember that
$\overline {\tau _0}$
is a superficial average rather than a phase average – a static problem with a uniform pressure within the fluid would generate the same uniform pressure within the solid, for which
$\overline {\tau _0} = -(1-V_f)p_0I$
.
Appendix B. Perturbations of the boundary
In section 2–3, we illustrated the theory with an example in which the boundaries of the inclusions were defined by the level sets of (2.9). With this definition of the boundary, the inclusions are approximately circular, but not exactly circular: as the slow coordinate varies over a given inclusion, the radius of the inclusion
$a({\textbf{x}})$
varies by an amount
$O(\delta)$
, so that these inclusions are actually elliptical, but with an eccentricity of
$O(\delta)$
. We would not expect such a small perturbation of the shape to affect the leading-order homogenised problem, which is why so many authors treat circular inclusions in this way. On the other hand, an
$O(\delta)$
perturbation to
$h$
would introduce new terms into
$\textbf{n}_1$
, and it is not obvious a priori that these would not appear in the homogenised model. In this appendix, we verify that
$O(\delta)$
corrections to
$h$
do not affect the leading-order homogenised problem.
We return to the problem of a dielectric composite considered in section 2, which we restate for convenience
\begin{equation} \nabla \cdot ( \varepsilon \nabla \phi) = -\rho, \end{equation}
with
\begin{align} \left [\textbf{n} \cdot ( \varepsilon \nabla \phi) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\begin{align} \left [ \phi \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}
We suppose that the boundary of the inclusion is given by the level set of
$h({\textbf{X}},{\textbf{x}})$
, but this time we suppose that
$h$
depends on
$\delta$
also, so that
\begin{equation*} h({\textbf{X}},{\textbf{x}}) = h_0({\textbf{X}},{\textbf{x}}) + \delta h_1({\textbf{X}},{\textbf{x}}) + \cdots.\end{equation*}
The expansion of the normal is then
\begin{equation} \textbf{n} = \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|} + \delta \left (\frac {\nabla _{{\textbf{X}}} h_1 + \nabla _{{\textbf{x}}} h_0}{|\nabla _{{\textbf{X}}} h_0|} - \nabla _{{\textbf{X}}} h_0 \cdot (\nabla _{{\textbf{X}}} h_1 + \nabla _{{\textbf{x}}} h_0) \frac {\nabla _{{\textbf{X}}} h_0}{|\nabla _{{\textbf{X}}} h_0 |^3} \right) + O(\delta ^2). \end{equation}
B.1 Expansion onto the leading-order domain
In addition to performing a multiple scales expansion, we must expand (B.2)–(B.3) onto the leading-order boundary position. We suppose that the boundary perturbation is such that positions
${\textbf{X}}(\mathbf {s})$
on the boundary may be expanded as
${\textbf{X}}(\mathbf {s}) = {\textbf{X}}_0(\mathbf {s}) + \delta \textbf{R}({\mathbf {s}})$
. To expand a quantity onto the leading-order boundary, we then write
\begin{equation} \begin{split} \textbf{Q}({\textbf{x}}, {\textbf{X}}_0 + \delta \textbf{R}) = \textbf{Q}_0({\textbf{x}}, & {\textbf{X}}_0) + \delta \left ( \textbf{Q}_1({\textbf{x}}, {\textbf{X}}_0) + (\textbf{R} \cdot \nabla _{{\textbf{X}}}) \textbf{Q}_0({\textbf{x}}, {\textbf{X}}_0) \right) + \cdots. \end{split} \end{equation}
Note that when we write the normal in terms of
$h$
as in (B.4), all these terms are evaluated on
${\textbf{X}}(\mathbf {s}) = {\textbf{X}}_0(\mathbf {s}) + \delta \textbf{R}$
. Expanding again to evaluate the leading-order boundary gives
\begin{equation} \textbf{n} = \textbf{n}_0 + \delta \textbf{n}_1 + \cdots, \end{equation}
with
\begin{align} \textbf{n}_0 & = \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|}, \end{align}
\begin{align} \textbf{n}_1 & = (\textbf{R} \cdot \nabla _{{\textbf{X}}}) \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|}+ \left (\frac {\nabla _{{\textbf{X}}} h_1 + \nabla _{{\textbf{x}}} h_0}{|\nabla _{{\textbf{X}}} h_0|} - \nabla _{{\textbf{X}}} h_0 \cdot (\nabla _{{\textbf{X}}} h_1 + \nabla _{{\textbf{x}}} h_0) \frac {\nabla _{{\textbf{X}}} h_0}{|\nabla _{{\textbf{X}}} h_0 |^3} \right), \end{align}
where all terms are now evaluated at
${\textbf{X}}_0$
. Finally, we need to relate
$\textbf{R}$
to
$h$
. Since
$h({\textbf{x}},{\textbf{X}}) = h({\textbf{x}},{\textbf{X}}_0+\delta \textbf{R})=0$
, expanding and equating coefficients at
$O(\delta)$
gives
$h_1 = - \textbf{R} \cdot \nabla _{{\textbf{X}}} h_0$
.
B.2 Homogenisation
Combining the multiple scales expansion (2.12) with the normal expansion (B.6) and boundary expansion, we find at leading-order
\begin{align} \nabla _{{\textbf{X}}} \cdotp (\varepsilon \nabla _{{\textbf{X}}} \phi _0) &= 0 \quad \text {in} \quad D, \end{align}
\begin{align} \left [\textbf{n}_0 \cdot (\varepsilon \nabla _{{\textbf{X}}} \phi _0) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\begin{align} \left [ \phi _0 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
with
$\phi _0$
$\boldsymbol {1}$
-periodic in
$\textbf{X}$
, where the jump is across the leading-order boundary. Thus, at leading-order, we have a potential which is independent of the fast-scale. At next order, we obtain
\begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} \phi _1 \right) &= 0 \quad \text {in} \quad D, \end{align}
\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}
\begin{align} \left [ \phi _1 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}
Writing
$ \phi _1 = \boldsymbol{\Psi }\cdot \nabla _{{\textbf{x}}} \phi _0 + \overline {\phi }_1$
for
$\overline {\phi }_1$
independent of
$\textbf{X}$
, we find that
$\Psi$
satisfies the cell problem given by (2.20)–(2.22). At
$O(\delta ^2)$
, we find
\begin{align} \nabla _{{\textbf{X}}} \cdot \textbf{Q}_1 + \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 &= -\rho \quad \text {in} \quad D, \end{align}
\begin{align} [\textbf{n}_0 \cdot \textbf{Q}_1 + \textbf{n}_1 \cdot \textbf{Q}_0 + \textbf{n}_0 \cdot ( \textbf{R} \cdot \nabla _{{\textbf{X}}})\textbf{Q}_0]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\begin{align} [\phi _2 + \textbf{R} \cdot \nabla _{{\textbf{X}}} \phi _1]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}
where
\begin{equation*} \textbf{Q}_0 = \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0, \qquad \textbf{Q}_1 =\nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1. \end{equation*}
Integrating (B.15) over the unit cell, we have
\begin{equation} \int _{D} \nabla _{{\textbf{X}}} \cdot \textbf{Q}_1 \, {\textrm {d}} \textbf{X} + \int _{D} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}, \end{equation}
where
$\rho _{\textrm {eff}}$
is given by (2.27). Using the divergence theorem in the first term, we obtain
\begin{equation} - \int _{\partial D_{\mathrm {ex}}} \textbf{Q}_1 \cdot \textbf{n}_0 \, {\textrm {d}} S_{\textbf{X}} + \int _{\partial D_{\mathrm {in}}} \textbf{Q}_1 \cdot \textbf{n}_0 \, {\textrm {d}} S_{\textbf{X}} + \int _{D} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{equation}
Using periodicity to eliminate the integral over the boundary of the unit cell in
${\partial } D_{\mathrm {ex}}$
, and substituting (B.16), we find
\begin{equation} \int _{\partial D_{\mathrm {in}}} \left [\textbf{n}_0 \cdot ( \textbf{R} \cdot \nabla _{{\textbf{X}}}) \textbf{Q}_0 \right]^{\mathrm {ex}}_{{\mathrm {in}}} \, {\textrm {d}} S_{\textbf{X}} + \int _{\partial D_{\mathrm {in}}} \left [\textbf{Q}_0 \cdot \textbf{n}_1\right]^{\mathrm {ex}}_{\mathrm {in}} {\textrm {d}} S_{\textbf{X}} + \int _{D} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{equation}
Applying the Reynolds transport theorem to the final integral, we obtain
\begin{equation} \int _{\partial D_{\mathrm {in}}} \left [\textbf{n}_0 \cdot ( \textbf{R} \cdot \nabla _{{\textbf{X}}}) \textbf{Q}_0 \right]^{\mathrm {ex}}_{{\mathrm {in}}} \, {\textrm {d}} S_{\textbf{X}} + \int _{\partial D_{\mathrm {in}}} \left [\textbf{Q}_0 \cdot \textbf{n}_1 + \textbf{Q}_0 \cdot \textbf{V}_0 \cdot \textbf{n}_0 \right]^{\mathrm {ex}}_{\mathrm {in}} {\textrm {d}} S_{\textbf{X}} + \nabla _{{\textbf{x}}} \cdot \int _{D} \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{equation}
The analogue of (2.29) is
\begin{equation} \textbf{V}_0 \cdot \textbf{n}_0 + \textbf{n}_1 = (\textbf{R} \cdot \nabla _{{\textbf{X}}}) \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|}+\frac {\nabla _{{\textbf{X}}} h_1}{|\nabla _{{\textbf{X}}} h_0|} - \frac {\nabla _{{\textbf{X}}} h_0 \cdotp (\nabla _{{\textbf{X}}} h_1 + \nabla _{{\textbf{x}}} h_0)}{|\nabla _{{\textbf{X}}} h_0|^2} \textbf{n}_0. \end{equation}
From (B.13),
$\textbf{Q}_0 \cdot \textbf{n}_0$
is continuous, which leaves
\begin{align} \int _{\partial D_{\mathrm {in}}}\left [\textbf{n}_0\cdot ( \textbf{R} \cdot \nabla _{{\textbf{X}}}) \textbf{Q}_0 + \textbf{Q}_0 \cdot (\textbf{R} \cdot \nabla _{{\textbf{X}}}) \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|}+\frac {\nabla _{{\textbf{X}}} h_1}{|\nabla _{{\textbf{X}}} h_0|} \cdot \textbf{Q}_0 \right]^{\mathrm {ex}}_{\mathrm {in}} {\textrm {d}} S_{\textbf{X}}+ \nabla _{{\textbf{x}}} \cdot \int _{D} \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{align}
Now
\begin{equation*} \nabla _{{\textbf{X}}} \times ( \textbf{R} \times \textbf{Q}_0) = \textbf{R} \nabla _{{\textbf{X}}} \cdot \textbf{Q}_0 - \textbf{Q}_0 \nabla _{{\textbf{X}}} \cdot \textbf{R} + ( \textbf{R} \cdot \nabla _{{\textbf{X}}})\textbf{Q}_0 - (\textbf{Q}_0\cdot \nabla _{{\textbf{X}}}) \textbf{R}.\end{equation*}
Thus, using that
$\textbf{Q}_0 \cdot \textbf{n}_0$
is continuous,
$\nabla \cdot \textbf{Q}_0=0$
, and Stokes’ theorem,
\begin{equation*} \int _{\partial D_{\mathrm {in}}}\left [\textbf{n}_0\cdot ( \textbf{R} \cdot \nabla _{{\textbf{X}}}) \textbf{Q}_0 \right]^{\mathrm {ex}}_{\mathrm {in}} {\textrm {d}} S_{\textbf{X}} =\int _{\partial D_{\mathrm {in}}}\left [\textbf{n}_0\cdot (\textbf{Q}_0 \cdot \nabla _{{\textbf{X}}}) \textbf{R} \right]^{\mathrm {ex}}_{\mathrm {in}} {\textrm {d}} S_{\textbf{X}} \end{equation*}
leaving, after substituting for
$h_1$
, and writing
$\textbf{n}_0$
in terms of
$h_0$
,
\begin{align*} \int _{\partial D_{\mathrm {in}}}\left [ \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|} \cdot (\textbf{Q}_0 \cdot \nabla _{{\textbf{X}}}) \textbf{R} + \textbf{Q}_0 \cdot (\textbf{R} \cdot \nabla _{{\textbf{X}}}) \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|} - \frac {\nabla _{{\textbf{X}}} ( \textbf{R}\cdot \nabla _{{\textbf{X}}} h_0)}{|\nabla _{{\textbf{X}}} h_0|} \cdot \textbf{Q}_0 \right]^{\mathrm {ex}}_{\mathrm {in}} {\textrm {d}} S_{\textbf{X}} + \nabla _{{\textbf{x}}} \cdot \int _{D} \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{align*}
In index notation, with
$\textbf{Q}_0 = (Q_1,Q_2,Q_3)$
,
$\textbf{R} = (R_1,R_2,R_3)$
,
\begin{align*} & {\frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|} \cdot (\textbf{Q}_0 \cdot \nabla _{{\textbf{X}}}) \textbf{R} + \textbf{Q}_0 \cdot (\textbf{R} \cdot \nabla _{{\textbf{X}}}) \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|} - \frac {\nabla _{{\textbf{X}}} ( \textbf{R}\cdot \nabla _{{\textbf{X}}} h_0)}{|\nabla _{{\textbf{X}}} h_0|} \cdot \textbf{Q}_0 } \\ &\quad = \frac {1}{|\nabla _{{\textbf{X}}} h_0|} \mathchoice {\frac {{\partial }h_0}{{\partial }X_i}}{{\partial }h_0/{\partial }X_i}{{\partial }h_0/{\partial }X_i}{{\partial }h_0/{\partial }X_i} \left (Q_j \mathchoice {\frac {{\partial }R_i}{{\partial }X_j}}{{\partial }R_i/{\partial }X_j}{{\partial }R_i/{\partial }X_j}{{\partial }R_i/{\partial }X_j}\right) + Q_j R_i \mathchoice {\frac {{\partial }}{{\partial }X_i}}{{\partial }/{\partial }X_i}{{\partial }/{\partial }X_i}{{\partial }/{\partial }X_i} \left ( \mathchoice {\frac {{\partial }h_0}{{\partial }X_j}}{{\partial }h_0/{\partial }X_j}{{\partial }h_0/{\partial }X_j}{{\partial }h_0/{\partial }X_j}\frac {1 }{|\nabla _{{\textbf{X}}} h_0|}\right) - \frac {1}{|\nabla _{{\textbf{X}}} h_0|} Q_j \mathchoice {\frac {{\partial }}{{\partial }X_j}}{{\partial }/{\partial }X_j}{{\partial }/{\partial }X_j}{{\partial }/{\partial }X_j} \left (R_i\mathchoice {\frac {{\partial }h_0}{{\partial }X_i}}{{\partial }h_0/{\partial }X_i}{{\partial }h_0/{\partial }X_i}{{\partial }h_0/{\partial }X_i}\right) \\ &\quad = Q_j R_i \mathchoice {\frac {{\partial }h_0}{{\partial }X_j}}{{\partial }h_0/{\partial }X_j}{{\partial }h_0/{\partial }X_j}{{\partial }h_0/{\partial }X_j}\mathchoice {\frac {{\partial }}{{\partial }X_i}}{{\partial }/{\partial }X_i}{{\partial }/{\partial }X_i}{{\partial }/{\partial }X_i} \left (\frac {1 }{|\nabla _{{\textbf{X}}} h_0|}\right) = -(\textbf{Q}_0 \cdot \textbf{n}_0) (\textbf{R} \cdot \nabla _{{\textbf{X}}})|\nabla _{{\textbf{X}}} h_0|. \end{align*}
Since
$\textbf{Q}_0 \cdot \textbf{n}_0$
is continuous, we are left with
\begin{equation*}\nabla _{{\textbf{x}}} \cdot \int _{D} \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}\end{equation*}
as expected.
B.3 Homogenisation with
$\textbf{n}$
independent of
$\delta$
The calculation in section B.2 can be used to show that the alternative approach of fixing the normal
$\textbf{n}$
to be independent of
$\delta$
gives the same homogenised problem providing the geometry (through
$h$
or
$\textbf{r}$
) is expanded in powers of
$\delta$
.
To see this explicitly, suppose now that
$h_1$
is such that
$\textbf{n}_1 = \textbf{0}$
. As usual, we find at leading-order
\begin{align} \nabla _{{\textbf{X}}} \cdotp (\varepsilon \nabla _{{\textbf{X}}} \phi _0) &= 0 \quad \text {in} \quad D, \end{align}
\begin{align} \left [\textbf{n}_0 \cdot (\varepsilon \nabla _{{\textbf{X}}} \phi _0) \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\begin{align} \left [ \phi _0 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
with
$\phi _0$
$\boldsymbol {1}$
-periodic in
$\textbf{X}$
, where the jump is across the leading-order boundary, and we have continued to use the symbol
$\textbf{n}_0$
even though
$\textbf{n}$
is now independent of
$\delta$
. Thus, as usual,
$\phi _0$
is independent of
$\textbf{X}$
. At next order, we obtain
\begin{align} \nabla _{{\textbf{X}}} \cdot \left ( \varepsilon \nabla _{{\textbf{X}}} \phi _1 \right) &= 0 \quad \text {in} \quad D, \end{align}
\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}
\begin{align} \left [ \phi _1 \right]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}
Writing
$ \phi _1 = \boldsymbol{\Psi }\cdot \nabla _{{\textbf{x}}} \phi _0 + \overline {\phi }_1$
for
$\overline {\phi }_1$
independent of
$\textbf{X}$
, we find that
$\Psi$
satisfies the usual cell problem (2.20)–(2.22). At
$O(\delta ^2)$
, we find
\begin{align} \nabla _{{\textbf{X}}} \cdot \textbf{Q}_1 + \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 &= -\rho \quad \text {in} \quad D, \end{align}
\begin{align} [\textbf{n}_0 \cdot \textbf{Q}_1 + \textbf{n}_0 \cdot ( \textbf{R} \cdot \nabla _{{\textbf{X}}})\textbf{Q}_0]^{\mathrm {ex}}_{\mathrm {in}} &= 0, \end{align}
\begin{align} [\phi _2 + \textbf{R} \cdot \nabla _{{\textbf{X}}} \phi _1]^{\mathrm {ex}}_{\mathrm {in}} &= 0. \end{align}
where
\begin{equation*} \textbf{Q}_0 = \nabla _{{\textbf{X}}} \phi _1 + \nabla _{{\textbf{x}}} \phi _0, \qquad \textbf{Q}_1 =\nabla _{{\textbf{X}}} \phi _2 + \nabla _{{\textbf{x}}} \phi _1. \end{equation*}
The terms involving
$\textbf{R}$
here are usually (mistakenly) omitted in the literature [Reference Holmes12, Reference Penta, Ambrosi and Shipley16, Reference Penta, Ambrosi and Quarteroni17].
Integrating (B.30) over the unit cell, we have
\begin{equation} \int _{D} \nabla _{{\textbf{X}}} \cdot \textbf{Q}_1 \, {\textrm {d}} \textbf{X} + \int _{D} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}, \end{equation}
where
$\rho _{\textrm {eff}}$
is given by (2.27). Using the divergence theorem in the first term, we obtain
\begin{equation} - \int _{\partial D_{\mathrm {ex}}} \textbf{Q}_1 \cdot \textbf{n}_0 \, {\textrm {d}} S_{\textbf{X}} + \int _{\partial D_{\mathrm {in}}} \textbf{Q}_1 \cdot \textbf{n}_0 \, {\textrm {d}} S_{\textbf{X}} + \int _{D} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{equation}
Using periodicity to eliminate the integral over the boundary of the unit cell in
${\partial } D_{\mathrm {ex}}$
, and substituting (B.31), we find
\begin{equation} \int _{\partial D_{\mathrm {in}}} \left [\textbf{n}_0 \cdot ( \textbf{R} \cdot \nabla _{{\textbf{X}}}) \textbf{Q}_0 \right]^{\mathrm {ex}}_{{\mathrm {in}}} \, {\textrm {d}} S_{\textbf{X}} + \int _{D} \nabla _{{\textbf{x}}} \cdot \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{equation}
Applying the Reynolds transport theorem to the final integral, we obtain
\begin{equation} \int _{\partial D_{\mathrm {in}}} \left [\textbf{n}_0 \cdot ( \textbf{R} \cdot \nabla _{{\textbf{X}}}) \textbf{Q}_0 \right]^{\mathrm {ex}}_{{\mathrm {in}}} \, {\textrm {d}} S_{\textbf{X}} + \int _{\partial D_{\mathrm {in}}} \left [\textbf{Q}_0 \cdot \textbf{V}_0 \cdot \textbf{n}_0 \right]^{\mathrm {ex}}_{\mathrm {in}} {\textrm {d}} S_{\textbf{X}} + \nabla _{{\textbf{x}}} \cdot \int _{D} \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{equation}
Now, with
$\mathbf{n_1} = \textbf{0}$
, (B.8) implies
$\textbf{R}$
is such that
\begin{eqnarray} \textbf{V}_0 \cdot \textbf{n}_0 & = & -\frac {{\nabla }_{{\textbf{x}}} h_0}{|{\nabla }_{{\textbf{X}}} h_0|} = (\textbf{R} \cdot \nabla _{{\textbf{X}}}) \frac {\nabla _{{\textbf{X}}} h_0 }{|\nabla _{{\textbf{X}}} h_0|}+\frac {\nabla _{{\textbf{X}}} h_1}{|\nabla _{{\textbf{X}}} h_0|} - \frac {\nabla _{{\textbf{X}}} h_0 \cdotp (\nabla _{{\textbf{X}}} h_1 + \nabla _{{\textbf{x}}} h_0)}{|\nabla _{{\textbf{X}}} h_0|^2} \textbf{n}_0.\qquad \end{eqnarray}
The same arguments as in section B.2 then show that the surface integral vanishes, leaving
\begin{equation*}\nabla _{{\textbf{x}}} \cdot \int _{D} \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}\end{equation*}
as expected.
If we did not include the perturbation to the boundary, assuming that both
$h$
and
$\textbf{n}$
are independent of
$\delta$
, then (B.36) would become
\begin{equation} \int _{\partial D_{\mathrm {in}}} \left [\textbf{Q}_0 \cdot \textbf{V}_0 \cdot \textbf{n}_0 \right]^{\mathrm {ex}}_{\mathrm {in}} {\textrm {d}} S_{\textbf{X}} + \nabla _{{\textbf{x}}} \cdot \int _{D} \textbf{Q}_0 \, {\textrm {d}} \textbf{X} = -\rho _{\textrm {eff}}. \end{equation}
Since
\begin{eqnarray} \textbf{V}_0 \cdot \textbf{n}_0 & = & -\frac {{\nabla }_{{\textbf{x}}} h_0}{|{\nabla }_{{\textbf{X}}} h_0|}, \end{eqnarray}
the surface integral does not vanish (unless
${\nabla }_{{\textbf{X}}} h$
is parallel to
${\nabla }_{{\textbf{x}}} h$
) and is an extra term, which should not be present.
Appendix C. Derivation of (2.11)
Here, we explicitly derive (2.11) from (2.10). We have
\begin{eqnarray*} \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 |} =\frac {\nabla _{{\textbf{X}}} h + \delta \nabla _{{\textbf{x}}} h }{((\nabla _{{\textbf{X}}} h + \delta \nabla _{{\textbf{x}}} h)\cdot (\nabla _{{\textbf{X}}} h + \delta \nabla _{{\textbf{x}}} h))^{1/2}} \\ & = & \frac {\nabla _{{\textbf{X}}} h + \delta \nabla _{{\textbf{x}}} h }{(|\nabla _{{\textbf{X}}} h|^2 + 2\delta \nabla _{{\textbf{x}}} h\cdot \nabla _{{\textbf{X}}} h + \delta ^2 |\nabla _{{\textbf{x}}} h|^2)^{1/2}} \\ & = & \frac {\nabla _{{\textbf{X}}} h + \delta \nabla _{{\textbf{x}}} h }{|\nabla _{{\textbf{X}}} h|} \left (1 + \frac {2\delta \nabla _{{\textbf{x}}} h\cdot \nabla _{{\textbf{X}}} h+ \delta ^2 |\nabla _{{\textbf{x}}} h|^2}{|\nabla _{{\textbf{X}}} h|^2} \right)^{-1/2}\\ & = & \frac {\nabla _{{\textbf{X}}} h + \delta \nabla _{{\textbf{x}}} h }{|\nabla _{{\textbf{X}}} h|} \left (1 - \frac {\delta \nabla _{{\textbf{x}}} h\cdot \nabla _{{\textbf{X}}} h}{|\nabla _{{\textbf{X}}} h|^2} +O(\delta ^2)\right)\\ & = & \frac {\nabla _{{\textbf{X}}} h }{|\nabla _{{\textbf{X}}} h|} + \delta \left (\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} \right) +O(\delta ^2). \end{eqnarray*}
Open access
