## Abstract

We develop a method of functional equations to derive analytical approximate formulae for the effective conductivity tensor of the two-dimensional composites with elliptical inclusions. The sizes, the locations and the orientations of the ellipses can be arbitrary. The analytical formulae contain all the above geometrical parameters in symbolic form.

## 1. Introduction

Analysis concerning the transport properties of inhomogeneous materials is of fundamental theoretical interest. It is also important in engineering, for instance in optimal design problems. Analytical formulae for the macroscopic properties with physical and geometrical parameters in symbolic form are useful to predict the behaviour of composites. Using the terminology of electric or thermal conductivity, we study the transport properties of a two-dimensional, two-phase composite with parallel elliptical cylinders of conductivity *λ* embedded in a host material of the normalized unit conductivity.

The effective conductivity of random conducting materials has been studied since the nineteenth century by Maxwell and Rayleigh. Most of the work in literature is concerned with dilute composites with spherical and ellipsoidal inclusions (circular and elliptical inclusions in two-dimensional composites), because the corresponding boundary value problem for one inclusion in the whole space can be solved in closed form. The attention of this study is on two-dimensional composites for which few analytical formulae are known.

Let *ν* be the area fraction of the inclusions, and
1.1
denotes the contrast parameter introduced by Bergman (1978). The semi-axes of the ellipses are denoted by *r*(1+*α*) and *r*(1−*α*), where 0<*α*<1 and *r*>0. If all the ellipses are aligned with the coordinate axes, the effective conductivity tensor is diagonal. Its diagonal components were estimated by Galeener (1971*a*–*c*):
1.2
Cohen *et al.* (1973) discussed what they viewed as inconsistencies arising from formulae (1.2) and derived formulae based on an elliptical Lorentz cavity having the same depolarization factors as the inclusions:
1.3
The scalar macroscopic conductivity of randomly distributed ellipses were calculated by Hetherington & Thorpe (1992) and by Zimmerman (1996):
1.4
Having applied the differential method, Zimmerman (1996) deduced another approximate formula
1.5

The following formula was obtained by Hetherington & Thorpe (1992) in the case of holes whose shape is an *n*-sided regular polygon
1.6
where *Γ* is a gamma function. Garboczi & Douglas (1996) investigated the influence of the various shapes of inclusions on the effective conductivity up to *O*(*ν*^{2}). Ammari *et al.* (2005) expressed the effective conductivity tensor *Λ* via the dipole matrix *M*
1.7
where *I* is the identity matrix. A method of integral equations was applied to compute the matrix *M*. Movchan *et al.* (2002, 2003) explicitly calculated *M* for the shapes described by algebraic equations. Crowdy (2008) applied the Schottky–Klein prime function to effectively compute the local field of finite sets of discs. Having used the dipole matrix method, Rylko (2000) deduced formulae for the principal conductivities of a dilute composite with curvilinear inclusions unidirectionally located on the plane
1.8
where the shape of the inclusion is defined by the function *r*=1+*ϵ**f*(*θ*) in polar coordinates (*r*,*θ*), *ϵ* is a small positive parameter, and
General methods to deduce formulae (1.2)–(1.5) and discussion of other approaches were presented by Milton (2002).

The summarized works are based mainly on the solution to the boundary value problems with one inclusion in the plane. It should be emphasized that the above-discussed formulae are valid only to the first order of *ν*. In general, to evaluate higher-order terms, we need to discuss boundary value problems with an infinite number of inclusions and to investigate at least the first-order interactions between inclusions. Such problems in a class of doubly periodic functions for composites with circular inclusions were systematically studied by Berlyand & Mityushev (2001, 2005), Mityushev (2001, 2005) and Mityushev *et al.* (2008). In particular, an exact formula for the effective conductivity of the square array of cylinders was constructed. Approximate analytical formulae for the effective conductivity tensor for any distribution of the discs were deduced.

The problem of determining the conductivity of randomly distributed ellipses has not yet been dealt with rigorously for non-dilute composites. As has already been noted, Maxwell’s approach is based on a single inclusion and can produce formulae valid only to the first order of *ν*. Such approximations take into account the shape of inclusions, but not the geometry of their location. To the best of our knowledge, there are only four works devoted to the rigorous study of composites with many elliptical inclusions with the precision at least *O*(*ν*^{2}). McPhedran & Nicorovici (1997) and Nicorovici & McPhedran (1996) applied Rayleigh’s method to a rectangular array of elliptical cylinders and reduced the problem to an infinite system of linear algebraic equations. Analytical solution to a truncated system for sufficiently small *α* yields
1.9
where *S*_{2} is a lattice sum. More precise formulae for not small *α* were also deduced. The later ones differ from equation (1.9) by a coefficient on *ν* in the denominator. Formulae (1.9) and others presented by McPhedran & Nicorovici (1997) and Nicorovici & McPhedran (1996) are rigorous, because the applied Rayleigh’s method takes into account both the shape of the ellipses and the geometry of the array, by means of the lattice sum.

Yardley *et al.* (1999, 2001) calculated the effective conductivity for a structure consisting of a finite number of layers, or stack, of ellipses by Rayleigh’s method. It is interesting to extend the results by McPhedran & Nicorovici (1997), Nicorovici & McPhedran (1996) and Yardley *et al.* (1999, 2001) to random media by using the method of functional equations that effectively works for circular inclusions.

In this paper, the first step has been done in this direction. We propose an iterative method to compute the effective conductivity tensor *Λ* for an arbitrary distribution of non-overlapping ellipses. It is based on symbolic computations and yields analytical formulae for the local field and for *Λ* in the form of a series in *ν*. The final formulae (4.34)–(4.37) for *Λ* are derived up to *O*(*ν*^{3}). One can proceed the iterations and obtain higher-order formulae. However, it requires advanced symbolic computations.

This paper is organized as follows. A list of the nomenclature is given in table 1. The problem is presented in §2, and the -problem that governs the local complex potential is derived in §2*a*. The problem is reduced to integral equations in §2*c*. The main part of the paper, solution to the integral equations, is presented in §3. Applications to the calculation of the effective conductivity tensor are presented in §4. New analytical formulae up to *O*(*ν*^{3}) are summarized in §4*c*–*f*. They are universal for any distribution of the centres of inclusions, of their sizes and of orientations. This analytical dependence of *Λ* on locations of elliptical inclusions was not known before. Some concluding remarks are given in §5.

## 2. Statement of the problem and reduction to integral equations

### (a) Geometry

Non-overlapping elliptical inclusions *D*_{m} (*m*=1,2,…,*n*) are embedded in the plane (*x*,*y*). It is convenient to put the semi-axes equal to *r*_{m}(1+*α*) and *r*_{m}(1−*α*), respectively. The positive parameter *r*_{m} characterizes the size of the inclusion, and *α* the shape of the ellipse (0<*α*<1). The limit case when *α* is equal to 0 can be treated in the final formulae by calculating the corresponding limits. The limit case *α*=1 requires a separate investigation that is not presented in this paper. Let an inclusion *D*_{m} be centred at (*x*_{m},*y*_{m}) and the angle between the major semi-axis of the ellipse and the *x*-axis is equal to *θ*_{m}. Introduce the local coordinates (*X*,*Y*) for a fixed inclusion *D*_{m} as follows:
2.1
The local equation of the ellipse ∂*D*_{m}, i.e. the boundary of *D*_{m}, has the form
2.2
The foci of the ellipse ∂*D*_{m} in the local coordinates are located at .

Introduce the local and global complex coordinates *Z*=*X*+i*Y* and *z*=*x*+i*y*, where i denotes the imaginary unit. The Joukowsky conformal mapping
2.3
transforms the annulus onto *D*_{m}−*Γ*_{m}, where *Γ*_{m} is the slit along the *X*-axis. The inverse mapping to equation (2.3) has the form
2.4
where the branch of the square root is chosen in such a way that
2.5
Formulae (2.3) and (2.4) in the global coordinates become
2.6
and
2.7
where *a*_{m}=*x*_{m}+i*y*_{m} and *s*_{m}=*r*_{m}e^{iθm}.

### (b) Conjugation conditions

Let *D* denote the complement of the closures of all domains *D*_{m} to the plane. It is convenient to assume that infinity belongs to *D* as it is usually done in complex analysis. We study the conductivity of the two-dimensional composite, when the domains *D* and *D*_{m} are occupied by materials of unit and *λ* conductivity, respectively, where . Then, the potentials *u*(*z*) and *u*_{m}(*z*) are harmonic in *D* and *D*_{m} (*m*=1,2,…,*n*) and satisfy the conjugation (transmission) conditions
2.8
where ∂/∂*n* is the outward normal derivative to the ellipses. For simplicity, it is assumed that the potential *u*(*z*) has singularities only in the domain *D* described by a function Re*f*(*z*), where *f*(*z*) is analytic in all the inclusions *D*_{k} and Re stands for the real part of a complex number. For instance, if *u*(*z*)∼*x*=Re*z* at infinity, i.e.
2.9
the external flux is applied parallel to the *x*-axis.

Following Mityushev & Rogosin (2000), introduce complex potentials *φ*(*z*) and *φ*_{m}(*z*) analytic in *D* and *D*_{m}, respectively, in such a way that *u*(*z*) and *u*_{m}(*z*) are related to the complex potentials, respectively, by
2.10
Then the conditions (2.8) can be reduced to the following -linear problem
2.11
where the overbar stands for the complex conjugation.

### (c) Integral equations

In this section, we reduce the problem (2.11) to integral equations. First, recall the Sokhotskij–Plemelj formulae. The curve divides the complex plane onto two domains *D*^{+} and *D*. Here, each curve ∂*D*_{k} is orientated in the clockwise direction. Let μ(*t*) be a Hölder continuous function on *L*. Then the function
2.12
is analytic and continuous on the complex plane except *L* where its limit boundary values and satisfy the jump condition
2.13
The relation (2.11) can be written in the form (2.13) with *Φ*^{+}(*t*)=*φ*_{k}(*t*)−*f*(*t*), *Φ*^{−}(*t*)=*φ*(*t*) and . Then equation (2.11) yields
2.14
and
2.15
One can consider equation (2.14) as a system of integral equations with respect to *φ*_{k}(*z*) analytic in *D*_{k} and continuously differentiable in its closure. When *φ*_{k}(*z*) are found, *φ*(*z*) is calculated by equation (2.15). It is worth noting that equations (2.14) are not the classic integral equations of the potential theory. They correspond to integral equations that can be deduced from the Schwarz alternating method due to Mityushev & Rogosin (2000).

Following Mityushev & Rogosin (2000), one can prove that the system of integral equations (2.14) has a unique solution up to an additive constant when |*ρ*|<1 that corresponds to the physical meaning of the problem (see equation (1.1)). This solution can be found by the method of successive approximations uniformly convergent in any compact subset of *D*^{+}. This theoretical result will be proved in detail in a separate paper. In this paper, attention is paid to finding a constructive solution to equations (2.14).

## 3. Solution to integral equations

Integral equations (2.14) are deduced for general shape of the inclusions *D*_{k} and can be numerically solved, for instance, by the boundary element method due to Lifanov (1996). This method and other numerical methods are effective for numerically given geometries and data. However, if one wishes to preserve the contrast parameter *ρ*, the concentration *ν* and parameters of the shape *D*_{k} in symbolic form, then he should apply symbolic computations.

In this section, equations (2.14) are solved for the elliptical shape of inclusions with the parameters *ρ*, *α*, *θ*_{k}, *a*_{k} and *n* in symbolic form. First, the integral equations (2.14) on the complex plane *z* are transformed to equations on the complex plane *w*. The latter equations are reduced to functional equations. A simple constructive algorithm to solve the functional equations in symbolic form is described.

### (a) Transformation of integral equations

Let *k* be fixed in equations (2.14). The doubly connected domain *D*_{k}−*Γ*_{k} is mapped onto the annulus by the conformal mapping (2.7); *D*_{k} is transformed onto the unit circle |*w*|=1, and *Γ*_{k} onto the circle . Introduce the functions
3.1
analytic in and continuous in . Substitute equation (3.1) in equations (2.14) and change the variables in the integrals as follows:
3.2
Then equations (2.14) become
3.3
where *F*(*w*)=*f*(*z*).

Moreover, it follows from the continuity of *φ*_{k}(*z*) when *z* passes the slit *Γ*_{k} that
3.4
Equation (3.4) implies that *Φ*_{k}(*w*) is represented in the form
3.5
where *ϕ*_{k}(*w*) is analytic in the unit disc |*w*|<1. Equation (3.5) follows from the representation of *Φ*_{k}(*w*) in the form of the Laurent series in the annulus and form the application of equation (3.4). The same arguments yield the representation of *F*(*w*) in the form *F*(*w*)=*g*_{k}(*w*)+*g*_{k}(*α*/*w*), where *g*_{k}(*w*) is analytic in the unit disc. Substitution of equation (3.5) into equation (3.3) yields
3.6
where
3.7
and
3.8
Here, the relation on the unit circle is used.

### (b) Calculation of the integrals

In this section, the integrals (3.7) and (3.8) are analytically calculated by residua. To study singularities of the integrands, we investigate the roots of the denominator. We have the quadratic equation with respect to τ,
3.9
The cases of equal and non-equal *k* and *m* have to be separately investigated.

Let

*k*=*m*. Then equation (3.9) becomes 3.10 Its two solutions have the form 3.11 One can see τ_{1}and τ_{2}lie in the unit disc, as .Let

*k*≠*m*. To avoid confusion with equation (3.11), the roots of equation (3.9) in this case are denoted by*w*_{1}and*w*_{2}: 3.12 The branch of the square root is chosen in accordance with equation (2.5). We will prove that |*w*_{1}|<1 and |*w*_{2}|>1. Consider equation (3.9) in the domain |τ|≥1. Let denote the complement of*D*_{m}to the extended complex plane*z*. Changing the variable τ by*t*(see equation (3.2)), we obtain 3.13 with respect to , where the parameter 3.14 belongs to*D*_{k}. Therefore, equation (3.13) has the unique solution , as (*k*≠*m*). Hence, equation (3.9) has a unique solution in the domain |τ|≥1. The second root*w*_{2}from equation (3.12) in the variables τ and*w*corresponds to the root*t*=*z*in the variables*t*and*z*. Therefore, |*w*_{2}|≥1. It follows from Viète’s formulae that*w*_{1}*w*_{2}=*α*. This implies that |*w*_{1}|<1, as*α*<1.

We now proceed to analytically calculate the integral (3.7). For *m*≠*k*, it can be written in the form
3.15
where *w*_{1} and *w*_{2} have the form (3.12). Calculate *P*_{km}(*w*) by residua using the analyticity of in |τ|>1,
3.16
where the relation *w*_{1}*w*_{2}=*α* is used. Similar arguments yield formulae
3.17
Here, the following formula for *Φ*^{−}(*w*) analytic in |*w*|>1 is used due to Gakhov (1966):

### (c) Functional equations

The integrals (3.7) and (3.8) have been exactly calculated in the previous section. Substituting the results (3.16) and (3.17) into equation (3.6), we transform the integral equations (3.6) to the following functional equations:
3.18
Here, for convenience the root *w*_{1} is written as the function of *w*,
3.19
where
3.20

The left-hand part of equation (3.18) consists of the functions *ϕ*_{k}(*w*) and *ϕ*_{k}(*α*/*w*) analytic in |*w*|<1 and , respectively. Denote by *P*^{+} the project operator that transforms a function analytic in to its part analytic in the unit disc. This operator can be considered as taking the regular part of the Laurent series or as the integral operator with |ζ|<1. Application of *P*^{+} to equation (3.18) yields
3.21
Here, the relation is used.

If the external field is applied parallel to the *x*-axis, then equation (2.9) implies that *F*(*w*)=*z*=*s*_{k}(*w*+*α*/*w*)+*a*_{k}, hence
3.22
Hereafter, equations (3.21) are considered with the function (3.22).

One can consider equations (3.21) as a system of functional equations with respect to *ϕ*_{k}(*w*) analytic in the unit disc and continuous in its closure. Despite equations (3.21) involving the operator *P*^{+}, they are more convenient in symbolic computations than equations (3.18). As stated at the end of §2*c*, the system (3.21) has a unique solution. This solution can be found by the method of successive approximations. Equations (3.21) can be considered as iterative functional equations with shift into domain (Kuczma *et al.* 1990), as |*β*_{km}(*w*)|<1. The later inequality is equivalent to the relation |*w*_{1}|<1 justified in §3*b*.

The special case of *α*=0 can be investigated by the limit *α*→0 when the ellipses become circles. In this case we arrive at the functional equations discussed in Mityushev (1993). Applications of such equations and their modifications to determining the effective conductivity tensor of a composite containing circular inclusions were presented by Berlyand & Mityushev (2001, 2005) and Mityushev *et al.* (2008). Discussion of dilute regime and regular locations of circular inclusions was presented by Milton (2002), McPhedran (1986) and McPhedran & Milton (1987) by application of Rayleigh’s method.

## 4. Effective conductivity tensor *Λ*

### (a) General formulae for *Λ*

The effective conductivity tensor
4.1
can be calculated by the Maxwell approach. Following Mityushev (1999, 2005), one can obtain formula for dilute composites
4.2
where *I*_{k} denotes the double integral
4.3
where *φ*_{k}(*z*) is a solution of equation (2.11) with *f*(*z*)=*z*. The component *λ*^{y} can be calculated by the formula due to Mityushev (2001):
4.4
where the double integral *I*_{k} is calculated with *φ*_{k}(*z*), a solution of equation (2.11) with *f*(*z*)=−i*z*. The later corresponds to the external field applied in the *y*-direction.

The integral *I*_{k} can be transformed by Green’s formula written in complex form
4.5
where is the area of *D*_{k}. The integral (4.5) can be calculated by the change of variables (see equation (3.2)). Then
Using the relation on the unit circle, we obtain
Substituting equations (3.1) and (3.5) into equation (4.5) and calculating the integrals by the residua (at zero and at infinity), we arrive at the simple formula
4.6
Therefore, equation (4.2) becomes
4.7
The value *ϕ*′_{k}(0) can be computed after solution to functional equations (3.21). It is noted at the end of §3*b* that these equations can be solved by the method of successive approximations. Differentiating equations (3.21), multiplying the result by and introducing the function , then it implies that
4.8
Here, the commutativity of the operators *P*^{+} and *d*/d*w* is used. Formula (4.7) becomes
4.9
Therefore, to determine *λ*^{x} and *λ*^{xy}, one can solve the functional equation (4.8), compute *ψ*_{k}(0) and substitute it in equation (4.9).

### (b) Algorithm to compute *Λ* in symbolic form

Introducing the dimensionless parameter characterizes the ratio of the line sizes of inclusions to the distances between them. To describe the structure of the functional equation (4.8), we shortly write it in the form
4.10
where *A*_{km} is a bounded linear operator. It will be clear later why the power ℓ^{2} is taken as a coefficient. Equation (4.10) can be solved by the following two iterative schemes. First, the direct iterations can be applied to it to obtain
4.11
where is the *p*th approximation of *ψ*_{k}(*w*). As it is noted at the end of §3*b*, the iterations (4.11) always converge. The rate of the convergence is equal to |*ρ*|. Even in the case |*ρ*|=1, the method (4.11) can be modified to be convergent as it was done by Mityushev & Rogosin (2000).

The second iterative scheme is constructed on the basis of the parameter ℓ^{2}:
4.12
Therefore, at each step the functional equation has to be solved
4.13
Here, the function is expressed through computed at the previous step
Mityushev & Rogosin (2000) proved that equation (4.13) has a unique solution represented in the form of the absolutely convergent series
4.14
The iterative scheme (4.12) corresponds to the generalized alternating Schwarz method discussed by Mityushev & Rogosin (2000).

The zero approximation of equation (4.12) can be easily obtained from equation (4.14): 4.15 Substitution of equation (4.15) into equation (4.9) yields the Maxwell approximation 4.16 Averaging over the orientation of the inclusions yields formula (1.4).

Higher approximations can be obtained by numerical and symbolic computations. It is worth noting that the realization of the iterative schemes (4.11) and (4.12) involves computations of the compositions of functions and of the regular parts of the Laurent series. There is not any integral that requires numerical computations that rather could not be computed in symbolic form.

### (c) First-order approximation for general distribution of inclusions

We now proceed to determine via given by equation (4.15) from the functional equation
4.17
where
4.18
The value can be immediately determined by substitution of *w*=0 into equation (4.17):
4.19
Therefore, to determine , it is sufficient to calculate (*P*^{+}*β*′_{km})(0). It is convenient to calculate (*α*+1/*α*)(*P*^{+}*β*′_{km})(0) in the form of the expansion on *α*,
4.20
This formula explicitly shows that the value is of order ℓ^{2}. Further, only the terms to *O*(ℓ^{2}) are kept in equation (4.20) and in the following computations. Substitution of equation (4.20) into equation (4.18) and of equation (4.19) into equation (4.9) yield
4.21
where
4.22
Formula (4.21) is valid up to *O*(*ν*^{2}) and takes into account first-order interactions between the inclusions.

### (d) Parallel inclusions

Consider equal elliptical inclusions with the parallel semi-axes. Then *θ*_{k}=*θ* and *r*_{k}=*r* for all *k*=1,2,…,*n* and equations (4.22) and (4.21) are simplified, respectively, to
4.23
and
4.24

Consider the case *θ*=0. Introduce the sum
4.25
Then the components of the effective conductivity tensor take the form
4.26
where Re and Im stand for the real and the imaginary parts of complex values, respectively. To calculate *λ*^{y}, it is sufficient to replace *a*_{k} by i*a*_{k}. Then in accordance with equation (4.4), we obtain
4.27

### (e) Randomly distributed inclusions

Consider equal elliptical inclusions randomly distributed on the plane. More precisely, it is assumed that all *r*_{k}=*r* are fixed and each angle *θ*_{k} is a random value uniformly distributed on [0,*π*]. Moreover, the centres *a*_{k} are independent identically distributed random variables such that the random vector (*a*_{1},*a*_{2},…,*a*_{n}) obeys a non-overlapping distribution for ellipses. The distribution of the centres and of the angles are mutually independent. Then the averages over *θ*_{k} from equation (4.21) vanish
4.28
and equation (4.21) becomes
4.29

Formula (4.29) is a generalization of the Maxwell formula (1.4). The term describes the first-order interactions of the inclusions.

### (f) Infinite number of inclusions

The most interesting begins when *n* tends to infinity. Then the first-order term in ℓ^{2} can be transformed into the first-order approximation in the concentration *ν*. The following limit of equation (4.25) plays the main role in the investigation
4.30
The formal pass to an infinite sum transforms equation (4.30) to a conditionally convergent series. To treat this series, fix an order of |*a*_{k}|. Let 0<|*a*_{1}|≤|*a*_{2}|≤⋯≤|*a*_{k}|≤…. The limit (4.30) was written via the absolutely convergent sums by Mityushev (1999)
4.31
where *m*≠*k* in the sum . The second limit from equation (4.31) is well defined, if it does not depend on the choice of point *z** and on the curve along which .

Mityushev (1999) proved that the first limit from equation (4.31) exists, if the area fraction
4.32
exists. Here *A*(*R*) is the area of inclusions located in a disc of radius *R*. The limit (4.32) has to be independent on the centre of the large disc. Hence, the first term from equation (4.32) is defined as the average of the absolutely convergent sums . It is forbidden to interchange the limit and the sum , because, in general, this produces the conditionally convergent series . Therefore, formula (4.31) cannot be simplified in general cases. The limits from equation (4.31) exist under the above described conditions on the distribution of *a*_{k}. This is not surprising because the composite discussed may not be homogenized.

The above theoretical discussion concerning is simplified, if we consider a doubly periodic distribution of *a*_{k}. Then the considered composite is homogenized and the limit (4.31) always exists. Let *N* points *a*_{k} (*k*=1,2,…,*N*) are placed in a rectangular periodicity cell *Q* with the sides γ and γ^{−1}. The area of the cell holds 1. All other points *a*_{k} are obtained from the first *N* points by translations on the vectors γ and iγ^{−1} expressed in the form of complex numbers. Further, it is convenient to apply to the limit (4.30) the Eisenstein summation (for detail, see Mityushev *et al.* 2008). Then equation (4.30) becomes
4.33
where *E*_{2}(*z*) denotes the Eisenstein function, coincides with *E*_{2}(*z*) for *z*≠0 and . Here, ζ is the Weierstrass ζ-function corresponding to the periodicity cell *Q*. In the case *N*=1, we have , where *S*_{2} stands for the standard lattice sum when all points *a*_{k} are produced from *z*=0 by translations on the fundamental lattice vectors. In the case of the square cell, .

Substitution of equation (4.33) instead of into equation (4.26) and simple transformations up to *O*(*ν*^{3}) yield
4.34
Using relations between the sums calculated in the perpendicular direction, we obtain
4.35
In the case *N*=1, we have . The value *S*_{2}=2γ^{−1}ζ(γ/2) is a real number for a rectangular array (see Rylko (2000)). For small *α*, formulae (4.34) and (4.35) are coincided to equations (1.9).

Along similar lines, equation (4.29) implies that
4.36
for randomly distributed ellipses with major semi-axes parallel to the *x*-axis. For the square array, *N*=1 and *S*_{2}=*π*. Then equation (4.36) becomes
4.37

A discrepancy between equations (1.9) and (4.34) with *N*=1 can be explained in terms of the iterative schemes (4.11) and (4.12). Rayleigh’s method is based on representations of the unknown functions and on the addition theorems that produce an infinite system of linear algebraic equations. Analysis of Rayleigh’s infinite system for circular inclusions made by Rylko (2000) shows that Rayleigh’s method corresponds to the iterative scheme (4.11) with *α*=0, i.e., Rayleigh’s method presented in continuous form via functional equations yields a series on *ρ*. In this section, we use the iterative scheme (4.12) that produces, in general, other approximate formulae.

## 5. Conclusion and discussion

Exact and approximate formulae for the effective conductivity tensor *Λ* of a two-dimensional medium containing elliptical inclusions are important in many applications. In this paper, we specialize the terminology of electric or thermal conductivity to make the discussion more concrete.

The Maxwell approach concerning the dilute regime yields various approximate formulae for *Λ* presented in §1. These formulae are valid up to *O*(*ν*^{2}), where each inclusion independently influences the macroscopic properties of the medium. According to the homogenization theory, higher approximations require the study of a medium with infinitely many inclusions, in particular, periodically arranged. Such higher approximations would involve interactions between inclusions written via the lattice sums, the correlation functions, etc.

In this paper, the method of the functional equation is developed to boundary value problems with elliptical inclusions. The investigation is confined to the case when the first-order interactions between inclusions are taken into account. New analytical approximate formulae for *Λ* are derived. These formulae take into account the shape, the sizes and the location of the inclusions. The influence of the geometry is expressed by the sum (4.25), by the limit (4.30) and by the sum (4.33) for doubly periodic composites. The final formulae (4.26), (4.27) and (4.34)–(4.37) for *Λ* are valid up to *O*(*ν*^{3}). Formulae (4.21) and (4.22) are universal for any distribution of the centres of inclusions. One can directly consider them from the statistical point of view as is done in §4*e*. Formulae (4.34)–(4.37) can also be treated as a statistical result. Theoretically, this procedure is equivalent to estimate the pair correlation function *P*(*z*|0) discussed by Milton (2002). From a practical point of view, it is convenient instead of *P*(*z*|0) (or of given by equation (4.30)) to estimate by equation (4.33) for sufficiently large *N*. It is also possible to investigate the expression in the right-hand side of equation (4.33) and formulae involving this expression for any *N* in symbolic form. After, one can put *N* infinitely large as it was done by Berlyand & Mityushev (2001, 2005) for circular inclusions.

Formulae (4.34)–(4.37) explicitly express the dependence of the effective conductivity tensor on the centres *a*_{k} of inclusions. If *a*_{k} are located at the sites of a regular array, from equation (4.33) becomes the standard lattice sum *S*_{2} and we arrive at formulae by McPhedran & Nicorovici (1997), Nicorovici & McPhedran (1996) and Yardley *et al.* (1999, 2001) deduced for special regular locations of ellipses. It is worth noting that the later papers contain higher-order approximations in *ν* for a regular array of ellipses. The method proposed in this paper can also be applied to higher-order approximations for an arbitrary distribution of inclusions. This can be done for numerically fixed *a*_{k}. However, direct computations with *a*_{k} in symbolic form yield enormous formulae that should be written in a simple form via the generalized lattice sums (*k*≥2). This requires advanced symbolic computations. Such computations were performed for circular inclusions (*α*=0) by Mityushev (2001), Mityushev *et al.* (2008) and Berlyand & Mityushev (2001, 2005) with a practically arbitrary precision in *ν*. The obtained analytical formulae are consistent with the results presented by Milton (2002), McPhedran (1986) and McPhedran & Milton (1987). One can believe that such computations could be developed to the general case of elliptical inclusions to derive higher-order formulae for the tensor *Λ*.

## Acknowledgements

I am grateful to referees for remarks improving the presentation of the results.

## Footnotes

- Received April 25, 2009.
- Accepted June 12, 2009.

- © 2009 The Royal Society