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

HAL
open science

# Homogenization of fiber materials with equivalent inclusion method

Antoine Martin, Sébastien Brisard

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

## To cite this version:

Antoine Martin, Sébastien Brisard. Homogenization of fiber materials with equivalent inclusion method. 25e Congrès Français de Mécanique 2022, Association Française de Mécanique, Aug 2022, Nantes, France. 10 p. hal-04089803

HAL Id: hal-04089803

https://hal.science/hal-04089803

Submitted on 5 May 2023

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.

Public Domain
# Homogenization of fiber materials with equivalent inclusion method

A. MARTIN<sup>a,b</sup>, S. BRISARD<sup>b</sup>

a. Spie batignolles génie civil, 30 avenue du Général Gallieni, CS 10192, 92023 Nanterre CEDEX –
antoine.martin@spiebatignolles.fr

b. Navier, Ecole des Ponts, Univ Gustave Eiffel, IFSTTAR, CNRS, Marne-la-Vallée, France –
sebastien.brisard@univ-eiffel.fr

## Résumé :

Grâce à un coût de calcul relativement bas, la méthode de l'inclusion équivalente [8] est une alternative intéressante aux méthodes traditionnelles de calcul en champ complet, pour des microstructures de type matrice/inclusion. Brisard et al. [2] ont introduit une formulation variationnelle de cette méthode, bien adaptée à l'homogénéisation. Dans cet article, cette méthode variationnelle est mise en oeuvre dans le cas des fibres, dont l'élancement permet de décrire le champ de polarisation avec peu de degrés de liberté, en assimilant les fibres à des éléments linéiques. Les résultats sont présentés ici sur quelques fibres, pour des problèmes de conductivité. Ils peuvent être étendus à des problèmes d'élasticité.

## Abstract :

Due to its low computational cost, the equivalent inclusion method [8] is an attractive alternative to traditional full-field computations of heterogeneous materials formed of simple inhomogeneities embedded in a homogeneous matrix. Brisard et al. [2] introduced a variational form of this method, well-suited to homogenization. In this paper, we implement this variational method in the case of fibres, whose aspect ratio allows to describe the internal fields with few degrees of freedom. The accuracy of the method is proved for 3D conductivity problems with a few fibres.

## Mots clefs : Homogenization, Conductivity, Equivalent inclusion method, Lippmann-Schwinger equation

## 1 Introduction

Recently, Spie batignolles and EuroMC have developed an electromagnetic shielding process for confidential facilities, based on a concrete reinforced by conducting fibers [7]. Within the framework of homogenization, we study the effective electric conductivity of such heterogeneous materials.

Our goal is to propose an estimation of the homogenized conductivity, that improves on 'mean field' models, while being less costly computationally than 'full field' methods.
In order to get around this difficulty, we have chosen an intermediate technique, the equivalent inclusion method, based on Eshelby's work [5] and extended to several inclusions by Moschovidis et al. [8]. We use a variational form of the equivalent inclusion method, introduced by Brisard et al. [2], whose benefits are the following: (i) the resulting linear system is well-posed, (ii) the numerical solution converges to the exact solution, (iii) the method can provide rigorous bounds on the apparent properties of the statistical volume element. This method is called EIM in the following. The application of such method can be found in electrostatics (homogenized conductivity) but also in mechanics (homogenized stiffness), thermal science...

In this paper, we present the implementation of the EIM in the case of fibres in 3D, which is a case which would require too much computational resources for a full field simulation. We take advantage of fibre aspect ratio to reduce the number of DOF. We assess the accuracy of the method through 3D conductivity problems involving a limited number of cylinders.

The paper is organized as follows. Section 2 is a presentation of the variational form of the equivalent inclusion method [2] and of its implementation in the case of fibres. In Section 3, we present a first application on one fibre and another application with two fibres.

## 2 Equivalent inclusion method for fibres

In this section, we implement the variational form of the equivalent inclusion method introduced by Brisard et al. [2] in the case of linear conductivity, for cylindrical inclusions.

The modified Lippmann-Schwinger equation is first introduced in Sec. 2.1 and its Galerkin discretization (for cylinders) is presented in Sec. 2.2.

Following Ostoja-Starzewski [9], we use the terminology statistical volume element to denote a large (but finite) realization of the random heterogeneous material under consideration.

### 2.1 The modified Lippmann-Schwinger equation

Within the framework of electric conductivity, we consider a statistical volume element $\Omega$ of a random heterogeneous material. The conductivity at $x \in \Omega$ is $\sigma(x)$ (symmetric, positive definite, second-order tensor); besides, $E(x)$, $\phi(x)$ and $j(x)$ denote the electric field, electric potential and volumic current, respectively, at point $x$.

The apparent conductivity of the SVE $\Omega$, $\sigma^{app}$, is found from the solution to the following boundary-value problem

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

$$ \Omega: \quad j = \sigma \cdot E, \qquad (2) $$

$$ \Omega: \quad E = \mathrm{grad} \, \phi, \qquad (3) $$

with boundary conditions that must be specified. The boundary conditions chosen here, in order for the equivalent inclusion method to be applicable, are mixed boundary conditions, previously discussed within the framework of linear elasticity [3, 2].

In order to apply these boundary conditions, the SVE $\Omega$ must be embedded in a homogeneous, infinite matrix with conductivity $\sigma_0$; $\chi$ denotes the indicator function of $\Omega$ ($\chi(x) \in \{0, 1\}$ and $\chi(x) = 1$ iff $x \in \Omega$). It was shown in Brisard et al. [3] that the introduced mixed boundary conditions (not stated in
the present paper for the sake of conciseness) allow to reformulate problem (1)–(3) as a unique integral equation of the Lippmann–Schwinger type [12, 14] :

$$ (\sigma - \sigma_0)^{-1} \cdot \tau + \Gamma_0^{\infty}(\tau - \chi\langle\tau\rangle) = \bar{E}, \qquad (4) $$

where the macroscopic electric field $\bar{E}$ is prescribed at infinity :

$$ \|\mathbf{x}\| \to +\infty : \quad \mathbf{E}(\mathbf{x}) \to \bar{\mathbf{E}}, \qquad (5) $$

and the main unknown is the polarization $\tau$

$$ \tau = (\sigma - \sigma_0) \cdot \mathbf{E}. \qquad (6) $$

In the above modified Lippmann–Schwinger equation, $\Gamma_0^\infty$ denotes the Green operator, associated with the conductivity $\sigma_0$. It maps the field $\tau$ onto the field $\Gamma_0^\infty(\tau)$. In the remainder of this paper, it will always be assumed that the reference material is isotropic : $\sigma = \sigma_0 \mathbf{1}$. Then, we have the following expression of the Green operator [4, 11]

$$ \Gamma_0^\infty(\tau)(\mathbf{x}) = \mathbf{P}_0 \cdot \tau(\mathbf{x}) + \mathrm{PV}_{\mathbf{x}} \int_{\mathbf{y} \in \Omega} \mathbf{G}_0(\mathbf{y} - \mathbf{x}) \cdot \tau(\mathbf{y}) \, d^3\mathbf{y}, \qquad (7) $$

where

$$ \mathbf{P}_0 = \frac{1}{3\sigma_0} \quad \text{and} \quad \mathbf{G}_0(\mathbf{r}) = \frac{1 - 3\mathbf{n} \otimes \mathbf{n}}{4\pi\sigma_0 r^3} \qquad (8) $$

$(r = \|\mathbf{r}\|$ and $\mathbf{n} = \mathbf{r}/r)$. In the above equation, “PV$_x$” refers to the principal value at $x$ for spherical excluded regions

$$ \mathrm{PV}_{\mathbf{x}} \int_{\mathbf{y} \in \Omega} f(\mathbf{y}) \, d^3\mathbf{y} = \lim_{\delta \to 0} \int_{\substack{\mathbf{y} \in \Omega \\ \|\mathbf{y} - \mathbf{x}\| \ge \delta}} f(\mathbf{y}) \, d^3\mathbf{y} \qquad (9) $$

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

$$ \sigma^{\mathrm{app}} \cdot \bar{\mathbf{E}} = \sigma_0 \cdot \langle \mathbf{E} \rangle + \langle \tau \rangle = \sigma_0 \cdot \bar{\mathbf{E}} + \langle \tau \rangle \qquad (10) $$

To close this section, it is emphasized that the equivalence between, on the one hand problem (1)–(3) with mixed boundary conditions and, on the other hand, the modified Lippmann–Schwinger equation (4) holds for SVEs $\Omega$ of ellipsoidal shape; of course, the limit for large SVEs of the apparent conductivity does not depend on the radii of the ellipsoid. In the remainder of this paper, the microstructures under consideration are statistically isotropic. It is therefore natural to use spherical SVEs.

## 2.2 Discretization of the Lippmann–Schwinger equation

The Lippmann–Schwinger equation (4) is first reformulated as a variational problem as follows

$$ \text{Find } \tau \in \mathcal{V} \text{ such that, for all } \varpi \in \mathcal{V}: \qquad \langle \varpi \cdot (\sigma - \sigma_0)^{-1} \cdot \tau \rangle + \langle \varpi \cdot \Gamma_0^\infty(\tau - \chi\langle\tau\rangle) \rangle = \bar{\mathbf{E}} \cdot \langle \varpi \rangle \qquad (11) $$
Following a Galerkin approach, we then introduce the finite-dimensional discretization subspace $V^p$ of $V$, ($p \in \mathbb{N}$ is a discretization parameter to be defined), and search for a solution of (11) in $V^p$, denoted $\tau^p$.

In the following, we consider a spherical SVE $\Omega \subset \mathbb{R}^3$ that hosts $N$ cylindrical inhomogeneities $\Omega_1, \dots, \Omega_N$ embedded in a homogeneous matrix $\Omega_0$. We introduce the common radius $R$ and common total length $2L$ of the fibres. Then, $\Omega_\alpha$ is a circular cylinder centered at $\mathbf{x}_\alpha \in \Omega$, oriented by the unit-vector $\mathbf{n}_\alpha$ (see Fig. 2). The aspect ratio $e$ of the fibers is defined as : $e = L/R$. The inhomogeneities do not overlap. The corresponding indicator functions $\chi_\alpha$ ($\alpha = 0, \dots, N$) are

$$ \chi_{\alpha}(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \in \Omega_{\alpha}, \\ 0 & \text{otherwise.} \end{cases} \qquad (12) $$

It will be convenient to introduce volume averages over inhomogeneity $\Omega_\alpha$

$$ \langle \bullet \rangle_{\alpha} = \frac{1}{V_{\alpha}} \int_{\Omega_{\alpha}} \bullet(\mathbf{x}) \, d^3\mathbf{x}, \qquad (13) $$

where $V_\alpha$ denotes the volume of $\Omega_\alpha$. Finally, $f_\alpha = V_\alpha/V$ denotes the volume fraction occupied by inhomogeneity $\alpha$ within the SVE.

Within the framework of the EIM, the reference material must coincide with the matrix. In other words, the SVE $\Omega$ is embedded in an infinite, homogeneous material with same conductivity $\sigma_0$ as the matrix. Similarly, each inhomogeneity $\alpha = 1, \dots, N$ has homogeneous conductivity $\sigma_\alpha$.

For the applications considered in this paper, it is assumed that the matrix is isotropic : $\sigma_0 = \sigma_0 \mathbf{1}$. The discretization space $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} \Psi_{\alpha}^k(\mathbf{x}) \tau_{\alpha}^k, \qquad (14) $$

where $K_p$ is the number of scalar shape functions $\Psi_\alpha^k$ that are supported in $\Omega_\alpha$ and $\tau_\alpha^k$ are constant, unknown vectors.

Typically, for slender cylindrical inhomogeneities, the longitudinal coordinate $z_\alpha$ defined as follows

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

plays a specific role. This suggests to consider polarizations that are polynomials of high-order $p$ of the longitudinal coordinate, and polynomials of low-order $q$ of the two other, transverse, coordinates. Here, only the simplest case $q = 0$ was investigated. In other words, the discretization space $V^p$ is now the space of tensor fields that are, over each inhomogeneity $\Omega_\alpha$, polynomial of the local longitudinal coordinate $z_\alpha$. We then have $K_p = p$ and

$$ \Psi_{\alpha}^{k}(\mathbf{x}) = \chi_{\alpha}(\mathbf{x}) z_{\alpha}^{k}. \qquad (16) $$

Plugging this decomposition into the discrete variational problem, and testing with test functions $\tau^p \in$
<figure><img src="image_3.png" /><figcaption>FIGURE 1 – Problem on a RVE</figcaption></figure>

V$^p$ decomposed similarly, the following linear system is derived

$$ \sum_l [M_{\alpha}^{k+l}(\sigma_{\alpha} - \sigma_0)^{-1} + S_{\alpha}^{kl} - f_{\alpha} M_{\alpha}^k M_{\alpha}^l \mathbf{P}_0] \cdot \mathbf{\tau}_{\alpha}^l + \sum_{\beta \neq \alpha, l} (\mathbf{T}_{\alpha\beta}^{kl} - f_{\beta} M_{\alpha}^k M_{\beta}^l \mathbf{P}_0) \cdot \mathbf{\tau}_{\beta}^l = M_{\alpha}^k \bar{\mathbf{E}}, \quad (17) $$

where

$$ M_{\alpha}^{k} = \frac{1 + (-1)^{k}}{2(k + 1)} L^{k}, \qquad (18) $$

while $T_{\alpha\beta}^{kl}$ denotes the linear mapping over the space of 3d vectors defined below

$$ T_{\alpha\beta}^{kl} \cdot \tau = \langle \Psi_{\alpha}^{k} \Gamma_{0}^{\infty} (\Psi_{\beta}^{l} \tau) \rangle_{\alpha}. \qquad (19) $$

Eq. (19) shows that $T_{\alpha\beta}^{kl} \cdot \tau$ is the $\Psi_{\alpha}^k$-weighted average of the (opposite of the) electric field induced over inclusion $\alpha$ by the non-uniform polarization $\Psi_{\beta}^l \tau$ supported in inclusion $\beta$. As such, $T_{\alpha\beta}^{kl}$ will be called influence tensors.

For $\alpha = \beta$, these tensors measure the electric field induced on inclusion $\alpha$ by a polarization applied to the same inclusion $\alpha$. The $T_{\alpha\alpha}^{kl}$ will be called self-influence tensors, and denoted $S_{\alpha}^{kl}$.

The evaluation of these tensors $T_{\alpha\beta}^{kl}$ and $S_{\alpha}^{kl}$ will be discussed in the next Section.

Finally, Eq. (10) is applied to the solution to the above linear system 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,l} f_{\alpha} M_{\alpha}^{k} \tau_{\alpha}^{k}. \qquad (20) $$
## 2.3 Integral expression of the influence and self-influence tensors

### 2.3.1 Influence tensors

Plugging Eq. (16) in Eq. (19), and using Eqs. (7) and (8) (which give the expression of the Green operator for isotropic reference material), we obtain :

$$ T_{12}^{kl} = \sigma_0^{-1} (U_{12}^{kl} \mathbf{1} - \mathbf{V}_{12}^{kl}), \qquad (21) $$

with

$$ U_{12}^{kl} = \frac{1}{4\pi V_1} \int_{\mathbf{y}_1 \in \Omega_1} \int_{\mathbf{y}_2 \in \Omega_2} \frac{z_1^k z_2^l}{\|\mathbf{y}_2 - \mathbf{y}_1\|^3} d^3\mathbf{y}_2 d^3\mathbf{y}_1 \qquad (22) $$

and

$$ \mathbf{V}_{12}^{kl} = \frac{3}{4\pi V_1} \int_{\mathbf{y}_1 \in \Omega_1} \int_{\mathbf{y}_2 \in \Omega_2} z_1^k z_2^l \frac{(\mathbf{y}_2 - \mathbf{y}_1) \otimes (\mathbf{y}_2 - \mathbf{y}_1)}{\|\mathbf{y}_2 - \mathbf{y}_1\|^5} d^3\mathbf{y}_2 d^3\mathbf{y}_1, \qquad (23) $$

with

$$ z_{\alpha} = (\mathbf{y}_{\alpha} - \mathbf{x}_{\alpha}) \cdot \mathbf{n}_{\alpha} \quad (\alpha = 1, 2). \qquad (24) $$

We then introduce the local cylindrical coordinates $(r_\alpha, \theta_\alpha, z_\alpha)$ ($\alpha = 1, 2$, see Fig. 2)

$$ \mathbf{y}_{\alpha} = \mathbf{x}_{\alpha} + z_{\alpha} \mathbf{n}_{\alpha} + r_{\alpha} \mathbf{e}_{r,\alpha}, \qquad (25) $$

where $\mathbf{e}_{r,\alpha}$ is the unit radial vector, while $\theta_\alpha$ is the polar angle with respect to a fixed (unspecified) direction. We also introduce the following dimensionless variables

$$ \zeta_1 = \frac{z_1}{L}, \quad \zeta_2 = \frac{z_2}{L}, \quad \rho_1 = \frac{r_1}{R}, \quad \text{and} \quad \rho_2 = \frac{r_2}{R}. \qquad (26) $$

Then

$$ \mathbf{y}_2 - \mathbf{y}_1 = L\zeta_2\mathbf{n}_2 - L\zeta_1\mathbf{n}_1 + \mathbf{x}_2 - \mathbf{x}_1 + R(\rho_2\mathbf{e}_{r,2} - \rho_1\mathbf{e}_{r,1}). \qquad (27) $$

We were not able to derive a closed-form expression of these 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 of, this assumption is certainly verified for most pairs of fibers. It is then assumed that $R \|\rho_2 \mathbf{e}_{r,2} - \rho_1 \mathbf{e}_{r,1}\| \ll L\zeta_2 \mathbf{n}_2 - L\zeta_1 \mathbf{n}_1 + \mathbf{x}_2 - \mathbf{x}_1$ in the above expression. Then $U_{12}^{kl}$ and $\mathbf{V}_{12}^{kl}$ read

$$ U_{12}^{kl} = \frac{L^{k+l+3}}{8\pi^2} \frac{R^2}{L^2} \int_{\substack{0 \le \rho_\alpha \le 1 \\ 0 \le \theta_\alpha \le 2\pi \\ -1 \le \zeta_\alpha \le 1}} \frac{\rho_1 \rho_2 \zeta_1^k \zeta_2^l}{\|\mathbf{w}_0\|^3} d\zeta_1 d\zeta_2 d\theta_1 d\theta_2 d\rho_1 d\rho_2, \qquad (28) $$

$$ \mathbf{V}_{12}^{kl} = \frac{3L^{k+l+3}}{8\pi^2} \frac{R^2}{L^2} \int_{\substack{0 \le \rho_\alpha \le 1 \\ 0 \le \theta_\alpha \le 2\pi \\ -1 \le \zeta_\alpha \le 1}} \frac{\rho_1 \rho_2 \zeta_1^k \zeta_2^l \mathbf{w}_0 \otimes \mathbf{w}_0}{\|\mathbf{w}_0\|^5} d\zeta_1 d\zeta_2 d\theta_1 d\theta_2 d\rho_1 d\rho_2, $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2 – Computation of the influence tensors</figcaption></figure>

where

$$ \mathbf{w}_0 = \mathbf{r}_{12} + L\zeta_2\mathbf{n}_2 - L\zeta_1\mathbf{n}_1. \qquad (29) $$

Upon integration with respect to $0 \le \rho_1, \rho_2 \le 1$ and $0 \le \theta_1, \theta_2 \le 2\pi$

$$ U_{12}^{kl} = \frac{L^{k+l+3}}{8} \frac{R^2}{L^2} \int_{-1}^{1} \int_{-1}^{1} \frac{\zeta_1^k \zeta_2^l}{\|\mathbf{w}_0\|^3} d\zeta_1 d\zeta_2, \qquad (30) $$

$$ V_{12}^{kl} = \frac{3L^{k+l+3}}{8} \frac{R^2}{L^2} \int_{-1}^{1} \int_{-1}^{1} \frac{\zeta_1^k \zeta_2^l \mathbf{w}_0 \otimes \mathbf{w}_0}{\|\mathbf{w}_0\|^5} d\zeta_1 d\zeta_2. \qquad (31) $$

In the above expressions, the first integration (with respect to $\zeta_1$) can be performed analytically, the second integral (with respect to $\zeta_2$) is evaluated numerically.

### 2.3.2 Self-influence tensors

Using Eqs (7) and (19), the self-influence tensor is given by

$$ S_{\alpha}^{kl} = M_{\alpha}^{k+l} P_0 + \frac{1}{V_{\alpha}} \int_{x \in \Omega_{\alpha}} \Psi_{\alpha}^{k}(\mathbf{x}) \mathrm{PV} \int_{y \in \Omega_{\beta}} \Psi_{\beta}^{l}(\mathbf{y}) G_0(\mathbf{y} - \mathbf{x}) \mathrm{d}^3\mathbf{y} \mathrm{d}^3\mathbf{x}. \qquad (32) $$

Eq. (32) shows that the singularity cannot be removed in the integral expression of the self-influence tensors, which makes their analytical evaluation difficult. Rather than attempting a numerical evaluation of these integrals, and observing that this tensor can be precomputed off-line prior to the full EIM calculation, we used a numerical approach based on a finite element analysis.

The resulting problem is formulated over the axisymmetric domain $\Omega_1$, although the loading is not axisymmetric. Still, it is possible to reduce the 3d problem to a finite set of 2d problems [13, 10]. These 2d variational problems were implemented within the FEniCS framework [6, 1].
<figure><img src="image_5.png" /><figcaption>FIGURE 3 – Eshelby’s problem with contrast $10^6$ for aspect ratio 50</figcaption></figure>

# 3 Applications with a few fibres

## 3.1 Problem with one fibre

The equivalent inclusion method introduced above is first applied to Eshelby's inhomogeneity problem [5] for a circular, cylindrical inhomogeneity. More precisely, we consider a single inhomogeneity $\Omega_1$ centered at the origin and embedded in a homogeneous, infinite matrix. The system is subjected to a uniform electric field $\bar{E} = \bar{E}_z e_z$ at infinity, where $\bar{E}$ is a scalar constant and $e_z$ denotes the axis of revolution of the inhomogeneity. Both matrix and inhomogeneity have isotropic conductivities $\sigma_0$ and $\sigma_1$, respectively. The contrast of conductivities is $\sigma_1/\sigma_0 = 10^6$, while the aspect ratio is $L/R = 50$.

We compared the results of EIM with finite element calculation. To present this finite element results, we compute the average field on the cylinder cross section, along the fibre longitudinal coordinate $z$. The EIM field, which is uniform in fibre section, is shown for every $z$. We can see on figure 3 a good agreement between finite element calculation and EIM method, for the longitudinal component of the electric field. Moreover, we can see that even for quadratic polynomials, the EIM estimation is close to the finite element solution.

## 3.2 Two parallel cylinders

Here we consider two parallel cylinders submitted to an external field $\bar{E} = e_z$ at infinity, where $e_z$ corresponds to the orientation of the cylinders. The distance between cylinders is $2R$. We work with an aspect ratio of $e = 50$ and a contrast of $10^6$. We can see on figure 4 that even for cylinders near from each other, their interaction is well captured and there is no need to search for high polynomial orders.
<figure><img src="image_6.png" /><figcaption>FIGURE 4 – 2 parallel cylinders with aspect ratio 50 and contrast $10^6$</figcaption></figure>

## 4 Conclusion

The method presented in this paper is an adaptation of the variational form of the equivalent inclusion method [8, 2], in the case of fibres. It is based on the same principles, but it takes advantage of the fibre aspect ratio. The tensors involved in the linear system to be assembled are pre-calculated via a combination of 2D finite element calculations, analytical calculations and numerical integrations. The examples presented in this paper with few cylinders show a good agreement between this approximate method and more accurate finite element calculations, even for a small number of degrees of freedom.

This work will now be extended to large assemblies of conducting fibres, where full field approaches become too expensive.

## Références

[1] Alnæs, Martin and Blechta, Jan and Hake, Johan and Johansson, August and Kehlet, Benjamin and Logg, Anders and Richardson, Chris and Ring, Johannes and Rognes, Marie E. and Wells, Garth N., The FEniCS Project Version 1.5, Archive of Numerical Software, 3 (2015) 2197-8263

[2] S. Brisard, L. Dormieux, K. Sab, A variational form of the equivalent inclusion method for numerical homogenization, International Journal of Solids and Structures, 51 (2014) 716-728.

[3] Brisard, S. and Sab, K. and Dormieux, L., New boundary conditions for the computation of the apparent stiffness of statistical volume elements, Journal of the Mechanics and Physics of Solids, 12, vol. 61 (2013)

[4] Solid Mixture Permittivities, Brown, William Fuller, The Journal of Chemical Physics, vol. 23 number 8 (1955), 1514-1517

[5] Eshelby, J. D., The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proceedings of the Royal Society of London, 241 (1957)
[6] Logg, Anders and Mardal, Kent-Andre and Wells, Garth N., Automated Solution of Differential Equations by the Finite Element Method, The FEniCS Book, Lecture Notes in Computational Science and Engineering, 84, Springer (2012)

[7] Mardiguian, M. and Caron-Fellens, J. P. The intelligent concrete : A new, economical technique for architectural shielding of buildings IEEE Electromagnetic Compatibility Magazine vol. 6, no 2 (2017)

[8] Moschovidis, Z. A. and Mura, T., Two-ellipsoidal inhomogeneities by the equivalent inclusion method, Journal of Applied Mechanics, 42 (1975)

[9] Ostoja-Starzewski, Martin, Material Spatial Randomness : From Statistical to Representative Volume Element, Probabilistic Engineering Mechanics, 21 (2006) 2 112-132

[10] Percy, John H. and Pian, Theodore H. H. and Klein, Stanley and Navaratna, Dhirendra R., Application of Matrix Displacement Method to Linear Elastic Analysis of Shells of Revolution, AIAA Journal, 3 (1965) 11 2138-2145

[11] Torquato, S., Effective Electrical Conductivity of Two-phase Disordered Composite Media, Journal of Applied Physics, vol. 58, number 10 (1985), 3790-3797

[12] Willis, J.R., Bounds and Self-Consistent Estimates for the Overall Properties of Anisotropic Composites, Journal of the Mechanics and Physics of Solids, volume 25 (1977) number 3 185-202

[13] Wilson, Edward L., Structural Analysis of Axisymmetric Solids, AIAA Journal, 3 (1965) 12 2269-2274

[14] Zeller, R. and Dederichs, P. H., Elastic Constants of Polycrystals, Physica Status Solidi (B), 55 (1973) 2 831-842