# Numerical homogenization of fibrous composites with equivalent inclusion method

Antoine Martin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Antoine Martin. Numerical homogenization of fibrous composites with equivalent inclusion method. 2024. hal-04436391

HAL Id: hal-04436391

https://hal.science/hal-04436391v1

Preprint submitted on 3 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical homogenization of fibrous composites with equivalent inclusion method

A. Martin<sup>a,*</sup>

<sup>a</sup>Spie batignolles génie civil, 30 avenue du Général Gallieni, CS 10192, 92023 Nanterre CEDEX

## Abstract

The equivalent inclusion method is proposed for the homogenization of fibrous composites, as an alternative to costly full-field methods and less accurate mean-field or effective-field approaches. General methodology was presented in our previous study [14], where a discretized polarization field in fibers was introduced, which took advantage of inhomogeneity slenderness. Here, a method for computing apparent 3d conductivities in periodic conditions is proposed and implemented for numerical homogenization of fibrous composites. A wide range of aspect ratios and volume fractions are investigated, and the data are used for providing analytical formulas which fit the curves.

**Keywords:** Homogenization, Conductivity, Equivalent inclusion method, Cylindrical inclusions, Periodic microstructure,

## 1. Introduction

Fibers, for their high aspect ratio, are often used in industry or construction to improve performance of materials. Various properties can then be enhanced, in fields of physics that are seemingly different but in fact obey to very similar mathematical laws, such as mechanics [4], thermics [3], or electromagnetism [12]. Fibers are most of the time added to a homogeneous matrix, and this is a challenging work to predict homogenized properties of such a composite. In particular, it is known that for high volume fractions, some high interactions exist between fibers, whereas these interactions are almost invisible at lower volume fractions.

On the one hand, mean-field or effective-field homogenization methods [16, 25, 20, 10] succeed to take into account slenderness of fibers, and give good estimates of homogenized properties in the limit of a low volume fraction. However, these methods become inaccurate when this fraction increases.

On the other hand, full-field simulations [23, 18, 11] are made difficult due to the fineness of the mesh, which is governed not only by the small transverse dimension of the inhomogeneities but also by the distance between them. The computation then makes intervene a large number of unknowns, when meshing operation is possible.

In a previous study [14], we present an intermediate technique between mean-field and full-field methods, within the framework of electrical conductivity, which turns out to be more accurate than mean-field methods, and requires less unknowns than finite element computation. This technique is based on the equivalent inclusion method, first introduced by Moschovidis and Mura [17] and applied by various authors [6, 2, 22, 9, 26]. More precisely, it comes from a Galerkin-based variational form of the equivalent inclusion method, introduced by Brisard et al. [5] for spherical inhomogeneities within the framework of linear elasticity.

The study shows that this technique (referred to as the EIM below) is accurate and efficient and then can be used for homogenization of fibrous composites. This homogenization procedure was lacking in our first study and it is the purpose of the present paper to carry out this procedure. We introduce here a new framework for computing EIM conductivities, in comparison to the boundary conditions introduced before [see [14]]. The gain in terms of unknowns is huge compared to finite element method (FEM), but the homogenized values are very close.

*Corresponding author
Email address: antoine.martin@spiebatignolles.fr (A. Martin)
This paper is organized as follows. In Sec. 2, EIM for slender, cylindrical inhomogeneities, which was presented in our first study [14], is recalled. The method for taking into account long distance interactions is presented and implemented in Sec. 3. In Sec. 4, the homogenization procedure is carried out for an aspect ratio of 50, and is shown to be both accurate and efficient. This procedure is repeated in Sec. 5 for aspect ratios between 40 and 400, providing useful analytical formulas for homogenization of isotropic fibrous composites with high contrast. Finally, Sec. 6 closes this paper with a few concluding remarks and a discussion of future research directions.

## 2. Background: the EIM for slender, cylindrical inhomogeneities

The EIM developed in Martin et al. [14] is an adaptation of the variational form of EIM, introduced by Brisard et al. [5], to slender, cylindrical inhomogeneities. It is based on a particular discretization of the Lippmann–Schwinger equation and is recalled in the following. In Martin et al. [14], the method was presented with mixed boundary conditions, but EIM can also be implemented with periodic boundary conditions, used in this paper.

### 2.1. Discretization of the Lippmann–Schwinger equation

We consider a periodic cell $\Omega$ of a random heterogeneous material. The conductivity at $x \in \Omega$ is $\sigma(x)$ (symmetric, positive definite, second-order tensor); $E(x)$, $\phi(x)$ and $j(x)$ denote the electric field, the electric potential and the volumic current, respectively, at point $x$. The apparent conductivity of the cell $\Omega$, $\sigma^{\text{app}}$, is found from the solution to the following problem

$$ \Omega: \quad \operatorname{div} \mathbf{j} = 0, \qquad (1) $$

$$ \Omega: \quad \mathbf{j} = \boldsymbol{\sigma} \cdot \mathbf{E}, \qquad (2) $$

$$ \Omega: \quad \mathbf{E} = \bar{\mathbf{E}} + \operatorname{grad} \phi^{\text{per}}, \qquad (3) $$

where $\phi^{\text{per}}$ is a $\Omega$-periodic field, $\mathbf{j}$ is a $\Omega$-antiperiodic field, and $\bar{\mathbf{E}}$ is a prescribed constant vector. Eq. (3) ensures that the electric field is curl-free, with $\bar{\mathbf{E}} = \langle \mathbf{E} \rangle$, where angle brackets denote volume averages over the cell $\Omega$

$$ \langle \bullet \rangle = \frac{1}{|\Omega|} \int_{\Omega} \bullet(\mathbf{x}) \, d^3\mathbf{x}. \qquad (4) $$

Problem (1)–(3) can then be replaced by a unique integral equation of the Lippmann–Schwinger type [27]

$$ (\boldsymbol{\sigma} - \boldsymbol{\sigma}_0)^{-1} \cdot \boldsymbol{\tau} + \boldsymbol{\Gamma}_0^{\text{per}}(\boldsymbol{\tau}) = \bar{\mathbf{E}}, \qquad (5) $$

where the main unknown is the polarization $\boldsymbol{\tau}$

$$ \boldsymbol{\tau} = (\boldsymbol{\sigma} - \boldsymbol{\sigma}_0) \cdot \mathbf{E}, \qquad (6) $$

and $\boldsymbol{\Gamma}_0^{\text{per}}$ denotes the periodic Green operator, associated with the conductivity $\boldsymbol{\sigma}_0$, a symmetric, constant, positive definite, second-order tensor.

From the volume average $\langle \boldsymbol{\tau} \rangle$ of the solution $\boldsymbol{\tau}$ to the Lippmann–Schwinger equation (5), the apparent conductivity is readily retrieved. Indeed, using Eq. (6)

$$ \boldsymbol{\sigma}^{\text{app}} \cdot \bar{\mathbf{E}} = \langle \boldsymbol{\sigma} \cdot \mathbf{E} \rangle = \langle \boldsymbol{\sigma}_0 \cdot \mathbf{E} + \boldsymbol{\tau} \rangle = \boldsymbol{\sigma}_0 \cdot \langle \mathbf{E} \rangle + \langle \boldsymbol{\tau} \rangle = \boldsymbol{\sigma}_0 \cdot \bar{\mathbf{E}} + \langle \boldsymbol{\tau} \rangle \qquad (7) $$

We then follow a Galerkin procedure to derive an approximate solution to the Lippmann–Schwinger equation (5). Contracting Eq. (5) with a test function $\boldsymbol{\sigma} \in \mathcal{V}$ and volume averaging over $\Omega$ delivers the following variational problem

$$ \text{Find } \boldsymbol{\tau} \in \mathcal{V} \text{ such that, for all } \boldsymbol{\sigma} \in \mathcal{V}: \quad a(\boldsymbol{\tau}, \boldsymbol{\sigma}) = \bar{\mathbf{E}} \cdot \langle \boldsymbol{\sigma} \rangle \qquad (8) $$

where $\mathcal{V}$ is the set of square integrable $\Omega$-periodic vector fields, and $a$ is the following bilinear form

$$ a(\boldsymbol{\tau}, \boldsymbol{\sigma}) = \langle \boldsymbol{\sigma} \cdot (\boldsymbol{\sigma} - \boldsymbol{\sigma}_0)^{-1} \cdot \boldsymbol{\tau} \rangle + \langle \boldsymbol{\sigma} \cdot \boldsymbol{\Gamma}_0^{\text{per}}(\boldsymbol{\tau}) \rangle. \qquad (9) $$
Discretization of the above variational problem then follows a standard Galerkin procedure. Introducing the discretization parameter $p \in \mathbb{N}$ (to be defined) and the finite-dimensional discretization subspace $\mathbb{V}^p$ of $\mathbb{V}$, we consider the following variational problem

$$ \text{Find } \tau^p \in \mathbb{V}^p \text{ such that, for all } \boldsymbol{\sigma}^p \in \mathbb{V}^p: \quad a(\tau^p, \boldsymbol{\sigma}^p) = \bar{\mathbf{E}} \cdot \langle \boldsymbol{\sigma}^p \rangle. \qquad (10) $$

which differs from problem (8) only by the space to which the trial and test functions belong. Since the dimension of $\mathbb{V}^p$ is finite, solving problem (10) amounts to solving a linear system. The resulting solution $\tau^p$ approximates the solution to the initial problem (8). The variational form of the EIM results from the specialization of the above procedure to matrix-inhomogeneity media, and a specific class of discretization spaces $\mathbb{V}^p$.

## 2.2. The discretization for slender, cylindrical inhomogeneities

We therefore consider a periodic cell $\Omega$ that hosts $N$ inhomogeneities $\Omega_1, \dots, \Omega_N$ embedded in a homogeneous matrix $\Omega_0$. The inhomogeneities do not overlap, and each of them has homogeneous conductivity $\sigma_\alpha$ ($\alpha = 1, \dots, N$). Conductivity $\sigma_0$ is chosen to be equal to the matrix conductivity. $\tau$ is a $\Omega$-periodic field and $\Omega_\alpha$ may be not connex (when an inhomogeneity crosses the cell boundary).

The discretization space $\mathbb{V}^p$ is generated by a finite number of linearly independent functions supported on the inhomogeneities. More precisely, we seek the following decomposition for the trial function $\tau^p$

$$ \tau^p(\mathbf{x}) = \sum_{\alpha=1}^{N} \sum_{k=0}^{K_p-1} \tau_{\alpha}^k \Psi_{\alpha}^k(\mathbf{x}), \qquad (11) $$

where $\Psi_{\alpha}^k$ is a vector shape function supported in $\Omega_\alpha$, and $\tau_{\alpha}^k$ are scalar unknowns in $\mathbb{R}$. The number of shape functions $K_p \in \mathbb{N}$ must be specified.

**Remark 1.** Unless otherwise noted, greek indices ($\alpha, \beta, \dots$) span the $1, \dots, N$ range in the remainder of this paper (note that the matrix $\alpha = 0$ is not included in the sum), while latin indices ($k, l, \dots$) span the $0, \dots, K_p - 1$ range.

For the sake of simplicity, *monodisperse* assemblies only will be considered in the derivations presented below, and we introduce the common radius $R$ and common total length $2L$ of the cylindrical inhomogeneities. The aspect ratio $e$ of the cylinders is defined as: $e = L/R$. Inhomogeneity $\alpha$ is centered at $\mathbf{x}_\alpha \in \Omega$ and oriented by the unit-vector $\mathbf{n}_\alpha$. For such inhomogeneities, the longitudinal coordinate $z_\alpha$ defined as follows

$$ z_{\alpha} = (\mathbf{x} - \mathbf{x}_{\alpha}) \cdot \mathbf{n}_{\alpha} \qquad (12) $$

clearly plays a specific role. In our previous paper [14], the chosen shape functions are composed of a vector field polynomial in $z_\alpha$, with one component uniform in the transverse section, and the other component radial and slightly dependent of the cylinder radial coordinate.

Plugging this decomposition into the discrete variational problem (10), and testing with test functions $\boldsymbol{\sigma}^p \in \mathbb{V}^p$ decomposed similarly, the following linear system is derived

$$ \sum_l R_{\alpha}^{kl} \tau_l^l + \sum_{\beta,l} T_{\alpha\beta}^{kl} \tau_{\beta}^l = \mathbf{M}_{\alpha}^k \cdot \bar{\mathbf{E}}, \qquad (13) $$

where

$$ R_{\alpha}^{kl} = \langle \Psi_{\alpha}^{k} \cdot (\sigma_{\alpha} - \sigma_{0})^{-1} \cdot \Psi_{\alpha}^{l} \rangle_{\alpha}, \qquad (14) $$

$$ T_{\alpha\beta}^{kl} = \langle \Psi_{\alpha}^{k} \cdot \Gamma_{0}^{\text{per}}(\Psi_{\beta}^{l}) \rangle_{\alpha}, \qquad (15) $$

$$ \mathbf{M}_{\alpha}^{k} = \langle \Psi_{\alpha}^{k} \rangle_{\alpha}, \qquad (16) $$

where the volume average over inhomogeneity $\Omega_\alpha$, $\langle \cdot \rangle_\alpha$, is defined as follows

$$ \langle \cdot \rangle_{\alpha} = \frac{1}{|\Omega_{\alpha}|} \int_{\Omega_{\alpha}} \bullet(\mathbf{x}) \, d^3\mathbf{x}. \qquad (17) $$
Finally, Eq. (7) is applied to the solution to the linear system (13) to derive the EIM estimate of the apparent conductivity, $\sigma^{\text{EIM}}$

$$ \sigma^{\text{EIM}} \cdot \bar{\mathbf{E}} = \sigma_0 \cdot \bar{\mathbf{E}} + \sum_{\alpha,k} f_{\alpha} \tau_{\alpha}^{k} \mathbf{M}_{\alpha}^{k}. \qquad (18) $$

where $f_{\alpha} = |\Omega_{\alpha}|/|\Omega|$ denotes the volume fraction occupied by inhomogeneity $\alpha$ within the cell $\Omega$.

It must be emphasized that computation of $T_{\alpha\beta}^{kl}$ depends on periodic Green operator $\Gamma_0^{\text{per}}$. In our previous work [14], this computation was facilitated by the use of Green operator on the infinite medium, $\Gamma_0^{\infty}$, which can be expressed with an explicit formula. For a periodic medium, it is possible to compute these coefficients in Fourier space, as it is down for example in To et al. [24]. Unfortunately, the involved summations usually converge slowly. We therefore need to look for a better strategy for computing $T_{\alpha\beta}^{kl}$.

### 3. Computation of coefficients $T_{\alpha\beta}^{kl}$

$T_{\alpha\beta}^{kl}$ represents the $\Psi_{\alpha}^k$-weighted average of the (opposite of the) electric field induced on inhomogeneity $\alpha$ by a polarization $\Psi_{\beta}^l$ situated on inhomogeneity $\beta$. However, in periodic conditions, we can consider that the polarization is situated on $\beta$ and all its periodic images. Therefore, it is tempting to sum the effects of each polarization to obtain the coefficient, but it fails to give the true effect of surrounding inhomogeneities, as it will be shown in the following. The "good" strategy is given by Reaction Field Method, used in molecular physics.

**Remark 2.** It must be noted that the case $\alpha = \beta$ is particular since the polarization is situated on $\alpha$ (and also all its periodic images). In this case, we note $T_{\alpha\beta}^{kl}=S_{\alpha}^{kl}$.

#### 3.1. Reaction Field Method: the "good" strategy for computing interactions

The technique proposed for computing periodic Green operator will be based on the same concepts as Reaction Field Method (RFM) in molecular simulation [7, 13, 8]. In this method, the operator $\Gamma_0^{\text{RFM}}$ is not computed within the framework of periodic boundary conditions, but rather on an infinite medium. The approach is intuitive and $\Gamma_0^{\text{RFM}}$ is defined on a polarization field which is not normalizable. Indeed, this polarization field is defined on the infinite medium, and may not tend to 0 at infinity.

Considering an inhomogeneity $\alpha$ centered on $\mathbf{x}_{\alpha}$, surrounded by an infinite distribution of inhomogeneities, the purpose is to compute the electric field $\Gamma_0^{\text{RFM}}(\tau)$ on $\alpha$, induced by a polarization field $\tau$ defined on all inhomogeneities (including $\alpha$). The first idea is given by the well-known cavity problem used by various authors [15], which consists in assimilate the polarization field defined on surrounding inhomogeneities to a uniform polarization field situated outside a cavity, centered on $\mathbf{x}_{\alpha}$. This first approximation is relevant, but could not take into account the effect of nearby inhomogeneities. Hence, the second idea is to compute this effect with Green operator $\Gamma_0^{\infty}$, defined on an infinite medium.

Hence, for each $\mathbf{x}$ on inhomogeneity $\alpha$, the electric field $\Gamma_0^{\text{RFM}}(\tau)$ is proposed to be:

$$ \Gamma_0^{\text{RFM}}(\tau)(\mathbf{x}) = \Gamma_0^{\infty}(\chi_{\alpha}^d \tau)(\mathbf{x}) - \mathbf{P}_0 \cdot \langle \tau \rangle \qquad (19) $$

with $\mathbf{P}_0 = 1/(3\sigma_0)$ and $\chi_{\alpha}^d$ is indicator function of a bowl centered on $\mathbf{x}_{\alpha}$, with a radius $d$ (to be defined), called here interaction distance. This bowl defines a spherical domain, in which polarization field $\tau$ is non zero on the inhomogeneities, whereas outside this domain, $\tau$ is considered as a uniform field filling the space and equal to $\langle \tau \rangle$. Hence, the (opposite of the) electric field induced by the neighbouring inhomogeneities is given by Green operator in an infinite medium $\Gamma_0^{\infty}$, whereas electric field induced by the polarization outside the cavity is simply given by $\mathbf{P}_0 \cdot \langle \tau \rangle$.

**Remark 3.** Here, inhomogeneities situated at the limit of the spherical domain can be split in two parts, which would necessitate to compute interactions between $\alpha$ and half part of inhomogeneities. However, for the sake of simplicity, we consider that an inhomogeneity is situated inside the spherical domain if and only if its center is strictly inside.

**Remark 4.** Here, definition of $\langle \tau \rangle$ is not clear because the domain $\Omega$ on which computing the average has not been defined. This problem will be resolved in the following considering a periodic cell.
<figure><img src="image_3.png" /><figcaption>Figure 1: Green operator in RFM, for a periodic distribution of fibers</figcaption></figure>

Finally, the field $\Gamma_0^{\text{RFM}}(\tau)$ defined on $\alpha$ can also be written as

$$ \Gamma_0^{\text{RFM}}(\tau)(\mathbf{x}) = \sum_{\beta \in I_\alpha} \Gamma_0^\infty(\chi_\beta \tau)(\mathbf{x}) - \mathbf{P}_0 \cdot \langle \tau \rangle \qquad (20) $$

where $I_\alpha$ is the set of indices $\beta$ such that $\mathbf{x}_\beta$ is strictly inside the cavity centered on $\mathbf{x}_\alpha$ (note that $\alpha \in I_\alpha$). Hence, it must be pointed out that terms involving $\Gamma_0^\infty$ in the sum can now be computed as it was already done in our previous work [14] (the basic principles are recalled in Appendix A).

We can now see that the RFM helps us to define our periodic Green operator $\Gamma_0^{\text{per}}$. Indeed, considering that $\tau$ defines a $\Omega$-periodic field $\tau^{\text{per}}$, we propose (see Fig.1) that for each $\mathbf{x} \in \Omega$ situated on fiber $\alpha$ (which may not be connex),

$$ \Gamma_0^{\text{per}}(\tau^{\text{per}})(\mathbf{x}) = \Gamma_0^{\text{RFM}}(\tau)(\mathbf{x}) \qquad (21) $$

And the average $\langle \tau \rangle$ is then performed on the basic cell $\Omega$.

Coefficients $T_{\alpha\beta}^{kl}$ and $S_\alpha^{kl}$ are then computed with expression (15), with $\Gamma_0^{\text{per}}$ defined above.

## 3.2. Applications and assessment of the new operator

Our new operator then defines a new EIM, and we can now apply it to a periodic cell. Particularly, we want to show that this operator is relevant, and that there is two important contributions that can not be omitted: first, contribution of fibers situated in the spherical domain (terms given by the operator $\Gamma_0^\infty$), and second, contribution of further fibers situated outside the domain (given by $\mathbf{P}_0$).

We then consider a cubic periodic cell aligned with the global cartesian basis ($\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z$), with a uniform distribution of fibers inside, whose angular distribution is uniform and isotropic. A single fiber, noted 1, is placed on the center of the cell, along $\mathbf{e}_z$. The imposed field is $\bar{\mathbf{E}} = \bar{\mathbf{E}}\mathbf{e}_z$. We compute the new EIM field, for a distance $d$ equal to the cell size $D = 4L$. We also compute a field with FEM with COMSOL Multiphysics®5.6¹. Results on fiber 1 are shown on Fig. 2. We can see an excellent agreement between the two fields.

<footnote n="1">¹COMSOL Multiphysics® v. 5.6. www.comsol.com. COMSOL AB, Stockholm, Sweden, last retrieved 2023-05-13.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2: The longitudinal electric field in fiber 1, with FEM and new EIM</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Same computation as Fig. 2, omitting the term related to $P_0$ in Green operator $\Gamma_0^{\text{per}}$</figcaption></figure>

Moreover, we also compute the field with the new operator, but omitting the term related to $P_0$ (the "wrong" strategy). As we can see on Figure 3, it seems that this latter field strongly underestimates the finite element field. Changing $d$ does not affect significatively the result. Intuitively, it seems that fibers situated within the spherical domain, tend to decrease the field in fiber 1, as it can be seen on Fig. 3, whereas, as expected, the uniform polarization $\langle \tau \rangle$ increases the field in fiber 1.

Finally, we consider the same case and compute EIM field, for different sizes $d$ of the spherical domain. As we can see on Fig. 4, EIM field in fiber 1 seems to converge, when $d$ increases, to a limit which is in agreement with FE field, but we will not investigate this point more in this work. In the following, we just set $d$ equal to the cell size $D$.

# 4. Homogenization with the new EIM

## 4.1. Computing apparent and effective conductivities with new EIM

Because the method presented previously was shown to be accurate, it is now possible to compute apparent conductivities on periodic cells, for different types of microstructures. We begin with isotropic microstructures with high contrast, knowing that other cases could be considered in a further work. In this Section and in the remainder, contrast in set to $10^6$, and microstructures are isotropic, generated by Random Sequential Addition (RSA).

**Remark 5.** To check that cylinders do not overlap, it would have been long, in term of computation time, to perform the tests on perfect cylinders. Here, we choose to perform the tests on sphero-cylinders, which is more efficient, and then compute the fields on perfect cylinders.
<figure><img src="image_6.png" /><figcaption>Figure 4: Same computation as Fig. 2, increasing the value of d</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Mean apparent conductivities for different sizes of cell. The volume fraction is $f = 0.01$</figcaption></figure>

Particularly in this Section, the aspect ratio is $e = 50$, and other cases are considered in Section 5. We begin by computing apparent conductivities on several microstructures whose volume fraction of fibers is 0.01. We recall that the limit of apparent conductivities when size of the periodic cell tends to infinity is called *effective conductivity*. Without loss of generality, we can also define this effective conductivity as

$$ \sigma^{\text{eff}} = \lim_{D \to \infty} \overline{\sigma^{\text{app}}(D)} \qquad (22) $$

where $\overline{\sigma^{\text{app}}(D)}$ is the average of several apparent conductivities computed on cells of the same size $D$. Hence, it is interesting to compute apparent conductivities for different values of $D$, in order to study behaviour of $\overline{\sigma^{\text{app}}(D)}$ when $D$ increases. We were able to perform this calculation for $D$ between 4L and 8L, as shown on Fig. 5. We also plotted the 95 percent confidence interval, calculated with the estimated standard deviation obtained with the microstructures generated for each size $D$ (the number of microstructures is between 11 and 16). We can see here that the apparent conductivity is not very sensitive to the value of $D$.

In the following, we thus set $D = 4L$ and compute the corresponding averages of apparent conductivities on a few microstructures, for different volume fractions. We also made the finite element computation for a few microstructures with COMSOL Multiphysics®5.6. Results are shown on Fig. 6. We can see an excellent agreement between FEM and our EIM. Moreover, dilute scheme (DS) and Mori–Tanaka scheme (for ellipsoids of same aspect ratio) underestimate FEM and EIM values, whereas Ponte-Castañeda & Willis scheme [21] gives highly overestimated conductivities.
<figure><img src="image_8.png" /><figcaption>Figure 6: Mean apparent conductivities along volume fraction, for EIM and FEM, with $e = 50$, contrast of $10^6$, and $D = 4L$</figcaption></figure>

<table><caption>Table 1: Comparison of number of unknowns of EIM and FEM</caption><thead><tr><td>$f$</td><td>$n$ FEM</td><td>$n$ EIM</td></tr></thead><tbody><tr><td>0.01</td><td>~ 8M</td><td>9 144</td></tr><tr><td>0.02</td><td>~ 15M</td><td>18 324</td></tr><tr><td>0.03</td><td>~ 28M</td><td>27 468</td></tr><tr><td>0.04</td><td>~ 49M</td><td>36 648</td></tr></tbody></table>

## 4.2. Discussion on the number of unknowns

In terms of number $n$ of unknowns, the gain between EIM and FEM is huge, as we can see on Table 1, for a cell size $D = 4L$. As a consequence, solving time for EIM is also very small compared to FEM. Indeed, for $f = 0.01$, solving the linear system only takes 3 s for EIM, and 13 min 16 s for FEM. Few remarks can be made here.

**Remark 6.** For a more precise evaluation of computation time, we here give the characteristics of the computer used for values given in this Section:

* RAM: 32 Go

* Processor base speed : 2.71 GHz

* Cores : 6

* Logical processors : 12

First, at fixed aspect ratio $e$, with a RSA algorithm for generation of microstructures, there is a limit $f_c$ on the volume fraction that can be reached, which is known to be proportional to $1/e$ [19, 1]. Moreover, the effect of interactions between fibers is supposed to be greater around this volume fraction, and thus it is more relevant to compute apparent conductivities around $f_c$. It must also be noted that for a fixed volume fraction, and a fixed size of cell, the number of fibers grows quadratically with $e$, so that this number of fibers grows linearly with $e$ when computing a conductivity around the volume fraction of interest $f_c$. As a consequence, we can expect that all numerical methods whose complexity depends on the number of fibers (like EIM) will encounter a limit when $e$ increases.

Second, it must be noted on the base of the above consideration that EIM will always have less number of unknowns than FEM. Indeed, number of unknowns with EIM is about $N \times K_p$ where $N$ is the number of fibers and $K_p$ does not depend on the aspect ratio. Hence, the number of unknowns around $f_c$ grows also linearly with $e$. On the contrary, for FEM, it is expected that the number of unknowns is at least proportional to $N \times e$, because the number of
<figure><img src="image_9.png" /><figcaption>Figure 7: Mean apparent conductivities along volume fraction with EIM, for different aspect ratios (contrast of $10^6$, and $D = 4L$). The aspect ratios are precised at the end of each curve.</figcaption></figure>

nods on the mesh of a unique fiber is proportional to $e$. Then, the number of unknowns for FEM is at least quadratic in $e$ around $f_c$. Nevertheless, it is a lower bound on the number of unknowns, because space between fibers also needs to be meshed. For a high volume fraction, the number of nods is not clear, but if we consider that each element in the cell is of the order of magnitude of the radius $R$, the number of unknowns becomes cubic in the aspect ratio $e$. All these considerations explain the huge gain between FEM and EIM in terms of number of unknowns.

However, for both FEM and EIM, time for assembling the linear system, must also be counted. One of the drawbacks of the EIM proposed here is the time required for computing interaction coefficients, which is quadratic in the number of unknowns $n$, whereas for FEM, it is linear with $n$. However, it can be noted that this operation is embarrassingly parallel<fnref n="2" />, and therefore depends on the characteristics of the processor used. As a consequence, EIM computation may turn to be long for complex microstructures. As an example, at $f = 0.01$, EIM requires 11 min for assembling the linear system, whereas FEM requires 5 min 40 s. However, FEM also requires time for computing geometry (55 s), and mesh (8 min 10 s).

It must be noted that meshing operation cannot sometimes be performed without intervention of an operator. Indeed, if distance between fibers is very small, or if the cell boundary splits a fiber in several parts, with a very small one, element size becomes small and meshing operation is more difficult. As a consequence, the number of unknowns for FEM can vary when volume fraction is high, whereas number of unknowns for EIM is always the same. It makes also FEM more instable (we do not know the number of unknowns and if the computation will work), whereas EIM is more robust.

Finally, even if computation time can be long for both methods, the number of unknowns very small for EIM is a huge advantage when little RAM is available. As a consequence, EIM homogenization was possible in this paper for an aspect ratio of 400, and for high volume fractions (6% for an aspect ratio of 50), while FEM could not be carried out for these parameters.

## 5. Useful analytical formulas for homogenization of isotropic fibrous composites

Previous procedure shows that computing an apparent conductivity on a periodic cell of size is $D = 4L$ gives a relatively good estimate of effective conductivity. Therefore, this homogenization procedure has been performed on cells of size $D = 4L$, for different aspect ratios (between 40 and 400) and different volume fractions (below 6%), for the same high contrast of $10^6$. Results are shown on Fig. 7, and provide useful abacuses for design of fiber materials.

<footnote n="2">²https://en.wikipedia.org/wiki/Embarrassingly_parallel, last retrieved 2023-05-13.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 8: $c_1(e)$, based on cylinder Eshelby problem, fitted with shape function (24)</figcaption></figure>

However, because these abacuses show only results for a few cases, computation must be carried out once again for other aspect ratios or volume fractions. Moreover, as it was explained previously, computation time may be long due to high aspect ratios involved here and depending on computational resources. Therefore, in the following, we propose to make a numerical fitting of previous curves, owing to provide analytical formulas that give an idea of the effective conductivities without carrying out the full EIM computation.

We propose a Taylor expansion for the mean apparent conductivity, in the volume fraction, for each aspect ratio:

$$ \overline{\sigma}^{\text{app}}(e, f)/\sigma_0 = 1 + c_1(e)f + c_2(e)f^2 + c_3(e)f^3 + \dots \qquad (23) $$

However, coefficients of degree 1, $c_1(e)$ are already known. Indeed, for very low volume fractions, effective conductivity may be computed by mean-field methods, based on Eshelby's problem of a single inhomogeneity in an infinite medium. Here, for each aspect ratio, $c_1(e)$ have been computed with COMSOL Multiphysics®5.6, as shown on Fig. 8 in a log-log scale.

We then propose to fit $c_1(e)$ obtained with COMSOL Multiphysics®5.6, with another shape function of the aspect ratio. To that extent, we inspire ourselves on the ellipsoidal case, for which solution of Eshelby's problem can be computed with analytical formulas. Indeed, in this case, effective conductivity is, for high aspect ratios and infinite contrast, equivalent to $\sigma_0 + f e^2/(3 \ln e) \sigma_0$. We then propose to fit $c_1(e)$ with the following shape function:

$$ f_1(e) = a e^b / (\ln e)^c \qquad (24) $$

where $a, b, c \in \mathbb{R}$. We perform this fitting with function `curve_fit` of library `scipy` of Python 3.8.3. We then find $a = 0.386$, $b = 1.91$ and $c = 0.925$. This fitting is also shown on Fig. 8.

In the following, we truncate expansion (23), until order 5 ($c_i = 0$ for $i \ge 5$) and propose to realize a fitting of datas obtained with EIM for aspect ratios between 40 and 400. However, we also propose to use values $c_1(e)$ obtained with COMSOL Multiphysics®5.6. In other words, we seek for $c_2(e), c_3(e)$ and $c_4(e)$ such that

$$ \overline{\sigma}^{\text{EIM}}(e, f)/\sigma_0 - 1 - c_{\text{COM}}(e)f = c_2(e)f^2 + c_3(e)f^3 + c_4(e)f^4 \qquad (25) $$

where $c_{\text{COM}}(e)$ are previous values of $c_1(e)$ obtained with COMSOL Multiphysics®5.6. We also assume the following form of the coefficients $c_2, c_3, c_4$:

$$ c_i(e) = a_i e^{b_i} \qquad (26) $$

where $a_i, b_i \in \mathbb{R}$ ($i = 2, 3, 4$). With `scipy` library, we find $a_2 = 4.49 \times 10^{-2}$, $a_3 = -3.17 \times 10^{-3}$, $a_4 = 2.20 \times 10^{-4}$, and $b_2 = 2.97$, $b_3 = 4.39$ and $b_4 = 5.60$. It then gives the following approximation of mean apparent conductivity as a function of aspect ratio $e$ and volumic fraction of fibers $f$:

$$ \overline{\sigma}^{\text{app}}(e, f)/\sigma_0 \approx 1 + (0.386 e^{1.91} / (\ln e)^{0.925}) f + 0.0449 e^{2.97} f^2 - 0.00317 e^{4.39} f^3 + 0.00022 e^{5.6} f^4 \qquad (27) $$
## 6. Conclusion and outlook

In this paper, a new framework for application of the equivalent inclusion method to fibrous composites is proposed. It is based on a Green operator often used in molecular simulations, and arises naturally from the need of an infinite medium computation. It is compared with FEM with periodic boundary conditions and results show that this new EIM is very accurate and efficient for homogenization. We then computed a large set of apparent conductivities for different volume fractions (below 10%) and aspect ratios (between 40 and 400), which probably could not have been possible with FEM. It allowed us to propose abacuses and analytical formulas which can be used as good estimates for homogenized conductivities of isotropic fibrous composites.

Nevertheless, we highlighted some limits of our EIM, like computation time of interaction coefficients, or dependence of the number of unknowns on the aspect ratio. Moreover, as it is explained in our previous study [14], computation of interactions is approximated considering that distance between fibers is large in front of fiber diameter. This means that our computation must be wrong for sufficiently higher volume fractions. In our paper, we showed that FEM and EIM were in agreement when FE computation was possible. However, when only EIM is possible, the estimation of the error made is not clear.

Finally, as an interesting perspective, the case of linear elasticity, which makes intervene very similar equations, may be studied in a further work.

## Acknowledgments

We acknowledge *Spie batignolles génie civil* and *Laboratoire Navier* for providing financial and technical resources to this work.

## Appendix A. Computation details for interaction and self-influence coefficients

For the practical implementation of the method, it is convenient to write the interaction and self-influence coefficients in a 4 × 4 matrix form:

$$ T_{\alpha\beta}^{mn} = \begin{pmatrix} T_{\alpha\beta}^{m,n} & T_{\alpha\beta}^{m,n+p+1} & T_{\alpha\beta}^{m,n+2(p+1)} & T_{\alpha\beta}^{m,n+3(p+1)} \\ T_{\alpha\beta}^{m+p+1,n} & T_{\alpha\beta}^{m+p+1,n+p+1} & T_{\alpha\beta}^{m+p+1,n+2(p+1)} & T_{\alpha\beta}^{m+p+1,n+3(p+1)} \\ T_{\alpha\beta}^{m+2(p+1),n} & T_{\alpha\beta}^{m+2(p+1),n+p+1} & T_{\alpha\beta}^{m+2(p+1),n+2(p+1)} & T_{\alpha\beta}^{m+2(p+1),n+3(p+1)} \\ T_{\alpha\beta}^{m+3(p+1),n} & T_{\alpha\beta}^{m+3(p+1),n+p+1} & T_{\alpha\beta}^{m+3(p+1),n+2(p+1)} & T_{\alpha\beta}^{m+3(p+1),n+3(p+1)} \end{pmatrix} \quad (A.1) $$

and

$$ S_{\alpha}^{mn} = \begin{pmatrix} S_{\alpha}^{m,n} & S_{\alpha}^{m,n+p+1} & S_{\alpha}^{m,n+2(p+1)} & S_{\alpha}^{m,n+3(p+1)} \\ S_{\alpha}^{m+p+1,n} & S_{\alpha}^{m+p+1,n+p+1} & S_{\alpha}^{m+p+1,n+2(p+1)} & S_{\alpha}^{m+p+1,n+3(p+1)} \\ S_{\alpha}^{m+2(p+1),n} & S_{\alpha}^{m+2(p+1),n+p+1} & S_{\alpha}^{m+2(p+1),n+2(p+1)} & S_{\alpha}^{m+2(p+1),n+3(p+1)} \\ S_{\alpha}^{m+3(p+1),n} & S_{\alpha}^{m+3(p+1),n+p+1} & S_{\alpha}^{m+3(p+1),n+2(p+1)} & S_{\alpha}^{m+3(p+1),n+3(p+1)} \end{pmatrix} \quad (A.2) $$

with $0 \le m, n \le p$ where $p$ is the order of polynomial in $z_{\alpha}$. Here, latin, superior indices to the matrices T and S are assumed to span 0, ..., $p$.

*Interaction coefficients.* For two distinct inhomogeneities $\alpha \neq \beta$, we have

$$ T_{\alpha\beta}^{mn} = \sigma_0^{-1} (U_{\alpha\beta}^{mn} - V_{\alpha\beta}^{mn}) \quad (A.3) $$

where $U_{\alpha\beta}^{mn}$ and $V_{\alpha\beta}^{mn}$ are 4 × 4 matrices. We have

$$ (U_{\alpha\beta}^{mn})_{ij} = \frac{1}{4\pi V_{\alpha}} \int_{y_{\alpha} \in \Omega_{\alpha}} \int_{y_{\beta} \in \Omega_{\beta}} \frac{z_{\alpha}^{m} z_{\beta}^{n} \delta_{ij}}{\|y_{\beta} - y_{\alpha}\|^3} d^3 y_{\beta} d^3 y_{\alpha}, \quad (A.4) $$
and

$$ (\mathbf{U}_{\alpha\beta}^{mn})_{4j} = \frac{1}{4\pi V_{\alpha}} \int_{\mathbf{y}_{\alpha} \in \Omega_{\alpha}} \int_{\mathbf{y}_{\beta} \in \Omega_{\beta}} \frac{z_{\alpha}^{m} z_{\beta}^{n}}{\|\mathbf{y}_{\beta} - \mathbf{y}_{\alpha}\|^3} r_{\alpha} \mathbf{e}_{r,\alpha} \cdot \mathbf{e}_{j} \, d^3\mathbf{y}_{\beta} \, d^3\mathbf{y}_{\alpha}, \quad (A.5) $$

for i, j = 1, 2, 3, and

$$ (\mathbf{U}_{\alpha\beta}^{mn})_{44} = \frac{1}{4\pi V_{\alpha}} \int_{\mathbf{y}_{\alpha} \in \Omega_{\alpha}} \int_{\mathbf{y}_{\beta} \in \Omega_{\beta}} \frac{z_{\alpha}^{m} z_{\beta}^{n}}{\|\mathbf{y}_{\beta} - \mathbf{y}_{\alpha}\|^3} r_{\alpha} r_{\beta} \mathbf{e}_{r,\alpha} \cdot \mathbf{e}_{r,\beta} \, d^3\mathbf{y}_{\beta} \, d^3\mathbf{y}_{\alpha}. \quad (A.6) $$

We also have

$$ (\mathbf{V}_{\alpha\beta}^{mn})_{ij} = \frac{3}{4\pi V_{\alpha}} \int_{\mathbf{y}_{\alpha} \in \Omega_{\alpha}} \int_{\mathbf{y}_{\beta} \in \Omega_{\beta}} z_{\alpha}^{m} z_{\beta}^{n} \frac{[(\mathbf{y}_{\beta} - \mathbf{y}_{\alpha}) \cdot \mathbf{e}_{i}] [(\mathbf{y}_{\beta} - \mathbf{y}_{\alpha}) \cdot \mathbf{e}_{j}]}{\|\mathbf{y}_{\beta} - \mathbf{y}_{\alpha}\|^5} \, d^3\mathbf{y}_{\beta} \, d^3\mathbf{y}_{\alpha}, \quad (A.7) $$

and

$$ (\mathbf{V}_{\alpha\beta}^{mn})_{4j} = \frac{3}{4\pi V_{\alpha}} \int_{\mathbf{y}_{\alpha} \in \Omega_{\alpha}} \int_{\mathbf{y}_{\beta} \in \Omega_{\beta}} z_{\alpha}^{m} z_{\beta}^{n} r_{\alpha} \frac{\mathbf{e}_{r,\alpha} \cdot (\mathbf{y}_{\beta} - \mathbf{y}_{\alpha})}{\|\mathbf{y}_{\beta} - \mathbf{y}_{\alpha}\|^5} (\mathbf{y}_{\beta} - \mathbf{y}_{\alpha})_j \, d^3\mathbf{y}_{\beta} \, d^3\mathbf{y}_{\alpha}, \quad (A.8) $$

for i, j = 1, 2, 3, and

$$ (\mathbf{V}_{\alpha\beta}^{mn})_{44} = \frac{3}{4\pi V_{\alpha}} \int_{\mathbf{y}_{\alpha} \in \Omega_{\alpha}} \int_{\mathbf{y}_{\beta} \in \Omega_{\beta}} z_{\alpha}^{m} z_{\beta}^{n} r_{\alpha} r_{\beta} \frac{\mathbf{e}_{r,\alpha} \cdot (\mathbf{y}_{\beta} - \mathbf{y}_{\alpha})}{\|\mathbf{y}_{\beta} - \mathbf{y}_{\alpha}\|^5} (\mathbf{y}_{\beta} - \mathbf{y}_{\alpha}) \cdot \mathbf{e}_{r,\beta} \, d^3\mathbf{y}_{\beta} \, d^3\mathbf{y}_{\alpha}. \quad (A.9) $$

We were not able to derive a closed-form expression of the above integrals for two cylinders. However, assuming that the radius $R$ is small compared to the smallest distance between the two cylinders, a multipole expansion can be produced. For low volume fractions, this assumption is certainly verified for most pairs of cylinders. We then propose the following approximation for $i, j = 1, 2, 3$,

$$ (\mathbf{U}_{\alpha\beta}^{mn})_{ij} = \frac{L^{m+n+3}}{8} \frac{R^2}{L^2} \int_{-1}^{1} \int_{-1}^{1} \frac{\zeta_{\alpha}^{m} \zeta_{\beta}^{n} \delta_{ij}}{\|\mathbf{w}^0\|^3} \, d\zeta_{\alpha} \, d\zeta_{\beta}, \quad (A.10) $$

$$ (\mathbf{V}_{\alpha\beta}^{mn})_{ij} = \frac{3L^{m+n+3}}{8} \frac{R^2}{L^2} \int_{-1}^{1} \int_{-1}^{1} \frac{\zeta_{\alpha}^{m} \zeta_{\beta}^{n} w_{i}^{0} w_{j}^{0}}{\|\mathbf{w}^0\|^5} \, d\zeta_{\alpha} \, d\zeta_{\beta}. \quad (A.11) $$

with

$$ \mathbf{w}^0 = \mathbf{r}_{\alpha\beta} + L\zeta_{\beta} \mathbf{n}_{\beta} - L\zeta_{\alpha} \mathbf{n}_{\alpha} \quad \text{and} \quad \mathbf{r}_{\alpha\beta} = \mathbf{x}_{\beta} - \mathbf{x}_{\alpha}. \quad (A.12) $$

And for $i = 4$ or $j = 4$, we can prove that the first term of the multipole expansion of $(\mathbf{U}_{\alpha\beta}^{mn})_{ij}$ and $(\mathbf{V}_{\alpha\beta}^{mn})_{ij}$ is null [see [14] for details].

In the above expressions, integrals can be evaluated numerically. However, this method becomes inefficient for a large number of inclusions, in which case it is more advantageous to evaluate the first integral (with respect to $\zeta_{\alpha}$) analytically [see [14] for details].

**Self-influence coefficients.** When $\alpha = \beta$, computing strategy is different because the principal value cannot be removed, which makes the analytical evaluation of $S_{\alpha}^{kl}$ difficult. Given that this coefficient can be precomputed off-line prior to the full EIM calculation, a numerical approach is used, based on a 2d finite element computation [see [14] for details].
References

[1] Balberg, I., Binenbaum, N., Wagner, N., 1984. Percolation Thresholds in the Three-Dimensional Sticks System. Physical Review Letters 52, 1465-1468. doi:10.1103/PhysRevLett.52.1465.

[2] Benedikt, B., Lewis, M., Rangaswamy, P., 2006. On elastic interactions between spherical inclusions by the equivalent inclusion method. Computational Materials Science 37, 380-392. doi:10.1016/j.commatsci.2005.10.002.

[3] Biercuk, M.J., Llaguno, M.C., Radosavljevic, M., Hyun, J.K., Johnson, A.T., Fischer, J.E., 2002. Carbon nanotube composites for thermal management. Applied Physics Letters 80, 2767-2769. doi:10.1063/1.1469696.

[4] Brandt, A.M., 2008. Fibre reinforced cement-based (FRC) composites after over 40 years of development in building and civil engineering. Composite Structures 86, 3-9. doi:10.1016/j.compstruct.2008.03.006.

[5] Brisard, S., Dormieux, L., Sab, K., 2014. A variational form of the equivalent inclusion method for numerical homogenization. International Journal of Solids and Structures 51, 716-728. doi:10.1016/j.ijsolstr.2013.10.037.

[6] Fond, C., Géhant, S., 2002. Effects of mechanical interactions on the hydrostatic stress in randomly distributed rubber particles in an amorphous polymer matrix. Polymer 43, 909-919. doi:10.1016/S0032-3861(01)00621-8.

[7] Fukuda, I., Nakamura, H., 2012. Non-Ewald methods: Theory and applications to molecular systems. Biophysical Reviews 4, 161-170. doi:10.1007/s12551-012-0089-4.

[8] Hansson, T., Oostenbrink, C., van Gunsteren, W.F., J. Molecular dynamics simulations.

[9] Hiroshi, H., Minoru, T., 1986. Equivalent inclusion method for steady state heat conduction in composites. International Journal of Engineering Science 24, 1159-1172. doi:10.1016/0020-7225(86)90011-X.

[10] Hori, M., Double-inclusion model and overall moduli of multi-phase composites.

[11] Kabel, M., Merkert, D., Schneider, M., 2015. Use of composite voxels in FFT-based homogenization. Computer Methods in Applied Mechanics and Engineering 294, 168-188. doi:10.1016/j.cma.2015.06.003.

[12] Kong, L.B., Li, Z.W., Liu, L., Huang, R., Abshinova, M., Yang, Z.H., Tang, C.B., Tan, P.K., Deng, C.R., Matitsine, S., 2013. Recent progress in some composite materials and structures for specific electromagnetic applications. International Materials Reviews 58, 203-259. doi:10.1179/1743280412Y.0000000011.

[13] Lee, F.S., Warshel, A., 1992. A local reaction field method for fast evaluation of long-range electrostatic interactions in molecular simulations. The Journal of Chemical Physics 97, 3100-3107. doi:10.1063/1.462997.

[14] Martin, A., Brisard, S., Deleville, S., Sab, K., 2023. Assessment of the equivalent inclusion method for the numerical homogenization of fibrous composites. Journal of Computational Physics 477, 111943. doi:10.1016/j.jcp.2023.111943.

[15] Maxwell Garnett, J.C., 1904. Colours in metal glasses and in metallic films. Philosophical Transactions of the Royal Society A 203, 385-420. doi:10.1098/rsta.1904.0024.

[16] Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica 21, 571-574. doi:10.1016/0001-6160(73)90064-3.

[17] Moschovidis, Z.A., Mura, T., 1975. Two-Ellipsoidal Inhomogeneities by the Equivalent Inclusion Method. Journal of Applied Mechanics 42, 847-852. doi:10.1115/1.3423718.

[18] Moulinec, H., Suquet, P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Computer Methods in Applied Mechanics and Engineering 157, 69-94. doi:10.1016/S0045-7825(97)00218-1, arXiv:2012.08962.

[19] Néda, Z., Florian, R., Brechet, Y., 1999. Reconsideration of continuum percolation of isotropically oriented sticks in three dimensions. Physical Review E 59, 3717-3719. doi:10.1103/PhysRevE.59.3717.

[20] Norris, A.N., 1985. A differential scheme for the effective moduli of composites. Mechanics of materials 4, 1-16.

[21] Ponte Castañeda, P., Willis, J.R., 1995. The effect of spatial distribution on the effective behavior of composite materials and cracked media. Journal of the Mechanics and Physics of Solids 43, 1919-1951. doi:10.1016/0022-5096(95)00058-Q.

[22] Rodin, G.J., Hwang, Y.L., 1991. On the problem of linear elasticity for an infinite region containing a finite number of non-intersecting spherical inhomogeneities. International Journal of Solids and Structures 27, 145-159. doi:10.1016/0020-7683(91)90225-5.

[23] Sab, K., 1992. On the homogenization and the simulation of random materials. European journal of mechanics. A. Solids 11, 585-607.

[24] To, Q.D., Bonnet, G., Hoang, D.H., 2016. Explicit effective elasticity tensors of two-phase periodic composites with spherical or ellipsoidal inclusions. International Journal of Solids and Structures 94-95, 100-111. doi:10.1016/j.ijsolstr.2016.05.005.

[25] Walpole, L.J., 1969. On the overall elastic moduli of composite materials. Journal of the Mechanics and Physics of Solids 17, 235-251. doi:10.1016/0022-5096(69)90014-3.

[26] Yin, H.M., Lee, PH., Liu, YJ., 2014. Equivalent inclusion method for the Stokes flow of drops moving in a viscous fluid. Journal of Applied Mechanics 81.

[27] Zeller, R., Dederichs, P.H., 1973. Elastic Constants of Polycrystals. physica status solidi (b) 55, 831-842. doi:10.1002/pssb.2220550241.