Random 2D Composites and the Generalized Method of Schwarz

. 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 space C 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 𝑛 -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.


Introduction
The effective properties of random heterogeneous materials and methods of their computation are of considerable interest [1][2][3][4][5][6][7][8][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][12][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 = ( 1 ,  2 ,  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  ,  = 1, 2, . ..} be given, where each D  has a fixed geometry.Let all the particles be randomly located in the space and each particle D  occupies a domain   without deformations.Thus, the deterministic elements D  are introduced independently but the joint set is introduced { 1 ,  2 , . ..} randomly.The diversity of random locations is expressed by joint probabilistic distributions of the nonoverlapping domains   .
Various methods were applied to such random composites.The most known theoretical approach is based on the -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 -point correlation functions for large .Let a distribution of inclusions be fixed.Then, the effective conductivity tensor Λ  is uniquely determined through D  and their physical properties.In order to estimate Λ  , 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 space C 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 -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 (] 20 ) in concentration ] (exactly written in Section 4.6 up to (] 8 ) and up to (] 20 ) [21]).Constructive formulae of the traditional approach [9] are deduced only up to (] 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, , is arbitrary taken, and the final formulae for the effective tensor contain  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 -particle problem in the infinite plane can be applied only to diluted distributions of the clusters consisting of  inclusions.
This paper is the first part of the work devoted to the GMS described in general form in Section 4.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.

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  is defined by the vertices ± 1 /2 and ± 2 /2.Without loss of generality the area of  can be normalized to one; hence, The points  1  1 +  2  2 ( 1 ,  2 ∈ Z) generate a doubly periodic lattice Q (torus), where Z stands for the set of integer numbers.
Consider  nonoverlapping simply connected domains   in the cell  with Lyapunov's boundaries   and the multiply connected domain  =  \ ⋃  =1 (  ∪   ), the complement of all the closures of   to .Each curve   leaves   to the left (see Figure 1).Let  denote a complex variable and a point   is arbitrary fixed in   ( = 1, 2, . . ., ).
We study conductivity of the doubly periodic composite when the host  +  1  1 +  2  2 and the inclusions   +  1  1 +  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 Here, x = ( where / denotes the outward normal derivative to   .The external field is modelled by the quasi periodicity conditions where  ℓ (ℓ = 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  satisfies E = −∇. ( The linear constitutive relation is fulfilled in  and   : Problem ( 3)-( 4) can be stated in the probabilistic context with randomly located inclusions.For definiteness, consider a set of the domains   as a realization of the random locations of a collection of particles D  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 (a) associated with finding a configuration of inclusions with position a := ( 1 ,  2 , . . .,   ).

Advances in Mathematical Physics
The point   determines the position of the particle D  in the cell .The set of all configurations is denoted by C. Let d  denote a small area element about point   and da = d 1 d 2 ⋅ ⋅ ⋅ d  .Then (a)da is equal to the probability of finding  1 in d 1 ,  2 in d 2 , . . .,   in d  .
The ensemble average is introduced as the expectation in a probabilistic space The effective conductivity tensor of deterministic composites can be found in terms of the solutions of two independent problems (3)- (4).For instance, one can take Re  2 and  1 = 0,  2 = Im  2 .In the first case, Here,  denotes the contrast parameter introduced by Bergman: The effective conductivity tensor Λ  of stochastic composites is determined by application of the operator (7) to the deterministic values (8).Let the probabilistic distribution of the domains   ( = 1, 2, . . ., ) in  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: where the deterministic value  11 (a) is calculated by (9).

Traditional Approach
As it is stressed in Introduction we apply the straightforward approach to determine   .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.

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 field E 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]): Here, M 1 (x; ) is the single-inclusion operator which accounts for the first order interactions over and above E 0 and inclusions periodically equivalent to the th inclusion, that is, to   +  1  1 +  2  2 for all integers  1 and  2 .
The operator M 1 (x; ) 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   in the infinite plane is solved and the single-inclusion operator K 1 (x; ) in the plane is constructed (see, e.g., formulae (17.4) and (19.6) from [9]).Next, M 1 (x; ) is constructed via K 1 (x; ) by a periodicity operator.The periodicity transformation P : K 1 (x; )  → M 1 (x; ) can be easily performed since K 1 (x;   ) is expressed via the dipole tensor.The periodicity operator transforms rational functions in  =  1 +  2 to elliptic functions [28]; for instance, P :  −1  → () and P :  −2  → ℘().The second method to construct M 1 (x; ) is based on the same scheme realized in the torus topology.The operator K 1 (x; ), hence M 1 (x; ), 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   in terms of the conformal mapping of   onto the unit disk.In general, K  (x; ,  1 , . . .,   ) and M  (x; ,  1 , . . .,   ) are the -inclusion operator which accounts for the th order interactions.
The effective conductivity can be found through the polarization field defined by (19.3)-(19.4)from [9]: Following the traditional approach [9] we introduce the "double" average operator over the unit cell  and over the configuration space C: where dx = d 1 d 2 and da = d 1 d 2 ⋅ ⋅ ⋅ d  .The effective tensor ( 8) can be defined by 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 where   (a) is a complicated functional of the -inclusion cluster operators M  (x; ,  1 , . . .,   ) and the probability density (a).It is worth noting that   (a) is derived through K  (x; ,  1 , . . .,   ) in [9].Actually, our periodic approach yields the same result by rearrangement of the terms in each   (a).

Contrast Expansions.
Exact contrast expansions can be presented by formula (20.1) from [9] which after some modifications becomes 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   ( = 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) where G(x, x  ) denotes the periodic Green's function.Substitution of ( 2) into (20) yields The method of successive approximations applied to (21) yields a series on the contrast parameter  which can be shortly written in the form 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 Λ  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: The -point coefficients  (1)   are exactly expressed in terms of the -point correlation functions ( = 2, 3, . . ., ) (see formulae (20.38)-(20.41)and (20.59)-(20.60)from [9]).Formula ( 23) can be written in the form (19).
The coefficients  (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).

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  simple domains, the cell  with one-inclusion   ( = 1, 2, . . ., ).In each step of the algorithm the boundary conditions for the th simple domain are corrected by influence of the th simple domains ( ̸ = ) 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 where  =  0 +  ext is the decomposition of the potential  onto the regular periodic part  0 and the external field  ext corresponding to the quasi periodicity conditions (4).
In the case of the space problem the external field  ext has a singularity at infinity.Introduce the function and the domains  + := ⋃  =1   ,  − := .It is convenient to represent the harmonic functions in different domains as one piecewise harmonic function  0 (x) introduced above in  =  − and equal to   (x) −  ext (x) in   ∪   ( = 1, 2, . . ., ).Then, (24) can be considered as the jump problem: where  + 0 and  − 0 correspond to the limit values of (x) when x tends to  from  + and  − , respectively.The jump problem (26) has a unique solution expressed in terms of the simple layer potential  for the curve  in the cell  in the torus topology.We have The operator  is considered as an operator in an appropriate functional space [29][30][31][32].In the classic statement,  should be replaced by the whole space.The operator  is decomposed onto the simple layer potentials   along curves   ( = 1, 2, . . ., ) as  = (1/( + 1)) ∑  =1   .Here, the multiplier 1/( + 1) is introduced for short form of the equations below.Then (27) implies that x ∈   ( = 1, 2, . . ., ) , Equations ( 28) can be considered as a system of integral equations on the potentials   (x) in the inclusions   ( = 1, 2, . . ., ).If it is solved, the potential (x) is calculated in  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): where w is a fixed point in .
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 in R  ( = 2, 3).
Let tensor fields U() in  and U  () in   satisfy the contact conditions where T and T  are boundary operators involving physical parameters of materials occupying the domains  and   , respectively.Let T  be presented in the form T  = T −   , where   denotes a contrast parameter tensor.Then (32) becomes Let P  denote the simple layer potential corresponding to the considered linear equation for the domains   and  −  separated by   .Application of ∑  =1 P  to (33) yields x ∈   ∪   ( = 1, 2, . . ., ) . ( 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 x ∈   ∪   ( = 1, 2, . . ., ) . ( One-inclusion problem for   ( = 1, 2, . . ., ) 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].

Integral Equations.
In the present section, we discuss the GMS for the double periodic problem (3)- (4).Following [34] we introduce the complex potentials () and   () in such a way that where x = ( 1 ,  2 ) and  =  1 +  2 .Two real equations (3) can be written as one R-linear condition: where the bar denotes the complex conjugation.The real part of (37) yields the first equation in (3).The imaginary part yields Im () = (1 + ) Im   ().After differentiation of the latter equality along   and application of the Cauchy-Riemann equations we arrive at the second equation in (3) (see for details [10]).
The unknown functions () and   () are analytic in  and   , respectively, and continuously differentiable in the closures of the domains considered.It follows from ( 4) that the function () is quasiperiodic: where the constants  1 and  2 are taken from (4) and  1 and  2 are undetermined real constants which should be found during solution to the problem.
The following formula was proved in [10]: Here, the differential operator is understood as / =    (/), where  = () is a parametric equation of   and    () is the boundary value of the complex derivative / = / 1 − (/ 2 ).Differentiation of (37) and application of (39) yield where () =   (),   () =    (), and   () stands for the unit normal outward vector to   expressed in terms of the complex values.The function () is doubly periodic.The Rlinear problem (40) holds also in the limit cases for perfectly conducting inclusions ( = 1) and for isolators ( = −1).
The R-linear problem (37)-( 38) can be reduced to a system of integral equation in the following way.The cell  in the torus topology is divided into two domains  + = ⋃  =1   (not connected) and  − =  (multiply connected).Let  = ⋃  =1   denote the boundary of  + and let () be a Hölder continuous function on .Let () denote the Weierstrass function for which [28]  ( +   ) −  () =   ( = 1, 2) , (41) where   = 2(  /2).Following [35] introduce the Eisenstein function where the Eisenstein summation is used [35].The Weierstrass and Eisenstein functions are related by formula  1 () = ()−  2 , where  2 = (2/ 1 )( 1 /2).It follows from Legendre's identity  1  2 −  2  1 = 2 (see [28,35]) and ( 41) that the jumps of  1 () have the form 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 where () is a Hölder continuous function on  and  and  0 are complex constants.The function Φ() is analytic in  ± and quasiperiodic; that is, it has constant jumps per periodicity cell calculated with the help of (43): The following Sochocki (Sokhotski-Plemelj) formulae take place on torus [36]: where Φ ± () denote the limit values of Φ() on  when  ∈  ± tends to  ∈ , respectively.The principal value of the integral in ( 46) is considered following [37].Formulae (46) imply that the function Φ() satisfies the jump condition Equation ( 47) can be considered as a C-linear conjugation problem [10] in a class of quasiperiodic functions with the given jump ().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 () =   () on , Φ() =   () in   ⊂  + , and Φ() = () in  =  − .Then, (44) can be written in  + as follows: 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: 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  is considered (0 <  ≤ 1).This is a Banach space endowed with the norm: Consider the closed subspace A ⊂ C (0,) consisting of functions analytically continued into all   .Therefore, every element of A is a function analytic in the domain ⋃  =1   and Hölder continuous in its closure.It is worth noting that the domain ⋃  =1   is not connected and convergence in A implies the uniform convergence of functions in ⋃  =1 (  ∪   ).
Equations (50) are written in the domains   and can be continued to the boundary by use of ( 46) 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  0 .The constants  1 +  1 ,  2 +  2 , and  0 disappear in the problem (40) after differentiation.The structure of the general solution of (40) has the form () =  1  (1) () +  2  (2) (), where the real constants  1 and  2 correspond to the external field (jumps of (x) per a cell).Hence, the Rlinear 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,  = 1 and  = .Therefore, the complex constant  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 .For simplicity, we put Then the external field corresponds to the complex potential  (ext) () =  and to the complex flux  (ext) () = 1.In the vector form, this external flux becomes (−1, 0).It follows from the first relation (45) that where in this relation and calculate The integral from the left hand part of (56) expresses the total flux through the low side of the parallelogram  parallel to the  1 -axis.It must be equal to zero because of the periodicity and that the external flux (−1, 0) vanishes in the  2 -direction.Therefore,  1 = 0; hence  = 1 by (55).

Convergence of the GMS.
This section is devoted to study of the integral equations ( 50), (53) in the space A.
Let Ψ() =   () for  ∈   ∪   and let  denote the operator from the right hand part of (48), Then, (48), Advances in Mathematical Physics 9 (53) can be shortly written as the following equation in the space A: (59) In the present section, we demonstrate that (59) has the unique solution for || < 1 given by the series We use the following general result of functional analysis.
Theorem 4 (see [38]).Let  be a linear bounded operator in a Banach space B. If for any element  ∈ B and for any complex number ] satisfying the inequality |]| ≤ 1 equation has a unique solution, then the unique solution of can be found by the method of successive approximations.The approximations converge in B to the solution Theorem 4 can be applied to (59) in the Banach space A.
Theorem 5. Let || < 1. Equation Ψ = Ψ + 1 has a unique solution.This solution can be found by the method of successive approximations convergent in the space A.
Uniform convergence of (76) for || = 1 can be justified by use of the arguments form [10].

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 ⋃  =1 (  ∪   ).The method of successive approximations can be also presented in the form of the following iterative scheme: When the functions   () ( = 1, 2, . . ., ) are found, the effective conductivity is calculated by (57).
In order to compare ( 21) and ( 48), (53) we introduce the complex intensity field () =  1 () −  2 () isomorphic to the vector field E(x) = ( 1 (x),  2 (x)).Then, the vector equation ( 21) can be written as a complex equation on ().For definiteness, we consider the complex potential   () in   .It follows from ( 5), (36) and   () =    () that Therefore, the complex equation on () in   can be written as (48) on   ().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.

Cluster
One can check that the operator   is singular and   ( ̸ = ) are compact in A  .Equation (79) can be shortly written in the form The (zero)th approximation in the concentration ] for (79) can be written as the following integral equation: Equation ( 82) corresponds to the R-linear problem (40) for one-inclusion   in the cell .Solution to the oneinclusion problem (82) defines the inverse operator ( +   ) −1 bounded in A  .Then, (81) can be written in the equivalent form As it is established in Section 4.4 the vector-function E(x) and the complex functions   () express the same vector field in   up to an isomorphism for E 0 = (−1, 0).Therefore, application of the successive approximations to (83) yields the same series as (12) in   .A detailed comparison can be performed.It shows, for instance, that the term E 0 + ∑  =1 M 1 (x; )⋅E 0 from (12) corresponds to (+  ) −1 1 from (83).
The solution of the problem (82) depends on the shape of   .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  (0)  () do not depend on .Then, we arrive at the formula valid up to (] 2 ): where the shape factor  is introduced following [23,40]: The Clausius-Mossotti formula equivalent to (84) up to (] 2 ) is used in applications [6]: 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 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   for arbitrary nonoverlapping locations of different disks [24,34].An example of equal disks is discussed below in Section 4.6.
Consider the integral operator from (48): where () is analytic in the disk   .We have  =  * () on   .Hence, the operator (89) can be written in the form where the function ( * () ) is analytically continued into | −   | > .
The function  1 () has poles of first order at the points  =  1  1 +  2  2 (see [35]).The function (/( −   )) 2 ( * () ) is analytic in the domain | −   | >  and the integral (  )() can be calculated by residues at | −   | > .Introduce the operator We have for  = and for  ̸ = 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 where  1 ,  2 run over integers in the sum ∑ * with the excluded term  1 =  2 = 0 for  = .The functional equations (94) do not contain any integral.Each term ( ()  1 , 2   )() 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  1 ,  2 can be overcome by application of Eisenstein functions   ().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 ] =  2 which is equivalent to the expansion in  2 [24]: Few first functions in the expansion (95) are given by the following exact formulae: where   ( −   ) is defined as Ẽ ( −   ) =   ( −   ) − ( −   ) − for  =  and   (0) := Ẽ (0).It is possible to proceed (96) and calculate the next approximations following the algorithm [24].
Let  be a natural number;   runs over 1 to ,   = 2, 3, . . . .Let C denote the operator of complex conjugation.Introduce the multiple convolution sums The integrals from (57) are simplified by the mean value theorem for harmonic functions Substitution of (96) and the next terms into (98) yields where The next coefficients   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 The infinite set {ê  1 ⋅⋅⋅  ,   = 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   ( ≤ 20) within a given distribution of   in few hours on an ordinary notebook when  = 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: where the relation ê2 =  from [41] is used.Another formula deduced in [41] can be useful in simulations and estimations of the third-order term of ( 102 One can find the numerical values of ê 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 ( 4 ): The latter formula gives the Torquato-Milton parameter  1 (see (20.59) and (20.66) from [9]): The three-and four-point contrast bounds on the effective conductivity are written by this parameter (see (21.33) 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.

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  = 1) was obtained in [42].The term "exact formula" should be explained precisely.
Here, the exact formula means that Λ  is written as an expression which contains the geometrical parameters (the fundamental translation vectors   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 Λ  can be written in the form of the series similar to (19) in which all the coefficients   were exactly written.Approximate analytical formulae for Λ  were obtained for arbitrary nonoverlapping locations of different disks in [34].Here, the term "approximate analytical formula" means that the effective conductivity   was written in the form of the series (19) in which the first few coefficients   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 inclusions   can have exotic shapes that complicate the criterion of the dilute composite.This difficulty can be overcome by a conformal mapping of  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 Section 4.1.

Figure 1 :
Figure 1: Random double periodic location of the nonoverlapping domains   (gray disks) bounded by   ;  is the complement of all the closed domains   ∪   to the periodicity cell  (dashed parallelogram).