Contents

This post collects every major analytical and semi-analytical model for the effective elastic moduli of two-phase composites in one place. It is intended as a reference hub: other posts on this site link here for derivations and DOI-linked primary sources. For numerical homogenization (FFT-based), see the FFT Homogenization post.

Overview

The table below summarises all models covered. They are grouped by the strength of their theoretical foundation.

Model Type Key assumption
Voigt Exact upper bound Uniform strain
Reuss Exact lower bound Uniform stress
Hashin-Shtrikman Tightest variational bounds Statistically isotropic, spherical inclusions
Dilute Estimate, $v_f \to 0$ No inclusion interactions
Mori-Tanaka Estimate (= HS lower) Mean-field; inclusions see matrix average strain
Self-consistent Estimate Inclusion embedded in unknown effective medium
Halpin-Tsai Semi-empirical Fitted reinforcement parameter $\xi$
Ponte Castañeda-Willis Higher-order bounds Spatial pair-distribution of inclusions

The Two-Phase RVE

The reference problem is a single spherical inclusion (fiber/particle) of radius $a$ embedded in an infinite matrix. The surrounding matrix is taken as the Representative Volume Element (RVE). This is the exact geometry for which the Eshelby, HS, MT and self-consistent solutions are derived.

Interactive 3D RVE — drag to rotate, scroll to zoom. Single spherical inclusion (blue, transparent) in a cubic matrix (grey wireframe). Load arrows show applied stress σ. Coordinate axes: x₁ (red), x₂ (green), x₃ (blue).

The volume fraction of the inclusion phase is

\[v_f = \frac{V_{\rm inclusion}}{V_{\rm RVE}}\]

and the matrix occupies $v_m = 1 - v_f$. All models express the effective elastic moduli as functions of $v_f$, the phase moduli $E_f, E_m$, and the phase Poisson ratios $\nu_f, \nu_m$.

Phase bulk and shear moduli:

\[K_i = \frac{E_i}{3(1-2\nu_i)}, \qquad G_i = \frac{E_i}{2(1+\nu_i)} \tag{1}\]

Voigt and Reuss Bounds

The Voigt and Reuss bounds are the widest possible rigorous bounds — valid for any microstructure geometry.

Longitudinal modulus E1

Under a strain applied along the fiber axis (direction 1), iso-strain (Voigt) is exact for an ideally bonded continuous fiber:

\[E_1^{\rm V} = v_f E_f + (1 - v_f) E_m \tag{2}\]

The iso-stress (Reuss) lower bound:

\[\frac{1}{E_1^{\rm R}} = \frac{v_f}{E_f} + \frac{1 - v_f}{E_m} \tag{3}\]

For continuous fibers, the rule of mixtures (Voigt) is the practically correct prediction for $E_1$.

Transverse modulus E2

Applied transversely (direction 2), the same mathematical expressions bound $E_2$:

\[E_2^{\rm V} = v_f E_f + (1 - v_f) E_m \tag{4}\] \[\frac{1}{E_2^{\rm R}} = \frac{v_f}{E_f} + \frac{1 - v_f}{E_m} \tag{5}\]

These bounds are very wide for typical fiber/matrix stiffness ratios. For a tight estimate of $E_2$, use the Hashin-Shtrikman bounds below.

References: Voigt [1]; Reuss [2]

See also: Interactive comparison chart — all models on one plot.

Hashin-Shtrikman Bounds

The Hashin-Shtrikman (HS) bounds [4] are the tightest possible bounds for a statistically isotropic microstructure of spherical inclusions, given only $v_f$ and the phase moduli. They are derived by a variational principle using a polarisation field.

For a 3D isotropic two-phase composite with spherical inclusions, the bounds on bulk modulus $K$ and shear modulus $G$ are:

HS lower bounds (assuming $K_f > K_m$, $G_f > G_m$):

\[K^{-} = K_m + \frac{v_f}{\dfrac{1}{K_f - K_m} + \dfrac{3(1-v_f)}{3K_m + 4G_m}} \tag{6}\] \[G^{-} = G_m + \frac{v_f}{\dfrac{1}{G_f - G_m} + \dfrac{6(1-v_f)(K_m + 2G_m)}{5G_m(3K_m + 4G_m)}} \tag{7}\]

HS upper bounds:

\[K^{+} = K_f + \frac{1-v_f}{\dfrac{1}{K_m - K_f} + \dfrac{3v_f}{3K_f + 4G_f}} \tag{8}\] \[G^{+} = G_f + \frac{1-v_f}{\dfrac{1}{G_m - G_f} + \dfrac{6v_f(K_f + 2G_f)}{5G_f(3K_f + 4G_f)}} \tag{9}\]

Effective Young’s modulus from $K^{\pm}$ and $G^{\pm}$:

\[E^{\pm} = \frac{9K^{\pm}G^{\pm}}{3K^{\pm} + G^{\pm}} \tag{10}\]

References: Hashin & Shtrikman [4]; Hill [5]

The Eshelby Tensor

Eshelby’s Inclusion Problem

Eshelby [3] showed that an ellipsoidal inclusion that would undergo a stress-free eigenstrain $\boldsymbol{\varepsilon}^*$ if unconstrained instead develops a uniform constrained strain when embedded in an infinite matrix:

\[\varepsilon^c_{ij} = S_{ijkl}\, \varepsilon^*_{kl} \tag{11}\]

The fourth-order tensor $\mathbb{S}$ depends only on the shape of the inclusion and the Poisson ratio $\nu_m$ of the surrounding matrix — not on the moduli contrast. This makes it the universal kernel in every inclusion-based homogenisation model.

Components for a Sphere

For a spherical inclusion embedded in an isotropic matrix with Poisson ratio $\nu_m$, all non-zero independent Eshelby tensor components have closed-form expressions:

\[S_{1111} = S_{2222} = S_{3333} = \frac{7 - 5\nu_m}{15(1-\nu_m)} \tag{12}\] \[S_{1122} = S_{2233} = S_{1133} = \frac{5\nu_m - 1}{15(1-\nu_m)} \tag{13}\] \[S_{1212} = S_{2323} = S_{1313} = \frac{4 - 5\nu_m}{15(1-\nu_m)} \tag{14}\]

All other components are zero. The hydrostatic and deviatoric projections are:

\[S_{\rm hyd} = S_{iikk}/9 = \frac{1+\nu_m}{3(1-\nu_m)} \tag{15}\] \[S_{\rm dev} = S_{1212} = \frac{4 - 5\nu_m}{15(1-\nu_m)} \tag{16}\]

From these components the HS scalar factors follow: $\beta_K = 3K_m/(3K_m + 4G_m)$ and $\beta_G = G_m(9K_m + 8G_m)/[6G_m(K_m + 2G_m)]$, which appear in the HS bound denominators (eqs. 6–9) and in the dilute-limit expressions (eqs. 19–20).

Full derivation for finite aspect ratio (prolate spheroid, $e = 1$ to $\infty$): see the Eshelby Tensor post.

References: Eshelby [3]; Mura [9]

Mori-Tanaka

The Mori-Tanaka (MT) model [6] places each inclusion in the average strain of the matrix rather than the far-field applied strain, accounting for inclusion-inclusion interactions in a mean-field sense.

For spherical inclusions with $K_f > K_m$, $G_f > G_m$:

\[K^{\rm MT} = K_m + \frac{v_f}{\dfrac{1}{K_f - K_m} + \dfrac{3(1-v_f)}{3K_m + 4G_m}} \tag{17}\] \[G^{\rm MT} = G_m + \frac{v_f}{\dfrac{1}{G_f - G_m} + \dfrac{6(1-v_f)(K_m + 2G_m)}{5G_m(3K_m + 4G_m)}} \tag{18}\]

MT = HS lower bound: Equations (17)–(18) are algebraically identical to the HS lower bounds (6)–(7). Benveniste [10] proved this rigorously: when inclusions are stiffer than the matrix, the Mori-Tanaka estimate coincides exactly with the HS lower bound.

References: Mori & Tanaka [6]; Benveniste [10]

Dilute Limit

When $v_f \to 0$, each inclusion is embedded in the unperturbed matrix strain field. The effective bulk and shear moduli are:

\[K^{\rm dil} = K_m + v_f \frac{(K_f - K_m)(3K_m + 4G_m)}{3K_f + 4G_m} \tag{19}\] \[G^{\rm dil} = G_m + v_f \frac{5G_m(G_f - G_m)(3K_m + 4G_m)}{G_m(9K_m + 8G_m) + 6G_f(K_m + 2G_m)} \tag{20}\]

These expressions are valid for $v_f \lesssim 0.1$. At higher volume fractions, inclusion interactions are significant and the dilute estimate overpredicts $E_{\rm eff}$.

Self-Consistent Scheme

Hill [7] proposed embedding the inclusion in the unknown effective medium rather than the matrix. The effective moduli satisfy the implicit equations:

\[K^{\rm SC} = \left[\frac{v_f K_f}{K_f + \frac{4}{3}G^{\rm SC}} + \frac{v_m K_m}{K_m + \frac{4}{3}G^{\rm SC}}\right] \left[\frac{v_f}{K_f + \frac{4}{3}G^{\rm SC}} + \frac{v_m}{K_m + \frac{4}{3}G^{\rm SC}}\right]^{-1} \tag{21}\] \[G^{\rm SC} = \left[\frac{v_f G_f}{G_f + \beta^{\rm SC}} + \frac{v_m G_m}{G_m + \beta^{\rm SC}}\right] \left[\frac{v_f}{G_f + \beta^{\rm SC}} + \frac{v_m}{G_m + \beta^{\rm SC}}\right]^{-1} \tag{22}\]

where

\[\beta^{\rm SC} = \frac{G^{\rm SC}(9K^{\rm SC} + 8G^{\rm SC})}{6(K^{\rm SC} + 2G^{\rm SC})} \tag{23}\]

These equations are solved iteratively (typically 20–100 iterations from an MT initial guess). The self-consistent scheme is symmetric in the two phases and is most appropriate for polycrystal-like microstructures where neither phase forms a continuous matrix. It can violate the HS bounds at extreme volume fractions or high phase contrast ($E_f/E_m \gtrsim 100$).

Reference: Hill [7]

Halpin-Tsai

The Halpin-Tsai equations [8] are semi-empirical, with a reinforcement parameter $\xi$ fitted to geometry:

\[E^{\rm HT} = E_m \frac{1 + \xi\,\eta\, v_f}{1 - \eta\, v_f}, \qquad \eta = \frac{E_f/E_m - 1}{E_f/E_m + \xi} \tag{24}\]

Recommended values of $\xi$:

Geometry Property $\xi$
Spherical inclusion $E$ $2$
Circular fiber (infinite) $E_2$ $2$
Short fiber (aspect ratio $l/d$) $E_2$ $2l/d$
Any fiber $G_{12}$ $1$ (use $G_f/G_m$)
Any fiber $E_1$ $\to\infty$ (recovers Voigt)

The Halpin-Tsai formula has no rigorous micromechanical derivation; it is accurate to ~5% for $v_f \lesssim 0.6$ at typical glass/epoxy contrast ($E_f/E_m \approx 20$).

Reference: Halpin & Kardos [8]

Higher-Order: Ponte Castaneda-Willis

The Ponte Castañeda-Willis (PCW) bounds [12] tighten the HS bounds by incorporating the spatial distribution of inclusions through a second morphological descriptor: the distribution tensor $\mathbb{P}^d$, which characterises the shape of the “security zone” — the excluded volume around each inclusion.

For a composite with spherical inclusions and a spherical spatial distribution (inclusions arranged with isotropic pair statistics, non-overlapping), $\mathbb{P}^d$ takes the same isotropic form as the Eshelby polarisation tensor evaluated at the matrix moduli. The PCW bounds on bulk and shear modulus are:

\[K^{\rm PCW,\pm} = K_{\rm ref}^{\pm} + \frac{v_f (K_f - K_{\rm ref}^{\pm}) (K_m - K_{\rm ref}^{\pm})}{v_f (K_m - K_{\rm ref}^{\pm}) + v_m (K_f - K_{\rm ref}^{\pm}) + \beta_K^{\pm}(K_f - K_{\rm ref}^{\pm})(K_m - K_{\rm ref}^{\pm}) / \delta_K^{\pm}} \tag{25}\]

where the reference moduli $K_{\rm ref}^{\pm}$ and $G_{\rm ref}^{\pm}$ are the HS reference moduli (stiffer phase for upper, softer for lower), and the scalar factors for spherical geometry are:

\[\beta_K = \frac{3K_{\rm ref} + 4G_{\rm ref}}{3(K_{\rm ref} + \tfrac{4}{3}G_{\rm ref})} = 1, \qquad \delta_K = K_{\rm ref} + \tfrac{4}{3}G_{\rm ref} \tag{26}\]

For spherical inclusions with a spherical spatial distribution, the PCW lower bound equals the HS lower bound (PCW lower = HS lower = MT). The PCW upper bound uses $K_{\rm ref} = K_f$, $G_{\rm ref} = G_f$, with “inclusions” = matrix phase at volume fraction $(1-v_f)$:

\[K^{\rm PCW,-} = K^{\rm HS,-} \tag{27}\] \[K^{\rm PCW,+} = K^{\rm HS,+} - \frac{(1-v_f)^2(K_m - K_f)^2 / \alpha_f}{\alpha_f + (1-v_f)(K_m - K_f)}, \qquad \alpha_f = K_f + \tfrac{4}{3}G_f \tag{28}\] \[G^{\rm PCW,+} = G^{\rm HS,+} - \frac{(1-v_f)^2(G_m - G_f)^2 \beta_f}{[1 + (1-v_f)(G_m - G_f)\beta_f]^2}, \qquad \beta_f = \frac{6(K_f+2G_f)}{5G_f(3K_f+4G_f)} \tag{29}\]

The PCW upper bound is strictly tighter than the HS upper bound for $0 < v_f < 1$:

\[K^{\rm PCW,+} \leq K^{\rm HS,+}, \qquad G^{\rm PCW,+} \leq G^{\rm HS,+} \tag{30}\]

The correction vanishes at $v_f = 0$ and $v_f = 1$, reaching its maximum at intermediate $v_f$ (roughly 0.3–0.6).

Reference: Ponte Castañeda & Willis [12]

Deep dive: How Microstructure Arrangement Affects Composite Bounds: A PCW Investigation — interactive 3D scenes showing isotropic, layered, columnar and periodic microstructures with live PCW bound charts for varying $w$.

Interactive Comparison

All models are computed in the browser. Adjust phase properties; charts update live. The self-consistent result is obtained by iterating eqs. (21)–(23) to convergence (tol = 10⁻⁸). The PCW upper bound (eqs. 28–29) is tighter than the HS upper; the PCW lower = HS lower (eq. 27).

Python Reference Implementation

import numpy as np

def phase_KG(E, nu):
    K = E / (3*(1 - 2*nu))
    G = E / (2*(1 + nu))
    return K, G

def E_from_KG(K, G):
    return 9*K*G / (3*K + G)

def hs_bounds(Ef, Em, nuf, num, vf):
    """Hashin-Shtrikman bounds for 3D isotropic composite (spherical inclusions)."""
    Kf, Gf = phase_KG(Ef, nuf)
    Km, Gm = phase_KG(Em, num)
    vm = 1 - vf
    K_lo = Km + vf / (1/(Kf-Km) + 3*vm/(3*Km+4*Gm))
    K_hi = Kf + vm / (1/(Km-Kf) + 3*vf/(3*Kf+4*Gf))
    G_lo = Gm + vf / (1/(Gf-Gm) + 6*vm*(Km+2*Gm)/(5*Gm*(3*Km+4*Gm)))
    G_hi = Gf + vm / (1/(Gm-Gf) + 6*vf*(Kf+2*Gf)/(5*Gf*(3*Kf+4*Gf)))
    return (E_from_KG(K_lo, G_lo), E_from_KG(K_hi, G_hi),
            K_lo, K_hi, G_lo, G_hi)

def dilute(Ef, Em, nuf, num, vf):
    """Dilute limit (valid for vf < 0.1)."""
    Kf, Gf = phase_KG(Ef, nuf)
    Km, Gm = phase_KG(Em, num)
    K = Km + vf*(Kf-Km)*(3*Km+4*Gm) / (3*Kf+4*Gm)
    G = Gm + vf*5*Gm*(Gf-Gm)*(3*Km+4*Gm) / (Gm*(9*Km+8*Gm) + 6*Gf*(Km+2*Gm))
    return E_from_KG(K, G)

def self_consistent(Ef, Em, nuf, num, vf, tol=1e-8, max_iter=300):
    """
    Self-consistent scheme (Hill 1965).
    Fixed-point iteration starting from Mori-Tanaka initial guess.
    """
    Kf, Gf = phase_KG(Ef, nuf)
    Km, Gm = phase_KG(Em, num)
    vm = 1 - vf
    # Initial guess: HS lower (= MT)
    K = Km + vf / (1/(Kf-Km) + 3*vm/(3*Km+4*Gm))
    G = Gm + vf / (1/(Gf-Gm) + 6*vm*(Km+2*Gm)/(5*Gm*(3*Km+4*Gm)))
    for _ in range(max_iter):
        alpha = 4*G/3
        beta  = G*(9*K + 8*G) / (6*(K + 2*G))
        K_new = (vf*Kf/(Kf+alpha) + vm*Km/(Km+alpha)) / (vf/(Kf+alpha) + vm/(Km+alpha))
        G_new = (vf*Gf/(Gf+beta)  + vm*Gm/(Gm+beta))  / (vf/(Gf+beta)  + vm/(Gm+beta))
        if abs(K_new - K) < tol and abs(G_new - G) < tol:
            K, G = K_new, G_new
            break
        K, G = K_new, G_new
    return E_from_KG(K, G)

def pcw_upper(Ef, Em, nuf, num, vf):
    """
    Ponte Castaneda-Willis upper bound (PCW 1995) for spherical inclusions
    with a spherical (isotropic) spatial distribution.
    PCW lower = HS lower (no correction); PCW upper is tighter than HS upper.
    """
    Kf, Gf = phase_KG(Ef, nuf)
    Km, Gm = phase_KG(Em, num)
    vm = 1 - vf

    # Bulk modulus upper bound (ref = fiber; eq. 28)
    K_hs_hi = Kf + vm / (1/(Km-Kf) + 3*vf/(3*Kf+4*Gf))
    alphaF  = Kf + 4*Gf/3
    num_cK  = vm*vm*(Km-Kf)**2 / alphaF
    den_cK  = alphaF + vm*(Km-Kf)
    K_pcw   = K_hs_hi - num_cK / den_cK

    # Shear modulus upper bound (ref = fiber; eq. 29)
    betaF   = 6*(Kf + 2*Gf) / (5*Gf*(3*Kf + 4*Gf))
    G_hs_hi = Gf + vm / (1/(Gm-Gf) + 6*vf*(Kf+2*Gf)/(5*Gf*(3*Kf+4*Gf)))
    dGmf    = Gm - Gf
    num_cG  = vm*vm*dGmf**2 * betaF
    den_cG  = (1 + vm*dGmf*betaF)**2 / betaF
    G_pcw   = G_hs_hi - num_cG / den_cG

    return E_from_KG(K_pcw, G_pcw), K_pcw, G_pcw

def halpin_tsai(Ef, Em, vf, xi=2):
    """Halpin-Tsai (xi=2 for sphere/circular fiber; xi=2*l/d for short fiber)."""
    eta = (Ef/Em - 1) / (Ef/Em + xi)
    return Em * (1 + xi*eta*vf) / (1 - eta*vf)


# Example: glass spheres (Ef=74 GPa, nuf=0.22) in epoxy (Em=3.5 GPa, num=0.35)
vf = np.linspace(0, 0.7, 200)
Elo, Ehi, *_ = hs_bounds(74e3, 3.5e3, 0.22, 0.35, vf)
Esc  = np.array([self_consistent(74e3, 3.5e3, 0.22, 0.35, v) for v in vf])
Eht  = halpin_tsai(74e3, 3.5e3, vf)

Epcw = np.array([pcw_upper(74e3, 3.5e3, 0.22, 0.35, v)[0] for v in vf])

idx = 60  # vf ~ 0.30
print(f"At vf=0.30:  HS  [{Elo[idx]:.0f}, {Ehi[idx]:.0f}] MPa")
print(f"             PCW [{Elo[idx]:.0f}, {Epcw[idx]:.0f}] MPa  (lower unchanged)")
print(f"             SC  = {Esc[idx]:.0f} MPa")
print(f"             HT  = {Eht[idx]:.0f} MPa")

Conclusion and Recommendations

Which method should you use?

Situation Recommended model Reason
Need a guaranteed bound, any microstructure Voigt + Reuss Widest but always valid
Statistically isotropic matrix-inclusion composite Hashin-Shtrikman lower bound Tightest rigorous bound; equals MT
Engineering estimate, moderate $v_f$ Halpin-Tsai ($\xi=2$) Simple, calibrated, accurate to ~5%
Symmetric two-phase (neither phase dominant) Self-consistent (Hill) Treats both phases equally
Very low $v_f$ ($< 0.10$) Dilute limit Exact in dilute regime
Tightest upper bound, isotropic distribution Ponte Castañeda-Willis upper Strictly tighter than HS upper; PCW lower = HS lower
Any geometry, full-field accuracy FFT homogenization Numerical, exact for periodic cell

Key results from the charts

  • HS lower = Mori-Tanaka when inclusions are stiffer than the matrix — a single curve covers both.
  • Self-consistent sits close to the HS lower bound at low $v_f$ but rises above it at moderate $v_f$, and can exit the HS band at very high volume fractions.
  • Halpin-Tsai closely tracks the HS band centre for typical glass/epoxy contrast ($E_f/E_m \approx 20$). At higher contrast it tends toward the HS upper bound.
  • PCW upper bound is visibly tighter than HS upper at intermediate $v_f$ (0.3–0.6). The PCW lower bound equals HS lower for an isotropic spatial distribution.
  • Dilute limit diverges from MT/HS above $v_f \approx 0.10$ — do not use it outside this range.
  • For a continuous fiber composite, $E_1$ is predicted almost exactly by Voigt (rule of mixtures) regardless of model. Model choice matters almost exclusively for transverse properties.

Practical guidance

Halpin-Tsai is the simplest engineering estimate; Mori-Tanaka / HS lower gives a tighter rigorous bound. FFT homogenization gives full-field accuracy for periodic RVEs. The self-consistent scheme is standard for polycrystal plasticity but rarely appropriate for matrix-inclusion composites.

References

[1] Voigt, W. (1889). Über die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper. Annalen der Physik, 274(12), 573–587.

[2] Reuss, A. (1929). Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle. ZAMM, 9(1), 49–58. doi:10.1002/zamm.19290090104

[3] Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proceedings of the Royal Society A, 241(1226), 376–396. doi:10.1098/rspa.1957.0133

[4] Hashin, Z. & Shtrikman, S. (1963). A variational approach to the theory of the elastic behaviour of multiphase materials. Journal of the Mechanics and Physics of Solids, 11(2), 127–140. doi:10.1016/0022-5096(63)90060-7

[5] Hill, R. (1963). Elastic properties of reinforced solids: some theoretical principles. Journal of the Mechanics and Physics of Solids, 11(5), 357–372. doi:10.1016/0022-5096(63)90036-X

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

[7] Hill, R. (1965). A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids, 13(4), 213–222. doi:10.1016/0022-5096(65)90011-6

[8] Halpin, J. C. & Kardos, J. L. (1976). The Halpin-Tsai equations: A review. Polymer Engineering & Science, 16(5), 344–352. doi:10.1002/pen.760160512

[9] Mura, T. (1987). Micromechanics of Defects in Solids (2nd ed.). Martinus Nijhoff.

[10] Benveniste, Y. (1987). A new approach to the application of Mori-Tanaka’s theory in composite materials. Mechanics of Materials, 6(2), 147–157. doi:10.1016/0167-6636(87)90005-6

[11] Christensen, R. M. (1979). Mechanics of Composite Materials. Wiley.

[12] 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(12), 1919–1951. doi:10.1016/0022-5096(95)00052-6

[13] Moulinec, H. & Suquet, P. (1994). A fast numerical method for computing the linear and nonlinear mechanical properties of composites. Comptes Rendus de l’Académie des Sciences, Série II, 318, 1417–1423.

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