Let $C(x)$ be the local stiffness tensor for a two phase material with stiffness tensors $C_A$ and $C_B$. The stiffness tensor at location $x$ can be represented at a perturbation from a reference stiffness tensor.
$$C(x) = C^R + C'(x)$$The strain field at location $(x)$ can also be defined in terms of a simular perturbation.
$$\varepsilon(x) = \bar{\varepsilon} + \varepsilon '(x)$$where $\bar{\varepsilon}$ is the average strain and $\varepsilon '(x)$ is the local strain perturbation from $\bar{\varepsilon}$.
The constitutive equation is therefore.
$$\sigma_{ij}(x) = \big(C^R_{ijlk} + C'_{ijlk}(x) \big ) \big (\bar{\varepsilon}_{lk} + \varepsilon'_{lk}(x) \big )$$The equilibrium condition is defined below.
$$\sigma_{ij,j}(x) = \Big [\big(C^R_{ijlk} + C'_{ijlk}(x) \big ) \big (\bar{\varepsilon}_{lk} + \varepsilon'_{lk}(x) \big )\Big ]_{,j} = 0$$$$\sigma_{ij,j}(x) = C^R_{ijlk}\varepsilon'_{lk,j}(x) + C'_{ijlk,j}(x)\bar{\varepsilon}_{lk} + \Big [C'_{ijlk}(x) \varepsilon'_{lk}(x)\Big ]_{,j} = 0$$Let
$$F_i(x) = C'_{ijlk,j}(x)\bar{\varepsilon}_{lk} + \Big [C'_{ijlk}(x) \varepsilon'_{lk}(x)\Big ]_{,j} $$Using the definitation of $F(x)$ above, the equilibrium equation above can be rearranged in the form of an inhomogenous differential equation.
$$C^R_{ijlk}\varepsilon'_{lk,j}(x) + F_i(x) = 0$$By using the relationship between strain and displacement, the equilibrium equation can be rewritten as follows.
$$ \varepsilon_{kl}(x) = \frac{\big (u_{k,l}(x) + u_{l,k}(x) \big)}{2} $$$$C^R_{ijkl} \frac{\big (u'_{k,lj}(x) + u'_{l,kj}(x) \big)}{2} + F_i(x) = 0$$The solution to the displacements can be found using green's functions.
$$C^R_{ijkl} G_{km,lj}(r) + \delta_{im}\delta(x-r) = 0$$$$u'_k(x) = \int_V G_{ik}(r) F_i (x-r)dr = \int_V G_{ik}(r) \Big [C'_{ijlk}(x-r)\bar{\varepsilon}_{lk} + \big [C'_{ijlk}(x-r)\varepsilon'_{lk}(x-r)\big ]\Big ]_{,j}dr$$and
$$u'_l(x) = \int_V G_{il}(r) F_i (x - r)dr = \int_V G_{ik}(r) \Big [C'_{ijlk}(x-r)\bar{\varepsilon}_{lk} + \big [C'_{ijlk}(x-r)\varepsilon'_{lk}(x-r)\big ]\Big ]_{,j}dr$$therefore the strain can also be found interns of green's functions.
$$\varepsilon'_{kl}(x) = \int_V \frac{\big (G_{ik,l}(r) + G_{il,k}(r) \big)}{2} F_i (x-r)dr = \int_V \frac{\big (G_{ik,l}(r) + G_{il,k}(r) \big)}{2} \Big [C'_{ijlk}(x-r)\bar{\varepsilon}_{lk} + \big [C'_{ijlk}(x-r)\varepsilon'_{lk}(x-r)\big ]\Big ]_{,j}dr$$Note that the $G(r)$ terms depend on the reference medium $C^R$.
The equation above can be recast with all of the derivatives on the green's functions by integrating by parts.
$$ \varepsilon'_{kl}(x) = \Bigg [ \int_S \frac{\big (G_{ik,l}(r) + G_{il,k}(r) \big)}{2} \Big [C'_{ijlk}(x-r)\bar{\varepsilon}_{lk} + \big [C'_{ijlk}(x-r)\varepsilon'_{lk}(x-r)\big ]\Big ] n_j dS\Bigg ]_{r \rightarrow 0}^{r \rightarrow \infty} - $$ $$ \int_V \frac{\big (G_{ik,lj}(r) + G_{il,kj}(r) \big)}{2} \Big [C'_{ijlk}(x-r)\bar{\varepsilon}_{lk} + \big [C'_{ijlk}(x-r)\varepsilon'_{lk}(x-r)\big ]\Big ]dr $$In the equation above, the surface term tending to zero is a principal value integral because of the singularity in the green's functions at $r = 0$. As a result, the integrand is not differentiable. Torquato shows that by excluding a sphere at the origin and using integration by parts and the divergence theorem we can arrive at the following equation [1].
$$\varepsilon'_{kl}(x) = I_{ikjl} - E_{ikjl} + \int_V \Phi_{ikjl}(r) \Big [C'_{ijlk}(x-r)\bar{\varepsilon}_{lk} + \big [C'_{ijlk}(x-r)\varepsilon'_{lk}(x-r)\big ]\Big ]dr $$where
$$\Phi_{ikjl}(r) = - \frac{\big (G_{ik,lj}(r) + G_{il,kj}(r) \big)}{2} $$is the green's function terms, and
$$I_{ikjl}^{\infty} = \lim_{r \rightarrow \infty} \int_S\frac{\big (G_{ik,l}(r) + G_{il,k}(r)\big)}{2} \Big [C'_{ijlk}(x-r)\bar{\varepsilon}_{lk} + \big [C'_{ijlk}(x-r)\varepsilon'_{lk}(x-r)\big ]\Big ]n_l dS $$$$E_{ikjl}(x) = \lim_{r \rightarrow 0} \int_S\frac{\big (G_{ik,l}(r) + G_{il,k}(r)\big)}{2} n_l dS $$are the contribution from the surface integrals at $\infty$ and from the singularity.
Finally let
$$\Gamma_{iklj}(r) = I_{ikjl}^{\infty}\delta(r)-E_{ikjl}\delta(r) + \Phi_{ikjl}(r)$$the strain can then be written in the following form.
$$\varepsilon'_{kl}(x) = \int_V \Gamma_{ikjl}(r) \Big [C'_{ijlk}(x-r)\bar{\varepsilon}_{lk} + \big [C'_{ijlk}(x-r)\varepsilon'_{lk}(x-r)\big ]\Big ]dr $$By recursively inserting $\varepsilon'(x)$ into the LHS of the equation, we get the following series.
$$ \varepsilon'(x) =\int_V \Gamma(r) C'(x-r) \bar{\varepsilon} dr +\int_V \int_V \Big[ \Gamma(r) C'(x-r)\bar{\varepsilon}\Big ]\Big [\Gamma(r') C'(x-r') \bar{\varepsilon}\Big] dr'dr + ...$$As long as
$$\Gamma(r) C'(x)\bar{\varepsilon} << 1$$the series can be truncated after a few terms and still provide resonable accuracy.
Let
$$ C'(x-r) = \int_H h m(h, x-r) dh$$where $m(h, r)$ is the microstructure function which is a probablity density that spans both the local state space $h$ and real space $x$. The expectation of local state variable for the microstructure function is the integral over the local state space $H$ and discribes the expected local state $h$ which is equal to $C'(r)$.
Also let
$$\alpha(h, r) = \Gamma(r) h \bar{\varepsilon} $$$$\alpha(h, h', r, r') = \Gamma(r) h \bar{\varepsilon} \Gamma(r') h' \bar{\varepsilon} $$ $$ etc... $$
where again $h$ is the local state variable.
Plugging these definitations into the Weak Contrast Expansion recasts the series in the following form.
The discrete version of this equation is the MKS localization.
[1] Torquato, S., 1997. Effective stiffness tensor of composite media. I. Exact series expansions. J. Mech. Phys. Solids 45, 1421–1448.
[2] Brent L.Adams, Surya Kalidindi, David T. Fullwood. Microstructure Sensitive Design for Performance Optimization.
[3] David T. Fullwood, Brent L.Adams, Surya Kalidindi. A strong contrast homogenization formulation for multi-phase anisotropic materials.