Effective conductivity of a random suspension of highly conducting spherical particles

TL;DR

This paper refines the analytical understanding of the effective conductivity in random suspensions of highly conducting spheres, correcting previous formulas and developing a new computational model to account for inclusion interactions.

Contribution

It corrects Jeffrey's formula for effective conductivity by properly analyzing the conditionally convergent sum and introduces a new symbolic computation model for higher-order terms.

Findings

Corrected Jeffrey's formula up to O(f^3)

Developed an algorithm for effective conductivity tensor up to O(f^{10/3})

Showed that previous universal formulas are limited to dilute or non-interacting composites.

Abstract

Randomly distributed non-overlapping perfectly conducting spheres are embedded in a conducting matrix with the concentration of inclusions $f$. Jeffrey (1973) suggested an analytical formula valid up to $O(f^3)$ for macroscopically isotropic random composites. A conditionally convergent sum arose in the spatial averaging. In the present paper, we apply a method of functional equations to random composites and correct Jeffrey's formula. The main revision concerns the proper investigation of the conditionally convergent sum and correction the $f^2$-term. A new model of symbolic computations is developed in order to compute the effective conductivity tensor. The corresponding algorithm is realized up to $O(f^{\frac{10}3})$. The obtained formulae explicitly demonstrate the dependence of the effective conductivity tensor on the deterministic and probabilistic distributions of inclusions in the $f^2$-term, and in the $f^3$-term. This leads to the conclusion that some previous formulae presented as universal, i.e., valid for all random composites, may be actually applied only to dilute or to special composites when interaction between inclusions do not matter.

Equations

$$\frac{{\lambda }_e}{\lambda }=\frac{1+2\beta f}{1-\beta f}+O(f^2).$$

$$\frac{{\mu }_e}{\mu }=1+\frac{5}{2}f+O(f^2),$$

$$\frac{{\lambda }_e}{\lambda }=1+3\beta f+3{\beta }^2{f}^2+3f^2{\beta }^3\left(\frac{3}{4}+\frac{9}{16}\frac{{\lambda }_1+2}{2{\lambda }_1+3}+\dots \right).$$

$$\mathbf{n}(\mathbf{x})=\frac{1}{r_k}(\mathbf{x}-{\mathbf{a}}_k),\mathbf{x}\in \partial D_k(k=1,2,\dots ,n).$$

$$u=u_k,\frac{\partial u}{\partial \mathbf{n}}={\lambda }_1\frac{\partial u_k}{\partial \mathbf{n}},|\mathbf{x}-{\mathbf{a}}_k|=r_k,k=1,2,\dots ,N.$$

$$[u{]}_1=1,[u{]}_2=0,[u{]}_3=0,$$

$$\frac{\partial }{\partial \mathbf{n}}=\frac{x_1-a_{k1}}{r_k}\frac{\partial }{\partial x_1}+\frac{x_2-a_{k2}}{r_k}\frac{\partial }{\partial x_2}+\frac{x_3-a_{k3}}{r_k}\frac{\partial }{\partial x_3}=\frac{\partial }{\partial r},r=r_k,$$

$$u=u_k,\frac{\partial u}{\partial r}={\lambda }_1\frac{\partial u_k}{\partial r},r=r_k,k=1,2,\dots ,N.$$

$$-⟨q_i⟩={\int }_D\frac{\partial u}{\partial x_i}d\mathbf{x}+{\lambda }_1\sum _{k=1}^N{\int }_{D_k}\frac{\partial u_k}{\partial x_i}d\mathbf{x}=$$

$${\int }_{\partial \mathcal{O}}u{n}_{i}ds+\sum _{k=1}^N{\int }_{\partial D_k}\left({\lambda }_1{u}_k-u\right)n_{i}ds,$$

$$-⟨q_i⟩={\delta }_{i1}+({\lambda }_1-1)\sum _{k=1}^N{\int }_{D_k}\frac{\partial u_k}{\partial x_i}d\mathbf{x}.$$

$$-⟨q_i⟩={\delta }_{i1}+({\lambda }_1-1)\sum _{k=1}^N\frac{4}{3}\pi r_k^3\frac{\partial u_k}{\partial x_i}({\mathbf{a}}_k).$$

$${\lambda }_{ij}={\delta }_{i1}+({\lambda }_1-1)\sum _{k=1}^N\frac{4}{3}\pi r_k^3\frac{\partial u_k}{\partial x_i}({\mathbf{a}}_k).$$

$${\lambda }_{ij}={\delta }_{ij}+({\lambda }_1-1)\frac{f}{N}\sum _{k=1}^N\frac{\partial u_k}{\partial x_i}({\mathbf{a}}_k)(i,j=1,2,3),$$

$$u=u^{(0)}+\frac{1}{{\lambda }_1}{u}^{(1)}+\frac{1}{{\lambda }_1^2}{u}^{(2)}+\dots ,u_k=u_k^{(0)}+\frac{1}{{\lambda }_1}{u}_k^{(1)}+\frac{1}{{\lambda }_1^2}{u}_k^{(2)}+\dots .$$

$$u^{(0)}=c_k,|\mathbf{x}-{\mathbf{a}}_k|=r_k(k=1,2,\dots ,N).$$

$${\lambda }_{ij}={\delta }_{ij}+\sum _{k=1}^N\frac{4}{3}\pi r_k^3{\int }_{\partial D_k}(x_i-a_{ki})\frac{\partial u^{(0)}}{\partial r}ds+O({\lambda }_1^{-1}).$$

$$u=c_k,\mathbf{x}\in \partial D_k(k=1,2,\dots ,n).$$

$$u(\mathbf{x})-x_j\text{mbox}tendsto0,as|\mathbf{x}|\to \infty .$$

$${\int }_{\partial D_k}\frac{\partial u}{\partial \mathbf{n}}ds=0⟺{\int }_{\partial D_k}\frac{\partial u}{\partial r}ds=0(k=1,2,\dots ,n).$$

$${\mathbf{x}}_{(k)}^{\ast }=\frac{r_k^2}{r^2}(\mathbf{x}-{\mathbf{a}}_k)+{\mathbf{a}}_k,$$

$${\mathcal{K}}_{k}w(\mathbf{x})=\frac{r_k}{r}w({\mathbf{x}}_{(k)}^{\ast }).$$

$$\frac{\partial u}{\partial \mathbf{n}}=\frac{u_k}{r_0}+2\frac{\partial u_k}{\partial \mathbf{n}},|\mathbf{x}-{\mathbf{a}}_k|=r_k,(k=1,2,\dots ,n)$$

$$u_k(\mathbf{x})=-\sum _{m\ne k}\frac{r_m}{|\mathbf{x}-{\mathbf{a}}_m|}{u}_m({\mathbf{x}}_{(m)}^{\ast })+x_j-c_k,|\mathbf{x}-{\mathbf{a}}_k|\le r_k,k=1,2,\dots ,n$$

$${\int }_{\partial D_k}{u}_{k}ds=0(k=1,2,\dots ,n).$$

$$u(\mathbf{x})=-\sum _{m=1}^n\frac{r_m}{|\mathbf{x}-{\mathbf{a}}_m|}{u}_m({\mathbf{x}}_{(m)}^{\ast })+x_j,\mathbf{x}\in \dot{D}.$$

$$x_{(k_m{k}_{m-1}\dots k_1)}^{\ast }:={\left(x_{(k_{m-1}{k}_{m-2}\dots k_1)}^{\ast }\right)}_{k_m}^{\ast }.$$

$$u_k(\mathbf{x})=(Ph)(\mathbf{x}),\mathbf{x}\in D_k\cup \partial D_k,$$

$$(P_{k}h)(\mathbf{x})=h(\mathbf{x})-$$

$$\sum _{m\ne k}\frac{r_m}{|\mathbf{x}-{\mathbf{a}}_m|}h({\mathbf{x}}_{(m)}^{\ast })+\sum _{m\ne k}\sum _{k_1\ne k}\frac{r_m}{|\mathbf{x}-{\mathbf{a}}_m|}\frac{r_{k_1}}{|{\mathbf{x}}_{(m)}^{\ast }-{\mathbf{a}}_{k_1}|}h({\mathbf{x}}_{(k_{1}m)}^{\ast })-$$

$$\sum _{m\ne k}\sum _{\begin{array}{c}{k}_1\ne m\\ k_2\ne k_1\end{array}}\frac{r_m}{|\mathbf{x}-{\mathbf{a}}_m|}\frac{r_{k_1}}{|{\mathbf{x}}_{(m)}^{\ast }-{\mathbf{a}}_{k_1}|}\frac{r_{k_2}}{|{\mathbf{x}}_{(k_{1}m)}^{\ast }-{\mathbf{a}}_{k_2}|}h({\mathbf{x}}_{(k_2{k}_{1}m)}^{\ast })+\dots ,$$

$$u_k(\mathbf{x})=(P{x}_j)(\mathbf{x})-(Pc)(\mathbf{x}),\mathbf{x}\in D_k\cup \partial D_k(k=1,2,\dots ,n).$$

Full text

Keywords: random composites; spherical inclusions; effective conductivity; Eisenstein summation; triple periodic functions ; symbolic computations

1 Introduction

The effective properties of composites with non-overlapping spherical inclusions are of considerable interest in a number of applied investigations. Due to universality of mathematical modeling we refer to conductivity and transport phenomena which arise in different physical problems, e.g., electrical conductivity, dielectrics, magnetism, heat conduction, diffusion, flow in porous media, and antiplane elasticity [36, 19]. Analytical approximate formulae for the macroscopic properties of random media are important in the fundamental sciences. They qualitatively describe the difference between regular and irregular composites, demonstrate precisely the dependence of their properties on geometrical structure. Such a dependence has the crucial importance in applications to polymer based composites and to nanocomposites.

Clear comprehension of the impact of random geometry onto the macroscopic properties is established in statistical mechanics [12, 16] and in the theory of bounds [19, 35, 4]. Analytical formulae for the effective properties of composites are used in material sciences to create composites with desired properties. However, engineers frequently derive and use analytical formulae valid for a narrow class of geometric structures, for instance dilute. What is more, they declare their universality. This mythology is shortly criticized in [8]. It is difficult to divide valuable and wrong analytical formulae because wrong formulae are not wrong in peculiar cases. We do not mention doubtful results in the present paper since such a mention has to be justified. An exeption concerns only the paper [7] for which such a justification is made. It is worth noting that without any doubt many wrong results are useful in investigations and lead to the proper one.

We now proceed to review the discussed results on the mathematical level. Regular structure are sufficiently well studied by a method of multipole expansions [10, 17, 18, 30, 31, 32], by Fourier series [38, 34], and by triply periodic functions [3]. Asymptotically equivalent expressions for the effective conductivity were discussed in [1] and [8, Chapter 8].

Various methods were developed to extend the above results to random composites, see an example of such a composite with $1000$ spheres per a periodicity cell is displayed in Fig.1. We do not discuss here pure numerical methods useful for special geometries [37] and pay attention to analytical approximate formulae.

Let equal spherical particles of conductivity $\lambda_{1}$ be embedded in a matrix of conductivity $\lambda$. Let $\beta=\frac{\lambda_{1}-\lambda}{\lambda_{1}+2\lambda}$ denote the contrast parameter. The first order approximation in $f$ for the effective conductivity is known as the Clausius-Mossotti approximation or Maxwell’s formula [19, 15, 35]

[TABLE]

In fluid mechanics, the effective viscosity of hard spherical particles can be found by the Einstein formula [15, 35]

[TABLE]

where $\mu$ stands for the viscosity of fluid.

Many attempts were applied to develop these formulae to high concentrations. Theoretical difficulties were related to a conditionally convergent integral (sum) arisen in the spatial averaging. This integral for the effective conductivity of the regular cubic array was determined by Rayleigh [10]. Batchelor [2] estimated a similar integral and the second order term of $\mu_{e}$ for random suspensions by rather physical intuition than a rigorous mathematical investigation. Jeffrey [7] modified Batchelor’s approach to conductivity problem having taken into account two-sphere interactions and derived the following formula for macroscopically isotropic composites

[TABLE]

This formula was compared with others in [35, p.493].

In the present paper, we revise Jeffrey’s formula (1.3) and find that the $f^{2}$-term is equal to $3f^{2}$ in the case $\beta=1$, i.e., for highly conducting inclusions. Jeffrey’s formula (1.3) in this case gives $4.51f^{2}$. We do not analyze where is the fallacy in [7]. We can say that it is certainly methodological, not a computational error, related to intuitive physical investigation of the conditionally convergent integral discussed in [7]. The critical review of the limit transition from finite to infinite number of inclusions in applications of self-consistent and cluster methods is presented in [8, 20, 21, 22] for 2D conductivity problems. In particular, it was demonstrated that self-consistent and cluster methods with various intuitive corrections are frequently applied within the accuracy $f$, but the obtained results are applied to high concentrations without any justification.

The conductivity problem for a finite number of spheres in $\mathbb{R}^{3}$ was solved in [23] in terms of the 3D Poincaré series by a method of functional equations. This method can be treated as the alternating method of Schwrarz in the form of contrast or cluster expansions [24]. The limit transition to infinite number of inclusions was performed in [8, Chapter 8] in the special 3D case when inclusions form the regular cubic lattice.

In the present paper, we develop the method of functional equations to the problem with $N$ spheres per unit cubic cell in Sec.2. First, the local field around a finite number $n$ of highly conducting spheres arbitrarily located without mutual intersections is exactly written in Sec.3. The averaged local field is calculated in Sec.4. This first part follows the scheme [8, Chapter 8]. However, we cannot fully address to [8], since it was ultimately applied only to the regular cubic lattice. Here, we modify the method in order to apply it to an arbitrary location of inclusions. The limit $n\to\infty$ of the averaged local field is calculated in Sec.5 by means of the Eisenstein summation following Rayleigh [10]. The justification of this approach was given in [17, Sec.2.4] by study of the shape-dependent sums. It was shown that the Eisenstein summation yielded the vanishing polarization charge on the exterior surface having tended to infinity (for details see [17, Fig.2]). Though this justification concerned a regular cubic lattice it can be extended verbatim to random composites. Application of the Eisenstein summation to random composites yields an analytical formula for the effective conductivity tensor. Its components are explicitly written in the form (5.21)-(5.22) up to $O(f^{\frac{10}{3}})$.

Numerical examples with $N=1000$ spheres per cell are presented in Sec.6. The randomness is considered by the straight-forward approach used in [8]. It can be shortly outlined as follows. First, a deterministic problem with an arbitrary location of non-overlapping spheres is solved and the effective conductivity tensor is explicitly written. Since the centers of spheres $\mathbf{a}_{k}$ and their number per cell $N$ are symbolically presented in the final formulae, we may consider them as the random variable $\{\mathbf{a}_{1},\mathbf{a}_{2},\ldots,\mathbf{a}_{N}\}$ which obeys a prescribed joint probabilistic distribution. We perform $10$ numerical experiments with $N=1000$ and obtain practically the same numerical values for the coefficients in the powers of $f$ in the expansion of the effective conductivity tensor. Sec.7 is devoted to application of the obtained formulae to a stir casting process used for the fabrication of composites. Appendix A contains computational formulae for the lattice sums, Appendix B a new model of the applied symbolic computations.

2 General formula for highly conducting spheres

Let the vectors $\bm{\omega}_{1}=(1,0,0)$, $\bm{\omega}_{2}=(0,1,0)$ and $\bm{\omega}_{3}=(0,0,1)$ form a cubic lattice. The fundamental periodicity cell (the $\mathbf{0}$-cell) is the cube $\mathcal{O}=\{\mathbf{x}\in\mathbb{R}^{3}:-\frac{1}{2}<x_{j}<\frac{1}{2}\;(j=1,2,3)\}$. Let the centers $\mathbf{a}_{k}$ of mutually disjoint balls $D_{k}=\{\mathbf{x}\in\mathbb{R}^{3}:|\mathbf{x}-\mathbf{a}_{k}|<r_{k}\}$ ($k=1,2,\ldots,N$) lie in $\mathcal{O}$ and $D$ denote the complement of all the balls $|\mathbf{x}-\mathbf{a}_{k}|\leq r_{k}$ to $\mathcal{O}$. One can consider the triply periodic set of balls $\{D_{k}+\sum_{i=1,2,3}k_{i}\bm{\omega}_{i}\}$ ($k_{i}\in\mathbb{Z}$) equivalent to $N$ balls in the triply periodic topology.

Let $\mathbf{n}=(n_{1},n_{2},n_{3})$ denote the unit outward normal vector to the sphere $\partial D_{k}$ and $\frac{\partial}{\partial\mathbf{n}}$ the corresponding normal derivative. The normal vector has the form

[TABLE]

We are looking for functions $u(\mathbf{x})$ harmonic in $D$ and $u_{k}(\mathbf{x})$ harmonic in $D_{k}$ ($k=1,2,\ldots,N$), respectively, and continuously differentiable in the closures of the considered domains with the conjugation (transmission) conditions

[TABLE]

Equations (2.2) express the perfect contact between materials of conductivity $\lambda_{1}$ occupying the balls and the host of the normalized unit conductivity. It is assumed that the functions $u(\mathbf{x})$ and $u_{k}(\mathbf{x})$ are quasi-periodic, namely,

[TABLE]

where $[u]_{j}:=u(\mathbf{x}+\omega_{j})-u(\mathbf{x})$ stands for the jump of $u(\mathbf{x})$ per cell along the axis $x_{j}$.

Instead of the normal derivative in (2.2) we will also consider the derivative $\frac{\partial}{\partial r}$ where $r=|\mathbf{x}-\mathbf{a}_{k}|$ is the radial local coordinate near $\partial D_{k}$. For a fixed $k$, we have

[TABLE]

where $\mathbf{a}_{k}=(a_{k1},a_{k2},a_{k3})$. Then, (2.2) becomes

[TABLE]

The averaged flux $\langle\mathbf{q}\rangle$ can be calculated by application of the Ostrogradsky - Gauss formula in $\mathcal{O}$. The flux component $\langle q_{i}\rangle$ along the axis $x_{i}$ becomes

[TABLE]

where $\mathbf{n}=(n_{1},n_{2},n_{3})$. The integral over $\partial\mathcal{O}=\partial D-\sum_{k=1}^{n}\partial D_{k}$ in (2.6) is equal to the Kronecker delta $\delta_{i1}$ because of (2.3). Using the first relation (2.2) and again the Ostrogradsky-Gauss formula, we obtain

[TABLE]

The mean value theorem for harmonic functions [25, p.294] yields

[TABLE]

Let $\Lambda=\{\lambda_{ij}\}$ denote the effective conductivity tensor defined by the relation $\mathbf{q}=-\Lambda\langle\overline{\nabla u}\rangle$. Here, $\overline{\nabla u}=([u]_{1},[u]_{2},[u]_{3})$ stands for the jump vector of the potential $u$ per cell. The normalized jump conditions (2.3) yield $\lambda_{ij}=-\langle q_{i}\rangle$, hence, (2.8) becomes

[TABLE]

The formula (2.9) is derived for $j=1$, i.e., for the external flux applied in the $x_{1}$-direction, see (2.3), and can be established analogously for $j=2,3$. Therefore, the functions $u_{k}(\mathbf{x})$ from (2.10) depend also on $j$.

In the case of equal radii we arrive at the formula

[TABLE]

where $f$ denote the concentration of inclusions.

Below, we study the case of highly conducting inclusions when $\lambda_{1}\gg 1$ and

[TABLE]

Substitute (2.11) into (2.2) and take the terms up to $O(\lambda_{1}^{-1})$. The zero-th coefficient yields the boundary value problem

[TABLE]

Equation (2.12) can be considered as the modified Dirichlet problem for the exterior of spheres with the undetermined constants $c_{k}$ on the boundary. It follows from [8, Chapter 8] that the leading part of the effective conductivity tensor is expressed by means of the function $u^{(0)}(\mathbf{x})$ by formula

[TABLE]

3 Modified Dirichlet problem

In the present section, we solve the modified Dirichlet problem (2.12) by the method of functional equations. This method was proposed in [23] for a finite number of inclusions $n$ in $\mathbb{R}^{3}$. The method was extended to triple periodic problems by the limit $n\to\infty$ in [8, Chapter 8] for $N=1$.

First, we shortly present the functional equations and their solution for a finite number of inclusions. The rigorous justification of the method of functional equations can be found in [8, 23]. It is worth noting that a system of linear algebraic equations on the undetermined constants $c_{k}$ was constructed but not solved in [8]. Next, we apply an asymptotic method to solve this system for equal radii and find the undetermined constants in the form of series. The obtained explicit formulae are used in the next sections to determine the effective conductivity tensor.

In the present section, we consider a problem for a finite domain $D(n)$ with $n$ balls in $\mathbb{R}^{3}$ and construct its solution $u(\mathbf{x},n)$. The required triple periodic function can be constructed by the limit transition $u^{(0)}(\mathbf{x})=\lim_{n\to\infty}u(\mathbf{x},n)$. It will be convenient to introduce auxiliary functions $u_{k}(\mathbf{x},n)$ harmonic in $D_{k}$ and express the effective conductivity in terms of $\frac{\partial u_{k}}{\partial x_{i}}(\mathbf{a}_{k},\infty)$. Formally, $u_{k}(\mathbf{x},n)$ differs from the potential $u_{k}(\mathbf{x})$ in $D_{k}$ for the general problem (2.2) with finite $\lambda_{1}$. The new function $u_{k}(\mathbf{x},n)$ can be considered as the potential in $D_{k}$ satisfying the special interface condition (3.6) discussed below. We write $D$ and $u(\mathbf{x})$ instead of $D(n)$ and $u(\mathbf{x},n)$ for shortness.

Consider mutually disjoint balls $D_{k}=\{\mathbf{x}\in\mathbb{R}^{3}:|\mathbf{x}-\mathbf{a}_{k}|<r_{k}\}$ ($k=1,2,\ldots,n$) and the domain $\dot{D}=\mathbb{R}^{3}\backslash\cup_{k=1}^{n}(D_{k}\cup\partial D_{k})$. It is convenient to add the infinite point and introduce the domain $D=\dot{D}\cup\{\infty\}$ lying in the one-point compactification of $\mathbb{R}^{3}$. We find a function $u(\mathbf{x})$ harmonic in $\dot{D}$ and continuously differentiable in $\dot{D}\cup_{k=1}^{n}\partial D_{k}$ with the boundary conditions

[TABLE]

Here, $c_{k}$ are undetermined constants which should be found during solution to the boundary value problem. Let the external flux at infinity be parallel to the $x_{j}$-axis, hence,

[TABLE]

We have

[TABLE]

The inversion with respect to a sphere $\partial D_{k}$ is introduced by formula

[TABLE]

where $r=|\mathbf{x}-\mathbf{a}_{k}|$ denotes the local spherical coordinate near the point $\mathbf{x}=\mathbf{a}_{k}$ . The Kelvin transform with respect to $\partial D_{k}$ has the form [36, 25]

[TABLE]

If a function $w(\mathbf{x})$ is harmonic in $|\mathbf{x}-\mathbf{a}_{k}|<r_{k}$, the function $\mathcal{K}_{k}w(\mathbf{x})$ is harmonic in $|\mathbf{x}-\mathbf{a}_{k}|>r_{k}$ and vanishes at infinity.

The method of functional equations was developed in [23, 8] to solve the modified Dirichlet problem (3.1). First, auxiliary functions $u_{k}(\mathbf{x})$ harmonic in the balls $|\mathbf{x}-\mathbf{a}_{k}|<r_{k}$, respectively, and continuously differentiable in their closures, satisfying the boundary conditions

[TABLE]

were introduced. Second, it was proved that these functions satisfy the system of functional equations

[TABLE]

and the conditions

[TABLE]

We now fix the constants $c_{k}$ in (3.7) and return to their determination at the end of this section. The functional equations (3.7) are solved by the method of successive approximations uniformly convergent in the union of the balls $|\mathbf{x}-\mathbf{a}_{k}|\leq r_{k}$, see [8] for details.

After solution to the system (3.7) the solution of the problem (3.2) is given up to undetermined constants by formula

[TABLE]

Introduce the space $\mathcal{C}(D_{k})$ of functions harmonic in all $D_{k}$ and continuous in their closures endowed with the norm $\|h\|_{k}=\max_{\mathbf{x}\in\partial D_{k}}|h(\mathbf{x})|$. Let $\mathcal{C}=\mathcal{C}\left(\cup_{k=1}^{n}D_{k}\right)$ denote the space of functions harmonic in all $D_{k}$ and continuous in their closures. The norm in $\mathcal{C}$ has the form $\|h\|=\max_{k}\|h\|_{k}=\max_{k}\max_{\mathbf{x}\in\partial D_{k}}|h(\mathbf{x})|$. The system of functional equations (3.7) can be considered as an equation in the space $\mathcal{C}$ on a function equal to $u_{k}(\mathbf{x})$ in each closed ball $D_{k}\cup\partial D_{k}$.

Let $k_{s}$ run over $1,2,\ldots,n$. Consider the sequence of inversions with respect to the spheres $\partial D_{k_{1}},\partial D_{k_{2}},\ldots,\partial D_{k_{m}}$ determined by the recurrence formula

[TABLE]

It is supposed that no equal neighbor numbers in the sequence $k_{1},k_{2},\ldots,k_{m}$. The transformations (3.10) for $m=1,2,\ldots$ with the identity map form the Schottky group $\mathcal{S}$ of maps acting in $\mathbb{R}^{3}$, see the 2D theory [26].

Introduce the function $h(\mathbf{x})$ in the space $\mathcal{C}$ defined by equations $h(\mathbf{x})=x_{j}-c_{k}$ in $|\mathbf{x}-\mathbf{a}_{k}|\leq r_{k}$ ($k=1,2,\ldots,n$). Straight-forward application of the successive approximations to (3.7) gives the exact formula

[TABLE]

where the operator $P$ acts in the space $\mathcal{C}$ and has the form $P=P_{k}$ in $D_{k}\cup\partial D_{k}$; the operator $P_{k}$ is defined in terms of the series

[TABLE]

where for instance $\sum_{\begin{subarray}{c}k_{1}\neq m\\ k_{2}\neq k_{1}\end{subarray}}:=\sum_{k_{1}\neq m}\sum_{k_{2}\neq k_{1}}$. Every sum $\sum_{k_{s}\neq k_{s-1}}$ contains terms with $k_{s}=1,2,\ldots,n$ except $k_{s}=k_{s-1}$. Following [27, Chapter 4] and [23] one can prove compactness of $P$ in $\mathcal{C}$.

The uniqueness of solution established in [23, 8] implies that the successive approximations can be applied separately to $h_{1}(\mathbf{x})$ and to $h_{2}(\mathbf{x})$ when $h(\mathbf{x})=h_{1}(\mathbf{x})+h_{2}(\mathbf{x})$. The unique result will be the same. Thus, $P$ can be applied separately to $x_{j}$ and to the piece-wise constant function $c(\mathbf{x})=c_{k}$ where $x\in D_{k}$ ($k=1,2,\ldots,n$). Then, the solution of the system (3.7) can be written in the form

[TABLE]

Substitution of (3.13) into (3.9) yields

[TABLE]

where the operator $P_{0}$ is introduced as follows

[TABLE]

[TABLE]

Consider a class of functions $\mathcal{R}$ harmonic in $\mathbb{R}^{3}$ except a finite set of isolated points located in $D$ where at most polynomial growth can take place. It follows from (3.12) that the operator $P:\mathcal{R}\to\mathcal{R}$ is properly defined. The series $(P_{0}h)(\mathbf{x})$ for $h\in\mathcal{R}$ converges uniformly in every compact subset of $\mathbf{x}\in\dot{D}\cup\partial D$. We call it by the 3D Poincaré $\theta$-series associated to the Schottky group $\mathcal{S}$ for the function $h(\mathbf{x})$.

One can see from (3.14) that $u(\mathbf{x})$ contains the undetermined constants $c_{k}$ only in the part

[TABLE]

[TABLE]

The function $(P_{0}c)(\mathbf{x})$ is represented in the form $(P_{0}c)(\mathbf{x})=\sum_{m=1}^{n}c_{m}p_{m}(\mathbf{x})$. Let $\delta_{km}$ denote the Kronecker symbol. The functions $p_{m}(\mathbf{x})$ can be obtained by Lemma 4.3 from [27, p.152] or by substitution $c(\mathbf{x})=\delta_{km}$, $\mathbf{x}\in D_{k}$ ($k=1,2,\ldots,n$) into (3.16)

[TABLE]

[TABLE]

In order to find the undetermined constants $c_{k}$ $(k=1,2,\ldots,n)$ we use equations (3.8). The integral from (3.8) is calculated by the mean value theorem for harmonic functions [25, p.294]. Then, (3.8) becomes

[TABLE]

The function $u_{k}$ is written exactly by (3.13)-(3.15). Consider the case of equal radii $r_{k}=r_{0}$ and write $u_{k}(\mathbf{x})$ explicitly up to $O(r_{0}^{7})$

[TABLE]

Introduce the designation $z_{k}:=a_{kj}-c_{k}$ and substitute $\mathbf{x}=\mathbf{a}_{k}$ into (3). Then, (3.18) becomes

[TABLE]

It can be easily obtained from (3.20) that

[TABLE]

Substitute (3.21), computed up to $O(r_{0}^{7})$, into (3)

[TABLE]

This is the main analytical approximate formula which will be used below to compute the effective conductivity tensor.

4 Averaged conductivity

The countable set of centers forms the lattice described in Sec.1 and can be ordered in the following way $\mathcal{A}=\{\mathbf{a}_{k}+\sum_{s=1,2,3}m_{s}\bm{\omega_{s}}\}$ . Here, $(m_{1},m_{2},m_{3})\in\mathbb{Z}^{3}$ is the number of the cell; the points $\mathbf{a}_{k}$ ($k=1,2,\ldots,N$) belong to the $\mathbf{0}$-cell.

First, we consider the finite number of inclusions $n=NM$ in the space $\mathbb{R}^{3}$ using the linear order numeration $\mathbf{a}_{k}$ ($k=1,2,\ldots,n$), i.e., $M$ is the number of cells and the fixed number $N$ is the number of inclusions per cell. We consider equal spheres, i.e., $r_{k}=r_{0}$ for simplicity. Introduce the value

[TABLE]

Then, the component of the effective conductivity tensor $\lambda_{ij}$ can be estimated by formula

[TABLE]

Following the scheme [6] we use the Eisenstein summation in the limit $n\to\infty$ $\Longleftrightarrow\{M\to\infty,N$ is fixed}. The subscript $j=1,2,3$ denotes the potential jump per cell along the axis $x_{j}$ in the triple periodic problem.

It follows from (3.6) that for every fixed $k$

[TABLE]

The first integral can be calculated by the Ostrogradsky-Gauss formula

[TABLE]

since the $i$th component of the unit outward normal vector to the sphere $\partial D_{k}$ is equal to $\frac{x_{i}-a_{ki}}{r_{k}}$. Let $|D_{k}|$ denote the volume of the ball $D_{k}$. Application of the mean value theorem for harmonic functions to (4.4) yields

[TABLE]

Green’s identity and the mean value theorem imply that

[TABLE]

Substitution of (4.4)-(4.6) into (4.1), (4.3) yields

[TABLE]

where the value $\frac{\partial u_{k}}{\partial x_{1}}(\mathbf{a}_{k})$ can be calculated by the approximation (3). When the external flux is applied in the $x_{1}$-direction, i.e., $j=1$, we have

[TABLE]

where the double sum $\sum_{m,s}:=\sum_{m\neq k}\sum_{s\neq m}\;.$

[TABLE]

[TABLE]

5 Effective conductivity

Consider a regular infinite array of spheres with $N=1$ and $M=n$, i.e., one sphere per periodicity cell. Rayleigh [10] calculated the limit (see formula after (64) in his paper [10]) having used the Eisenstein summation

[TABLE]

Substitution of (5.1) into (4.8) and further into the limit expression (4.7) yields the effective conductivity of the simple cubic array of sphere

[TABLE]

where $f=\frac{4}{3}\pi r_{0}^{3}$ denotes the concentration of spheres (one sphere per unit cubic cell).

For general $N$ the formula (5.2) is replaced by the following relation

[TABLE]

where

[TABLE]

[TABLE]

The vectors $\mathbf{R}=(R_{1},R_{2},R_{3})=\sum_{i=1,2,3}m_{i}\bm{\omega}_{i}$ $(m_{i}\in\mathbb{Z})$ generate the simple cubic lattice. The Eisenstein summation is used in (5.5) and it is assumed for shortness that $E_{11}(\mathbf{0}):=\frac{4}{3}\pi$ when $\mathbf{a}_{k}=\mathbf{a}_{m}$ is taken in (5.4).

The following relation holds $E_{11}=-\frac{\partial E_{1}}{\partial x_{1}}$, where $E_{1}(\mathbf{x})$ is the first coordinate of the function (8.4.77) introduced in [8]

[TABLE]

This function is related to Berdichevsky’s function $\wp_{11}(\mathbf{x})$ considered in [3]222The 3D $\wp$-functions are discussed only in the Russian edition [3]. (see also [8, Sec.3, Chapter 8])

[TABLE]

The function $\wp_{11}(\mathbf{x})$ was introduced in [3] through the absolutely convergent series

[TABLE]

where $\mathbf{R}$ runs over $\mathbb{Z}^{3}$ in the infinite sum $\sum\nolimits^{\prime}_{\mathbf{R}}$ except $\mathbf{R}=\mathbf{0}$. A power series of $E_{11}(\mathbf{x})$ on the coordinates of $\mathbf{x}$ can be applied to effectively compute the values $E_{11}(\mathbf{a}_{k}-\mathbf{a}_{m})$. Using the standard multidimensional Taylor expansion in the coordinates of $\mathbf{x}$ we obtain

[TABLE]

where, for instance,

[TABLE]

All the sums $\sum\nolimits^{\prime}_{\mathbf{R}}$ in (5.9) are at least of order $|\mathbf{R}|^{-4}$, hence, they converge absolutely [3]. In this case we can change the order of summation and use their symmetry in $R_{i}$. In particular, this implies that the coefficients in $x_{j}$ vanish.

The function $E_{12}=-\frac{\partial E_{1}}{\partial x_{2}}$ coincides with Berdichevsky’s function $\wp_{12}(\mathbf{x})$ [3] and can be presented as the Eisenstein series

[TABLE]

The function $\wp_{12}(\mathbf{x})$ was expanded in [3] into the absolutely convergent series

[TABLE]

The functions $E_{21}(\mathbf{x})$ and $\wp_{21}(\mathbf{x})$ are introduced analogously to (5.11) and (5.12) by the interchange of subscripts $1$ and $2$. It follows from the absolute convergence of (5.12) that

[TABLE]

The function $E_{13}=-\frac{\partial E_{1}}{\partial x_{3}}$ is similar to $E_{12}(\mathbf{x})$ and can be calculated analogously by replacement the subscript $2$ by $3$. The functions $E_{1j}$ determine the structural sums similar to (5.4)

[TABLE]

Introduce the absolutely convergent sums ($i_{1}+i_{2}+\ldots+i_{m}\leq\ell-4$) symmetric with respect to $R_{i}$

[TABLE]

Note that $e_{\ell}^{(i_{1},i_{2},i_{3})}=0$ when at least one superscript $i_{k}$ is odd. Then, (5.9) is reduced to

[TABLE]

The function $E_{12}(\mathbf{x})$ is calculated by the similar formula

[TABLE]

It is convenient for computations to replace combinations of $e_{\ell}^{(i_{1},i_{2},i_{3})}$ with Coulombic lattice sums [9] (see Appendix A)

[TABLE]

[TABLE]

The function $E_{13}(\mathbf{x})$ can be obtained from (5.19) by the interchange $x_{2}$ and $x_{3}$.

Introduce the discrete spatial convolution sums constructed by the sums (4.8)-(4.10)

[TABLE]

The values $e_{ij}$ and $e_{ij*pl}$ will called the structural sums. They generalize the classic lattice sums constructed for regular arrays. The approximation (5.3) is extended by application of (4.8)

[TABLE]

It follows from (4.9) that

[TABLE]

The component $\lambda_{13}$ can be derived from (4.10). Ultimately, it is written by the interchange of subscripts $2$ and $3$ in (5.22).

Remark 5.1*.*

In the case of the simple cubic array when $N=1$, the convolution (5.20) is simplified to $e_{11*11}=\left(\frac{4\pi}{3}\right)^{2}$ and $e_{12*12}=e_{13*13}=0$. The effective conductivity of the simple cubic array can be calculated by the Clausius-Mossotti approximation (1.1) for $\beta=1$

[TABLE]

what corresponds to the general formula (5.21). The most advanced formula for the simple cubic array was established by Berdichevsky [3, formula (11.80)]

[TABLE]

where $\lambda_{e}=\lambda_{11}$ is the scalar effective conductivity of the simple cubic array.

6 Numerical results for random composites

We now proceed to apply (5.20) to numerical estimation of the effective conductivity of random macroscopically isotropic composites. In samples generation we applied Random Sequential Adsorption (RSA) protocol, where consecutive objects are placed randomly in the cell, rejecting those that overlap previously absorbed ones. We generated 10 samples of the RSA distribution with $N=1000$ balls each. The concentration $f=0.3$ is fixed during the generation. Other concentrations about $f=0.3$ give practically the same result. Such a sample is displayed in Fig.1.

In calculations of (5.21) we used the following approximated values of convergent lattice sums

[TABLE]

computed as partial sums of (5.15) for $-250\leq R_{k}\leq 250$ ($k=1,2,3$) (see Appendix A). Obtained approximations agree to at least five significant digits with known numerical values [9]. Numerical results are presented in Table 1.

The structural sums (5.4) and (5.20) are special cases of discrete multidimensional convolution of functions defined in [29]. In our computations we applied efficient algorithms developed therein. The total time of calculations of results from Table 1, ran on a standard notebook equipped with Intel Core i7 Processor of 4th generation, was about 30 minutes.

One can see that $\overline{e_{11}}$ is close to $\frac{4\pi}{3}\approx 4.18879$. We suggest that $e_{11}=\frac{4\pi}{3}$ for ideally macroscopically isotropic composites. This hypothesis can be confirmed by the following observations. The expansion (5.18) of $E_{11}(\mathbf{x})$ is a linear combination of the even homogeneous polynomials $P_{l}(x_{1},x_{2},x_{3})$. One can directly check for few $l$ that

[TABLE]

We suggest that (6.2) holds for all $P_{l}(x_{1},x_{2},x_{3})$ in the expansion (5.18). Consider 24 points

[TABLE]

in the $\mathbf{0}$-cell generating a macroscopically isotropic structure. If (6.2) is true, $e_{11}=\frac{4\pi}{3}$ for these 24 points. Therefore, (5.21) for the considered macroscopically isotropic composites yields the scalar effective conductivity

[TABLE]

Remark 5.1 concerning the simple cubic array can be used to improve the polynomial formulae (5.21)-(5.22). These formulae are derived by the method of successive approximations applied to the system of functional equations (3.7) for $r_{m}=r_{0}$. We write exactly all the terms up to $O(r_{0}^{7})$ taking into account interactions among all the spheres. We shall increase this precision in a future work. But now we can exactly obtain some high order terms in the considered case of macroscopically isotropic composites. If we stay in the $p$th iteration applied to (3.7) only the term with $r_{0}^{p}$ we obtain a sequence of terms which leads to the Clausius-Mossotti approximation $1+3f+3f^{2}+3f^{3}+\ldots=\frac{1+2f}{1-f}$. Berdichevsky’s formula (5.24) takes into account interactions between spheres located in different periodicity cells. The same terms have to be in the formula for random case as in (5.24). The above argumentation makes possible to rewrite (6.3) in the asymptotically equivalent form

[TABLE]

7 Application to stir casting process

Properties of dispersed composites are strongly influenced by distributions of inclusions and interactions between them. It is observed in the experimental works [13, 14, 11] that composites with the same constitutes, shapes of inclusions and concentrations may have different macroscopic properties. The key of this diversity lies in different distributions of inclusions. Formulae (5.21)-(5.22) give us such a possibility to distinguish the geometric structures. This implies wide applications of formulae (5.21)-(5.22) to optimal design problems in fabrication of composites, to macroscopic description of polymers and to other problems where the geometry plays the crucial role.

In the present section, we develop a general effective scheme to distinguish macroscopically isotropic and anisotropic composites obtained in stir casting process. Stir casting is a technological process for the fabrication of Al-SiC and other composites [13, 14]. Below, we develop a method to optimize technological steps involved in 3D stir casting process. It is an extension of the 2D method proposed in [28].

It is established in Sec.5 that the normalized effective conductivity tensor for a general distribution of inclusions can be presented in the form

[TABLE]

where $I$ denotes the unit tensor. The symmetric tensor $\Lambda^{(2)}$ is explicitly written in terms of the structural sums

[TABLE]

where $e_{ij}$ are given by (5.4) and (5.14). Construction of the components $e_{11}^{*}$ and $e_{11}^{**}$ requires an explanation. The structural sum $e_{11}$ is calculated by equations (5.4)-(5.5) where the triple periodic function $E_{11}(\mathbf{x})$ is defined by the Eisenstein summation. The considered Eisenstein series is defined as a triple symmetric iterated series first summed along the $x_{1}$–axis and after along the $x_{2}$– and $x_{3}$– axes. According to the Rayleigh-McPhedran formalism the prime summation along the $x_{1}$–axis is justified by the applied external flux $(1,0,0)$. This method yields the component $e_{11}$.

The component $e_{11}^{*}$ in the tensor (7.2) corresponds to the conductivity in the $x_{2}$–direction when the external flux is applied along the $x_{2}$–axis. Hence, $e_{11}^{*}$ can be calculated analogously to $e_{11}$ by interchanging the $x_{1}$– and $x_{2}$– axes. In order to avoid ambiguity in the definition of Eisenstein summation we may introduce $e_{11}^{*}$ using the standard function $E_{11}(\mathbf{x})$ defined by (5.5)

[TABLE]

where $\mathbf{a}_{k}^{*}$ is obtained from $\mathbf{a}_{k}$ by interchanging the first and second coordinates, i.e. $\mathbf{a}_{k}^{*}=(a_{k2},a_{k1},a_{k3})$. The value $e_{11}^{**}$ is calculated analogously through the same Eisenstein series $E_{11}(\mathbf{x})$ by interchanging the first and third coordinates.

We now introduce a 3D anisotropy coefficient following the 2D method [28]. The deviatoric component of tensor $\Lambda^{(2)}$ has the form

[TABLE]

where

[TABLE]

We introduce the anisotropy scalar coefficient as the absolute value of the determinant of $\Lambda^{(2)}$,

[TABLE]

One can see that $\kappa\geq 0$ and $\kappa$ vanishes for macroscopically isotropic media.

The anisotropy coefficient (7.6) yields a powerful method to optimize parameters of stir casting processes for moderate concentrations. It is based on the following steps described in [13, 14] for 2D composites. First, a 3D sample has to be processed to determine the coordinates $\{\mathbf{a}_{1},\mathbf{a}_{2},\ldots\mathbf{a}_{N}\}$. Next, the coefficient $\kappa$ is calculated by (7.4)-(7.6) where $e_{ij}$ are exactly given by (5.4) and (5.14). If $\kappa$ is closed to zero, the considered composite is macroscopically isotropic. In the opposite case, the stir casting process should be continued to get a more uniform distribution of inclusions.

8 Discussion and conclusion

Derivation of analytical (approximate) formulae for the effective properties of composites requires a subtle mathematical study of the conditionally convergent sums arisen in the course of spatial averaging. The first order approximations (1.1) and (1.2) in $f$ were obtained as solutions to the single-inclusion problem.

For many years it was thought that Maxwell’s and Clausius-Mossotti approximations can be systematically and rigorously extended to higher orders in $f$ by taking into account interactions between pairs of spheres, triplets of spheres, and so on. However, it was recently demonstrated [21, 8] that the field around a finite cluster of inclusions can yield a correct formula for the effective conductivity only for non-interacting clusters. The higher order term can be properly found only after a subtle study of the conditionally convergent series. In the present paper, the Eisenstein summation is used following [10] for a regular cubic lattice. Justification of the Eisenstein summation was given in [17, Sec.2.4] by studying the shape-dependent sums. The present paper contains an algorithm to compute the effective conductivity tensor with an arbitrary precision in $f$. The symbolic computations are performed up to $O(f^{\frac{10}{3}})$. As a result, the analytical approximate formulae (5.21)-(5.22) are derived for an arbitrary locations of non-overlapping spheres, and (6.4) for a class of macroscopically isotropic composites. These formulae explicitly demonstrate the dependence of the effective conductivity tensor on the deterministic and probabilistic distributions of inclusions.

It is interesting to compare our analytical formulae (5.21)-(5.22) with others. We leave aside some numerical results for fixed locations of inclusions [37], since they do not contain the symbolic dependencies on the locations. Allthough such a dependence can be investigated by Monte Carlo methods, it requires a tedious, time-consuming computer simulations. It is rather impossible at the present time to get an accurate results. It is worth noting that our numerical results presented in Sec.6 are obtained for $10$ generated locations with $N=1000$ balls each, by using a standard laptop computer.

Some authors equate an approximate analytical formula with a model. Such an approach is misleading, since a mathematical modeling involves the fixed governed equations, interface and boundary conditions. Different approxmate formulae/solutions for the mathematical model hold under restrictions usually not discussed by authors. A serious methodological mistake may follow when intermediate manipulations are valid only within the precision $O(f)$, while the final formula is claimed to work with a higher precision, see an explicit example in [22]. In particular, it follows from our formulae (5.21)-(5.22) that it is impossible to write a universal higher order formula independent on locations of inclusions. Such a universal formula holds only for a limited class of composites with non-interacting inclusions, e.g., for a dilute composites and the Hashin-Shtrikman coated sphere assemblage [4].

It is exactly the case of misundersood precision which is confronted here. All the formulae (1.1) and (1.3) for $\beta=1$, and (6.4) are the same up to $O(f^{2})$. But not always formulae (1.1) and (1.3) can be extended to high concentrations. Moreover, there is a plenty of random distributions. Formula (6.4) holds up to $O(f^{\frac{10}{3}})$ for randomly distributed non-overlapping spheres realized with the RSA protocol. Another random distribution yields another values of the structural sums selected in Table 1 for the RSA protocol.

Let us address the Hashin-Shtrikman bounds [19]. The dependence of the effective conductivity of a random composite in $f$ corresponds to some monotonous curve drawn between the bounds. Such a curve can be sketched arbitrarily, and it will correspond to some unspecified distribution of inclusions. In the present paper, we deal with a uniform distribution corresponding to a stir-casting process described by the RSA model [13, 14]. The main theoretical requirement to the geometric model consists not only in writing a formula, but in a precise description of the geometrical conditions imposed on deterministic or random locations of inclusions. The rigorous statement and study of this theoretical problem is necessary for the proper approach to various applied problems, e.g., stir-casting process [13] discussed in Sec.7.

The comparison of (5.21)-(5.22) and (6.4) with numerical simulations and experimental results [11, 37] demonstrates their agreement for concentrations not exceeding $0.25$. For higher concentrations the results can be different. It is not surprising, since the formulae (5.21)-(5.22) explicitly demonstrate the dependence of the effective conductivity on the location of inclusions. The available experimental studies do not include the locations of inclusions, making comparison difficult.

As it is noted above, the formula (6.4) is derived for the mathematical expectation of the effective conductivity over the independent and identically distributed (i.i.d.) non-overlapping inclusions. Even the formula (6.4) for uniformly distributed non-overlapping balls is not universal, because the effective conductivity depends on the protocol of computer simulations or experimentral stirring, meaning that the very notion of randomness is non-universal [35, 8, 13, 33].

To summarize, the correct form of the $f^{2}$ term in (1.1) is determined. It is rigorously justified for a macroscopically isotropic media with highly conducting spherical inclusions. We explicitly demonstrate that the next $f^{3}$ term, depends on the deterministic and random locations of inclusions.

Appendix A Lattice sums

In three dimensions, Coulombic lattice sums are defined through the following formula [9]

[TABLE]

where

[TABLE]

As an example let us expand the sum ${\cal L}_{4}^{0}$

[TABLE]

Hence, by (5.15) we have

[TABLE]

Since the sums (5.15) are absolutely convergent, the permutation of indexes $i_{k}$ do not change their values. Therefore, the considered lattice sum simplifies to

[TABLE]

Other lattice sums can be simplified by the same method

[TABLE]

Appendix B Model of symbolic computations

Let us present a general model of symbolic computations of the solution of the system of functional equations (3.7) in form of the analytic approximation (3). The main difficulty in programming such a procedure is tackling nested summations and the fact that they are indefinite, i.e. we do not know in advance the number of functions $u_{k}$. Hence, we need a representation of the sum that maintains the summation index in a symbolic form.

B.1 Indefinite symbolic sums

Let $A_{1},A_{2},A_{3},\ldots,A_{n}$ be arbitrary sets and let $f(\mathbf{a})$ be arbitrary function of $\mathbf{a}=[a_{1},a_{2},a_{3},\ldots,a_{n}]\in A_{1}\times A_{2}\times A_{3}\times\ldots\times A_{n}$. Let us define the following sum:

[TABLE]

We do not discuss the existence of (B.1) and we assume that the information about the form of $A_{k}$ does not matter. Hereafter, all formulae are considered as symbolic expressions. Let us introduce the following indefinite symbolic sum as a representation of (B.1):

[TABLE]

where the operands are, correspondingly, a summand and a sequence of symbols over which the sum is taken. For example, using this representation, the sum $\sum_{{m,n=1}}^{N}|a_{m}-a_{m}|$ is expressed as follows:

[TABLE]

The key point of presented approach to computations using $sum$ is to prevent the expression, undergoing symbolic manipulations, from generating nested $sum$ subexpressions. This leads to the application of the following rules in the automatic simplification algorithms [5]:

Rule 1.

$sum(f(\mathbf{a})+g(\mathbf{a}),\mathbf{a})\rightarrow sum(f(\mathbf{a}),\mathbf{a})+sum(g(\mathbf{a}),\mathbf{a})$

Rule 2.

$f(\mathbf{a})sum(g(\mathbf{b}),\mathbf{b})\rightarrow sum(f(\mathbf{a})g(\mathbf{b}),\mathbf{b})$

Rule 3.

If $\mathbf{a}\cap\mathbf{b}=\emptyset$, then

         $sum(sum(f(\mathbf{a},\mathbf{b}),\mathbf{b}),\mathbf{a})\rightarrow sum(f(\mathbf{a},\mathbf{b}),[\mathbf{a},\mathbf{b}]),$

where $[\mathbf{a},\mathbf{b}]=[a_{1},\ldots,a_{n},b_{1},\ldots,b_{n},]$.

Rule 4.

$sum(0,\mathbf{a})\rightarrow 0$

B.2 Symbolic representation of successive approximations

The following system forms the generalized form of (3.7):

[TABLE]

where $u_{k}(\mathbf{x})$ ($k=1,2,\ldots,n$), $c_{k}$ are unknown and $\circ$ denotes function composition. The generalized form of the system can be applied in similar problems, for instance, in the effective conductivity of 2D composite materials [8, eq. (2.3.85)]. The system (B.3) can be rewritten in a symbolic form:

[TABLE]

assuming that summation excludes $a_{k}$. The application of the successive approximations yields following relations:

[TABLE]

where $u_{k,q}$ denotes the approximation of $u_{k}$ in $q$th iteration. Moreover, the remaining expressions indexed by $q$ correspond to expressions introduced to the solution by the $q$th iteration. For the sake of implementation, one can omit indexes $k$ and $m$. This yields relations:

[TABLE]

It is easy to prove, by induction and applying rules 1-4, that each approximation $u_{q}(\mathbf{x})$ can be rewritten in a form of a sum of non-nested $sum$ expressions and the term $f(\mathbf{x})+c_{q}$. Note that, the $Algebraic\_expand$ procedure in (B.2) keeps the summand expanded and let the rules work properly. The relation (B.2), as well as the rules 1-4 can be implemented in a Computer Algebra System (CAS). One can revert the symbolic expression $u_{q}(\mathbf{x})$ to an analytic formula, through two simple steps:

• •

symbols $c_{q}$, $a_{q}$ corresponds to $c_{k}$, $a_{k}$ in (B.3);

• •

symbols $c_{m}$, $a_{m}$ for $m<q$ transform to the summation, e.g. $\displaystyle sum(h(c_{3},c_{4},a_{3},a_{4}),[a_{3},a_{4}])\;\text{corresponds to}\;\sum_{\begin{subarray}{c}m\neq k\\ s\neq m\end{subarray}}h(c_{s},c_{m},a_{s},a_{m}).$

B.3 Application to the effective conductivity

Let us sketch the algorithm of computing the analytic approximation (3) of functions $u_{k}(\mathbf{x})$ explicitly up to $O(r_{0}^{q+1})$ using presented model. The procedure follows Sect. 3 and is expressed in the Mathematical Pseudo-language (MPL), described in details by Cohen [5], covering basic operations common for Computer Algebra Systems (CAS). The construction of the procedure is based on the application of following transformation rules, describing symbolic operations of substituting terms of a given form:

Rule 5.

$sum(f(\mathbf{a}),\mathbf{a})\rightarrow sum(Algebraic\_expand(f(\mathbf{a})),\mathbf{a})$

Rule 6.

$sum(f(\mathbf{a}),\mathbf{a})\rightarrow sum(Series(f(\mathbf{a}),r_{0},0,q),\mathbf{a})$

The $Series$ operator in the rule 6 produces a power series expansion for $f(\mathbf{a})$ with respect to $r_{0}$ about the point 0 to order $q$. Some CAS have the capability of implementation of transformation rules through the rule-based programming. Hence, the following procedure can be directly realized:

[TABLE]

Lines 8-11 find constants $c_{m}$ by calculating $z_{q}$ explicitly up to $O(r_{0}^{q+1})$ by successive approximations. Then, $c_{m}$ are substituted to the solution in lines 12-13. The $Reindex(expr,q,m)$ operator substitutes all symbols of the form $z_{m},a_{m}$ in $expr$ by $z_{2m-q},a_{2m-q}$ respectively, in order to prevent indices from duplicating in the resulting expression. Note that negative indices may appear, however it does not matter in the applied model of calculations. More details of the procedure $u(q)$, as well as an example implementation will be published in a separate paper.