# About the elastic energy related to an inclusion or an inhomogeneity in large or infinite domains and mixed boundary conditions

Christophe Fond

► To cite this version:

Christophe Fond. About the elastic energy related to an inclusion or an inhomogeneity in large or infinite domains and mixed boundary conditions. 2023. hal-03979614v2

HAL Id: hal-03979614

https://hal.science/hal-03979614v2

Preprint submitted on 15 Jun 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
About the elastic energy related to an inclusion or an
inhomogeneity in large or infinite domains and mixed
boundary conditions

Christophe Fond

Laboratoire ICube, Institut Universitaire de Technologie Robert Schuman - département
génie civil, 72 route du Rhin B.P. 315, F67411 Illkirch cedex

Abstract

The equivalent inclusion method allows to solve analytically and exactly problems for one inhomogeneity in an elastic domain in generally infinite spaces. It also allows to solve semi-analytically and approximately problems for several inhomogeneities in elastic domains. Knowing the additional elastic energy due to the presence of one or more inclusions / inhomogeneities is essential to estimate the elastic properties of the homogenized medium. The usual demonstrations proposed to calculate this energy supplement often make the assumption of stress free boundary condition. These demonstrations given in the literature are rigorous for solutions obtained under these conditions in finite domains. But analytical solutions in finite domains are extremely rare. Unfortunately, the assumption of stress free boundary condition is no longer valid when one uses analytical solutions established in infinite domains. For the latter solutions, which are mostly used, the stresses do actually tend towards zero when one moves away from the inclusion but it is shown here that the rate of decrease, compared to the rate of growth of the boundary surface, is not sufficient for one to consider them as zero. This is particularly apparent when calculating the work of external forces. Indeed, a correction - even a very small one - is essential to quantify precisely the strain energy when implementing numerically these calculation methods. I provide here a rigorous demonstration that does not make the usual stress-free boundary assumption by considering mixed boundary conditions when isolating a large finite domain in an infinite domain. This demonstration legitimizes the small correction to apply on the prescribed boundary conditions, enabling us to recover the correct result for calculating the work of
external forces, i. e. the elastic strain energy in the herein considered cases.
This allows credit for suitable numerical estimates of the strain energy in fi-
nite and infinite domains. Numerical illustrations are added to provide orders
of magnitude and finish convincing the readers.

*Keywords:* micromechanics, inhomogeneities, heterogeneities, inclusions, elasticity, equivalent inclusion method, homogenization

# 1. Introduction

The aim is to solve linear elasticity problems for heterogeneous media although plastic deformations can also be considered (Doghri et al., 2011). The estimation of the average elastic properties (Benveniste, 1987), (Brisard et al., 2014) is not the only objective since we are also interested in the local damage phenomena (Berthaud et al., 1994), (Fond et al., 1995), (Fond and Berthaud, 1995), (Fond et al., 1996), (Fond, 2001), (Shodja et al., 2003), (Benedikt et al., 2006), (Xu et al., 2023). We thus seek to know the localization tensors in the presence of mechanical interaction. In practice for our numerical applications, the elasticity is isotropic and the domains containing the inclusions or inhomogeneities are spherical. The calculation of elastic the strain energy is not immediate for the infinite domains considered. In particular, the remark often made in the literature concerning stresses that cancel at very distant boundaries<fnref n="1" /> leads to the assumption of a null work of these stresses in the boundary displacements. A closer look at this point shows that this is not generally the case for known analytic solutions since boundary conditions are mixed (Brisard et al., 2013). It is then necessary to reconsider this argument in the demonstration frequently provided in the literature - inexact when using analytical solutions calculated for an infinite domain - to analyze in which conditions one retrieves that this small energy supplement is indeed $\Delta W_p = \pm \frac{1}{2} V_p \underline{\underline{\sigma}}^\infty \underline{\underline{\beta}}^p$, ± according to displacement or stress boundary condition, where $\underline{\underline{\sigma}}^\infty$ denotes the boundary condition at infinity, $\underline{\underline{\beta}}^p$ the eigenstrain (Eshelby, 1957), (Eshelby, 1959), (Mura, 1993) and $V_p$ the volume of the inclusion or inhomogeneity. For the latter solutions, which are the majority of those used, these stresses do actually tend towards zero as one moves away from the inclusion. But their rate of decrease,

<footnote n="1">¹I. e. at infinity.</footnote>
compared to the rate of growth of the boundary surface, is not sufficient to
considered them equal to zero. Indeed, the calculation of the work of external
forces shows the product of a term that tends to zero by a term that tends
to infinity with inverse growth and decay rates so that this product tends
asymptotically to a constant and non-zero value which does not correspond
to the value of $\Delta W_p$ given above. The numerical implementation of this
computational method clearly shows that a numerical correction based on
the values at the far boundaries is essential to properly quantify the strain
energy. Even though this numerical correction tends to zero<fnref n="2" /> for larger and
larger domains, it remains numerically indispensable and forms the basis of
the rigorous demonstration provided in this paper.

## 2. Foreword

We will place ourselves under the small perturbations assumptions (S. P. A.), i. e. linear elasticity at small strains and small displacements<fnref n="3" />. The considered elasticity is linear so that it is possible to superimpose kinematically and statically admissible fields. Moreover, we find ourselves in a system where deformations, stresses and displacements are proportional to the solicitation so that:

- even if we assume, for the sake of simplicity, a strain close to the unit, we do not leave the framework of S. P. A. and it is enough to apply to the results a coefficient small in front of unity to obtain a realistic result corresponding to the small strains and displacements assumptions,

- since one deals with linear systems, a single calculation is sufficient, it is then enough to apply a coefficient of proportionality to adjust the result as one wants.

## 3. Illustration in a simple case

### 3.1. A sphere under pressure

Let's consider a sphere of external radius $R_{ext}$ made of a material which, for the sake of simplicity, Young's modulus $E$ is equal to 1 and Poisson's ratios $\nu$ equal to 0.3 so that its bulk modulus $K$ is equal to $\frac{1}{1.2}$. This sphere is submitted to a uniform equitriaxial tension equal to 1, i. e. the external

<footnote n="2">²As well as $\Delta W_p$ which is of the same order of magnitude.</footnote>

<footnote n="3">³One can assimilate the initial geometry to the deformed geometry.</footnote>
pressure $P_{ext}$ is equal to -1. The strain energy $W_0$ over the total volume $V_0$ of this sphere is then equal to $\frac{V_0 P_{ext}^2}{2K} = \frac{2\pi R_{ext}^3 P_{ext}^2}{3K}$.

## 3.2. A hollow sphere under pressure

In very rare cases, the analytical solution is known for a domain of finite size containing an inhomogeneity such as a hollow sphere under pressure. In the latter case, the problem is somehow uniaxial. The analytical solution is given in Appendix A. On the other hand, the analytical solution for a spherical cavity in a infinite domain under pressure is available (Eshelby, 1957), (Eshelby, 1959). In this case the eigenstrains $\beta_{ii}$ are non-zero and $\beta_{ij}, i \neq j$, are equal to zero.

## 3.3. Numerical results

Let's now consider, for the numerical values given in 3.1, a sphere of relatively large size compared to the size of an internal and centered spherical cavity, for instance typically $R_{ext} = 100R_{int}$. For the sake of simplicity, let's choose $R_{int} = 1m$ so that $R_{ext} = 100m$. The strain energy $W$ calculated from eq. A.3 is $W = 2.513280720 \times 10^6 J$ which leads to $W - W_0 = 6.597J$ since $W_0 = 2.513274123 \times 10^6 J$ in this case. Note that this corresponds to exact boundary condition prescribed in stress. The displacement at $R_{ext} = 100m$ is equal to $U_r(R_{ext}) \approx 40.000105m$. Imposing $P_{ext} = -1 \times \frac{40.0}{40.000105}$ in order to obtain exact prescribed boundary conditions in displacements, i. e. $U_r(R_{ext}) = 40.0m$ leads to $W - W_0 = -6.597J$.

For a domain of infinite size, with same material characteristics, containing a spherical cavity, the eigenstrains $\beta_{ii}$ are equal to 1.05 for $P_{ext} = -1Pa$ (Eshelby, 1957). This leads to $\Delta W_p = 6.597J$ for a spherical cavity of radius $R_{int} = 1m$. This shows that considering that $R_{ext} = 100R_{int}$ almost simulates an infinite domain and that $\Delta W_p = \pm \frac{1}{2}V_p \underline{\underline{\sigma}}_{\infty}^p \beta^p$ is valid for exact boundary conditions. This also shows, as obviously expected, that $\frac{\Delta W_p}{W_0} \ll 1$.

For homogenization problem, the strains or stresses are defined at infinity. In our case this would lead to $\epsilon_{ii}^\infty = 0.4$ and $\sigma_{ii}^\infty = 1^4$. It is commonly assumed that for infinite domains containing a finite number of inhomogeneities that the mean value of strain $\langle \epsilon \rangle$ is equal to $\underline{\underline{\epsilon}}^\infty$ and that the mean value of stress $\langle \underline{\underline{\sigma}} \rangle$ is equal to $\underline{\underline{\sigma}}^\infty$. In the case of Eshelby's solution for which $\beta_{ii} = 1.05$, at $R_{ext} = 100m$, $U_r(100m) = 40.00065m$ and $\sigma_{rr}(100m) = 0.9999990$.

<footnote n="4">⁴If $i \neq j$, $\epsilon_{ij}^\infty = 0$ and $\sigma_{ij}^\infty = 0$.</footnote>
Hence, calculating the strain energy from the work of external forces<fnref n="5" /> gives $W = 2.513275694 \times 10^6 J$ such that $W - W_0 \approx 1.571J$ at $R = 100m$. Even if they tend to, whatever the ratio $\frac{R}{R_{int}}$ is, neither the mean value of strain nor the mean value of stress corresponds to the exact desired mean state of strain or stress. Even for larger ratios of $\frac{R}{R_{int}}$, one always observe that $\langle \epsilon_{ii} \rangle > 0.4$ and $\langle \sigma_{ii} \rangle < 1$ and that the work of external forces reaches rapidly an asymptotic value equal to $W - W_0 \approx 1.571J$, as shown in Table 1. This value does not corresponds to $\Delta W_p$ even for infinitely large domains. Actually, if the "external forces"<fnref n="6" /> tend to zero for larger and larger domains, their work in the displacement converge to a non-zero finite value. The further use of variational principles instead of using the classical Equivalent Inclusion Method (E. I. M.), (Eshelby, 1957), (Eshelby, 1959), to numerically solve inhomogeneities problems, (Moschovidis and Mura, 1975), (Fond et al., 2002), must then manage properly the work of external forces. Indeed, the infinitesimal supplementary stress at infinity/very far from the inhomogeneity can never been considered as exactly null.

<table><caption>Table 1: Supplementary energy corresponding to the parameters given in 3.1 for an increasing radius $R$ of the integration contour, which surface is $4\pi R^2$, in an infinite domain in the case of the E. I. M. calculation. $U_{r0} = \epsilon_{ii}^\infty \times R$ and $\sigma_{rr0} = \sigma_{ii}^\infty$ denote the situation without any cavity.</caption><thead><tr><th>R<sub>int</sub>(m)</th><th>R(m)</th><th>(U<sub>r</sub>(R)-U<sub>r0</sub>(R))/U<sub>r0</sub>(R)(-)</th><th>(σ<sub>rr</sub>(R)-σ<sub>rr0</sub>(R))/σ<sub>rr</sub>(R)(-)</th><th>W - W<sub>0</sub>(J)</th></tr></thead><tbody><tr><td>1</td><td>10</td><td>+1.625 × 10<sup>-3</sup></td><td>-1.0 × 10<sup>-3</sup></td><td>1.56671</td></tr><tr><td>1</td><td>100</td><td>+1.625 × 10<sup>-6</sup></td><td>-1.0 × 10<sup>-6</sup></td><td>1.57079</td></tr><tr><td>1</td><td>1000</td><td>+1.625 × 10<sup>-9</sup></td><td>-1.0 × 10<sup>-9</sup></td><td>1.57081</td></tr></tbody></table>

# 4. Equivalent inclusion method (E. I. M.)

## 4.1. Notations

The E. I. M., see Appendix B, was initially proposed by (Eshelby, 1957). For several inclusions/inhomogeneities (Moschovidis and Mura, 1975), (Yu

<footnote n="5">⁵By isolating a spherical domain centered on the cavity in the infinite domain, the external forces being the action of the external matter on the imaginary surface of this sphere.</footnote>

<footnote n="6">⁶On the imaginary contour of integration.</footnote>
et al., 2021), (Yan et al., 2021), it uses as a functional basis eigenstrains<fnref n="7" /> generally distributed according to Taylor series in inclusions/inhomogeneities<fnref n="8" /> and is situated between a Galerkin method and a collocation method (Brisard et al., 2014). The additional elastic energy associated with the presence of an inclusion, compared to the elastic strain energy of the same elastic medium which would not contain an inclusion, was calculated by this same author. This calculation has been taken up by many authors (Mura, 1993). Let us therefore consider an exact solution for an eigenstrain $\underline{\underline{\beta}}^p$ in a matrix domain designated by $\mathcal{D}_p$ of volume $V_p$, see Fig. 1. This matrix domain $\mathcal{D}_p$ having undergone this eigenstrain is called inclusion. Let us denote by $\bar{U}^p \underline{\underline{u}}^p \underline{\underline{\beta}}^p$ and $\underline{\underline{\epsilon}}^p \underline{\underline{\beta}}^p = \underline{\underline{A}}^p \underline{\underline{\beta}}^p$ respectively the displacement and strain fields associated with the eigenstrain $\underline{\underline{\beta}}^p$. $\underline{\underline{A}}^p$, 4th-order tensor, denotes the strain localization tensor associated with the eigenstrain $\underline{\underline{\beta}}^p$. $\underline{\underline{U}}^p$, 3th-order tensor, denotes the displacement localization tensor associated with the eigenstrain $\underline{\underline{\beta}}^p$. The stress field inside the inclusion $p$ is given by $\underline{\underline{\sigma}}^p \underline{\underline{\beta}}^p = \underline{\underline{C}}^0 (\underline{\underline{\epsilon}}^p \underline{\underline{\beta}}^p - \underline{\underline{\beta}}^p)$ and outside the inclusion by $\underline{\underline{\sigma}}^p \underline{\underline{\beta}}^p = \underline{\underline{C}}^0 \underline{\underline{\epsilon}}^p \underline{\underline{\beta}}^p$ since $\underline{\underline{\beta}}^p = \underline{\underline{0}}$ outside the inclusion.

<footnote n="7">⁷I. e. a stress free strain.</footnote>

<footnote n="8">⁸Ellipsoidal in order to have analytical solutions.</footnote>

<footnote n="9">⁹Note that $\underline{\underline{A}}^p = \frac{1}{2}(\underline{\underline{grad}}[\underline{\underline{U}}^p] + \underline{\underline{grad}}^T[\underline{\underline{U}}^p])$.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: A *p* domain of a homogeneous matrix undergoes an eigenstrain **β**^p. An imaginary external action τ^*p. →n on its contour would restore it to the shape it had before undergoing this eigenstrain.</figcaption></figure>

4.2. Evolution of the fields away from inclusion

The exact solutions $\mathcal{A}^p$, in infinite domain, or semi-infinite spaces (Yang
et al., 2021), are known for ellipsoidal inclusions subjected to eigenstrains
generally distributed according to Taylor series (Ferrers, 1877), (Dyson et al.,
1891), (Eshelby, 1957), (Moschovidis and Mura, 1975), Fond et al. (2001),
Fond et al. (2002). For various shapes, it is necessary to use numerical tools
(Nakasone et al., 2000) to have a solution to this type of problem. In prac-
tice the domains are very large compared to the size of the inclusions but
not infinite. Most of the demonstrations concerning the calculation of the
energy supplement associated with inclusions/inhomogeneities assume that
the stress field associated with these solutions is self-equilibrated, i.e. that
there is no external force associated with this field. There would thus be
no work of external forces at external boundaries. But the solutions given
for infinite domains correspond to mixed boundary conditions for which nei-
ther the displacements nor the stresses ever exactly cancel out as one moves
away from the inclusion, although they tend to zero. Note therefore that,
for solution of eigenstrains in an infinite domain, displacements vary in $\frac{1}{L^2}$
and stresses and strains vary in $\frac{1}{L^3}$, where L denotes the distance to the
center of the inclusion. Under the S. P. A., it is relevant to superimpose on these fields uniform fields of strain or stress to, among other things, estimate the homogenized elastic properties<fnref n="10" />. Denoting by $\underline{\epsilon}^\infty$ any finite uniform imposed strain, the displacements at the boundaries $\overline{U^L}$ are proportional to the size $L$ of the considered domain. The surface integral, the latter being of size proportional to $L^2$, tends to a finite value $\Delta W_p$ when considering a medium of large size $L$ in front of the size $a$ of the inclusion, i. e. $\iint_{S_{\text{ext}}} [\underline{\sigma}^{A^p \beta^p} \cdot \overline{n}] \cdot \overline{U^L} dS \propto \frac{1}{L^3} L \times L^2 \neq 0$.

### 4.3. Satisfaction of boundary conditions on average

This explains that by superimposing the analytical solution - which we will call perturbation field, see Figs. C.3 and C.4 - to a uniform strain $\underline{\epsilon}^\infty$, one does not directly obtain the desired conditions nor those that cancel the work of the perturbation fields at the distant boundaries. In order to recover the desired boundary conditions from the superposition, it is necessary to subtract from $\underline{\epsilon}^\infty$ the average of the strains associated with the perturbation field, which we will denote by $\underline{\epsilon}^{<>} = \langle \underline{A^p} \cdot \underline{\beta^p} \rangle_D$, $\langle \rangle_D$ designating the average over the domain $D$. For a cubic domain, the result will not provide a flat deformed surface but an almost flat one, as shown in Fig. C.3, and all the more flat as the cube will be large compared to the size of the inclusion. Let's resume the calculations of Appendix B by imposing on the second member of the eq. B.2 $(\underline{\epsilon}^\infty - \underline{\epsilon}^{<>})$ so that we solve the equation $([\underline{C}^0 - \underline{C}^p]^{-1} \cdot [\underline{C}^0] - \underline{A^p}) \cdot \underline{\beta^p} = (\underline{\epsilon}^\infty - \underline{\epsilon}^{<>})$. It should be noted that $\underline{\epsilon}^{<>} \to 0$ with $L \to \infty$.

## 5. Demonstration

If one imposes $(\underline{\epsilon}^\infty - \underline{\epsilon}^{<>})$, the strain energy $W$ is:

$$W = \frac{1}{2} \iiint_D (\underline{\epsilon}^\infty - \underline{\epsilon}^{<>} + \underline{\epsilon}^{A^p \beta^p}) : \underline{C}^0 \cdot (\underline{\epsilon}^\infty - \underline{\epsilon}^{<>} + \underline{\epsilon}^{A^p \beta^p} - \underline{\beta}^p) \, dv \quad (1)$$

where $\underline{\epsilon}^{A^p \beta^p} = \underline{A^p} \cdot \underline{\beta^p}$. By developing we obtain:

<footnote n="10">¹⁰Since physically eigenstrains do not modify the elasticity of the inclusions.</footnote>
$$
\begin{equation}
\begin{aligned}
W ={}& \frac{1}{2} \iiint_D \underline{\epsilon}^{\infty} : \underline{C}^0 \underline{\epsilon}^{\infty} dv - \frac{1}{2} \iiint_D \underline{\epsilon}^{\infty} : \underline{C}^0 \underline{\beta}^p dv \\
& + \frac{1}{2} \iiint_D \underline{\epsilon}^{\infty} : \underline{C}^0 \underline{\epsilon}^{A^p \beta^p} dv - \frac{1}{2} \iiint_D \underline{\epsilon}^{\infty} : \underline{C}^0 \underline{\epsilon}^{<>} dv \\
& - \frac{1}{2} \iiint_D \underline{\epsilon}^{<>} : \underline{C}^0 \underline{\epsilon}^{\infty} dv + \frac{1}{2} \iiint_D \underline{\epsilon}^{<>} : \underline{C}^0 \underline{\epsilon}^{<>} dv \\
& - \frac{1}{2} \iiint_D \underline{\epsilon}^{<>} : \underline{C}^0 \underline{\epsilon}^{A^p \beta^p} dv + \frac{1}{2} \iiint_D \underline{\epsilon}^{<>} : \underline{C}^0 \underline{\beta}^p dv \\
& + \frac{1}{2} \iiint_D \underline{\epsilon}^{A^p \beta^p} : \underline{C}^0 \underline{\epsilon}^{\infty} dv - \frac{1}{2} \iiint_D \underline{\epsilon}^{A^p \beta^p} : \underline{C}^0 \underline{\epsilon}^{<>} dv \\
& + \frac{1}{2} \iiint_D \underline{\epsilon}^{A^p \beta^p} : \underline{C}^0 \underline{\epsilon}^{A^p \beta^p} dv - \frac{1}{2} \iiint_D \underline{\epsilon}^{A^p \beta^p} : \underline{C}^0 \underline{\beta}^p dv
\end{aligned}
\tag{2}
\end{equation}
$$

The first line gives $W_0 - \Delta W_p$ since the eigenstrain is zero outside the inclusion, i. e. $\iiint_D \underline{\beta}^p dv = \iiint_{D_p} \underline{\beta}^p dv$. Then we notice that the uniform fields can be taken out of the integral so that:

$$ \iiint_D \underline{\epsilon}^{<>} : \underline{C}^0 \underline{\epsilon}^{A^p \beta^p} dv = \underline{\epsilon}^{<>} : \underline{C}^0 \underline{\epsilon}^{A^p \beta^p} dv = V_0 \underline{\epsilon}^{<>} : \underline{C}^0 \underline{\epsilon}^{<>} dv \quad (3) $$

since $\frac{1}{V_0} \iiint_D \underline{\epsilon}^{A^p \beta^p} dv = \underline{\epsilon}^{<>}$. In the same way:

$$ \iiint_D \underline{\epsilon}^{\infty} : \underline{C}^0 \underline{\epsilon}^{A^p \beta^p} dv = V_0 \underline{\epsilon}^{\infty} : \underline{C}^0 \underline{\epsilon}^{<>} dv \quad (4) $$

Due to tensor symmetries<fnref n="11" />, eq. 1 becomes:

$$
\begin{align}
W ={}& W_0 - \Delta W_p \nonumber \\
& + \frac{1}{2} \iiint_D \underline{\epsilon}^{<>} : \underline{C}^0 \underline{\beta}^p dv - \frac{1}{2} \iiint_D \underline{\epsilon}^{A^p \beta^p} : \underline{C}^0 \underline{\beta}^p dv \tag{5} \\
& + \frac{1}{2} \iiint_D \underline{\epsilon}^{A^p \beta^p} : \underline{C}^0 \underline{\epsilon}^{A^p \beta^p} dv - \frac{1}{2} \iiint_D \underline{\epsilon}^{A^p \beta^p} : \underline{C}^0 \underline{\epsilon}^{<>} dv \nonumber
\end{align}
$$

<footnote n="11">¹¹ $\underline{\epsilon}^1 : \underline{C} \underline{\epsilon}^2 = \underline{\epsilon}^2 : \underline{C} \underline{\epsilon}^1$</footnote>
Recalling that $\underline{\underline{\tau}}^{*p} = -\underline{\underline{C}}^0 \underline{\underline{\beta}}^p$, since the work of external forces is equal to the work of internal efforts, it comes:

$$ - \iiint_D \underline{\underline{\epsilon}}^{A^p \beta^p} : \underline{\underline{C}}^0 \underline{\underline{\beta}}^p dv = \iiint_{\mathcal{D}_p} \underline{\underline{\epsilon}}^{A^p \beta^p} : \underline{\underline{\tau}}^{*p} dv = \iint_{\partial \mathcal{D}_p} \overline{U^{C^p \beta^p}} \cdot (\underline{\underline{\tau}}^{*p} \cdot \overline{\underline{n}}) ds \quad (6) $$

where $\partial \mathcal{D}_p$ denote the surface of the domain $\mathcal{D}_p$. Also by virtue of the equality of the work of internal and external efforts, see Fig. 1, it comes:

$$ \begin{aligned} & \iiint_D \underline{\underline{\epsilon}}^{A^p \beta^p} : \underline{\underline{C}}^0 \underline{\underline{\epsilon}}^{A^p \beta^p} dv \\ &= \iint_{\partial \mathcal{D}_p} \overline{U^{C^p \beta^p}} \cdot (\underline{\underline{\tau}}^{*p} \cdot \overline{\underline{n}}) ds + \iint_{\partial \mathcal{D}} \overline{U^{C^p \beta^p}} \cdot (\underline{\underline{C}}^0 \underline{\underline{\epsilon}}^{A^p \beta^p} \cdot \overline{\underline{n}}) ds \quad (7) \\ &= \iint_{\partial \mathcal{D}_p} \overline{U^{C^p \beta^p}} \cdot (\underline{\underline{\tau}}^{*p} \cdot \overline{\underline{n}}) ds + \iint_{\partial \mathcal{D}} \overline{U^{C^p \beta^p}} \cdot (\underline{\underline{\sigma}}^{A^p \beta^p} \cdot \overline{\underline{n}}) ds \end{aligned} $$

Since $\underline{\underline{\beta}}^p = \underline{\underline{0}}$ outside of $\mathcal{D}_p$:

$$ \iiint_D \underline{\underline{\epsilon}}^{<>} : \underline{\underline{C}}^0 \underline{\underline{\beta}}^p dv = \underline{\underline{\epsilon}}^{<>} : \iiint_{\mathcal{D}_p} \underline{\underline{C}}^0 \underline{\underline{\beta}}^p dv \quad (8) $$

then finally:

$$ \iiint_D \underline{\underline{\epsilon}}^{A^p \beta^p} : \underline{\underline{C}}^0 \underline{\underline{\epsilon}}^{<>} dv = V_0 \underline{\underline{\epsilon}}^{<>} : \underline{\underline{C}}^0 \underline{\underline{\epsilon}}^{<>} \quad (9) $$

This leads to:

$$ \begin{aligned} & W = W_0 - \Delta W_p \\ & + \frac{1}{2} \underline{\underline{\epsilon}}^{<>} : \iiint_{\mathcal{D}_p} \underline{\underline{C}}^0 \underline{\underline{\beta}}^p dv \\ & - \frac{1}{2} V_0 \underline{\underline{\epsilon}}^{<>} : \underline{\underline{C}}^0 \underline{\underline{\epsilon}}^{<>} \\ & + \frac{1}{2} \iint_{\partial \mathcal{D}} \overline{U^{C^p \beta^p}} \cdot (\underline{\underline{\sigma}}^{A^p \beta^p} \cdot \overline{\underline{n}}) ds \end{aligned} \quad (10) $$
For a finite domain of large size L in front of the size of the inclusion, given that far from the inclusion $\overline{U^{C^p \beta^p}} \propto \frac{1}{L^2}$, $\underline{\underline{\epsilon}} \propto \frac{1}{L^3}$, $\underline{\underline{\sigma}}^{A^p \beta^p} \propto \frac{1}{L^3}$, $V_0 \propto L^3$, $\partial \mathcal{D} \propto L^2$ and that $\iint_{\mathcal{D}_p} \beta^p dv$ is a constant value, it can be said that:

$$ 
\begin{align} 
V_0 \underline{\underline{\epsilon}} & : \underline{\underline{C}}^0 \underline{\underline{\epsilon}} \rightarrow 0 \quad \text{with } L \rightarrow \infty \nonumber \\ 
\underline{\underline{\epsilon}} & : \iint_{\mathcal{D}_p} \underline{\underline{C}}^0 \underline{\underline{\beta}}^p dv \rightarrow 0 \quad \text{with } L \rightarrow \infty \tag{11} \\ 
\iint_{\partial \mathcal{D}} \overline{U^{C^p \beta^p}} \cdot (\underline{\underline{\sigma}}^{A^p \beta^p} \cdot \vec{n}) ds & \rightarrow 0 \quad \text{with } L \rightarrow \infty \nonumber 
\end{align} 
$$

From which we can finally deduce that:

$$ W \approx W_0 - \Delta W_p \quad \text{with } L \to \infty \tag{12} $$

taking the precaution for numerical computations concerning finite domains of large size in front of that of the inclusion, to consider that the superposed fields are $\underline{\underline{\epsilon}}^{A^p \beta^p}$ and $(\underline{\underline{\epsilon}}^{\infty} - \underline{\underline{\epsilon}}^{<>)}$ and not only $\underline{\underline{\epsilon}}^{\infty}$.

# 6. Numerical illustration

## 6.1. Classical calculation

Let us consider for simplicity a cube of dimension $2L \times 2L \times 2L$ centered on a spherical inclusion of radius $a$ as illustrated in Fig. 2. Let us consider any case of isotropic linear elasticity under the S. P. A. For example $E^0 = 1$ a. u.$^{12}$ and $\nu^0 = 0.3$ for the matrix and $E^p = 0.5 \times E^0$ and $\nu^p = 0.1$ for the inhomogeneity $p$. Let us calculate the work of the external forces for a uniaxial traction of axis $\vec{z}$ $\epsilon_{33}^{\infty} = 1$ centered on the cube$^{13}$. The eigenstrain corresponding to the solution of the equivalent inclusion method is $\beta_{11} = \beta_{22} = 0.0335$ and $\beta_{33} = 0.62687$, the other terms being zero. Strains and stresses in the equivalent inclusion are uniform and are respectively $\epsilon_{11}^{A^p \beta^p} = \epsilon_{22}^{A^p \beta^p} = 0.0107$,

<footnote>$^{12}$a. u. for arbitrary unit.</footnote>

<footnote>$^{13}$All the results are proportional to $\epsilon^{\infty}$ so that for a small strain it is enough to apply a low value coefficient small in front of the unit.</footnote>
$\epsilon_{33}^{A^p \beta^p} = 0.325$, $\sigma_{11}^{A^p \beta^p} = \sigma_{22}^{A^p \beta^p} = 0.0208$ a. u., $\sigma_{33}^{A^p \beta^p} = 0.4498$ a. u.. For the inhomogeneity we obtain $\sigma_{11}^{A^p \beta^p} = \sigma_{22}^{A^p \beta^p} = 0.024$ a. u., $\sigma_{33}^{A^p \beta^p} = 0.1675$ a. u., the strains being the same as for the equivalent inclusion. The strain energy is given in Table 2 as a function of $L$. It is computed from the analytical solution in infinite domain and by numerical integration of 500 × 500 points per face of the cube. We denote by $W_0 = \frac{1}{2}V_0\epsilon^\infty : \sigma^\infty = \frac{1}{2}V_0\epsilon^\infty : C^0\epsilon^\infty$, where $V_0$ denotes the total volume of the domain, the elastic strain energy of the same domain that would not contain an inclusion. Let us now cal-

<table><caption>Table 2: Deformation energy and energy supplement as a function of the dimension of the cube containing the inhomogeneity.</caption><thead><tr><th>L/a (-)</th><th>W (a. u.)</th><th>W - W<sub>0</sub> (a. u.)</th></tr></thead><tbody><tr><td>3</td><td>108.216328100416</td><td>0.216328100416</td></tr><tr><td>10</td><td>4000.24628599888</td><td>0.24628599888</td></tr><tr><td>100</td><td>4000000.24851483</td><td>0.24851483</td></tr><tr><td>1000</td><td>40000000000.25054</td><td>0.25054</td></tr></tbody></table>

culate the additional energy with the classical formula $\Delta W_p = \pm \frac{1}{2} V_p \sigma^\infty \beta^p$ (Eshelby, 1957). $V_p = \pm \frac{4}{3} \pi a^3$ and $\sigma_{33}^\infty = 1$, the other stresses being zero in uniaxial tension. We obtain $\Delta W_p = -1.3213$ a. u. for boundary conditions defined in displacements and $\Delta W_p = +1.3213$ a. u. for boundary conditions defined in forces. This does not correspond to the values in Table 2 which estimates the value at +0.25 a. u. instead of -1.32 a. u..
<figure><img src="image_3.png" /></figure>

49103 noeuds &262191 elements

<figure><img src="image_4.png" /><figcaption>Figure 2: Left: spatial discretization of the finite element model used for a cube of dimensions 20 times the radius of a spherical inclusion, centered on this inclusion. All symmetries are exploited so that only $1/8^{th}$ of the volume constitutes the geometrical model for a unixial traction. Right: schematic representation of a section of a cube of dimensions $L$ centered on a spherical inclusion of radius $a$ and oriented along the axes of the reference frame.</figcaption></figure>

## 6.2. Adjustment of boundary conditions

Given that $\epsilon_{33}^{<>} \approx 2 \cdot 10^{-4}$, let us impose $\epsilon_{33}^{\infty} = 0.9998$ instead of unity, the strain energy calculated as in the previous calculation gives $W = 3998.7 \text{ a. u.}$. The additional energy is $\Delta W_p = \pm \frac{1}{2} V_p (C^{\underline{0}} \cdot \epsilon^{\infty}) \underline{\underline{\underline{\beta}}}^p = -1.3 \text{ a. u.}$. In this case, we find the theoretical prediction for a finite domain of large size in front of that of the inclusion, although the constraints are not zero at the boundaries and they work in the imposed displacements. It is therefore appropriate to try to estimate the accuracy that can be expected in finite domains by using exact analytical solutions valid for infinite domains. The term $\frac{1}{2} \iiint_{D_p} \epsilon^{A\underline{p}\beta\underline{p}} : C^{\underline{0}} \cdot \underline{\underline{\underline{\beta}}}^p dv$ is not negligible and is worth $0.554 \text{ a. u.}$. Similarly, the term $\frac{1}{2} \iiint_{D_p} \epsilon^{<>} : \underline{\underline{\underline{\sigma}}}^{\infty} ds$ which is equal to $0.08 \text{ a. u.}$ and the term $\frac{1}{2} \iiint_{D_p} \epsilon^{\infty} : \underline{\underline{\underline{\sigma}}}^{<>} ds$ which is equal to $-0.14 \text{ a. u.}$ may mar the accuracy of results.

## 6.3. Finite element model

The finite element method allows easy comparisons of numerical results (Novák et al., 2012). In order to complete the numerical verifications, a finite element model with CAST3M software (CAST3M, 2022) is implemented for the same geometric dimensions and elastic characteristics of the matrix and
inhomogeneity. The spatial discretization is shown in Fig. 2. It uses 262191 4-nodes tetrahedra for 147309 degrees of freedom. The boundary conditions are imposed in displacement so that the faces of the cube remain flat. The strain energy W is equal to 3998.69 a. u.. The term $\frac{1}{2} \iiint_D \beta^p : C^0 \cdot \epsilon^\infty dv$ is 1.30904 a. u. and the term $\frac{1}{2} \iiint_D \epsilon^{A^p\beta^p} : C^0 \cdot \epsilon^{A^p\beta^p} dv$ is 0.5546 a. u.. These results are as expected and validate the previous theoretical calculations. Other finite element results are compared in Appendix C to the exact solution in infinite domains.

## 7. Case of elastic inhomogeneities and mechanical interactions between inclusions/inhomogeneities

The equivalence equation is given in appendix Appendix B. In the case of N elastic inhomogeneities, there is interaction between inhomogeneities. Indeed, for general cases each inhomogeneity "feels" the presence of a softer or stiffer part of the surrounding medium. The exact solutions used are established for a homogeneous elastic matrix<fnref n="14" />. The simple superposition of their fields can no longer provide an exact solution. Indeed, the strain field induced by the presence of an inhomogeneity outside this inhomogeneity generates stress discontinuities at the interfaces of other inhomogeneities.

$$ \Delta_{\partial D_q} \underline{\underline{\underline{\sigma^q}}} = (\underline{\underline{\underline{C^q}}} - \underline{\underline{\underline{C^0}}}) \cdot \underline{\underline{\underline{\epsilon^{A^p\beta^p}}}} \quad (13) $$

where $\partial D_q$ denotes the $q^{th}$-inhomogeneity/matrix interface, i. e. the boundary of the $D_q$ domain. These "parasitic" forces $t_q^*$ are given by :

$$ \vec{t}_q^* = \Delta_{\partial D_q} \underline{\underline{\underline{\sigma^q}}} \cdot \vec{n}_q \quad (14) $$

These "parasitic" forces can be seen as external forces acting on the contours of the inhomogeneities (Benedikt et al., 2006), (Brisard et al., 2014). It has been proposed in Fond et al. (2002) to evaluate the quality of the solution provided by the equivalent inclusion method proposed by (Moschovidis and

<footnote n="14">¹⁴Linked to the Green's functions.</footnote>
Mura, 1975) by making these forces work in the displacements of the inhomogeneities contours. This estimator provides an additional strain energy $W^*$ which we wish to be as small as possible in front of the total additional strain energy.

$$W^* = \sum_{q=1}^{N} = \sum_{q=1}^{N} W_q^* = \iint \partial D_q \overrightarrow{t_q^*} \overrightarrow{U} ds \quad (15)$$

where $\overrightarrow{U}$ denotes the total displacement field, i.e. $\overrightarrow{U} = \overrightarrow{U}^\infty + \sum_{p=1}^{N} \overrightarrow{U}^{A^p \beta^p}$ for $N$ inhomogeneities. The solution will be considered globally suitable if $W^* \ll |W - W_0|$ and locally suitable if $|W_q^*| \ll |\Delta W^q|^{15}$. On distant boundaries, i. e. on $\partial D$, there is no problem of static admissibility but simply imperfect boundary conditions ($\overline{(\sum_{p=1}^{N} \sigma^{A^p \beta^p})} \cdot \overline{n_{\partial D}}$ and $\overline{(\sum_{p=1}^{N} \overline{U}^{A^p \beta^p})}$) with respect to what one wishes to impose.

The solution provided for several interacting inhomogeneities is therefore, except in special cases, not statically admissible. However, one can consider the provided solution as exact and fully admissible with imperfect boundary conditions, the $\partial D_p$ contours of the inhomogeneities being considered as boundaries as well as the contours $\partial D$ of the $D$ domain which contains them.

## 8. Discussion

I have shown here that it is necessary to cancel the average stress supplement at the distant boundaries to recover the results of the classical demonstration which assumes them to be zero. An analogous demonstration can be made to try to satisfy on average the conditions in stress at the distant boundaries by trying to cancel the average displacement supplement. Indeed, for a cube of volume $V_0 = 80000(\frac{L}{a})^3$ embedded in an infinite domain, thus an elastic energy of strain of this cube worth $W_0 = 40000 a \cdot u$. and by resuming the values of the case previously studied, one obtains, by imposing $\epsilon_{33} = 1$ in B.1, on the faces perpendicular to the tensile axis for the mean

<footnote n="15">¹⁵The absolute values are necessary because, depending on the boundary conditions and stiffness, the energy supplements can be positive or negative.</footnote>
displacements 1.00001952 * L and the mean stress $\sigma_{33}$ 0.99998463 a. u.. In order to adjust the boundary conditions on average to obtain $< \epsilon_{33} >$= 1 it is necessary to multiply $\epsilon^{\infty}$ by a coefficient worth 0.9999805. Similarly, in order to obtain $< \sigma_{33} >$= 1 a. u. it is necessary to multiply $\epsilon^{\infty}$ by a coefficient worth 1.00001537. The Table 3 shows that indeed we tend to the theoretical values in infinite domain by adjusting properly the boundary conditions.

<table><caption>Table 3: Deformation energy and energy supplement as a function of multiplier coefficient for a cube containing inhomogeneity.</caption><thead><tr><td>L/a (-)</td><td>coefficient</td><td>W (a. u.)</td><td>W - W<sub>0</sub> (a. u.)</td></tr></thead><tbody><tr><td>21.5443469</td><td>0.9999805</td><td>39998.69</td><td>-1.31</td></tr><tr><td>21.5443469</td><td>1.0000154</td><td>40001.48</td><td>+1.48</td></tr></tbody></table>

## 9. Conclusion

It is therefore clear that the numerical implementation of solutions adapted to infinite domains must be treated with vigilance when considering problems in very large - but finite - domains in front of the size of the inhomogeneities. Indeed, the very slight correction proposed here to be performed on the loading, even if it may seem negligible, is essential to obtain the desired results. It is shown here that the $\Delta W_p/W_0$ terms are of the same order of magnitude as the corrective terms of $\epsilon^{<}/\epsilon^{\infty}$ whatever the size of a finite domain. Moreover, once these aspects are well mastered, it is possible to consider the development of an alternative variational method to the E. I. M. proposed by (Moschovidis and Mura, 1975), (Hashin and Shtrikman, 1962), (Hill, 1963). Indeed, while using the same functional basis from exact analytical solutions, one can hope to obtain better numerical results by trying to minimize the strain energy in order to minimize the stress discontinuities at the interfaces of the inhomogeneities.

## Appendix A. An hollow sphere under pressure

Considering small strains and displacements and linear elasticity, the analytical solution for an hollow sphere of internal radius $R_{int}$ and external radius $R_{ext}$ under an internal pressure $p_{int}$ and an external pressure $p_{ext}$ is
given by:

$$ \sigma_{rr}(r) = \frac{p_{int}R_{int}^3 - p_{ext}R_{ext}^3}{R_{ext}^3 - R_{int}^3} - \frac{(p_{int} - p_{ext})R_{int}^3 R_{ext}^3}{(R_{ext}^3 - R_{int}^3)r^3} \quad (A.1) $$

where r is the distance from the centre of the sphere and $\sigma_{rr}$ is the radial stress. The displacement $U_r$ at a distance r is given by:

$$ U_r(r) = \frac{r}{2\mu + 3\lambda} \frac{p_{int}R_{int}^3 - p_{ext}R_{ext}^3}{R_{ext}^3 - R_{int}^3} + \frac{1}{4\mu} \frac{(p_{int} - p_{ext})R_{int}^3 R_{ext}^3}{(R_{ext}^3 - R_{int}^3)r^2} \quad (A.2) $$

where $\lambda$ and $\mu$ are the Lamé coefficients. The strain energy calculated from the work of external forces is given by:

$$ W = \frac{1}{2} p_{int} U_r(R_{int}) + \frac{1}{2} p_{ext} U_r(R_{ext}) \quad (A.3) $$

## Appendix B. Equivalence equation

Recall the equation of equivalence (Eshelby, 1957) :

$$ \underline{\underline{C}}^0(\underline{\underline{\epsilon}}^\infty + \underline{\underline{A}}^p \cdot \underline{\underline{\beta}}^p - \underline{\underline{\beta}}^p) = \underline{\underline{\underline{C}}}^p(\underline{\underline{\epsilon}}^\infty + \underline{\underline{\underline{A}}}^p \cdot \underline{\underline{\beta}}^p) \quad (B.1) $$

which can be rewritten as:

$$ ([\underline{\underline{C}}^0 - \underline{\underline{C}}^p]^{-1} \cdot [\underline{\underline{C}}^0] - \underline{\underline{A}}^p) \cdot \underline{\underline{\beta}}^p = \underline{\underline{\epsilon}}^\infty \quad (B.2) $$

## Appendix C. Example of fields far from the inclusion: mixed boundary conditions

For the case of the example cited in section 6, Fig. C.3 and Fig. C.4 provide respectively the perturbations, on the face of the cube corresponding to $x_3 = 10 a$, of the displacement fields $U_3$ and of the $\sigma_{33}$ component stress field, related to the presence of an elastic inhomogeneity. These perturbations correspond to the analytical solution calculated for infinite domains.
<figure><img src="image_5.png" /><figcaption>Figure C.3: Perturbations at $x_3 = 10a$ of the displacement field $U_3$ due to the presence of an elastic inhomogeneity of radius $a$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure C.4: Perturbations at $x_3 = 10a$ of the displacement field $\sigma_{33}$ due to the presence of an elastic inhomogeneity of radius $a$.</figcaption></figure>
Figs. C.5 and C.6 compare the analytical solution obtained for an infinite domain with that of the finite element model. The results of the analytical solution are calculated to satisfy average strain conditions on the finite domain $\mathcal{D}$, i. e. $\langle \underline{\epsilon} \rangle_{\mathcal{D}} = \underline{\epsilon}^{\infty}$ or to satisfy these conditions on an infinite domain, i. e. $\langle \underline{\epsilon} \rangle_{\text{infinite domain}} = \underline{\epsilon}^{\infty}$.

<figure><img src="image_7.png" /><figcaption>Figure C.5: Comparison of the results concerning the displacements $U_3$ along the $x_3$ axis, for $x_1 = x_2 = 0$, of the analytical solution and that of the finite element model. $\Delta U_3 = U_3 - \langle \epsilon_{33} \rangle > x_3$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure C.6: Comparison of the results concerning the stress component $\sigma_{33}$ along the $x_3$ axis, for $x_1 = x_2 = 0$, of the analytical solution and that of the finite element model.</figcaption></figure>
References

Benedikt, B., Lewis, M., Rangaswamy, P., 2006. On elastic interactions between spherical inclusions by the equivalent inclusion method. Computational Materials Science 37, 380–392.

Benveniste, Y., 1987. A new approach to the application of mori-tanaka's theory in composite materials. Mechanics of Materials 6, 147–157.

Berthaud, Y., Fond, C., Brun, P., 1994. Interactions between cracks and circular holes in bidimensionnal elastic media. Mechanics Research Communications 21, 525–533.

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.

Brisard, S., Sab, K., Dormieux, L., 2013. New boundary conditions for the computation of the apparent stiffness of statistical volume elements. Journal of the Mechanics and Physics of Solids 61, 2638–2658.

CAST3M, C.E.A., 2022. Logiciel de calculs aux éléments finis, www-cast3m.cea.fr.

Doghri, I., Brassart, L., Adam, L., Gérard, J.S., 2011. A second-moment incremental formulation for the mean-field homogenization of elasto-plastic composites. International Journal of Plasticity 27, 352–371.

Dyson, F., et al., 1891. The potentials of ellipsoids of variable densities. QJ Pure Appl. Math 25, 259–288.

Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proceedings of the royal society of London. Series A. Mathematical and physical sciences 241, 376–396.

Eshelby, J.D., 1959. The elastic field outside an ellipsoidal inclusion. Proceedings of the royal society of London. Series A. Mathematical and physical sciences 252, 561–569.

Ferrers, N., 1877. On the potentials of ellipsoids, ellipsoidal shells, elliptic laminae and elliptic rings of variable densities. QJ Pure Appl. Math 14, 1–22.
Fond, C., 2001. Cavitation criterion for rubber materials: a review of void-growth models. Journal of Polymer Science, part B - Polymer Physics 39, 2081-2096.

Fond, C., Berthaud, Y., 1995. Extensions of the pseudo tractions technique for friction in cracks, circular cavities and external boundaries; effect of the interactions on the homogenised stiffness. International Journal of Fracture 74, 1-28.

Fond, C., Flejou, J.L., Berthaud, Y., 1995. Interactions between cracks and circular holes in bidimensionnal elastic media. European Journal of Mechanics A-solids 14, 73-96.

Fond, C., Gehant, S., Schirrer, R., 2002. Effects of mechanical interactions on the hydrostatic stress in randomly distributed rubber particles in an amorphous polymer matrix. Polymer 43, 909-919.

Fond, C., Lobbrecht, A., Schirrer, R., 1996. Polymers toughened with rubber microspheres; an analytical solution for stresses and strains in the rubber particles at equilibrium and rupture. International Journal of fracture 77, 141-159.

Fond, C., Riccardi, A., Schirrer, R., Montheillet, F., 2001. Mechanical interaction between spherical inhomogeneities: an assessment of a method based on the equivalent inclusion. European Journal of Mechanics A-solids 20, 59-75.

Hashin, Z.a., Shtrikman, S., 1962. On some variational principles in anisotropic and nonhomogeneous elasticity. Journal of the Mechanics and Physics of Solids 10, 335-342.

Hill, R., 1963. New derivation of some elastic extremum principles. Progress in applied mechanics, The Prager anniversary volume, 99-106.

Moschovidis, Z.A., Mura, T., 1975. Two-ellipsoidal inhomogeneities by the equivalent inclusion method. Journal of Applied Mechanics 42, 847-852.

Mura, T., 1993. Micromechanics of Defects in Solids. Kluwer Academic Publishers, Dordrecht/Boston/London. revised edition, reprinted edition.
Nakasone, Y., Nishiyama, H., Nojiri, T., 2000. Numerical equivalent inclusion method: a new computational method for analyzing stress fields in and around inclusions of various shapes. Materials Science and Engineering: A 285, 229-238.

Novák, J., Łukasz Kaczmarczyk, Grassl, P., Zeman, J., Pearce, C.J., 2012. A micromechanics-enhanced finite element formulation for modelling heterogeneous materials. Computer Methods in Applied Mechanics and Engineering 201-204, 53-64.

Shodja, H., Rad, I., Soheilifard, R., 2003. Interacting cracks and ellipsoidal inhomogeneities by the equivalent inclusion method. Journal of the Mechanics and Physics of Solids 51, 945-960.

Xu, Y., Xiang, M., Yu, J., Tang, S., Chen, J., 2023. A variational fracture method based on eshelby transformation. European Journal of Mechanics - A/Solids 97, 104846.

Yan, Z., Zhang, M., Jiang, S., 2021. Numerical methods for solving the equivalent inclusion equation in semi-analytical models. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 235, 2455-2473.

Yang, W., Zhou, Q., Wang, J., Khoo, B.C., Phan-Thien, N., 2021. Equivalent inclusion method for arbitrary cavities or cracks in an elastic infinite/semi-infinite space. International Journal of Mechanical Sciences 195, 106259.

Yu, X.W., Wang, Z.W., Wang, H., 2021. Multiple ellipsoidal/elliptical inhomogeneities embedded in infinite matrix by equivalent inhomogeneous inclusion method. Mathematics and Mechanics of Solids, 108128652110071.