## Abstract

We study transport properties of equal infinite unidirectional circular cylinders that are arbitrarily distributed in a uniform host. The problem is reduced to a conjugation problem for the two-dimensional Laplace equation. The flux is written exactly in the form of a power series in Bergman's contrast parameter. Asymptotic formulae for the flux are deduced for densely packed inclusions. These formulae involve two small parameters: the ratio of the distance between the close discs to their radius and a dimensionless parameter that characterizes the degree of perfect conductor/isolator. Validity of the structural approximation is discussed.

## 1. Introduction

The prediction of the behaviour of composites by the use of transport properties of constitutes and by microstructure is an important problem of micromechanics. It can be formulated as a computational problem for the field exerted by an external source and disturbed by inclusions. The present paper is addressed to the study of the unidirectional infinite circular cylinders of conductivity *λ* embedded in a host matrix of the normalized unit conductivity. It is assumed that the section of any cylinder perpendicular to its axis is a disc of radius *r*. Then we arrive at the two-dimensional problem consisting of the determination of the potential satisfying the Laplace equation with appropriate boundary or conjugation conditions on the boundaries of the discs.

Theoretical results already obtained can be divided into pure numerical and analytical. Extensive reviews of modern numerical methods are given by Helsing (Helsing & Wadbro 2005), Kolodziej (Kolodziej & Strek 2001), Mityushev & Rogosin (1999) and Nishimura (2002; see also works cited therein).

The effective conductivity for a composite medium consisting of a dense square array of identical, perfectly conducting discs was estimated by Keller (1963)(1.1)for , where *ν* is the concentration of inclusions.

McPhedran *et al.* (McPhedran 1986; McPhedran & Milton 1987; McPhedran *et al*. 1988) calculated the effective conductivity of periodic arrays of cylinders, having used multipole expansions for the fields in and around each inclusion and generalized an identity deduced by Rayleigh (1892) to provide a set of linear equations for the multipole coefficients. A pair of nearly touching cylinders of high conductivity was considered, and truncation errors occurring in the matrix solution of the corresponding transport problem for the square array of cylinders were estimated. These investigations were developed and summarized by Movchan *et al*. (2002).

Having used a method of functional equations, Mityushev (1993) solved, in analytic form, a two-dimensional problem for finite numbers of non-overlapping discs. The locations of discs, their number and radii were arbitrary. The solution of the problem was obtained in the form of power series on a contrast parameter introduced by Bergman (1978, 1985)(1.2)One can consider the results of Mityushev (1993) as a modification of the alternating Schwarz method discussed by Mikhlin (1964), Mityushev & Rogosin (1999) and Helsing (Helsing & Wadbro 2005) in order to obtain a convergent algorithm for any geometry. This result is briefly presented in §2. Furthermore, it was extended to periodic problems by Rylko (2000), Berlyand & Mityushev (2001) and Mityushev (2001, 2005, 2007). Numerical computations were performed by Golubeva & Mityushev (1994) by the method of functional equations for |*ρ*| not closed to unity. Movchan *et al*. (2002) used the dipole matrix introduced by Pólya & Szëgo (1951) to evaluate the energy change associated with the presence of inclusions in dilute composites.

The structural approximation of composites in which a discrete network corresponds to the location of inclusions in composites was discussed by Berlyand & Kolpakov (2001, in preparation) and papers cited therein. This theory was applied to identical perfectly conducting discs randomly distributed in the representative cell. An asymptotic formula for the effective conductivity similar to Keller's formula (1.1) was deduced (see also formula (4.24) for the specific flux).

To summarize the previous results, it is worth noting that the method of functional equations yields constructive analytical formulae when the distances between the discs are not small or |*ρ*| is not closed to unity. The structural approximation yields reasonable results when the distances are small and *λ*=∞ (*ρ*=1). Therefore, it is interesting to investigate composites when some of the discs are close to others and some are far away from others. It is also interesting to consider the case *ρ*=−1 that was not discussed by Berlyand & Kolpakov (2001, in preparation).

In the present paper, we combine the method of functional equations and the structural approximation to investigate the field in two-dimensional composites with circular inclusions. In §2, the method of functional equations is described and its computational efficiency is verified. Section 3 is devoted to computations of the dipole matrix. Following Berlyand & Kolpakov (2001, in preparation) in §4, we introduce the small parameter *δ*/*r*, where *δ* is the minimum gap between discs and *r* is their radius. Another dimensionless small parameter *ϵ* is introduced to characterize the degree of perfect conductor/isolator(1.3)Different scales of *δ*/*r* and *ϵ* and their influence on the field are discussed in §§4 and 5. As a consequence, we obtain criteria of the validity of the structural approximation and new analytical formulae for the field in various scales of the distances between the inclusions. Using these formulae and the structural approximation of Berlyand & Kolpakov (2001, in preparation), we develop the method of functional equations to compute the flux between densely packed discs.

## 2. Method of functional equations

In the present section, we follow Mityushev (1993) with the following modifications: instead of the complex potential, its derivative is found and the final form solution is presented, not in the form of a series as given by Mityushev (1993), but through successive approximations that are more effective in symbolic computations.

It is convenient to use the complex variable *z*=*x*+i*y*. Consider non-overlapping discs (*k*=1, 2, …, *n*) on the complex plane . Let *D* denote the complement of the closures of *D*_{k} to the extended complex plane (see figures 1 and 4).

Let the domains *D* and *D*_{k} be occupied by materials of the normalized conductivities 1 and *λ*, respectively. The potential *u*(*x*, *y*)≡*u*(*z*) satisfies the Laplace equation in the discs *D*_{k} and in the matrix *D* with the conjugation conditions (transmission conditions)(2.1)where is the outward normal derivative to the circle Let the external field be applied in the *x*-direction and the flux at infinity be equal to (−1, 0). Then *u*(*z*)+*x* is bounded at infinity.

Introduce the complex potentials *ψ*(*z*) and *ψ*_{k}(*z*) that are analytic in *D* and *D*_{k}, respectively, and which are continuous in their closures in such a way that(2.2)and *ψ*(∞)=0. The two real conditions of (2.1) can be written as one complex one known as the -linear condition (details are explained by Mityushev (1993) and Mityushev & Rogosin (1999))(2.3)where the bar stands for the complex conjugation.

The -linear problem (2.3) can be reduced to the system of functional equations of Mityushev (1993) and Mityushev & Rogosin (1999)(2.4)wherestands for the inversion with respect to the circle .

Let *ψ*_{k}(*z*) be a solution of (2.4), then(2.5)It is convenient to write the flux as the complex value . Then in the force of (2.2), it can be written in terms of the complex potentials as follows:(2.6)

The above formulae have the following physical interpretation. Formula (2.4) means that the flux at the *k*th inclusion is cancelled by the fluxes at all other inclusions (, ) and by the external flux. Here, the external flux is defined by the vector (−1, 0). It is expressed by the complex number −1+i0. Formulae (2.5) and (2.6) show that the flux at the matrix is cancelled by the flux at all inclusions and by the external flux (compare to Rayleigh's approach of the same facts in discrete form presented by Rayleigh (1892), McPhedran *et al*. (McPhedran 1986; McPhedran & Milton 1987) and Rylko (2000)).

(*Mityushev 1993*; *Mityushev & Rogosin 1999*). *The system of functional equations* *(2.4)* *has a unique solution. This solution can be found by the method of successive approximations*(2.7)*where*(2.8)*The series* *(2.7)* *converges absolutely and uniformly in* *for any fixed* .

Formulae (2.7) and (2.8) can be used to obtain the flux (2.6) in analytic form. One can expect that the partial sums of the series (2.7) give good results for separated discs. Computation of the flux can be expensive for *ρ*=±1 between the discs when the discs are close enough to be touching. This suggestion is confirmed by the results of the example presented in figure 1. Here, the centres of three discs are located at the points *a*_{1}=−1, *a*_{2}=0 and *a*_{3}=1. The flux is calculated between the second and the third discs at the middle point *z*=0.5. One can see in figure 2 that six iterations for *ρ*=0.7 give reasonable results for *r* closed to 0.5. For *ρ*=1, the iterations are slowly convergent for *r* closed to 0.5, as demonstrated in figure 3.

## 3. Dipole matrix

Let *D* be an arbitrary multiply connected domain bounded by smooth curves ∂*D*_{k}, (*k*=1, 2, …, *n*). Let *n*(*t*) denote the unit outward normal vector to ∂*D*_{k} at the point *t* written as a complex value. Consider the -linear condition(3.1)where *ξ*=1 or −*i*. The functions and are related to the flux by (2.6), which is exerted by the dipole at infinite *ξ*. The complex number *ξ*=1 corresponds to the flux at infinity (−1, 0); *ξ*=−*i* corresponds to (0, 1). The function at infinity has the following asymptotic behaviour:(3.2)The matrix(3.3)is called the dipole matrix. This definition is in accordance with the definition of Movchan *et al*. (2002) that was given in terms of harmonic functions. The matrix is symmetric. Constructive formulae for are given by Movchan *et al*. (2002) for simply connected domains.

Making use of theorem 2.1, one can compute for circular multiply connected domains. It follows from (3.2) and (2.5) that(3.4)where satisfies the functional equations (2.4) and can be computed by (2.7) and (2.8). The value *m*^{(−i)} has the form(3.5)where satisfies the functional equations(3.6)

Consider two numerical examples with *ρ*=1.

Consider seven discs with the centres at

*a*_{1}=−1,*a*_{2}=1,*a*_{3}=−3+i,*a*_{4}=2.3−0.4i,*a*_{5}=3+2i,*a*_{6}=2+i and*a*_{7}=5+1.3i (see figure 4). The traces for the numbers of iterations beginning at 3 coincide, even for*r*closed to 0.68007 when the discs touch (see figure 5).Consider

*n*non-overlapping discs uniformly distributed in the unit square. Let*n*and*r*change but the concentration*ϕ*of the discs in the square is fixed,*ϕ*=*nπr*^{2}=0.3. Let denote the dipole matrix of*n*discs computed by*k*iterations and tr_{k}(*n*) be its trace. The number of iterations*k*≥3 gives the same numerical results. The dipole matrix and the corresponding traces tr_{k}(*n*) of randomly distributed*n*discs are presented belowDespite the fact that the inclusions are randomly generated, the theoretical isotropy is not reached in the matrices. It is related to slow convergence in*n*(not in the number of iterations). It is surprising that the isotropy is not reached, even for 5000 iterations (*n*=50).

The dipole matrix can be approximated by analytical formulae involving the centres *a*_{j} in symbolic forms. However, the final formulae are too long. They are used in the above examples even though they are hardly observed. Consider a simple example with three discs when the centres *a*_{1}, *a*_{2} and *a*_{3} are located at the real axis. Without loss of generality, it is assumed that *a*_{3}=0. For instance, the trace calculated by three iterations up to *O*(*r*^{8}) has the form(3.7)

In this section, we have justified that the flux between closed discs does not impact essentially on the dipole matrix of the group of discs. Therefore, the dipole matrix can be efficiently calculated by the method of functional equations with a small number of iterations. Movchan *et al*. (2002) described a relation between the dipole matrix and the polarizability tensor that is directly related to the effective conductivity of dilute composites by the Maxwell approximation formula (see Milton 2002). Here, by dilute composites we mean such composites in which groups of the discs are diluted in the matrix. Even though discs can be locally dense in each group, they do not generate a cluster in the infinite plane.

## 4. Flux between closed discs

Let be the distance between the discs *D*_{k} and *D*_{m}; . In this section, we investigate the problem when the dimensionless parameter *δ*/*r* is small, that is, some discs are almost touching.

Consider two discs for which the minimum *δ* is attached. For definiteness, we take *δ=δ*_{12}. Our goal is to estimate the flux (2.6) between the discs *D*_{1} and *D*_{2} at the point *z*=(*a*_{1}+*a*_{2})/2. Pair interactions between these discs produce the main part of the local field, but the influence of all *n* discs on the field between the marked discs *D*_{1} and *D*_{2} is also taken into account.

In order to investigate the functions *ψ*_{1} and *ψ*_{2} satisfying the system (2.4), we introduce the composition of two inversions with respect to the circles and (4.1)The Möbius mapping (4.1) has two fixed points that can be found from the square equation *α*(*z*)=*z*(4.2)where *a*=*a*_{2}−*a*_{1}. Assuming that *δ*/*r* is a small parameter we get(4.3)where *θ*=arg *a* is the argument of the complex number *a*.

Using (2.5) and (2.6), we estimate the flux between the discs at the point *z*=(*a*_{1}+*a*_{2})/2. At the beginning, we consider the value from (2.5)(4.4)The term with *m*=1 in (4.4) is equal to , where . The distance between the points *w* and *z*_{1} is of order , since(4.5)It follows from (4.5) that(4.6)The term with *m*=2 in (4.4) is estimated by the same method(4.7)

The terms with *m*≠1, 2 from (4.4) are of order *O*(1)(4.8)where *C* is a constant not depending on *δ*. In the force of (4.6)–(4.8), we obtain from (4.4)(4.9)

We now proceed to estimate *ψ*_{1}(*z*_{1}) and *ψ*_{2}(*z*_{2}). Substitute *z*=*z*_{1} in the first equation (*k*=1) of (2.4)(4.10)Using (4.3), one can check the relations(4.11)and(4.12)since . It follows from (4.12) that(4.13)

Substitution of (4.11) and (4.12) into (4.10) yields(4.14)where *c*_{1} is a constant not depending on *δ*. Here, the following relation is used:Application of the same arguments to the second functional equation (2.4) yields(4.15)

In order to estimate (4.9), we add the relations (4.14) and (4.15)(4.16)where . Solving equation (4.16), we have(4.17)Using the relation and the definition of *X*, we rewrite (4.9) in the form(4.18)Formulae (4.17) and (4.18) determine the behaviour of for small *δ*. Below we discuss the limit cases of (4.17) and (4.18) when *ρ*=1 (highly conducting inclusions) and *ρ*=−1 (poorly conducting inclusions). The more complicated cases when *ρ*=±(1−*ϵ*), as *ϵ* tends to zero, are also discussed.

Let *ρ*=1, then (4.17) and (4.18) and (2.6) imply that the flux at the point (*a*_{1}+*a*_{2})/2 has the form(4.19)where Re stands for the real part. Let *ρ*=−1. Then(4.20)where Im stands for the imaginary part. Formulae (4.19) and (4.20) show that the fluxes *q*_{−1} and *q*_{1} have, in general, a singularity of order as *δ* tends to zero. The coefficients and are of order *O*(1) and can vanish. The constant *c* can be determined by matching the fluxes computed by Berlyand–Kolpakov's method and by (4.19) (see §5 for details). Another consequence of formulae (4.19) and (4.20) is that the main terms of the vectors *q*_{−1} and *q*_{1} are perpendicular.

We now consider the asymptotic behaviour of the flux when not only *δ*/*r* tends to zero, but *ϵ* also tends to zero where *ρ*=1−*ϵ*. Then (2.6), (4.17) and (4.18) imply the following formula:(4.21)where , as *δ*/*r* and *ϵ* tends to zero. The flux *q*_{1} depends on two small parameters *δ*/*r* and *ϵ* and on one constant *c*, which is determined by the field induced by the discs *D*_{k}, *k*=3, 4, …, *n* (see §5).

The following formula takes place for poorly conducting inclusions:(4.22)when *ρ*=−1+*ϵ*, *ϵ* is a small positive parameter.

Formulae (4.21) and (4.22) can be simplified for various scales of the small parameters and *ϵ*. Consider an example when , then (4.21) and (4.22), respectively, become(4.23)

The asymptotic formulae (4.17)–(4.23) are in accordance with Berlyand–Kolpakov's theory (Berlyand & Kolpakov 2001, in preparation). Recall one of their results was devoted to densely packed discs of high conductivity (*ρ*=1). The continuous composite is approximated by a network connecting the centres of the discs chosen in accordance with the Voronoi tessellation. The specific flux between the discs with the numbers *i* and *j* is given by the following asymptotic formula (actually, a more precise formula was deduced):(4.24)Below, we compare this formula with the results obtained in §4.

It is worth noting that (4.17) and (4.18) imply that the flux between closed discs has not any singularity in *δ*, if |*ρ*| is separated from unity. Therefore, the structural approximation cannot be applied in this case.

The asymptotic formula (4.19) confirms Berlyand–Kolpakov's formula (4.24) obtained for *ρ*=1 (see §5 for details). It follows from (4.20) that the structural approximation can be applied for *ρ*=−1. This fact can be also deduced from the following general consideration. The flux *q*(*z*) has the form (2.6), where the complex potentials *ψ*(*z*) and *ψ*_{k}(*z*) satisfy the -linear problem (3.1) with *ξ*=1. Let us multiply (3.1) with *ξ*=1 by −i(4.25)One can see that and satisfy the -linear problem (4.25) that describes the flux in another composite with the same geometry, but the external field is rotated by 90° and the ratio of the conductivities of inclusions and matrix, *λ* is replaced by 1/*λ* (*ρ* is replaced by −*ρ*). As a result of these changes, the flux in the second composite is equal to the flux in the original one, but it is rotated by 90°. Similar arguments were used by Keller (1964) to deduce his famous duality relation for the effective conductivity. Therefore, the structural approximation can be applied to poorly conducting inclusions (*ρ*=−1) for two-dimensional composites.

Making use of formulae (4.21)–(4.23), one can investigate the behaviour of the flux between densely packed discs when |*ρ*| is closed to 1 via the small parameter *ϵ* introduced by (1.3) and its relation to . If *ϵ* is of the order of , the asymptotic behaviour of the flux is described by (4.21) and (4.22). If , the flux is estimated only via *ϵ*. As it follows from (4.23), the flux between the discs is proportional to 1/*ϵ*. Thus, the flux is sensitive to deviations of *ρ* near the limit points *ρ*=±1, when *ϵ* exceeds . Hence, the structural approximation can be also applied, but with the asymptotics (4.23) depending on *ϵ*.

## 5. Combination of functional equations and structural approximation

Theoretical results presented in the previous sections can be used not only to separate out the cases when the method of functional equations or structural approximation can be applied to dilute or densely packed composites, respectively. They are also useful to combine both approaches to composites with all scales of distances between the discs.

The flux in the matrix is given by formulae (2.5) and (2.6), where *ψ*_{k} is calculated in theorem 2.1. These calculations can be expensive for the flux in the gaps between close discs (see numerical computations at the end of §2). In order to estimate the flux between close discs when *ρ*=1, we apply the structural approximation. Formula (4.24) would not do, since it was deduced as a total specific flux between two discs (see (5.17)). More precisely, formula (4.24) provides a discrete approximation of the flux in the whole exterior of the discs. But the flux far away from the gaps is effectively computed by the method of functional equations and we need a local approximation between the discs.

We estimate the flux at the middle point between the discs. Introduce local coordinates near two discs chosen in such a way that the centres of the discs −*a*/2, *a*/2 lie on the real axis. Introduce the positive constant(5.1)The fixed point *z*_{2} defined by (4.2) takes the form(5.2)The potential(5.3)is equal to on the circles . According to (2.6), the flux can be expressed as the complex value(5.4)The flux at the point *z*=0 is equal to(5.5)Assuming that tends to zero and using (5.1) and (5.2), we obtain from (5.5)(5.6)

One can see that *q*(0) depends on the difference of the potentials on the discs Δ*u* and the specific flux *δ*^{−1} at the middle point, which differs from the total specific flux (4.24). In order to estimate Δ*u* in the external field, we consider the following problem:(5.7)where the function *u*(*z*) is harmonic in except infinity where . The constant Δ*u* in (5.7) is unknown and has to be determined by a solution to the problem.

The conformal mapping (compare to (5.3))(5.8)maps *D* onto the annulus , where *R* is given by (5.1) and *z*_{2} by (5.2). Then the problem (5.7) becomes(5.9)where is harmonic in *D*′. The boundary value problem (5.9) can be solved by the method of functional equations. Introduce the complex potentials , and that are analytic in *D*′, and , respectively, satisfying the -linear problem(5.10)Mityushev & Rogosin (1999) have established that the problems (5.9) and (5.10) are equivalent and . The functions and satisfy the functional equations(5.11)and(5.12)It is possible to solve equations (5.11) and (5.12). But for our purposes, it is sufficient to determine Δ*u*. It is easy to get(5.13)by substituting *w*=0 in (5.11) and *w*=∞ in (5.12).

The flux at the middle point *q*(0) from (5.5) becomes(5.14)Assuming that *δ*/*r* tends to zero, we obtain from (5.14) and (5.5)(5.15)This is the main formula to estimate the flux between close discs. One can see from (5.13) and (4.3) that the difference of the potentials between the discs is(5.16)

The total flux between the discs is calculated by substitution of *z*=i*y* in (5.4) and by calculation of the integral(5.17)Hence, the total specific flux is equal to . This coincides with (4.24) and confirms that asymptotic formulae deduced by functional equations coincide with formulae of the structural approximation.

The difference between the results obtained by the two methods, and a possibility to match these results are demonstrated in the following example. Consider two perfectly conducting discs with centres *a*_{1}=−1 and 1. The flux *q* is calculated at the middle point *z*=0. It follows from figure 6 that 5 iterations in the method of functional equations (see theorem 2.1) are not enough for *r*∼0.99. The results obtained by 200 iterations coincide with the results obtained by formula (5.15). However, one can see from figure 7 that 200 iterations are not enough for *r*∼0.9999. Hence, the flux can be precisely computed by functional equations exterior to a sufficiently small domain . The value of *η* is chosen in such a way that the flux is matched with the flux given by (5.15) in along the circle . It follows from figures 6 and 7 that *η* can be taken as 0.05 for 5 iterations and as 0.00005 for 200 iterations.

Based on Berlyand–Kolpakov's constructive method and on the asymptotic formulae (4.19)–(4.23), it is possible to apply the following method to determine the local flux for arbitrarily fixed locations of the discs. First, using the method of functional equations with a moderate number of iterations yields an analytic approximate formula for the flux out of the gaps between closed discs. Second, Berlyand–Kolpakov's method yields exact expressions for the difference of the potentials Δ*u*_{ij} between the closed discs with the numbers *i* and *j* for any fixed location of the discs. Then, one can determine the constants from (4.19) by the relation and obtain the flux (4.19) in the gaps between the discs. Note that (5.16) implies that Δ*u*_{ij} is of order .

The same constants *c*_{ij} are used in (4.20) for *ρ*=−1 and in (4.23) for .

## 6. Conclusion

The transport properties of equal infinite unidirectional circular cylinders arbitrarily distributed in a uniform host are systematically investigated by the method of functional equations and by a structural approximation. Asymptotic formulae for the flux are deduced for densely packed inclusions. These formulae involve two small parameters: the ratio of the distance between the discs to their radius *δ*/*r* and the degree of perfect conductor/isolator *ϵ*. Using these formulae, we amplify the method of functional equations by structural approximation to compute the flux between densely packed discs. Therefore, a combination of the methods yields analytical formulae for the flux in all scales of the distances between the discs.

We have also verified in §4 that the dipole matrix of the set of discs can be calculated by the method of functional equations with a small number of iterations.

## Footnotes

- Received July 3, 2007.
- Accepted November 2, 2007.

- © 2007 The Royal Society