# Random 2D Composites and the Generalized Method of Schwarz

Authors: Vladimir Mityushev
Journal: Advances in Mathematical Physics
 (2015)

---

## Abstract

Two-phase composites with nonoverlapping inclusions randomly embedded in matrix are investigated. A straightforward approach is applied to estimate the effective properties of random 2D composites. First, deterministic boundary value problems are solved for all locations of inclusions, that is, for all events of the considered probabilistic spaceC by the generalized method of Schwarz. Second, the effective properties are calculated in analytical form and averaged over C. This method is related to the traditional method based on the average probabilistic values involving the n-point correlation functions. However, we avoid computation of the correlation functions and compute their weighted moments of high orders by an indirect method which does not address the correlation functions. The effective properties are exactly expressed through these moments. It is proved that the generalized method of Schwarz converges for an arbitrary multiply connected doubly periodic domain and for an arbitrary contrast parameter. The proposed method yields an algorithm which can be applied with symbolic computations. The Torquato-Milton parameter ζ
1 is exactly written for circular inclusions.


---

## Body

## 1. Introduction

The effective properties of random heterogeneous materials and methods of their computation are of considerable interest [1–9]. Randomness in such problems is revealed through random tensor-functions locally describing the physical properties of medium. Despite the considerable progress made in the theory of disordered media, the main tool for studying such systems remains numerical simulations. Frequently, it is just asserted that it is impossible to get general analytical formulae for the effective properties except dilute composites when interactions among inclusions are neglected and except regular composites having simple geometric structures. This opinion sustained by unlimited belief in numerics has to be questioned by the recent pure mathematical investigations devoted to explicit solution to the Riemann-Hilbert problem for multiply connected domains [10] and by significant progress in symbolic computations [11–13]. In the present paper, we demonstrate that the theoretical results [10] can be effectively implemented in symbolic form that yields analytical formulae for random composites. First applications of the theory were performed in [14, 15].For conductivity problems governed by Laplace’s equation, the local conductivity tensorλ
(
x
) can be considered as a random function of the spatial variable x
=
(
x
1
,
x
2
,
x
3
). In the present paper, we restrict ourselves to two-phase composites with nonoverlapping inclusions when a collection of particles with fixed shapes and sizes is embedded in matrix (see Figure 1). More precisely, let a set of hard particles {
D
k
,
k
=
1,2
,
…
} be given, where each D
k has a fixed geometry. Let all the particles be randomly located in the space and each particle D
k occupies a domain D
k without deformations. Thus, the deterministic elements D
k are introduced independently but the joint set is introduced {
D
1
,
D
2
,
…
} randomly. The diversity of random locations is expressed by joint probabilistic distributions of the nonoverlapping domains D
k.Figure 1
Random double periodic location of the nonoverlapping domainsD
k (gray disks) bounded by L
k; D is the complement of all the closed domains D
k
∪
L
k to the periodicity cell Q (dashed parallelogram).Various methods were applied to such random composites. The most known theoretical approach is based on then-point correlation functions presented by Torquato [9] and summarized below in Section 3. This theory yields analytic formulae and bounds for the effective properties. The main shortage of this theory is computational difficulties to calculate the n-point correlation functions for large n.Let a distribution of inclusions be fixed. Then, the effective conductivity tensorΛ
e is uniquely determined through D
k and their physical properties. In order to estimate Λ
e, various statistical approaches were applied to construct the representative volume element (RVE) [16]. These approaches are based on the straightforward computations of the effective properties for various locations of inclusions and on the statistical investigations of the obtained numerical results. The main lack of this physical theory is that the notation of the RVE is not correctly defined and the results do not have precise meaning. This theory is rather a statistical analysis of the computational experiments and physical measures.The main purpose of this paper is to work out constructive analytical-numerical methods to calculate the effective constants of composites with nonoverlapping inclusions and to develop the corresponding RVE theory. We apply the straightforward approach to estimate the effective conductivity of random 2D composites. First, deterministic boundary value problems are solved for all locations of inclusions, that is, for all events in the considered probabilistic spaceC by the generalized method of Schwarz (GMS). This method was proposed by Mikhlin [17] as a generalization of the classical alternating method to multiply connected domains; convergence of the GMS was established in [10, 18]. The GMS is referred to decomposition methods [19] frequently used in numerical computations and realized in the form of alternating methods. A method of integral equations closely related to the GMS was developed in [20].After solution to the deterministic problem by the GMS the ensemble average for the effective conductivity (mathematical expectation) is calculated in accordance with the given distribution. This method is related to the traditional method [9] based on the average probabilistic values involving the n-point correlation functions. In our approach following the GMS, we avoid computation of the correlation functions and compute their weighted moments. The effective properties are exactly written through these moments. Analytical formulae and numerical simulations demonstrate advantages of our approach in the case of nonoverlapping inclusions. For instance, for circular inclusions, our method yields formulae of order O
(
ν
20
) in concentration ν (exactly written in Section 4.6 up to O
(
ν
8
) and up to O
(
ν
20
) [21]). Constructive formulae of the traditional approach [9] are deduced only up to O
(
ν
5
) for nonhigh contrast parameters. In particular, the Torquato-Milton parameter ζ
1 [6, 9] is exactly written for circular inclusions.Stochastic 2D problems are posed and solved in doubly periodic statement in the plane. Theoretically, doubly periodic problems constitute the special class of problems in the plane with infinite number of inclusions [22]. However, the number of inclusions per a periodicity cell, N, is arbitrary taken, and the final formulae for the effective tensor contain N in symbolic form. Similar nonperiodic statements can be used following [23]. Such an approach explains, for instance, the divergent sum (integral) arising in applications of self-consistent methods when the divergent sum does really diverge for some distributions. This implies that the corresponding composites cannot be homogenised over the whole plane even if the concentration is properly defined. The homogenization theory (see [8] and works cited therein) justifies existence of the effective properties for statistically homogeneous random fields which constitute a subclass of heterogeneous fields discussed in [9, 22, 23]. The divergence and other similar effects do not arise for doubly periodic composites when homogenization is always valid. A periodicity cell can be treated as a representative cell and vice versa. From practical point of view, each sample is finite; hence it can be considered as a representative cell with many inclusions. Application of the new RVE theory (see the next paragraph) can reduce the number of inclusions per periodicity cell but not always significantly, because reduction to a small number of inclusions per cell can distort the effective properties of the original sample. It was noted in [24] that periodic arrange of inclusions decreases the effective conductivity when conductivity of inclusions is greater than conductivity of host. Moreover, the theory of elliptic functions can be applied to constructive symbolic computations for periodic problems which yields exact and approximate analytical formulae for the effective tensor.The proposed method leads to construction of the rigorous RVE theory for nonoverlapping inclusions [25, 26]. The effective tensor can be written in the form of expansion in the moments of the correlation functions which can be considered as “basic elements” depending only on locations of inclusions. The RVE is defined as the minimal size periodicity cell corresponding to the set of basic elements calculated for the composite. A simple fast algorithm to determine the representative cell for a given composite is based on reduction of the number of inclusions per periodicity cell having the same basic elements as the given composite. It follows from simulations [11] that for the uniform nonoverlapping distribution of disks in order to reach high accuracy in the effective conductivity one needs to solve the corresponding boundary value problems at least for 64 inclusions per cell repeated at least 1500 times. These technical parameters essentially exceed possibilities of the traditional statistical approach to the RVE.Here, we are restricted by composites with nonoverlapping inclusions. Other types of composites and porous media, for instance, whenλ
(
x
) obeys the Gaussian distribution, the correlation functions, and reconstructions of media, give excellent practical results (see [1] and works cited therein).The goal of this paper is to show that the results obtained by the GMS [18, 24] can be considered as a constructive approach to resolve some problems of the classic theory of random composites [9] for nonoverlapping inclusions. The GMS allows us to overcome myths frequently met in the theory of random composites:(1)
It is asserted that only application of the correlation functions is the proper method to treat random composites. Hence, restrictions in this field are related to computational difficulties arising in numerical computations of the correlation functions. As it is noted above, we obviate the need for the correlation functions calculating the mathematical expectation of the effective tensor after complete solution to the deterministic problems.
(2)
It is thought that the contrast parameter expansions are valid only for sufficiently small contrast parameters. However, it was first justified in [27] that such expansions for plane conductivity problems are valid for all possible values of the contrast parameters.
(3)
Self-consistent methods for random composites are considered as a simple and constructive alternative to the method of correlation functions. However, applications of the GMS demonstrate essential limitations of these methods [23]. More precisely, the self-consistent methods do not capture interactions among inclusions. Even solution to the n-particle problem in the infinite plane can be applied only to diluted distributions of the clusters consisting of n inclusions.This paper is the first part of the work devoted to the GMS described in general form in Section4.1. The special attention is paid to conductivity problems. Integral equations corresponding to the GMS are described in Section 4.2. These equations can be solved by different iterative schemes. The first scheme is based on the contrast expansions (the traditional approach is in Section 3.2 and the GMS form in Section 4.4). The second method uses the cluster expansions (the traditional approach is in Section 3.1 and the GMS form is in Section 4.5). Applications of the GMS to elasticity problems will be described in a separate paper.


## 2. Deterministic and Stochastic Approaches

Letω
1 and ω
2 be the fundamental pair of periods on the complex plane C such that ω
1
>
0 and Im
⁡
ω
2
>
0, where Im stands for the imaginary part. The fundamental parallelogram Q is defined by the vertices ±
ω
1
/
2 and ±
ω
2
/
2. Without loss of generality the area of Q can be normalized to one; hence, (1)
ω
1
Im
⁡
ω
2
=
1
.The points m
1
ω
1
+
m
2
ω
2 (m
1
,
m
2
∈
Z) generate a doubly periodic lattice Q (torus), where Z stands for the set of integer numbers.ConsiderN nonoverlapping simply connected domains D
k in the cell Q with Lyapunov’s boundaries L
k and the multiply connected domain D
=
Q
∖
⋃
k
=
1
N
(
D
k
∪
L
k
), the complement of all the closures of D
k to Q. Each curve L
k leaves D
k to the left (see Figure 1). Let z denote a complex variable and a point a
k is arbitrary fixed in D
k (k
=
1,2
,
…
,
N).We study conductivity of the doubly periodic composite when the hostD
+
m
1
ω
1
+
m
2
ω
2 and the inclusions D
k
+
m
1
ω
1
+
m
2
ω
2 are occupied by conducting materials. Without loss of generality the conductivity of matrix can be normalized to one. The respective conductivity of inclusions is denoted by λ. Introduce the local conductivity as the function (2)
λ
x
=
1
,
x
∈
D
,
λ
,
x
∈
D
k
,
k
=
1,2
,
…
,
N
.Here, x
=
(
x
1
,
x
2
) is related to the above introduced complex variable z by formula z
=
x
1
+
i
x
2, where i denotes the imaginary unit. The potentials u
(
x
) and u
k
(
x
) are harmonic in D and D
k (k
=
1,2
,
…
,
N) and continuously differentiable in closures of the considered domains. The conjugation conditions are fulfilled on the contact interface (3)
u
=
u
k
,
∂
u
∂
n
=
λ
∂
u
k
∂
n
on
L
k
,
k
=
1,2
,
…
,
N
,where ∂
/
∂
n denotes the outward normal derivative to L
k. The external field is modelled by the quasi periodicity conditions (4)
u
z
+
ω
1
=
u
z
+
ξ
1
,
u
z
+
ω
2
=
u
z
+
ξ
2
,where ξ
l (l
=
1,2) are constants. The external field applied to the considered doubly periodic composites is given by the vector E
0
=
-
(
ξ
1
,
ξ
2
).In the book [9], problem (3)-(4) is written in terms of the electric conduction by use of the local electric current J and the intensity field E. The potential u satisfies (5)
E
=
-
∇
u
.The linear constitutive relation is fulfilled in D and D
k: (6)
J
x
=
λ
x
E
x
.Problem (3)-(4) can be stated in the probabilistic context with randomly located inclusions. For definiteness, consider a set of the domains D
k as a realization of the random locations of a collection of particles D
k with fixed shapes and sizes. It corresponds to location of hard particles in a vessel. Following Definition (
3.4
) from [9] we introduce the doubly periodic specific probability density P
(
a
) associated with finding a configuration of inclusions with position a
:
=
(
a
1
,
a
2
,
…
,
a
N
). The point a
k determines the position of the particle D
k in the cell Q. The set of all configurations is denoted by C. Let d
a
k denote a small area element about point a
k and d
a
=
d
a
1
d
a
2
⋯
d
a
N. Then P
(
a
)
d
a is equal to the probability of finding a
1 in d
a
1, a
2 in d
a
2
,
…
,
a
N in d
a
N.The ensemble average is introduced as the expectation in a probabilistic space(7)
f
^
=
∫
C
f
a
P
a
d
a
.The effective conductivity tensor of deterministic composites (8)
Λ
e
=
λ
11
λ
12
λ
21
λ
22can be found in terms of the solutions of two independent problems (3)-(4). For instance, one can take ξ
1
=
ω
1, ξ
2
=
R
e
ω
2 and ξ
1
=
0, ξ
2
=
I
m
ω
2. In the first case, (9)
λ
11
=
1
+
2
ρ
∑
k
=
1
N
∫
D
k
∂
u
k
∂
x
1
d
x
1
d
x
2
.Here, ρ denotes the contrast parameter introduced by Bergman: (10)
ρ
=
λ
-
1
λ
+
1
.The effective conductivity tensorΛ
e of stochastic composites is determined by application of the operator (7) to the deterministic values (8). Let the probabilistic distribution of the domains D
k (k
=
1,2
,
…
,
N) in Q be such that the corresponding two-dimensional composite is isotropic in macroscale; that is, the expected effective conductivity of the composite is expressed by a scalar. It can be calculated as follows: (11)
λ
e
=
λ
11
^
=
∫
C
λ
11
a
P
a
d
a
,where the deterministic value λ
11
(
a
) is calculated by (9).


## 3. Traditional Approach

As it is stressed in Introduction we apply the straightforward approach to determineλ
e. First, the deterministic problem (3)-(4) is solved for every configuration a
∈
C by the GMS. Further, the probabilistic average (11) is calculated. In order to compare the results obtained by this method and by others we follow the traditional approach to random composites described and summarized by Torquato [9]. We follow the book [9] and do not present formulae from earlier classic works (see references and historical notes therein and also in [6]) in order to avoid multiple notations. We stress that only two-dimensional two-phase dispersed composites with one phase embedded in the connected matrix are considered in the framework of the general theory. Torquato [9] and others distinguish two main methods: cluster expansion and contrast expansion discussed below.


### 3.1. Cluster Expansions

This method is presented as a correction of the single-inclusion approach by taking into account interactions between pairs of particle, triplets, and so forth. It is worth noting that only infinite number of inclusions on plane gives a correct result [9, 23]. Study of the finite number of inclusions on plane yields the effective properties only of dilute composites as it is demonstrated in [23].The intensity fieldE is considered as a function of x with the configuration parameter a and denoted by E
(
x
;
a
). Let the constant external field E
0 be fixed. Torquato describes the following cluster expansion (see (
19.8
) from [9]): (12)
E
x
;
a
=
E
0
+
∑
k
=
1
N
M
1
x
;
k
·
E
0
+
∑
k
=
1
N
∑
k
1
≠
k
N
M
2
x
;
k
,
k
1
·
E
0
+
⋯
.Here, M
1
(
x
;
k
) is the single-inclusion operator which accounts for the first order interactions over and above E
0 and inclusions periodically equivalent to the kth inclusion, that is, to D
k
+
m
1
ω
1
+
m
2
ω
2 for all integers m
1 and m
2. The operator M
1
(
x
;
k
) can be derived by the following two methods. The first method corresponds to the Maxwell approach applied to dilute composites. It consists of two steps. At the first step, a boundary value problem for one-inclusion D
k in the infinite plane is solved and the single-inclusion operator K
1
(
x
;
k
) in the plane is constructed (see, e.g., formulae (
17.4
) and (
19.6
) from [9]). Next, M
1
(
x
;
k
) is constructed via K
1
(
x
;
k
) by a periodicity operator. The periodicity transformation P
:
K
1
(
x
;
k
)
↦
M
1
(
x
;
k
) can be easily performed since K
1
(
x
;
a
k
) is expressed via the dipole tensor. The periodicity operator transforms rational functions in z
=
x
1
+
i
x
2 to elliptic functions [28]; for instance, P
:
z
-
1
↦
ζ
(
z
) and P
:
z
-
2
↦
℘
(
z
). The second method to construct M
1
(
x
;
k
) is based on the same scheme realized in the torus topology.The operatorK
1
(
x
;
k
), hence M
1
(
x
;
k
), is written in closed form for disk, ellipse, and many other domains (see [9] and works cited therein). It can be constructed for any shape D
k in terms of the conformal mapping of D
k onto the unit disk. In general, K
p
(
x
;
k
,
k
1
,
…
,
k
p
) and M
p
(
x
;
k
,
k
1
,
…
,
k
p
) are the p-inclusion operator which accounts for the pth order interactions.The effective conductivity can be found through the polarization field defined by(
19.3
)-(
19.4
) from [9]: (13)
P
x
=
0
,
z
∈
D
,
λ
-
1
E
x
,
z
∈
D
k
,
k
=
1,2
,
…
,
N
.Following the traditional approach [9] we introduce the “double” average operator over the unit cell Q and over the configuration space C: (14)
g
=
∫
Q
∫
C
g
x
;
a
d
a
d
x
,where d
x
=
d
x
1
d
x
2 and d
a
=
d
a
1
d
a
2
⋯
d
a
N. The effective tensor (8) can be defined by (15)
P
=
Λ
e
-
I
·
E
,where I is the unit tensor. Having used (15) Torquato (see (
19.33
) from [9] and earlier investigations cited therein) deduced formula for macroscopically isotropic composites (16)
λ
e
=
1
+
∑
n
=
1
∞
1
n
!
∫
C
W
n
a
d
a
,where W
n
(
a
) is a complicated functional of the n-inclusion cluster operators M
n
(
x
;
k
,
k
1
,
…
,
k
n
) and the probability density P
(
a
). It is worth noting that W
n
(
a
) is derived through K
n
(
x
;
k
,
k
1
,
…
,
k
n
) in [9]. Actually, our periodic approach yields the same result by rearrangement of the terms in each W
n
(
a
).LetD
k denote the area of D
k. The concentration of inclusions is defined in the following way: (17)
ν
=
∑
k
=
1
N
D
k
.It is noted in [9] that (16) yields the volume-fraction expansion (see (
19.34
) from [9]): (18)
λ
e
=
1
+
∑
n
=
1
∞
B
n
ν
n
,where the coefficients B
n are multidimensional integrals on M
p
(
x
;
k
,
k
1
,
…
,
k
p
) for p
=
1,2
,
…
,
n. Direct applications of the correlation functions yield analytical formulae and numerical results for B
n, n
=
2,3
,
4 (see [9] and a discussion in [11]).


### 3.2. Contrast Expansions

Exact contrast expansions can be presented by formula(
20.1
) from [9] which after some modifications becomes (19)
λ
e
=
1
+
∑
n
=
1
∞
C
n
ρ
n
,where ρ is the contrast parameter (10). Many authors have been thinking that the series (19) converges only for sufficiently small ρ, though convergence of the series for local fields had been proved in [27] for circular inclusions and in [10] (see Section 4.9
.
2) for general shapes for all admissible |
ρ
|
≤
1. This of course implies the convergence of (19).In order to deduce the general representation (19) and estimate the first coefficients C
n (n
=
1,2
,
3,4) Torquato [9] derives an integral equation on E
(
x
) which can be written in the form (compare to (
20.17
) from [9] for the cavity intensity field) (20)
E
x
=
E
0
+
∫
Q
d
x
′
G
x
,
x
′
·
λ
x
-
1
E
x
′
,where G
(
x
,
x
′
) denotes the periodic Green’s function. Substitution of (2) into (20) yields (21)
E
x
=
E
0
+
ρ
λ
+
1
∑
k
=
1
N
∫
D
k
d
x
′
G
x
,
x
′
·
E
x
′
.The method of successive approximations applied to (21) yields a series on the contrast parameter ρ which can be shortly written in the form (22)
E
x
=
∑
n
=
0
∞
ρ
n
H
n
·
E
0
,where the operator H is defined by the right hand part of (21). Use of (13) yields an analogous series for P
(
x
). The tensor Λ
e satisfies (15) which can be written in the form (20.37) from [9]. The latter formula for macroscopically isotropic media becomes (20.57) from [9]. In our designations it reads as follows: (23)
ρ
2
ν
2
λ
e
-
1
-
1
λ
e
+
1
=
ν
ρ
-
∑
n
=
3
∞
A
n
1
ρ
n
.The n-point coefficients A
n
(
1
) are exactly expressed in terms of the p-point correlation functions (p
=
2,3
,
…
,
n) (see formulae (20.38)–(20.41) and (20.59)-(20.60) from [9]). Formula (23) can be written in the form (19).The coefficientsA
n
(
1
) were calculated in [9] as sums of multiple integrals containing the correlation functions which are presented also through multiple integrals. Calculation of these integrals is the main difficulty for numerical application of (23).


## 3.1. Cluster Expansions

This method is presented as a correction of the single-inclusion approach by taking into account interactions between pairs of particle, triplets, and so forth. It is worth noting that only infinite number of inclusions on plane gives a correct result [9, 23]. Study of the finite number of inclusions on plane yields the effective properties only of dilute composites as it is demonstrated in [23].The intensity fieldE is considered as a function of x with the configuration parameter a and denoted by E
(
x
;
a
). Let the constant external field E
0 be fixed. Torquato describes the following cluster expansion (see (
19.8
) from [9]): (12)
E
x
;
a
=
E
0
+
∑
k
=
1
N
M
1
x
;
k
·
E
0
+
∑
k
=
1
N
∑
k
1
≠
k
N
M
2
x
;
k
,
k
1
·
E
0
+
⋯
.Here, M
1
(
x
;
k
) is the single-inclusion operator which accounts for the first order interactions over and above E
0 and inclusions periodically equivalent to the kth inclusion, that is, to D
k
+
m
1
ω
1
+
m
2
ω
2 for all integers m
1 and m
2. The operator M
1
(
x
;
k
) can be derived by the following two methods. The first method corresponds to the Maxwell approach applied to dilute composites. It consists of two steps. At the first step, a boundary value problem for one-inclusion D
k in the infinite plane is solved and the single-inclusion operator K
1
(
x
;
k
) in the plane is constructed (see, e.g., formulae (
17.4
) and (
19.6
) from [9]). Next, M
1
(
x
;
k
) is constructed via K
1
(
x
;
k
) by a periodicity operator. The periodicity transformation P
:
K
1
(
x
;
k
)
↦
M
1
(
x
;
k
) can be easily performed since K
1
(
x
;
a
k
) is expressed via the dipole tensor. The periodicity operator transforms rational functions in z
=
x
1
+
i
x
2 to elliptic functions [28]; for instance, P
:
z
-
1
↦
ζ
(
z
) and P
:
z
-
2
↦
℘
(
z
). The second method to construct M
1
(
x
;
k
) is based on the same scheme realized in the torus topology.The operatorK
1
(
x
;
k
), hence M
1
(
x
;
k
), is written in closed form for disk, ellipse, and many other domains (see [9] and works cited therein). It can be constructed for any shape D
k in terms of the conformal mapping of D
k onto the unit disk. In general, K
p
(
x
;
k
,
k
1
,
…
,
k
p
) and M
p
(
x
;
k
,
k
1
,
…
,
k
p
) are the p-inclusion operator which accounts for the pth order interactions.The effective conductivity can be found through the polarization field defined by(
19.3
)-(
19.4
) from [9]: (13)
P
x
=
0
,
z
∈
D
,
λ
-
1
E
x
,
z
∈
D
k
,
k
=
1,2
,
…
,
N
.Following the traditional approach [9] we introduce the “double” average operator over the unit cell Q and over the configuration space C: (14)
g
=
∫
Q
∫
C
g
x
;
a
d
a
d
x
,where d
x
=
d
x
1
d
x
2 and d
a
=
d
a
1
d
a
2
⋯
d
a
N. The effective tensor (8) can be defined by (15)
P
=
Λ
e
-
I
·
E
,where I is the unit tensor. Having used (15) Torquato (see (
19.33
) from [9] and earlier investigations cited therein) deduced formula for macroscopically isotropic composites (16)
λ
e
=
1
+
∑
n
=
1
∞
1
n
!
∫
C
W
n
a
d
a
,where W
n
(
a
) is a complicated functional of the n-inclusion cluster operators M
n
(
x
;
k
,
k
1
,
…
,
k
n
) and the probability density P
(
a
). It is worth noting that W
n
(
a
) is derived through K
n
(
x
;
k
,
k
1
,
…
,
k
n
) in [9]. Actually, our periodic approach yields the same result by rearrangement of the terms in each W
n
(
a
).LetD
k denote the area of D
k. The concentration of inclusions is defined in the following way: (17)
ν
=
∑
k
=
1
N
D
k
.It is noted in [9] that (16) yields the volume-fraction expansion (see (
19.34
) from [9]): (18)
λ
e
=
1
+
∑
n
=
1
∞
B
n
ν
n
,where the coefficients B
n are multidimensional integrals on M
p
(
x
;
k
,
k
1
,
…
,
k
p
) for p
=
1,2
,
…
,
n. Direct applications of the correlation functions yield analytical formulae and numerical results for B
n, n
=
2,3
,
4 (see [9] and a discussion in [11]).


## 3.2. Contrast Expansions

Exact contrast expansions can be presented by formula(
20.1
) from [9] which after some modifications becomes (19)
λ
e
=
1
+
∑
n
=
1
∞
C
n
ρ
n
,where ρ is the contrast parameter (10). Many authors have been thinking that the series (19) converges only for sufficiently small ρ, though convergence of the series for local fields had been proved in [27] for circular inclusions and in [10] (see Section 4.9
.
2) for general shapes for all admissible |
ρ
|
≤
1. This of course implies the convergence of (19).In order to deduce the general representation (19) and estimate the first coefficients C
n (n
=
1,2
,
3,4) Torquato [9] derives an integral equation on E
(
x
) which can be written in the form (compare to (
20.17
) from [9] for the cavity intensity field) (20)
E
x
=
E
0
+
∫
Q
d
x
′
G
x
,
x
′
·
λ
x
-
1
E
x
′
,where G
(
x
,
x
′
) denotes the periodic Green’s function. Substitution of (2) into (20) yields (21)
E
x
=
E
0
+
ρ
λ
+
1
∑
k
=
1
N
∫
D
k
d
x
′
G
x
,
x
′
·
E
x
′
.The method of successive approximations applied to (21) yields a series on the contrast parameter ρ which can be shortly written in the form (22)
E
x
=
∑
n
=
0
∞
ρ
n
H
n
·
E
0
,where the operator H is defined by the right hand part of (21). Use of (13) yields an analogous series for P
(
x
). The tensor Λ
e satisfies (15) which can be written in the form (20.37) from [9]. The latter formula for macroscopically isotropic media becomes (20.57) from [9]. In our designations it reads as follows: (23)
ρ
2
ν
2
λ
e
-
1
-
1
λ
e
+
1
=
ν
ρ
-
∑
n
=
3
∞
A
n
1
ρ
n
.The n-point coefficients A
n
(
1
) are exactly expressed in terms of the p-point correlation functions (p
=
2,3
,
…
,
n) (see formulae (20.38)–(20.41) and (20.59)-(20.60) from [9]). Formula (23) can be written in the form (19).The coefficientsA
n
(
1
) were calculated in [9] as sums of multiple integrals containing the correlation functions which are presented also through multiple integrals. Calculation of these integrals is the main difficulty for numerical application of (23).


## 4. The Generalized Method of Schwarz

### 4.1. General

The generalized method of Schwarz (GMS) was proposed by Mikhlin [17] as a generalization of the classical alternating method of Schwarz for finitely connected domains. The GMS is based on the decomposition of the considered domain with complex geometry onto simple domains and subsequent solution to boundary value problems for simple domains. In our case, we have N simple domains, the cell Q with one-inclusion D
k (k
=
1,2
,
…
,
N). In each step of the algorithm the boundary conditions for the kth simple domain are corrected by influence of the mth simple domains (m
≠
k) computed at the precedent steps.We now shortly present the GMS, first for the conductivity problem (3)-(4) and further for general linear problems. The ideal contact condition (3) can be rewritten in the form (24)
u
0
=
u
k
-
u
e
x
t
,
∂
u
0
∂
n
=
∂
u
k
-
u
e
x
t
∂
n
+
ρ
λ
+
1
∂
u
k
∂
n
on
L
k
,
k
=
1,2
,
…
,
N
,where u
=
u
0
+
u
e
x
t is the decomposition of the potential u onto the regular periodic part u
0 and the external field u
e
x
t corresponding to the quasi periodicity conditions (4). In the case of the space problem the external field u
e
x
t has a singularity at infinity.Introduce the function(25)
f
t
=
λ
+
1
∂
u
k
∂
n
t
,
t
∈
L
k
k
=
1,2
,
…
,
Nand the domains D
+
:
=
⋃
k
=
1
n
D
k, D
-
:
=
D. It is convenient to represent the harmonic functions in different domains as one piecewise harmonic function u
0
(
x
) introduced above in D
=
D
- and equal to u
k
(
x
)
-
u
e
x
t
(
x
) in D
k
∪
L
k (k
=
1,2
,
…
,
N). Then, (24) can be considered as the jump problem: (26)
u
0
+
=
u
0
-
,
∂
u
0
+
∂
n
=
∂
u
0
-
∂
n
+
ρ
f
t
on
L
=
⋃
k
=
1
n
L
k
,where u
0
+ and u
0
- correspond to the limit values of u
(
x
) when x tends to t from D
+ and D
-, respectively. The jump problem (26) has a unique solution expressed in terms of the simple layer potential P for the curve L in the cell Q in the torus topology. We have (27)
u
0
z
=
ρ
P
f
x
,
x
∈
D
+
∪
D
-
.The operator P is considered as an operator in an appropriate functional space [29–32]. In the classic statement, Q should be replaced by the whole space. The operator P is decomposed onto the simple layer potentials P
k along curves L
k (k
=
1,2
,
…
,
N) as P
=
(
1
/
(
λ
+
1
)
)
∑
k
=
1
N
P
k. Here, the multiplier 1
/
(
λ
+
1
) is introduced for short form of the equations below. Then (27) implies that (28)
u
k
x
=
ρ
∑
m
=
1
N
P
m
∂
u
m
+
∂
n
x
+
u
e
x
t
x
,
x
∈
D
k
k
=
1,2
,
…
,
N
,
(29)
u
x
=
ρ
∑
m
=
1
N
P
m
∂
u
m
+
∂
n
x
+
u
e
x
t
x
,
x
∈
D
.Equations (28) can be considered as a system of integral equations on the potentials u
k
(
x
) in the inclusions D
k
(
k
=
1,2
,
…
,
N
). If it is solved, the potential u
(
x
) is calculated in D by (29).Remark 1.
Equations (28) correspond to the GMS in Mikhlin’s form [17] when the method of successive approximations can diverge for closely spaced inclusions. The following slight modification yields a convergent algorithm for any location of inclusions and every contrast parameter (|
ρ
|
≤
1): (30)
u
k
x
=
ρ
∑
m
=
1
N
P
m
∂
u
m
+
∂
n
x
-
P
m
∂
u
m
+
∂
n
w
+
u
e
x
t
x
+
ρ
∑
m
=
1
N
P
m
∂
u
m
+
∂
n
w
,
x
∈
D
k
∪
L
k
k
=
1,2
,
…
,
N
,where w is a fixed point in D.Equations (27)–(30) are written for potentials. Their differentiation yields equations on the intensity field: (31)
E
k
x
=
ρ
∑
m
=
1
N
Q
m
E
m
x
+
E
0
z
,
x
∈
D
k
∪
L
k
k
=
1,2
,
…
,
N
,where E
k
(
x
)
=
E
(
x
) in D
k and Q
m are appropriate operators (similar to double layer potentials). Equation (31) is similar to (20).This GMS for harmonic functions can be extended to general problems governed by linear equations describing electrical and heat conduction, flow in porous media, viscous flow, elasticity, and coupled fields inR
d (d
=
2,3).Let tensor fieldsU
(
z
) in D and U
k
(
z
) in D
k satisfy the contact conditions (32)
T
·
U
=
T
k
·
U
k
,
on
⁡
L
k
k
=
1,2
,
…
,
N
,where T and T
k are boundary operators involving physical parameters of materials occupying the domains D and D
k, respectively. Let T
k be presented in the form T
k
=
T
-
ρ
k, where ρ
k denotes a contrast parameter tensor. Then (32) becomes (33)
U
=
U
k
-
T
-
1
·
ρ
k
·
U
k
,
on
L
k
k
=
1,2
,
…
,
N
.Let P
k denote the simple layer potential corresponding to the considered linear equation for the domains D
k and D
k
- separated by L
k. Application of ∑
k
=
1
N
P
k to (33) yields (34)
U
k
x
=
∑
m
=
1
N
P
m
·
T
-
1
·
ρ
m
·
U
m
x
+
U
e
x
t
x
,
x
∈
D
k
∪
L
k
k
=
1,2
,
…
,
N
.There are two different methods to solve (34). The first method is based on the direct iterations and corresponds to the contrast expansion in ρ (compare to Section 3.2). The second method is based on implicit iterations applied to the same equations (34) written in the form (35)
U
k
x
-
P
k
·
T
-
1
·
ρ
k
·
U
k
x
=
∑
m
≠
k
N
P
m
·
T
-
1
·
ρ
m
·
U
m
x
+
U
e
x
t
x
,
x
∈
D
k
∪
L
k
k
=
1,2
,
…
,
N
.One-inclusion problem forD
k
(
k
=
1,2
,
…
,
N
) is solved at each step of iterations. This scheme corresponds to the cluster expansions discussed in Section 3.1.Remark 2.
Equations (27) and (35) are not reduced to the classic integral equations constructed in the framework of the potentials theory [32]. It is another type of equations.Remark 3.
Equations (27)–(29) and (34)-(35) can be derived in the limit cases of inclusions (soft and hard inclusions in terminology [33]) by introduction of fictive potentials [10].


### 4.2. Integral Equations

In the present section, we discuss the GMS for the double periodic problem (3)-(4). Following [34] we introduce the complex potentials ϕ
(
z
) and ϕ
k
(
z
) in such a way that (36)
u
x
=
Re
⁡
ϕ
z
,
u
k
x
=
2
λ
+
1
Re
⁡
ϕ
k
z
,where x
=
(
x
1
,
x
2
) and z
=
x
1
+
i
x
2. Two real equations (3) can be written as one R-linear condition: (37)
ϕ
t
=
ϕ
k
t
-
ρ
ϕ
k
t
¯
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
,where the bar denotes the complex conjugation. The real part of (37) yields the first equation in (3). The imaginary part yields Im
⁡
ϕ
t
=
1
+
ρ
Im
⁡
ϕ
k
t. After differentiation of the latter equality along L
k and application of the Cauchy-Riemann equations we arrive at the second equation in (3) (see for details [10]).The unknown functionsϕ
z and ϕ
k
z are analytic in D and D
k, respectively, and continuously differentiable in the closures of the domains considered. It follows from (4) that the function ϕ
(
z
) is quasiperiodic: (38)
ϕ
z
+
ω
j
-
ϕ
z
=
ξ
j
+
i
d
j
j
=
1,2
,where the constants ξ
1 and ξ
2 are taken from (4) and d
1 and d
2 are undetermined real constants which should be found during solution to the problem.The following formula was proved in [10]: (39)
d
d
t
ϕ
k
t
=
-
n
2
t
¯
ϕ
k
′
t
¯
,
t
∈
L
k
.Here, the differential operator is understood as d
/
d
t
=
t
s
′
(
d
/
d
s
), where t
=
t
(
s
) is a parametric equation of L
k and ϕ
k
′
(
t
) is the boundary value of the complex derivative ∂
/
∂
z
=
∂
/
∂
x
1
-
i
(
∂
/
∂
x
2
). Differentiation of (37) and application of (39) yield (40)
ψ
t
=
ψ
k
t
+
ρ
n
k
2
t
¯
ψ
k
t
¯
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
,where ψ
(
z
)
=
ϕ
′
(
z
), ψ
k
(
z
)
=
ϕ
k
′
(
z
), and n
k
(
t
) stands for the unit normal outward vector to L
k expressed in terms of the complex values. The function ψ
(
z
) is doubly periodic. The R-linear problem (40) holds also in the limit cases for perfectly conducting inclusions (ρ
=
1) and for isolators (ρ
=
-
1).TheR-linear problem (37)-(38) can be reduced to a system of integral equation in the following way. The cell Q in the torus topology is divided into two domains D
+
=
⋃
k
=
1
n
D
k (not connected) and D
-
=
D (multiply connected). Let L
=
⋃
k
=
1
n
L
k denote the boundary of D
+ and let μ
(
t
) be a Hölder continuous function on L. Let ζ
(
z
) denote the Weierstrass ζ-function for which [28] (41)
ζ
z
+
ω
j
-
ζ
z
=
δ
j
j
=
1,2
,where δ
j
=
2
ζ
(
ω
j
/
2
). Following [35] introduce the Eisenstein function (42)
E
1
z
=
∑
m
1
,
m
2
∈
Z
1
z
+
m
1
ω
1
+
m
2
ω
2
,where the Eisenstein summation is used [35]. The Weierstrass and Eisenstein functions are related by formula E
1
(
z
)
=
ζ
(
z
)
-
S
2
z, where S
2
=
(
2
/
ω
1
)
ζ
(
ω
1
/
2
). It follows from Legendre’s identity δ
1
ω
2
-
δ
2
ω
1
=
2
π
i (see [28, 35]) and (41) that the jumps of E
1
(
z
) have the form (43)
E
1
z
+
ω
1
-
E
1
z
=
0
,
E
1
z
+
ω
2
-
E
1
z
=
-
2
π
i
ω
1
.The Cauchy-type integral on torus was introduced with theζ-function in the kernel [36]. We introduce it in the slight different form using the Eisenstein function (44)
Φ
z
=
1
2
π
i
∫
L
μ
t
E
1
t
-
z
d
t
+
C
z
+
C
0
,
z
∈
D
±
,where μ
(
t
) is a Hölder continuous function on L and C and C
0 are complex constants. The function Φ
(
z
) is analytic in D
± and quasiperiodic; that is, it has constant jumps per periodicity cell calculated with the help of (43): (45)
Φ
z
+
ω
1
-
Φ
z
=
C
ω
1
,
Φ
z
+
ω
2
-
Φ
z
=
-
1
ω
1
∫
L
μ
t
d
t
+
C
ω
2
.The following Sochocki (Sokhotski-Plemelj) formulae take place on torus [36]: (46)
Φ
±
t
=
±
1
2
μ
t
+
1
2
π
i
∫
L
μ
τ
E
1
τ
-
t
d
τ
+
C
t
+
C
0
,
t
∈
L
,where Φ
±
(
t
) denote the limit values of Φ
(
z
) on L when z
∈
D
± tends to t
∈
L, respectively. The principal value of the integral in (46) is considered following [37]. Formulae (46) imply that the function Φ
(
z
) satisfies the jump condition (47)
Φ
+
t
-
Φ
-
t
=
μ
t
,
t
∈
L
.Equation (47) can be considered as a C-linear conjugation problem [10] in a class of quasiperiodic functions with the given jump μ
(
t
). It follows from [36] that the general solution of (47) up to an additive arbitrary constant has the form (44).Consider (37) as the problem (47) with μ
(
t
)
=
ρ
ϕ
k
(
t
)
¯ on L, Φ
(
z
)
=
ϕ
k
(
z
) in D
k
⊂
D
+, and Φ
(
z
)
=
ϕ
(
z
) in D
=
D
-. Then, (44) can be written in D
+ as follows: (48)
ϕ
k
z
=
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
ϕ
m
t
¯
E
1
t
-
z
d
t
+
C
z
+
C
0
,
z
∈
D
k
k
=
1,2
,
…
,
N
.Equation (44) in D
- becomes(49)
ϕ
z
=
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
ϕ
m
t
¯
E
1
t
-
z
d
t
+
C
z
+
C
0
,
z
∈
D
.Equations (48)-(49) are deduced from the R-linear problem (3)-(4). Using the same arguments one can derive analogous equations from the problem (40) on the complex flux: (50)
ψ
k
z
=
-
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
n
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
C
,
z
∈
D
k
k
=
1,2
,
…
,
N
,
(51)
ψ
z
=
-
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
n
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
C
,
z
∈
D
.It is also possible to obtain these equations by differentiation of (48)-(49) using integration by parts and formula (39).The integral equations (50) are considered as an equation in the following Banach space. First, the space C
(
0
,
α
) of functions Hölder continuous on L is considered (0
<
α
≤
1). This is a Banach space endowed with the norm: (52)
ω
=
sup
t
∈
L
⁡
ω
t
+
sup
t
1,2
∈
L
t
1
≠
t
2
⁡
ω
t
1
-
ω
t
2
t
1
-
t
2
α
.Consider the closed subspace A
⊂
C
(
0
,
α
) consisting of functions analytically continued into all D
k. Therefore, every element of A is a function analytic in the domain ⋃
k
=
1
N
D
k and Hölder continuous in its closure. It is worth noting that the domain ⋃
k
=
1
N
D
k is not connected and convergence in A implies the uniform convergence of functions in ⋃
k
=
1
N
(
D
k
∪
L
k
).Equations (50) are written in the domains D
k and can be continued to the boundary by use of (46) (53)
ψ
k
t
=
-
ρ
∑
m
=
1
N
1
2
n
m
2
t
¯
ψ
m
t
¯
+
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
C
,
t
∈
L
k
k
=
1,2
,
…
,
N
.Therefore, the integral equations in the space A have the form (50), (53). It is proved in Section 4.3 that (50), (53) have a unique solution. Solution to these equations corresponds to the GMS which can be realized in different forms as it is noted in Section 4.1.The problem (37)-(38) for the complex potentials has a unique solution up to an arbitrary additive constant C
0. The constants ξ
1
+
i
d
1, ξ
2
+
i
d
2, and C
0 disappear in the problem (40) after differentiation. The structure of the general solution of (40) has the form ψ
(
z
)
=
ξ
1
ψ
(
1
)
(
z
)
+
ξ
2
ψ
(
2
)
(
z
), where the real constants ξ
1 and ξ
2 correspond to the external field (jumps of u
(
x
) per a cell). Hence, the R-linear problem (40) has two R-linear independent solutions. From the other side, these two independent solutions can be constructed by (50)-(51) taking, for instance, C
=
1 and C
=
i. Therefore, the complex constant C can be expressed through the real constants ξ
1 and ξ
2. In order to calculate the effective properties of macroscopically isotropic composites it is sufficient to fix ξ
1, ξ
2 and find C. For simplicity, we put (54)
ξ
1
=
ω
1
,
ξ
2
=
Re
⁡
ω
2
.Then the external field corresponds to the complex potential ϕ
(
e
x
t
)
(
z
)
=
z and to the complex flux ψ
(
e
x
t
)
(
z
)
=
1. In the vector form, this external flux becomes (
-
1,0
). It follows from the first relation (45) that (55)
C
=
1
+
i
d
1
ω
1
,where d
1
=
Im
⁡
[
ϕ
(
z
+
ω
1
)
-
ϕ
(
z
)
]. Put z
=
-
(
1
/
2
)
(
ω
1
+
ω
2
) in this relation and calculate(56)
d
1
=
Im
⁡
∫
-
1
/
2
ω
1
+
ω
2
1
/
2
ω
1
-
ω
2
ψ
t
d
t
=
-
∫
-
1
/
2
ω
1
+
R
e
ω
2
1
/
2
ω
1
-
R
e
ω
2
∂
u
∂
x
2
x
1
-
i
Im
⁡
ω
2
2
d
x
1
.The integral from the left hand part of (56) expresses the total flux through the low side of the parallelogram Q parallel to the x
1-axis. It must be equal to zero because of the periodicity and that the external flux (
-
1,0
) vanishes in the x
2-direction. Therefore, d
1
=
0; hence C
=
1 by (55).The effective conductivity tensor of deterministic composites can be calculated through the complex potentials in the inclusions (see [22, 34]) (57)
λ
11
-
i
λ
12
=
1
+
2
ρ
∑
k
=
1
N
∫
D
k
ψ
k
x
1
+
i
x
2
d
x
1
d
x
2
,where ψ
k
(
z
) (k
=
1,2
,
…
,
N) satisfy (50)-(51), (53) with C
=
1. The value (57) depends on the locations of inclusions, that is, on a
∈
C. The effective conductivity of macroscopically isotropic random composites is defined through the operator (7) applied to (57):(58)
λ
e
=
λ
11
^
,where the macroscopic isotropy implies that λ
12
^ vanishes and λ
11
^
=
λ
22
^.


### 4.3. Convergence of the GMS

This section is devoted to study of the integral equations (50), (53) in the space A.LetΨ
(
z
)
=
ψ
k
(
z
) for z
∈
D
k
∪
L
k and let A denote the operator from the right hand part of (48), (53). Then, (48), (53) can be shortly written as the following equation in the space A: (59)
Ψ
=
ρ
A
Ψ
+
1
.In the present section, we demonstrate that (59) has the unique solution for |
ρ
|
<
1 given by the series (60)
Ψ
=
∑
n
=
0
∞
ρ
n
A
n
1
.We use the following general result of functional analysis.Theorem 4 (see [38]).
LetB be a linear bounded operator in a Banach space B. If for any element f
∈
B and for any complex number ν satisfying the inequality |
ν
|
≤
1 equation (61)
Ψ
=
ν
B
Ψ
+
fhas a unique solution, then the unique solution of (62)
Ψ
=
B
Ψ
+
fcan be found by the method of successive approximations. The approximations converge in B to the solution (63)
Ψ
=
∑
k
=
0
∞
B
k
f
.Theorem4 can be applied to (59) in the Banach space A.Theorem 5.
Let|
ρ
|
<
1. Equation Ψ
=
ρ
A
Ψ
+
1 has a unique solution. This solution can be found by the method of successive approximations convergent in the space A.Proof.
Let|
ν
|
≤
1. Consider (64)
ψ
k
z
=
-
ν
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
f
k
z
,
z
∈
D
k
k
=
1,2
,
…
,
N
,
ψ
k
t
=
-
ν
ρ
1
2
n
k
2
t
¯
ψ
k
t
¯
+
∑
m
=
1
N
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
f
k
t
,
t
∈
L
k
k
=
1,2
,
…
,
N
,where f
k
∈
A
k.
Letψ
k
(
z
) be a solution of (53). Introduce the function ψ
(
z
) analytic in D and Hölder continuous in its closure as follows: (65)
ψ
z
=
-
ν
ρ
∑
m
=
1
n
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
,
z
∈
D
.The expression in the right hand part of (65) can be considered as Cauchy’s integral on torus (66)
Φ
z
=
1
2
π
i
∫
L
μ
t
E
1
t
-
z
d
tfor z
∈
Q
∖
L with the density μ
(
t
)
=
-
ν
ρ
n
m
2
(
t
)
¯
ψ
m
(
t
)
¯ on L
m. Using the property (47) of (66) and (64) on L
k we arrive at the R-linear conjugation relation on each fixed curve L
k: (67)
ψ
k
t
-
f
k
t
-
ψ
t
=
-
ν
ρ
n
k
2
t
¯
ψ
k
t
¯
,
t
∈
L
k
.Here, Φ
+
(
t
)
=
ψ
k
(
t
)
-
f
k
(
t
) and Φ
-
(
t
)
=
ψ
(
t
). Integration of (67) along L
k yields (68)
ϕ
t
=
ϕ
k
t
-
ν
ρ
ϕ
k
t
¯
-
g
k
t
+
γ
k
,
t
∈
L
k
k
=
1,2
,
…
,
N
,where the functions ϕ
(
z
), ϕ
k
(
z
) are analytic in D, D
k and belong to C
(
1
,
α
) in the closures of the domains considered, g
k
(
z
) is a primitive function of f
k
(
z
), and γ
k is an arbitrary constant. Moreover, the function ϕ
(
z
) is quasiperiodic; hence (69)
ϕ
z
=
φ
z
+
c
z
,where φ
(
z
) is periodic and c is a constant. The constant γ
k can be included into φ
k
(
z
) by the representation γ
k
=
α
k
-
ν
ρ
α
k
¯. Then, (68) becomes (70)
φ
t
=
φ
k
t
-
ν
ρ
φ
k
t
¯
-
g
k
t
-
c
t
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
.It follows from [20] and Appendix B of [22] that the general solution of the problem (70) has the following structure: (71)
φ
z
=
φ
1
z
+
φ
2
z
,
φ
k
z
=
φ
k
1
z
+
φ
k
2
z
,where (72)
φ
1
t
=
φ
k
1
t
-
ν
ρ
φ
k
1
t
¯
-
g
k
t
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
,
(73)
φ
2
t
=
φ
k
2
t
-
ν
ρ
φ
k
2
t
¯
-
c
t
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
.Each of the problems (72) and (73) (with fixed c) has a unique solution since |
ν
ρ
|
<
1 (the elliptic case in Mikhailov’s terminology [39]). The uniqueness also follows from [20]. This unique solution of the problem (73) can be easily found as φ
(
2
)
(
z
)
=
c
z and φ
k
(
2
)
(
z
)
=
0. Then, differentiation of (71)–(73) implies that the general solution of the problem (67) has the form (74)
ψ
z
=
ψ
1
z
+
c
,where the function ψ
k
(
1
)
(
z
) is uniquely determined from the problem (75)
ψ
1
t
=
ψ
k
1
t
+
ν
ρ
n
k
2
t
¯
ψ
k
1
t
¯
-
f
k
t
,
t
∈
L
k
k
=
1,2
,
…
,
N
.Comparison of (67) and (75) proves the theorem.The series(76)
Ψ
=
∑
k
=
0
∞
ρ
k
A
k
1converges absolutely for |
ρ
|
<
1, since it can be represented in the form Ψ
=
∑
k
=
0
∞
s
k
ρ
0
k
A
k
1 with s
ρ
0
=
ρ and |
ρ
|
<
|
ρ
0
|
<
1. The latter series converges absolutely with the rate |
s
| because the series ∑
k
=
0
∞
ρ
0
k
A
k
1 converges uniformly; hence, its general term tends to zero as k
→
∞.Uniform convergence of (76) for |
ρ
|
=
1 can be justified by use of the arguments form [10].


### 4.4. Contrast Expansions for the GMS

It is proved in the previous section that the unique solution of (50), (53) can be found by a method of successive approximations uniformly convergent in ⋃
k
=
1
N
(
D
k
∪
L
k
). The method of successive approximations can be also presented in the form of the following iterative scheme: (77)
ψ
k
0
z
=
1
,
ψ
k
p
z
=
-
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
p
-
1
t
¯
E
1
t
-
z
d
t
+
1
,
p
=
1,2
,
…
,
z
∈
D
k
k
=
1,2
,
…
,
N
.When the functions ψ
k
(
z
) (k
=
1,2
,
…
,
N) are found, the effective conductivity is calculated by (57).In order to compare (21) and (48), (53) we introduce the complex intensity field E
(
z
)
=
E
1
(
z
)
-
i
E
2
(
z
) isomorphic to the vector field E
(
x
)
=
(
E
1
(
x
)
,
E
2
(
x
)
). Then, the vector equation (21) can be written as a complex equation on E
(
z
). For definiteness, we consider the complex potential ψ
k
(
z
) in D
k. It follows from (5), (36) and ψ
k
(
z
)
=
ϕ
k
′
(
z
) that (78)
E
z
=
∂
u
k
∂
x
1
-
i
∂
u
k
∂
x
2
=
2
λ
+
1
ψ
k
z
.Therefore, the complex equation on E
(
z
) in D
k can be written as (48) on ψ
k
(
z
). We do not directly prove that (21) and (59) are the same. But it is easy to show that the forms of their solutions (22) and (60) coincide. Since the solutions E and Ψ coincide (more precisely, isomorphic) for E
0
=
(
-
1,0
), the power series (22) and (60) in ρ can be considered as the presentations of the same function for sufficiently small |
ρ
|. Therefore, the coefficients of (22) and (60) coincide. In particular, this fact implies that substitution of (60) into (57) yields a series which coincides with the series (19) obtained by the traditional approach. It follows from [10, 27] that the contrast expansion (60) converges for all admissible |
ρ
|
≤
1.


### 4.5. Cluster Expansions for the GMS

The GMS can be presented by other iterative schemes corresponding to the traditional cluster expansions summarized in Section3.1. First, we rewrite (48) in the form (79)
ψ
k
z
+
ρ
2
π
i
∫
L
k
n
k
2
t
¯
ψ
k
t
¯
E
1
t
-
z
d
t
=
-
ρ
∑
m
≠
k
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
1
,
z
∈
D
k
k
=
1,2
,
…
,
n
.Following Section4.2 for every fixed m
=
1,2
,
…
,
N introduce the Banach space A
m of functions analytic in D
m and Hölder continuous in D
m
∪
L
m. Let f
m
(
z
)
∈
A
m. For each fixed k
=
1,2
,
…
,
N introduce the operators (80)
P
k
m
f
m
z
=
1
2
π
i
∫
L
m
n
m
2
t
¯
f
m
t
¯
E
1
t
-
z
d
t
,
z
∈
D
k
m
=
1,2
,
…
,
n
.One can check that the operator P
k
k is singular and P
k
m (m
≠
k) are compact in A
k. Equation (79) can be shortly written in the form (81)
I
+
ρ
P
k
k
ψ
k
=
-
ρ
∑
m
≠
k
P
k
m
ψ
m
+
1
,
k
=
1,2
,
…
,
n
.The (zero)th approximation in the concentration ν for (79) can be written as the following integral equation: (82)
ψ
k
0
z
+
ρ
2
π
i
∫
L
k
n
k
2
t
¯
ψ
k
0
t
¯
E
1
t
-
z
d
t
=
1
,
z
∈
D
k
.Equation (82) corresponds to the R-linear problem (40) for one-inclusion D
k in the cell Q. Solution to the one-inclusion problem (82) defines the inverse operator (
I
+
ρ
P
k
k
)
-
1 bounded in A
k. Then, (81) can be written in the equivalent form (83)
ψ
k
=
-
ρ
∑
m
≠
k
I
+
ρ
P
k
k
-
1
P
k
m
ψ
m
+
I
+
ρ
P
k
k
-
1
1
,
k
=
1,2
,
…
,
n
.As it is established in Section 4.4 the vector-function E
(
x
) and the complex functions ψ
k
(
z
) express the same vector field in D
k up to an isomorphism for E
0
=
(
-
1,0
). Therefore, application of the successive approximations to (83) yields the same series as (12) in D
k. A detailed comparison can be performed. It shows, for instance, that the term E
0
+
∑
k
=
1
N
M
1
(
x
;
k
)
·
E
0 from (12) corresponds to (
I
+
ρ
P
k
k
)
-
1
1 from (83).The solution of the problem (82) depends on the shape of D
k. Therefore, the first order approximation in ν for the effective conductivity also depends on the shape and can be calculated by (57). If all the inclusions have the same form, the functions ψ
k
(
0
)
(
z
) do not depend on k. Then, we arrive at the formula valid up to O
(
ν
2
): (84)
λ
e
≈
1
+
2
ρ
ν
α
,where the shape factor α is introduced following [23, 40]: (85)
α
=
1
D
1
∫
D
1
ψ
1
0
x
1
+
i
x
2
d
x
1
d
x
2
.The Clausius-Mossotti formula equivalent to (84) up to O
(
ν
2
) is used in applications [6]: (86)
λ
e
≈
1
+
ρ
ν
α
1
-
ρ
ν
α
.One can find limitations of the formulae (84), (86) and a comparison of their accuracy in [23].The higher order iterations in the cluster expansion are based on the integral equations onψ
k
(
p
)
(
z
): (87)
ψ
k
p
z
+
ρ
2
π
i
∫
L
k
n
k
2
t
¯
ψ
k
p
t
¯
E
1
t
-
z
d
t
=
-
ρ
∑
m
≠
k
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
p
-
1
t
¯
E
1
t
-
z
d
t
+
1
,
z
∈
D
k
k
=
1,2
,
…
,
n
;
p
=
1,2
,
…
.The uniform convergence of the iterative scheme (87) was proved in the case ρ
=
1 (perfectly conducting inclusions) in Section 4.9
.
2 of [10]. The higher order iterations (87) yield approximate analytical formulae for λ
e for arbitrary nonoverlapping locations of different disks [24, 34]. An example of equal disks is discussed below in Section 4.6.


### 4.6. Circular Inclusions

Let the inclusionsD
k be disks of radius r centered at the points a
k. Then, the normal unit vector n
k
(
t
) to the circle |
t
-
a
k
|
=
r has the form n
k
(
t
)
=
(
t
-
a
k
)
/
r. The inversion with respect to |
t
-
a
k
|
=
r is defined by formula (88)
z
k
∗
=
r
2
z
-
a
k
¯
+
a
k
.Consider the integral operator from (48): (89)
J
k
m
f
z
=
1
2
π
i
∫
L
m
n
m
2
t
¯
f
t
¯
E
1
t
-
z
d
t
,
z
∈
D
k
,where f
(
z
) is analytic in the disk D
m. We have t
=
t
(
m
)
∗ on L
m. Hence, the operator (89) can be written in the form (90)
J
k
m
f
z
=
1
2
π
i
∫
L
m
r
t
-
a
m
2
f
t
m
∗
¯
E
1
t
-
z
d
t
,
z
∈
D
k
,where the function f
(
t
(
m
)
∗
)
¯ is analytically continued into |
z
-
a
k
|
>
r.The functionE
1
(
z
) has poles of first order at the points z
=
m
1
ω
1
+
m
2
ω
2 (see [35]). The function r
/
(
z
-
a
m
)
2
f
(
z
(
m
)
∗
)
¯ is analytic in the domain |
z
-
a
m
|
>
r and the integral (
J
k
m
f
)
(
z
) can be calculated by residues at |
z
-
a
m
|
>
r. Introduce the operator (91)
W
m
1
,
m
2
m
f
z
=
r
z
+
m
1
ω
1
+
m
2
ω
2
2
f
a
m
+
r
2
z
-
a
m
+
m
1
ω
1
+
m
2
ω
2
¯
¯
.We have for m
=
k
(92)
J
k
k
f
z
=
-
∑
m
1
,
m
2
∈
Z
∗
W
m
1
,
m
2
k
f
zand for m
≠
k
(93)
J
k
m
f
z
=
-
∑
m
1
,
m
2
∈
Z
W
m
1
,
m
2
m
f
z
.The double series (91)–(93) are defined by the Eisenstein summation that correctly determines their convergence [24].Substitution of (91)–(93) into (48) yields the system of functional equations (94)
ψ
k
z
=
ρ
∑
m
=
1
N
∑
m
1
,
m
2
∈
Z
∗
W
m
1
,
m
2
m
ψ
m
z
+
1
,
z
-
a
k
≤
r
k
=
1,2
,
…
,
N
,where m
1
,
m
2 run over integers in the sum ∑
∗ with the excluded term m
1
=
m
2
=
0 for m
=
k. The functional equations (94) do not contain any integral. Each term (
W
m
1
,
m
2
(
m
)
ψ
m
)
(
z
) can be treated as the shift operator written in the form of functional compositions. Such equations can be easily solved by use of symbolic computations. The difficulty related to the infinite double summations in m
1
,
m
2 can be overcome by application of Eisenstein functions E
n
(
z
). Equations (94) can be solved by the contrast expansions in ρ (see Section 4.4) and by the cluster expansions in ν (see Section 4.5). We refer to [24, 34] where this approach is explained in detail. The computational efficiency of the method can be demonstrated by few analytical formulae below. Consider the cluster expansion in ν
=
N
π
r
2 which is equivalent to the expansion in r
2 [24]: (95)
ψ
m
z
=
∑
n
=
0
∞
ψ
m
n
z
r
2
n
,
z
-
a
k
≤
r
k
=
1,2
,
…
,
N
.Few first functions in the expansion (95) are given by the following exact formulae: (96)
ψ
m
0
z
=
1
,
ψ
m
1
z
=
ρ
∑
k
=
1
N
E
2
z
-
a
k
,
ψ
m
2
z
=
ρ
2
∑
k
,
k
1
=
1
N
E
2
a
k
-
a
k
1
¯
E
2
z
-
a
k
,
ψ
m
3
z
=
ρ
3
∑
k
,
k
1
,
k
2
=
1
N
E
2
a
k
-
a
k
1
E
2
a
k
1
-
a
k
2
¯
E
2
z
-
a
k
2
-
2
ρ
2
∑
k
,
k
1
=
1
N
E
3
a
k
-
a
k
1
¯
E
3
z
-
a
k
,where E
n
(
z
-
a
k
) is defined as E
~
n
(
z
-
a
k
)
=
E
n
(
z
-
a
k
)
-
(
z
-
a
k
)
-
p for k
=
m and E
n
(
0
)
:
=
E
~
n
(
0
). It is possible to proceed (96) and calculate the next approximations following the algorithm [24].Letq be a natural number; k
s runs over 1 to N, m
q
=
2,3
,
…
. Let C denote the operator of complex conjugation. Introduce the multiple convolution sums (97)
e
m
1
⋯
m
q
=
1
N
1
+
1
/
2
m
1
+
⋯
+
m
q
∑
k
0
k
1
⋯
k
q
E
m
1
a
k
0
-
a
k
1
E
m
2
a
k
1
-
a
k
2
¯
⋯
C
q
E
m
q
a
k
q
-
1
-
a
k
q
.The integrals from (57) are simplified by the mean value theorem for harmonic functions (98)
λ
11
-
i
λ
12
=
1
+
2
ρ
ν
1
N
∑
k
=
1
N
ψ
k
a
k
.Substitution of (96) and the next terms into (98) yields (99)
λ
11
-
i
λ
12
=
1
+
2
ρ
ν
1
+
A
1
ν
+
A
2
ν
2
+
⋯
,where (100)
A
1
=
ρ
π
e
2
,
A
2
=
ρ
2
π
2
e
22
,
A
3
=
1
π
3
-
2
ρ
2
e
33
+
ρ
3
e
222
,
A
4
=
1
π
4
3
ρ
2
e
44
-
2
ρ
3
e
332
+
e
233
+
ρ
4
e
2222
,
A
5
=
1
π
5
-
4
ρ
2
e
55
+
ρ
3
3
e
442
+
6
e
343
+
3
e
244
-
2
ρ
4
e
3322
+
e
2332
+
e
2233
+
ρ
5
e
22222
,
A
6
=
1
π
6
5
ρ
2
e
66
-
4
ρ
3
e
255
+
3
e
354
+
3
e
453
+
e
552
+
ρ
4
3
e
2244
+
6
e
2343
+
4
e
3333
+
3
e
2442
+
6
e
3432
+
3
e
4422
-
2
ρ
5
e
22233
+
e
22332
+
e
23322
+
e
33222
+
ρ
6
e
222222
.The next coefficients A
n can be written in closed form by application of the algorithm presented in [18]. The deterministic values (97) provide the benchmark to compute the effective conductivity of random composites. Following the Monte Carlo method one can take representative sets of a
∈
C to statistically calculate the expectations (101)
e
^
m
1
⋯
m
q
=
∫
C
e
m
1
⋯
m
q
a
P
a
d
a
.The infinite set {
e
^
m
1
⋯
m
q
,
m
j
=
2,3
,
…
} completely determine the random geometric structure of the considered class of composites and can be taken as the basic set in the RVE theory [25].A fast algorithm to compute the convolution sums (97) is described in [11]. It allows us to compute A
n (n
≤
20) within a given distribution of a
k in few hours on an ordinary notebook when N
=
64 and with the number of computational experiments 1500 in the Monte Carlo method. It follows from [11] that six terms (100) are sufficient to deduce an analytic formula for a uniform nonoverlapping distribution of the disks.Formulae (99)-(100) correspond to the cluster expansion. Application of the explicit iterations to the functional equations (94) yields the contrast expansion as a series in ρ. One can see that the following third-order formula in ρ takes place: (102)
λ
11
-
i
λ
12
=
1
+
2
ρ
ν
+
2
ρ
2
ν
2
e
2
π
+
2
ρ
3
ν
3
∑
n
=
2
∞
-
1
n
n
-
1
e
n
n
ν
n
-
2
π
n
+
O
ρ
4
.For macroscopically isotropic composites, the expectation of (102) has the form (103)
λ
e
=
1
+
2
ρ
ν
+
2
ρ
2
ν
2
+
2
ρ
3
ν
3
∑
n
=
2
∞
-
1
n
n
-
1
e
^
n
n
ν
n
-
2
π
n
+
O
ρ
4
,where the relation e
^
2
=
π from [41] is used. Another formula deduced in [41] can be useful in simulations and estimations of the third-order term of (102)-(103): (104)
e
n
n
=
-
1
n
N
n
+
1
∑
m
=
1
N
∑
k
=
1
N
E
n
a
m
-
a
k
2
.One can find the numerical values of e
^
n
n for uniform nonoverlapping distributions in [11].In order to compare the contrast expansion formula (103) and formula (23) of the traditional theory, we write (23) up to O
(
ρ
4
): (105)
ρ
2
ν
2
λ
e
-
1
-
1
λ
e
+
1
=
ν
ρ
-
A
3
1
ρ
3
+
O
ρ
4
.Substitution of (103) into (105) yields (106)
A
3
1
=
ν
3
∑
n
=
2
∞
-
1
n
n
-
1
e
^
n
n
ν
n
-
2
π
n
-
1
.The latter formula gives the Torquato-Milton parameter ζ
1 (see (
20.59
) and (
20.66
) from [9]): (107)
ζ
1
=
ν
2
1
-
ν
∑
n
=
2
∞
-
1
n
n
-
1
e
^
n
n
ν
n
-
2
π
n
-
1
.The three- and four-point contrast bounds on the effective conductivity are written by this parameter (see (21.33)–(21.35) and (21.42)–(21.44) from [9]).Remark 6.
Formula (107) is valid for arbitrary distributed circular inclusions. Any other shape of inclusions can be approximated by a disk packing.Remark 7.
In accordance with (77) computation of the Torquato-Milton parameter ζ
1 for general shapes of inclusions can be reduced to a sum of the triple integrals (108)
∫
L
k
∫
L
k
1
∫
L
k
2
t
2
E
1
t
2
-
t
1
¯
E
1
t
1
-
t
d
t
¯
d
t
1
d
t
2
¯that should be easier than computation of the corresponding triple integral (20.66) from [9] on the 3-point correlation function which is computed by another triple integral. The integral (108) can be calculated by residues for algebraic curves.


## 4.1. General

The generalized method of Schwarz (GMS) was proposed by Mikhlin [17] as a generalization of the classical alternating method of Schwarz for finitely connected domains. The GMS is based on the decomposition of the considered domain with complex geometry onto simple domains and subsequent solution to boundary value problems for simple domains. In our case, we have N simple domains, the cell Q with one-inclusion D
k (k
=
1,2
,
…
,
N). In each step of the algorithm the boundary conditions for the kth simple domain are corrected by influence of the mth simple domains (m
≠
k) computed at the precedent steps.We now shortly present the GMS, first for the conductivity problem (3)-(4) and further for general linear problems. The ideal contact condition (3) can be rewritten in the form (24)
u
0
=
u
k
-
u
e
x
t
,
∂
u
0
∂
n
=
∂
u
k
-
u
e
x
t
∂
n
+
ρ
λ
+
1
∂
u
k
∂
n
on
L
k
,
k
=
1,2
,
…
,
N
,where u
=
u
0
+
u
e
x
t is the decomposition of the potential u onto the regular periodic part u
0 and the external field u
e
x
t corresponding to the quasi periodicity conditions (4). In the case of the space problem the external field u
e
x
t has a singularity at infinity.Introduce the function(25)
f
t
=
λ
+
1
∂
u
k
∂
n
t
,
t
∈
L
k
k
=
1,2
,
…
,
Nand the domains D
+
:
=
⋃
k
=
1
n
D
k, D
-
:
=
D. It is convenient to represent the harmonic functions in different domains as one piecewise harmonic function u
0
(
x
) introduced above in D
=
D
- and equal to u
k
(
x
)
-
u
e
x
t
(
x
) in D
k
∪
L
k (k
=
1,2
,
…
,
N). Then, (24) can be considered as the jump problem: (26)
u
0
+
=
u
0
-
,
∂
u
0
+
∂
n
=
∂
u
0
-
∂
n
+
ρ
f
t
on
L
=
⋃
k
=
1
n
L
k
,where u
0
+ and u
0
- correspond to the limit values of u
(
x
) when x tends to t from D
+ and D
-, respectively. The jump problem (26) has a unique solution expressed in terms of the simple layer potential P for the curve L in the cell Q in the torus topology. We have (27)
u
0
z
=
ρ
P
f
x
,
x
∈
D
+
∪
D
-
.The operator P is considered as an operator in an appropriate functional space [29–32]. In the classic statement, Q should be replaced by the whole space. The operator P is decomposed onto the simple layer potentials P
k along curves L
k (k
=
1,2
,
…
,
N) as P
=
(
1
/
(
λ
+
1
)
)
∑
k
=
1
N
P
k. Here, the multiplier 1
/
(
λ
+
1
) is introduced for short form of the equations below. Then (27) implies that (28)
u
k
x
=
ρ
∑
m
=
1
N
P
m
∂
u
m
+
∂
n
x
+
u
e
x
t
x
,
x
∈
D
k
k
=
1,2
,
…
,
N
,
(29)
u
x
=
ρ
∑
m
=
1
N
P
m
∂
u
m
+
∂
n
x
+
u
e
x
t
x
,
x
∈
D
.Equations (28) can be considered as a system of integral equations on the potentials u
k
(
x
) in the inclusions D
k
(
k
=
1,2
,
…
,
N
). If it is solved, the potential u
(
x
) is calculated in D by (29).Remark 1.
Equations (28) correspond to the GMS in Mikhlin’s form [17] when the method of successive approximations can diverge for closely spaced inclusions. The following slight modification yields a convergent algorithm for any location of inclusions and every contrast parameter (|
ρ
|
≤
1): (30)
u
k
x
=
ρ
∑
m
=
1
N
P
m
∂
u
m
+
∂
n
x
-
P
m
∂
u
m
+
∂
n
w
+
u
e
x
t
x
+
ρ
∑
m
=
1
N
P
m
∂
u
m
+
∂
n
w
,
x
∈
D
k
∪
L
k
k
=
1,2
,
…
,
N
,where w is a fixed point in D.Equations (27)–(30) are written for potentials. Their differentiation yields equations on the intensity field: (31)
E
k
x
=
ρ
∑
m
=
1
N
Q
m
E
m
x
+
E
0
z
,
x
∈
D
k
∪
L
k
k
=
1,2
,
…
,
N
,where E
k
(
x
)
=
E
(
x
) in D
k and Q
m are appropriate operators (similar to double layer potentials). Equation (31) is similar to (20).This GMS for harmonic functions can be extended to general problems governed by linear equations describing electrical and heat conduction, flow in porous media, viscous flow, elasticity, and coupled fields inR
d (d
=
2,3).Let tensor fieldsU
(
z
) in D and U
k
(
z
) in D
k satisfy the contact conditions (32)
T
·
U
=
T
k
·
U
k
,
on
⁡
L
k
k
=
1,2
,
…
,
N
,where T and T
k are boundary operators involving physical parameters of materials occupying the domains D and D
k, respectively. Let T
k be presented in the form T
k
=
T
-
ρ
k, where ρ
k denotes a contrast parameter tensor. Then (32) becomes (33)
U
=
U
k
-
T
-
1
·
ρ
k
·
U
k
,
on
L
k
k
=
1,2
,
…
,
N
.Let P
k denote the simple layer potential corresponding to the considered linear equation for the domains D
k and D
k
- separated by L
k. Application of ∑
k
=
1
N
P
k to (33) yields (34)
U
k
x
=
∑
m
=
1
N
P
m
·
T
-
1
·
ρ
m
·
U
m
x
+
U
e
x
t
x
,
x
∈
D
k
∪
L
k
k
=
1,2
,
…
,
N
.There are two different methods to solve (34). The first method is based on the direct iterations and corresponds to the contrast expansion in ρ (compare to Section 3.2). The second method is based on implicit iterations applied to the same equations (34) written in the form (35)
U
k
x
-
P
k
·
T
-
1
·
ρ
k
·
U
k
x
=
∑
m
≠
k
N
P
m
·
T
-
1
·
ρ
m
·
U
m
x
+
U
e
x
t
x
,
x
∈
D
k
∪
L
k
k
=
1,2
,
…
,
N
.One-inclusion problem forD
k
(
k
=
1,2
,
…
,
N
) is solved at each step of iterations. This scheme corresponds to the cluster expansions discussed in Section 3.1.Remark 2.
Equations (27) and (35) are not reduced to the classic integral equations constructed in the framework of the potentials theory [32]. It is another type of equations.Remark 3.
Equations (27)–(29) and (34)-(35) can be derived in the limit cases of inclusions (soft and hard inclusions in terminology [33]) by introduction of fictive potentials [10].


## 4.2. Integral Equations

In the present section, we discuss the GMS for the double periodic problem (3)-(4). Following [34] we introduce the complex potentials ϕ
(
z
) and ϕ
k
(
z
) in such a way that (36)
u
x
=
Re
⁡
ϕ
z
,
u
k
x
=
2
λ
+
1
Re
⁡
ϕ
k
z
,where x
=
(
x
1
,
x
2
) and z
=
x
1
+
i
x
2. Two real equations (3) can be written as one R-linear condition: (37)
ϕ
t
=
ϕ
k
t
-
ρ
ϕ
k
t
¯
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
,where the bar denotes the complex conjugation. The real part of (37) yields the first equation in (3). The imaginary part yields Im
⁡
ϕ
t
=
1
+
ρ
Im
⁡
ϕ
k
t. After differentiation of the latter equality along L
k and application of the Cauchy-Riemann equations we arrive at the second equation in (3) (see for details [10]).The unknown functionsϕ
z and ϕ
k
z are analytic in D and D
k, respectively, and continuously differentiable in the closures of the domains considered. It follows from (4) that the function ϕ
(
z
) is quasiperiodic: (38)
ϕ
z
+
ω
j
-
ϕ
z
=
ξ
j
+
i
d
j
j
=
1,2
,where the constants ξ
1 and ξ
2 are taken from (4) and d
1 and d
2 are undetermined real constants which should be found during solution to the problem.The following formula was proved in [10]: (39)
d
d
t
ϕ
k
t
=
-
n
2
t
¯
ϕ
k
′
t
¯
,
t
∈
L
k
.Here, the differential operator is understood as d
/
d
t
=
t
s
′
(
d
/
d
s
), where t
=
t
(
s
) is a parametric equation of L
k and ϕ
k
′
(
t
) is the boundary value of the complex derivative ∂
/
∂
z
=
∂
/
∂
x
1
-
i
(
∂
/
∂
x
2
). Differentiation of (37) and application of (39) yield (40)
ψ
t
=
ψ
k
t
+
ρ
n
k
2
t
¯
ψ
k
t
¯
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
,where ψ
(
z
)
=
ϕ
′
(
z
), ψ
k
(
z
)
=
ϕ
k
′
(
z
), and n
k
(
t
) stands for the unit normal outward vector to L
k expressed in terms of the complex values. The function ψ
(
z
) is doubly periodic. The R-linear problem (40) holds also in the limit cases for perfectly conducting inclusions (ρ
=
1) and for isolators (ρ
=
-
1).TheR-linear problem (37)-(38) can be reduced to a system of integral equation in the following way. The cell Q in the torus topology is divided into two domains D
+
=
⋃
k
=
1
n
D
k (not connected) and D
-
=
D (multiply connected). Let L
=
⋃
k
=
1
n
L
k denote the boundary of D
+ and let μ
(
t
) be a Hölder continuous function on L. Let ζ
(
z
) denote the Weierstrass ζ-function for which [28] (41)
ζ
z
+
ω
j
-
ζ
z
=
δ
j
j
=
1,2
,where δ
j
=
2
ζ
(
ω
j
/
2
). Following [35] introduce the Eisenstein function (42)
E
1
z
=
∑
m
1
,
m
2
∈
Z
1
z
+
m
1
ω
1
+
m
2
ω
2
,where the Eisenstein summation is used [35]. The Weierstrass and Eisenstein functions are related by formula E
1
(
z
)
=
ζ
(
z
)
-
S
2
z, where S
2
=
(
2
/
ω
1
)
ζ
(
ω
1
/
2
). It follows from Legendre’s identity δ
1
ω
2
-
δ
2
ω
1
=
2
π
i (see [28, 35]) and (41) that the jumps of E
1
(
z
) have the form (43)
E
1
z
+
ω
1
-
E
1
z
=
0
,
E
1
z
+
ω
2
-
E
1
z
=
-
2
π
i
ω
1
.The Cauchy-type integral on torus was introduced with theζ-function in the kernel [36]. We introduce it in the slight different form using the Eisenstein function (44)
Φ
z
=
1
2
π
i
∫
L
μ
t
E
1
t
-
z
d
t
+
C
z
+
C
0
,
z
∈
D
±
,where μ
(
t
) is a Hölder continuous function on L and C and C
0 are complex constants. The function Φ
(
z
) is analytic in D
± and quasiperiodic; that is, it has constant jumps per periodicity cell calculated with the help of (43): (45)
Φ
z
+
ω
1
-
Φ
z
=
C
ω
1
,
Φ
z
+
ω
2
-
Φ
z
=
-
1
ω
1
∫
L
μ
t
d
t
+
C
ω
2
.The following Sochocki (Sokhotski-Plemelj) formulae take place on torus [36]: (46)
Φ
±
t
=
±
1
2
μ
t
+
1
2
π
i
∫
L
μ
τ
E
1
τ
-
t
d
τ
+
C
t
+
C
0
,
t
∈
L
,where Φ
±
(
t
) denote the limit values of Φ
(
z
) on L when z
∈
D
± tends to t
∈
L, respectively. The principal value of the integral in (46) is considered following [37]. Formulae (46) imply that the function Φ
(
z
) satisfies the jump condition (47)
Φ
+
t
-
Φ
-
t
=
μ
t
,
t
∈
L
.Equation (47) can be considered as a C-linear conjugation problem [10] in a class of quasiperiodic functions with the given jump μ
(
t
). It follows from [36] that the general solution of (47) up to an additive arbitrary constant has the form (44).Consider (37) as the problem (47) with μ
(
t
)
=
ρ
ϕ
k
(
t
)
¯ on L, Φ
(
z
)
=
ϕ
k
(
z
) in D
k
⊂
D
+, and Φ
(
z
)
=
ϕ
(
z
) in D
=
D
-. Then, (44) can be written in D
+ as follows: (48)
ϕ
k
z
=
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
ϕ
m
t
¯
E
1
t
-
z
d
t
+
C
z
+
C
0
,
z
∈
D
k
k
=
1,2
,
…
,
N
.Equation (44) in D
- becomes(49)
ϕ
z
=
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
ϕ
m
t
¯
E
1
t
-
z
d
t
+
C
z
+
C
0
,
z
∈
D
.Equations (48)-(49) are deduced from the R-linear problem (3)-(4). Using the same arguments one can derive analogous equations from the problem (40) on the complex flux: (50)
ψ
k
z
=
-
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
n
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
C
,
z
∈
D
k
k
=
1,2
,
…
,
N
,
(51)
ψ
z
=
-
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
n
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
C
,
z
∈
D
.It is also possible to obtain these equations by differentiation of (48)-(49) using integration by parts and formula (39).The integral equations (50) are considered as an equation in the following Banach space. First, the space C
(
0
,
α
) of functions Hölder continuous on L is considered (0
<
α
≤
1). This is a Banach space endowed with the norm: (52)
ω
=
sup
t
∈
L
⁡
ω
t
+
sup
t
1,2
∈
L
t
1
≠
t
2
⁡
ω
t
1
-
ω
t
2
t
1
-
t
2
α
.Consider the closed subspace A
⊂
C
(
0
,
α
) consisting of functions analytically continued into all D
k. Therefore, every element of A is a function analytic in the domain ⋃
k
=
1
N
D
k and Hölder continuous in its closure. It is worth noting that the domain ⋃
k
=
1
N
D
k is not connected and convergence in A implies the uniform convergence of functions in ⋃
k
=
1
N
(
D
k
∪
L
k
).Equations (50) are written in the domains D
k and can be continued to the boundary by use of (46) (53)
ψ
k
t
=
-
ρ
∑
m
=
1
N
1
2
n
m
2
t
¯
ψ
m
t
¯
+
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
C
,
t
∈
L
k
k
=
1,2
,
…
,
N
.Therefore, the integral equations in the space A have the form (50), (53). It is proved in Section 4.3 that (50), (53) have a unique solution. Solution to these equations corresponds to the GMS which can be realized in different forms as it is noted in Section 4.1.The problem (37)-(38) for the complex potentials has a unique solution up to an arbitrary additive constant C
0. The constants ξ
1
+
i
d
1, ξ
2
+
i
d
2, and C
0 disappear in the problem (40) after differentiation. The structure of the general solution of (40) has the form ψ
(
z
)
=
ξ
1
ψ
(
1
)
(
z
)
+
ξ
2
ψ
(
2
)
(
z
), where the real constants ξ
1 and ξ
2 correspond to the external field (jumps of u
(
x
) per a cell). Hence, the R-linear problem (40) has two R-linear independent solutions. From the other side, these two independent solutions can be constructed by (50)-(51) taking, for instance, C
=
1 and C
=
i. Therefore, the complex constant C can be expressed through the real constants ξ
1 and ξ
2. In order to calculate the effective properties of macroscopically isotropic composites it is sufficient to fix ξ
1, ξ
2 and find C. For simplicity, we put (54)
ξ
1
=
ω
1
,
ξ
2
=
Re
⁡
ω
2
.Then the external field corresponds to the complex potential ϕ
(
e
x
t
)
(
z
)
=
z and to the complex flux ψ
(
e
x
t
)
(
z
)
=
1. In the vector form, this external flux becomes (
-
1,0
). It follows from the first relation (45) that (55)
C
=
1
+
i
d
1
ω
1
,where d
1
=
Im
⁡
[
ϕ
(
z
+
ω
1
)
-
ϕ
(
z
)
]. Put z
=
-
(
1
/
2
)
(
ω
1
+
ω
2
) in this relation and calculate(56)
d
1
=
Im
⁡
∫
-
1
/
2
ω
1
+
ω
2
1
/
2
ω
1
-
ω
2
ψ
t
d
t
=
-
∫
-
1
/
2
ω
1
+
R
e
ω
2
1
/
2
ω
1
-
R
e
ω
2
∂
u
∂
x
2
x
1
-
i
Im
⁡
ω
2
2
d
x
1
.The integral from the left hand part of (56) expresses the total flux through the low side of the parallelogram Q parallel to the x
1-axis. It must be equal to zero because of the periodicity and that the external flux (
-
1,0
) vanishes in the x
2-direction. Therefore, d
1
=
0; hence C
=
1 by (55).The effective conductivity tensor of deterministic composites can be calculated through the complex potentials in the inclusions (see [22, 34]) (57)
λ
11
-
i
λ
12
=
1
+
2
ρ
∑
k
=
1
N
∫
D
k
ψ
k
x
1
+
i
x
2
d
x
1
d
x
2
,where ψ
k
(
z
) (k
=
1,2
,
…
,
N) satisfy (50)-(51), (53) with C
=
1. The value (57) depends on the locations of inclusions, that is, on a
∈
C. The effective conductivity of macroscopically isotropic random composites is defined through the operator (7) applied to (57):(58)
λ
e
=
λ
11
^
,where the macroscopic isotropy implies that λ
12
^ vanishes and λ
11
^
=
λ
22
^.


## 4.3. Convergence of the GMS

This section is devoted to study of the integral equations (50), (53) in the space A.LetΨ
(
z
)
=
ψ
k
(
z
) for z
∈
D
k
∪
L
k and let A denote the operator from the right hand part of (48), (53). Then, (48), (53) can be shortly written as the following equation in the space A: (59)
Ψ
=
ρ
A
Ψ
+
1
.In the present section, we demonstrate that (59) has the unique solution for |
ρ
|
<
1 given by the series (60)
Ψ
=
∑
n
=
0
∞
ρ
n
A
n
1
.We use the following general result of functional analysis.Theorem 4 (see [38]).
LetB be a linear bounded operator in a Banach space B. If for any element f
∈
B and for any complex number ν satisfying the inequality |
ν
|
≤
1 equation (61)
Ψ
=
ν
B
Ψ
+
fhas a unique solution, then the unique solution of (62)
Ψ
=
B
Ψ
+
fcan be found by the method of successive approximations. The approximations converge in B to the solution (63)
Ψ
=
∑
k
=
0
∞
B
k
f
.Theorem4 can be applied to (59) in the Banach space A.Theorem 5.
Let|
ρ
|
<
1. Equation Ψ
=
ρ
A
Ψ
+
1 has a unique solution. This solution can be found by the method of successive approximations convergent in the space A.Proof.
Let|
ν
|
≤
1. Consider (64)
ψ
k
z
=
-
ν
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
f
k
z
,
z
∈
D
k
k
=
1,2
,
…
,
N
,
ψ
k
t
=
-
ν
ρ
1
2
n
k
2
t
¯
ψ
k
t
¯
+
∑
m
=
1
N
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
f
k
t
,
t
∈
L
k
k
=
1,2
,
…
,
N
,where f
k
∈
A
k.
Letψ
k
(
z
) be a solution of (53). Introduce the function ψ
(
z
) analytic in D and Hölder continuous in its closure as follows: (65)
ψ
z
=
-
ν
ρ
∑
m
=
1
n
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
,
z
∈
D
.The expression in the right hand part of (65) can be considered as Cauchy’s integral on torus (66)
Φ
z
=
1
2
π
i
∫
L
μ
t
E
1
t
-
z
d
tfor z
∈
Q
∖
L with the density μ
(
t
)
=
-
ν
ρ
n
m
2
(
t
)
¯
ψ
m
(
t
)
¯ on L
m. Using the property (47) of (66) and (64) on L
k we arrive at the R-linear conjugation relation on each fixed curve L
k: (67)
ψ
k
t
-
f
k
t
-
ψ
t
=
-
ν
ρ
n
k
2
t
¯
ψ
k
t
¯
,
t
∈
L
k
.Here, Φ
+
(
t
)
=
ψ
k
(
t
)
-
f
k
(
t
) and Φ
-
(
t
)
=
ψ
(
t
). Integration of (67) along L
k yields (68)
ϕ
t
=
ϕ
k
t
-
ν
ρ
ϕ
k
t
¯
-
g
k
t
+
γ
k
,
t
∈
L
k
k
=
1,2
,
…
,
N
,where the functions ϕ
(
z
), ϕ
k
(
z
) are analytic in D, D
k and belong to C
(
1
,
α
) in the closures of the domains considered, g
k
(
z
) is a primitive function of f
k
(
z
), and γ
k is an arbitrary constant. Moreover, the function ϕ
(
z
) is quasiperiodic; hence (69)
ϕ
z
=
φ
z
+
c
z
,where φ
(
z
) is periodic and c is a constant. The constant γ
k can be included into φ
k
(
z
) by the representation γ
k
=
α
k
-
ν
ρ
α
k
¯. Then, (68) becomes (70)
φ
t
=
φ
k
t
-
ν
ρ
φ
k
t
¯
-
g
k
t
-
c
t
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
.It follows from [20] and Appendix B of [22] that the general solution of the problem (70) has the following structure: (71)
φ
z
=
φ
1
z
+
φ
2
z
,
φ
k
z
=
φ
k
1
z
+
φ
k
2
z
,where (72)
φ
1
t
=
φ
k
1
t
-
ν
ρ
φ
k
1
t
¯
-
g
k
t
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
,
(73)
φ
2
t
=
φ
k
2
t
-
ν
ρ
φ
k
2
t
¯
-
c
t
,
t
∈
L
k
,
k
=
1,2
,
…
,
N
.Each of the problems (72) and (73) (with fixed c) has a unique solution since |
ν
ρ
|
<
1 (the elliptic case in Mikhailov’s terminology [39]). The uniqueness also follows from [20]. This unique solution of the problem (73) can be easily found as φ
(
2
)
(
z
)
=
c
z and φ
k
(
2
)
(
z
)
=
0. Then, differentiation of (71)–(73) implies that the general solution of the problem (67) has the form (74)
ψ
z
=
ψ
1
z
+
c
,where the function ψ
k
(
1
)
(
z
) is uniquely determined from the problem (75)
ψ
1
t
=
ψ
k
1
t
+
ν
ρ
n
k
2
t
¯
ψ
k
1
t
¯
-
f
k
t
,
t
∈
L
k
k
=
1,2
,
…
,
N
.Comparison of (67) and (75) proves the theorem.The series(76)
Ψ
=
∑
k
=
0
∞
ρ
k
A
k
1converges absolutely for |
ρ
|
<
1, since it can be represented in the form Ψ
=
∑
k
=
0
∞
s
k
ρ
0
k
A
k
1 with s
ρ
0
=
ρ and |
ρ
|
<
|
ρ
0
|
<
1. The latter series converges absolutely with the rate |
s
| because the series ∑
k
=
0
∞
ρ
0
k
A
k
1 converges uniformly; hence, its general term tends to zero as k
→
∞.Uniform convergence of (76) for |
ρ
|
=
1 can be justified by use of the arguments form [10].


## 4.4. Contrast Expansions for the GMS

It is proved in the previous section that the unique solution of (50), (53) can be found by a method of successive approximations uniformly convergent in ⋃
k
=
1
N
(
D
k
∪
L
k
). The method of successive approximations can be also presented in the form of the following iterative scheme: (77)
ψ
k
0
z
=
1
,
ψ
k
p
z
=
-
ρ
∑
m
=
1
N
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
p
-
1
t
¯
E
1
t
-
z
d
t
+
1
,
p
=
1,2
,
…
,
z
∈
D
k
k
=
1,2
,
…
,
N
.When the functions ψ
k
(
z
) (k
=
1,2
,
…
,
N) are found, the effective conductivity is calculated by (57).In order to compare (21) and (48), (53) we introduce the complex intensity field E
(
z
)
=
E
1
(
z
)
-
i
E
2
(
z
) isomorphic to the vector field E
(
x
)
=
(
E
1
(
x
)
,
E
2
(
x
)
). Then, the vector equation (21) can be written as a complex equation on E
(
z
). For definiteness, we consider the complex potential ψ
k
(
z
) in D
k. It follows from (5), (36) and ψ
k
(
z
)
=
ϕ
k
′
(
z
) that (78)
E
z
=
∂
u
k
∂
x
1
-
i
∂
u
k
∂
x
2
=
2
λ
+
1
ψ
k
z
.Therefore, the complex equation on E
(
z
) in D
k can be written as (48) on ψ
k
(
z
). We do not directly prove that (21) and (59) are the same. But it is easy to show that the forms of their solutions (22) and (60) coincide. Since the solutions E and Ψ coincide (more precisely, isomorphic) for E
0
=
(
-
1,0
), the power series (22) and (60) in ρ can be considered as the presentations of the same function for sufficiently small |
ρ
|. Therefore, the coefficients of (22) and (60) coincide. In particular, this fact implies that substitution of (60) into (57) yields a series which coincides with the series (19) obtained by the traditional approach. It follows from [10, 27] that the contrast expansion (60) converges for all admissible |
ρ
|
≤
1.


## 4.5. Cluster Expansions for the GMS

The GMS can be presented by other iterative schemes corresponding to the traditional cluster expansions summarized in Section3.1. First, we rewrite (48) in the form (79)
ψ
k
z
+
ρ
2
π
i
∫
L
k
n
k
2
t
¯
ψ
k
t
¯
E
1
t
-
z
d
t
=
-
ρ
∑
m
≠
k
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
t
¯
E
1
t
-
z
d
t
+
1
,
z
∈
D
k
k
=
1,2
,
…
,
n
.Following Section4.2 for every fixed m
=
1,2
,
…
,
N introduce the Banach space A
m of functions analytic in D
m and Hölder continuous in D
m
∪
L
m. Let f
m
(
z
)
∈
A
m. For each fixed k
=
1,2
,
…
,
N introduce the operators (80)
P
k
m
f
m
z
=
1
2
π
i
∫
L
m
n
m
2
t
¯
f
m
t
¯
E
1
t
-
z
d
t
,
z
∈
D
k
m
=
1,2
,
…
,
n
.One can check that the operator P
k
k is singular and P
k
m (m
≠
k) are compact in A
k. Equation (79) can be shortly written in the form (81)
I
+
ρ
P
k
k
ψ
k
=
-
ρ
∑
m
≠
k
P
k
m
ψ
m
+
1
,
k
=
1,2
,
…
,
n
.The (zero)th approximation in the concentration ν for (79) can be written as the following integral equation: (82)
ψ
k
0
z
+
ρ
2
π
i
∫
L
k
n
k
2
t
¯
ψ
k
0
t
¯
E
1
t
-
z
d
t
=
1
,
z
∈
D
k
.Equation (82) corresponds to the R-linear problem (40) for one-inclusion D
k in the cell Q. Solution to the one-inclusion problem (82) defines the inverse operator (
I
+
ρ
P
k
k
)
-
1 bounded in A
k. Then, (81) can be written in the equivalent form (83)
ψ
k
=
-
ρ
∑
m
≠
k
I
+
ρ
P
k
k
-
1
P
k
m
ψ
m
+
I
+
ρ
P
k
k
-
1
1
,
k
=
1,2
,
…
,
n
.As it is established in Section 4.4 the vector-function E
(
x
) and the complex functions ψ
k
(
z
) express the same vector field in D
k up to an isomorphism for E
0
=
(
-
1,0
). Therefore, application of the successive approximations to (83) yields the same series as (12) in D
k. A detailed comparison can be performed. It shows, for instance, that the term E
0
+
∑
k
=
1
N
M
1
(
x
;
k
)
·
E
0 from (12) corresponds to (
I
+
ρ
P
k
k
)
-
1
1 from (83).The solution of the problem (82) depends on the shape of D
k. Therefore, the first order approximation in ν for the effective conductivity also depends on the shape and can be calculated by (57). If all the inclusions have the same form, the functions ψ
k
(
0
)
(
z
) do not depend on k. Then, we arrive at the formula valid up to O
(
ν
2
): (84)
λ
e
≈
1
+
2
ρ
ν
α
,where the shape factor α is introduced following [23, 40]: (85)
α
=
1
D
1
∫
D
1
ψ
1
0
x
1
+
i
x
2
d
x
1
d
x
2
.The Clausius-Mossotti formula equivalent to (84) up to O
(
ν
2
) is used in applications [6]: (86)
λ
e
≈
1
+
ρ
ν
α
1
-
ρ
ν
α
.One can find limitations of the formulae (84), (86) and a comparison of their accuracy in [23].The higher order iterations in the cluster expansion are based on the integral equations onψ
k
(
p
)
(
z
): (87)
ψ
k
p
z
+
ρ
2
π
i
∫
L
k
n
k
2
t
¯
ψ
k
p
t
¯
E
1
t
-
z
d
t
=
-
ρ
∑
m
≠
k
1
2
π
i
∫
L
m
n
m
2
t
¯
ψ
m
p
-
1
t
¯
E
1
t
-
z
d
t
+
1
,
z
∈
D
k
k
=
1,2
,
…
,
n
;
p
=
1,2
,
…
.The uniform convergence of the iterative scheme (87) was proved in the case ρ
=
1 (perfectly conducting inclusions) in Section 4.9
.
2 of [10]. The higher order iterations (87) yield approximate analytical formulae for λ
e for arbitrary nonoverlapping locations of different disks [24, 34]. An example of equal disks is discussed below in Section 4.6.


## 4.6. Circular Inclusions

Let the inclusionsD
k be disks of radius r centered at the points a
k. Then, the normal unit vector n
k
(
t
) to the circle |
t
-
a
k
|
=
r has the form n
k
(
t
)
=
(
t
-
a
k
)
/
r. The inversion with respect to |
t
-
a
k
|
=
r is defined by formula (88)
z
k
∗
=
r
2
z
-
a
k
¯
+
a
k
.Consider the integral operator from (48): (89)
J
k
m
f
z
=
1
2
π
i
∫
L
m
n
m
2
t
¯
f
t
¯
E
1
t
-
z
d
t
,
z
∈
D
k
,where f
(
z
) is analytic in the disk D
m. We have t
=
t
(
m
)
∗ on L
m. Hence, the operator (89) can be written in the form (90)
J
k
m
f
z
=
1
2
π
i
∫
L
m
r
t
-
a
m
2
f
t
m
∗
¯
E
1
t
-
z
d
t
,
z
∈
D
k
,where the function f
(
t
(
m
)
∗
)
¯ is analytically continued into |
z
-
a
k
|
>
r.The functionE
1
(
z
) has poles of first order at the points z
=
m
1
ω
1
+
m
2
ω
2 (see [35]). The function r
/
(
z
-
a
m
)
2
f
(
z
(
m
)
∗
)
¯ is analytic in the domain |
z
-
a
m
|
>
r and the integral (
J
k
m
f
)
(
z
) can be calculated by residues at |
z
-
a
m
|
>
r. Introduce the operator (91)
W
m
1
,
m
2
m
f
z
=
r
z
+
m
1
ω
1
+
m
2
ω
2
2
f
a
m
+
r
2
z
-
a
m
+
m
1
ω
1
+
m
2
ω
2
¯
¯
.We have for m
=
k
(92)
J
k
k
f
z
=
-
∑
m
1
,
m
2
∈
Z
∗
W
m
1
,
m
2
k
f
zand for m
≠
k
(93)
J
k
m
f
z
=
-
∑
m
1
,
m
2
∈
Z
W
m
1
,
m
2
m
f
z
.The double series (91)–(93) are defined by the Eisenstein summation that correctly determines their convergence [24].Substitution of (91)–(93) into (48) yields the system of functional equations (94)
ψ
k
z
=
ρ
∑
m
=
1
N
∑
m
1
,
m
2
∈
Z
∗
W
m
1
,
m
2
m
ψ
m
z
+
1
,
z
-
a
k
≤
r
k
=
1,2
,
…
,
N
,where m
1
,
m
2 run over integers in the sum ∑
∗ with the excluded term m
1
=
m
2
=
0 for m
=
k. The functional equations (94) do not contain any integral. Each term (
W
m
1
,
m
2
(
m
)
ψ
m
)
(
z
) can be treated as the shift operator written in the form of functional compositions. Such equations can be easily solved by use of symbolic computations. The difficulty related to the infinite double summations in m
1
,
m
2 can be overcome by application of Eisenstein functions E
n
(
z
). Equations (94) can be solved by the contrast expansions in ρ (see Section 4.4) and by the cluster expansions in ν (see Section 4.5). We refer to [24, 34] where this approach is explained in detail. The computational efficiency of the method can be demonstrated by few analytical formulae below. Consider the cluster expansion in ν
=
N
π
r
2 which is equivalent to the expansion in r
2 [24]: (95)
ψ
m
z
=
∑
n
=
0
∞
ψ
m
n
z
r
2
n
,
z
-
a
k
≤
r
k
=
1,2
,
…
,
N
.Few first functions in the expansion (95) are given by the following exact formulae: (96)
ψ
m
0
z
=
1
,
ψ
m
1
z
=
ρ
∑
k
=
1
N
E
2
z
-
a
k
,
ψ
m
2
z
=
ρ
2
∑
k
,
k
1
=
1
N
E
2
a
k
-
a
k
1
¯
E
2
z
-
a
k
,
ψ
m
3
z
=
ρ
3
∑
k
,
k
1
,
k
2
=
1
N
E
2
a
k
-
a
k
1
E
2
a
k
1
-
a
k
2
¯
E
2
z
-
a
k
2
-
2
ρ
2
∑
k
,
k
1
=
1
N
E
3
a
k
-
a
k
1
¯
E
3
z
-
a
k
,where E
n
(
z
-
a
k
) is defined as E
~
n
(
z
-
a
k
)
=
E
n
(
z
-
a
k
)
-
(
z
-
a
k
)
-
p for k
=
m and E
n
(
0
)
:
=
E
~
n
(
0
). It is possible to proceed (96) and calculate the next approximations following the algorithm [24].Letq be a natural number; k
s runs over 1 to N, m
q
=
2,3
,
…
. Let C denote the operator of complex conjugation. Introduce the multiple convolution sums (97)
e
m
1
⋯
m
q
=
1
N
1
+
1
/
2
m
1
+
⋯
+
m
q
∑
k
0
k
1
⋯
k
q
E
m
1
a
k
0
-
a
k
1
E
m
2
a
k
1
-
a
k
2
¯
⋯
C
q
E
m
q
a
k
q
-
1
-
a
k
q
.The integrals from (57) are simplified by the mean value theorem for harmonic functions (98)
λ
11
-
i
λ
12
=
1
+
2
ρ
ν
1
N
∑
k
=
1
N
ψ
k
a
k
.Substitution of (96) and the next terms into (98) yields (99)
λ
11
-
i
λ
12
=
1
+
2
ρ
ν
1
+
A
1
ν
+
A
2
ν
2
+
⋯
,where (100)
A
1
=
ρ
π
e
2
,
A
2
=
ρ
2
π
2
e
22
,
A
3
=
1
π
3
-
2
ρ
2
e
33
+
ρ
3
e
222
,
A
4
=
1
π
4
3
ρ
2
e
44
-
2
ρ
3
e
332
+
e
233
+
ρ
4
e
2222
,
A
5
=
1
π
5
-
4
ρ
2
e
55
+
ρ
3
3
e
442
+
6
e
343
+
3
e
244
-
2
ρ
4
e
3322
+
e
2332
+
e
2233
+
ρ
5
e
22222
,
A
6
=
1
π
6
5
ρ
2
e
66
-
4
ρ
3
e
255
+
3
e
354
+
3
e
453
+
e
552
+
ρ
4
3
e
2244
+
6
e
2343
+
4
e
3333
+
3
e
2442
+
6
e
3432
+
3
e
4422
-
2
ρ
5
e
22233
+
e
22332
+
e
23322
+
e
33222
+
ρ
6
e
222222
.The next coefficients A
n can be written in closed form by application of the algorithm presented in [18]. The deterministic values (97) provide the benchmark to compute the effective conductivity of random composites. Following the Monte Carlo method one can take representative sets of a
∈
C to statistically calculate the expectations (101)
e
^
m
1
⋯
m
q
=
∫
C
e
m
1
⋯
m
q
a
P
a
d
a
.The infinite set {
e
^
m
1
⋯
m
q
,
m
j
=
2,3
,
…
} completely determine the random geometric structure of the considered class of composites and can be taken as the basic set in the RVE theory [25].A fast algorithm to compute the convolution sums (97) is described in [11]. It allows us to compute A
n (n
≤
20) within a given distribution of a
k in few hours on an ordinary notebook when N
=
64 and with the number of computational experiments 1500 in the Monte Carlo method. It follows from [11] that six terms (100) are sufficient to deduce an analytic formula for a uniform nonoverlapping distribution of the disks.Formulae (99)-(100) correspond to the cluster expansion. Application of the explicit iterations to the functional equations (94) yields the contrast expansion as a series in ρ. One can see that the following third-order formula in ρ takes place: (102)
λ
11
-
i
λ
12
=
1
+
2
ρ
ν
+
2
ρ
2
ν
2
e
2
π
+
2
ρ
3
ν
3
∑
n
=
2
∞
-
1
n
n
-
1
e
n
n
ν
n
-
2
π
n
+
O
ρ
4
.For macroscopically isotropic composites, the expectation of (102) has the form (103)
λ
e
=
1
+
2
ρ
ν
+
2
ρ
2
ν
2
+
2
ρ
3
ν
3
∑
n
=
2
∞
-
1
n
n
-
1
e
^
n
n
ν
n
-
2
π
n
+
O
ρ
4
,where the relation e
^
2
=
π from [41] is used. Another formula deduced in [41] can be useful in simulations and estimations of the third-order term of (102)-(103): (104)
e
n
n
=
-
1
n
N
n
+
1
∑
m
=
1
N
∑
k
=
1
N
E
n
a
m
-
a
k
2
.One can find the numerical values of e
^
n
n for uniform nonoverlapping distributions in [11].In order to compare the contrast expansion formula (103) and formula (23) of the traditional theory, we write (23) up to O
(
ρ
4
): (105)
ρ
2
ν
2
λ
e
-
1
-
1
λ
e
+
1
=
ν
ρ
-
A
3
1
ρ
3
+
O
ρ
4
.Substitution of (103) into (105) yields (106)
A
3
1
=
ν
3
∑
n
=
2
∞
-
1
n
n
-
1
e
^
n
n
ν
n
-
2
π
n
-
1
.The latter formula gives the Torquato-Milton parameter ζ
1 (see (
20.59
) and (
20.66
) from [9]): (107)
ζ
1
=
ν
2
1
-
ν
∑
n
=
2
∞
-
1
n
n
-
1
e
^
n
n
ν
n
-
2
π
n
-
1
.The three- and four-point contrast bounds on the effective conductivity are written by this parameter (see (21.33)–(21.35) and (21.42)–(21.44) from [9]).Remark 6.
Formula (107) is valid for arbitrary distributed circular inclusions. Any other shape of inclusions can be approximated by a disk packing.Remark 7.
In accordance with (77) computation of the Torquato-Milton parameter ζ
1 for general shapes of inclusions can be reduced to a sum of the triple integrals (108)
∫
L
k
∫
L
k
1
∫
L
k
2
t
2
E
1
t
2
-
t
1
¯
E
1
t
1
-
t
d
t
¯
d
t
1
d
t
2
¯that should be easier than computation of the corresponding triple integral (20.66) from [9] on the 3-point correlation function which is computed by another triple integral. The integral (108) can be calculated by residues for algebraic curves.


## 5. Conclusion

The iterative scheme (77) was constructively applied to special inclusions. An exact formula for regular array of disks (and solution of the problem (40) for N
=
1) was obtained in [42]. The term “exact formula” should be explained precisely. Here, the exact formula means that Λ
e is written as an expression which contains the geometrical parameters (the fundamental translation vectors ω
j of the lattice Q, the radius of the disk) and the physical parameter ρ explicitly in symbolic form. The exact formula does not contain any parameter calculated by a numerical method (by an integral equation method or by truncation of an infinite system of equations). The exact formulae for Λ
e can be written in the form of the series similar to (19) in which all the coefficients C
n were exactly written.Approximate analytical formulae forΛ
e were obtained for arbitrary nonoverlapping locations of different disks in [34]. Here, the term “approximate analytical formula” means that the effective conductivity λ
e was written in the form of the series (19) in which the first few coefficients C
n were exactly written. Other results obtained by the GMS are discussed also in Section 4.5.In the previous works, the absolute convergence of the method was proved under geometrical restrictions [17]. Roughly speaking, these restrictions mean that each inclusion is far away from others (see [9] and many papers cited therein). However, it was proved in [10] that a modified method always uniformly converges; that is, these geometrical restrictions are redundant. This is an interesting example of the difference between absolute and uniform convergence which shows that estimations on the absolute values or on the norm are too strong in comparison to the study of the uniform convergence.The inclusionsD
k can have exotic shapes that complicate the criterion of the dilute composite. This difficulty can be overcome by a conformal mapping of D onto a circular multiply connected domain. The capacity [43] (effective conductivity) is invariant under conformal mappings. Next, application of the functional equations for circular inclusions of different radii [34] can be applied to estimate their interactions.The GMS is effective in symbolic and numerical computations for nonoverlapping inclusions of simple shapes (disks and ellipses) in the deterministic statement. When an approximate analytical formula is deduced, the ensemble average is calculated by the straightforward computations. It is interesting to combine numerical methods effective for one-inclusion problems with the analytically described interactions between inclusions by the GSM in Section4.1.



---

*Source: 535128-2015-12-21.xml*