Basem Rajjoub
Eshelby's Inclusion: Interactive 3D Visualization
Contents
The Eshelby tensor $\mathbb{S}$ maps the eigenstrain an ellipsoidal inclusion would freely undergo to the constrained strain it actually develops when embedded in an elastic matrix:
\[\boldsymbol{\varepsilon}^c = \mathbb{S}{:}\boldsymbol{\varepsilon}^*\]$\mathbb{S}$ depends only on the inclusion shape and the matrix Poisson ratio $\nu_m$, not on the moduli contrast. This makes it the kernel of every inclusion-based homogenization scheme.
See also: Eshelby tensor components vs. aspect ratio · Analytical bounds and homogenization models · Tensor visualization (rank 2, 3, 4) · FFT homogenization
RVE Homogenization — Mori-Tanaka
Enter matrix and fiber isotropic properties, fiber volume fraction, and aspect ratio. The Mori-Tanaka scheme uses the Eshelby tensor to compute the full effective stiffness $\mathbb{C}^{\rm eff}$ for a unidirectional composite (fiber along $x_1$). The result is a transversely isotropic Voigt 6×6 matrix plus the five engineering constants $E_1, E_2, G_{12}, \nu_{12}, \nu_{23}$.
RVE with ellipsoidal inclusion (blue). Fiber along x₁ (red). Drag to orbit.
Theory
Eigenstrain and the inclusion problem
An eigenstrain $\boldsymbol{\varepsilon}^{\ast}$ is the strain an inclusion would undergo free of its surroundings: a thermal expansion, transformation strain, or plastic misfit. When embedded in an elastic matrix, the matrix resists this deformation and the inclusion develops a constrained strain $\boldsymbol{\varepsilon}^c < \boldsymbol{\varepsilon}^{\ast}$.
Eshelby (1957) [1] showed that for an ellipsoidal inclusion this constrained strain is uniform throughout the interior, a result that does not hold for any other shape. The constrained strain is related to the eigenstrain by:
\[\varepsilon^c_{ij} = S_{ijkl}\,\varepsilon^*_{kl} \tag{1}\]$\mathbb{S}$ depends only on the shape of the ellipsoid and the matrix Poisson ratio $\nu_m$, not on the moduli of the inclusion.
Eshelby conjecture. Only ellipsoidal inclusions produce a uniform interior strain field. Proven for 2D by Sendeckyj (1970) and for 3D by Liu (2008) [3].
Stress field
Inside the inclusion the elastic strain is $\boldsymbol{\varepsilon}^{\rm el} = \boldsymbol{\varepsilon}^c - \boldsymbol{\varepsilon}^{\ast}$ and the stress is
\[\boldsymbol{\sigma} = \mathbb{C}_m : (\boldsymbol{\varepsilon}^c - \boldsymbol{\varepsilon}^{\ast}) \tag{2}\]Both are uniform inside. Outside, stress and strain decay as $r^{-3}$.
The equivalent inclusion method (inhomogeneity problem)
The inclusion problem assumes equal moduli. For an inhomogeneity with $\mathbb{C}_f \neq \mathbb{C}_m$ under far-field strain $\bar{\boldsymbol{\varepsilon}}$, replace it with a fictitious homogeneous inclusion carrying an eigenstrain $\boldsymbol{\varepsilon}^{\ast}$ that produces the same interior strain. Matching stress states gives the equivalent inclusion condition:
\[\mathbb{C}_f : (\bar{\boldsymbol{\varepsilon}} + \boldsymbol{\varepsilon}^c) = \mathbb{C}_m : (\bar{\boldsymbol{\varepsilon}} + \boldsymbol{\varepsilon}^c - \boldsymbol{\varepsilon}^{\ast}) \tag{3}\]Substituting $\boldsymbol{\varepsilon}^c = \mathbb{S}:\boldsymbol{\varepsilon}^{\ast}$ and solving gives the dilute strain concentration tensor
\[\mathbb{A}_{\rm dil} = \bigl[\mathbb{I} + \mathbb{S}:\mathbb{C}_m^{-1}:(\mathbb{C}_f - \mathbb{C}_m)\bigr]^{-1} \tag{4}\]so that $\boldsymbol{\varepsilon}^{\rm inc} = \mathbb{A}_{\rm dil}:\bar{\boldsymbol{\varepsilon}}$.
Mori-Tanaka homogenization
The dilute estimate ($v_f \lesssim 0.05$) treats each inclusion as isolated. Mori-Tanaka [5] accounts for interactions by replacing the far-field with the average matrix strain $\bar{\boldsymbol{\varepsilon}}_m$. The effective stiffness is:
\[\mathbb{C}^{\rm eff} = \mathbb{C}_m + v_f\,(\mathbb{C}_f - \mathbb{C}_m):\mathbb{A}_{\rm dil}:\bigl[(1-v_f)\mathbb{I} + v_f\,\mathbb{A}_{\rm dil}\bigr]^{-1} \tag{5}\]This is what the RVE panel computes. MT coincides with the Hashin-Shtrikman lower bound for stiffer spherical inclusions and stays accurate to $v_f \approx 0.4$. See the analytical bounds post for a full comparison.
Component limits and physical meaning
For a prolate spheroid (semi-axes $a_1 = e \geq 1$, $a_2 = a_3 = 1$, fiber along $x_1$) there are six independent non-zero components. Unlike the stiffness tensor, $\mathbb{S}$ has no major symmetry: $S_{ijkl} \neq S_{klij}$ in general, so $S_{2211} \neq S_{1122}$ for $e \neq 1$. Transverse isotropy of the 2–3 plane requires $S_{2323} = (S_{2222}-S_{2233})/2$.
| Component | Physical meaning | Sphere ($e=1$) | ∞-cylinder |
|---|---|---|---|
| $S_{1111}$ | Axial eigenstrain → axial constrained strain | $\dfrac{7-5\nu_m}{15(1-\nu_m)}$ | $0$ |
| $S_{2222}$ | Transverse eigenstrain → transverse constrained strain | $\dfrac{7-5\nu_m}{15(1-\nu_m)}$ | $\dfrac{5-4\nu_m}{8(1-\nu_m)}$ |
| $S_{1122}$ | Transverse eigenstrain → axial constrained strain | $\dfrac{5\nu_m-1}{15(1-\nu_m)}$ | $0$ |
| $S_{2211}$ | Axial eigenstrain → transverse constrained strain | $\dfrac{5\nu_m-1}{15(1-\nu_m)}$ | $\dfrac{\nu_m}{2(1-\nu_m)}$ |
| $S_{2233}$ | Transverse coupling within the 2–3 plane | $\dfrac{5\nu_m-1}{15(1-\nu_m)}$ | $\dfrac{4\nu_m-1}{8(1-\nu_m)}$ |
| $S_{1212}$ | Shear coupling in the axial-transverse plane | $\dfrac{4-5\nu_m}{15(1-\nu_m)}$ | $\dfrac{1}{4}$ |
$S_{1111} \to 0$ for the infinite cylinder: the matrix cannot constrain axial deformation of an infinitely long fiber, reducing the problem to plane strain in the 2–3 plane. $S_{1122} \to 0$ for the same reason, while $S_{2211} \to \nu_m/(2(1-\nu_m))$ remains non-zero — axial eigenstrain still Poisson-couples into transverse strain.
The Eshelby tensor in closed form
The components of $\mathbb{S}$ for a prolate spheroid ($a_1 = e$, $a_2 = a_3 = 1$) are expressed via five shape integrals $I_1, I_2, I_{11}, I_{12}, I_{22}$ (Eshelby 1957 [1], 1959 [7]; notation following Mura [2]).
Shape integrals
\[\begin{aligned} I_1 &= 2\pi e \!\int_0^{\infty}\! \frac{ds}{(e^2\!+\!s)^{3/2}(1\!+\!s)}, & I_2 = I_3 &= 2\pi e \!\int_0^{\infty}\! \frac{ds}{(1\!+\!s)^{2}\sqrt{e^2\!+\!s}}, \\[6pt] I_{11} &= 2\pi e \!\int_0^{\infty}\! \frac{ds}{(e^2\!+\!s)^{5/2}(1\!+\!s)}, & I_{12} = I_{13} &= 2\pi e \!\int_0^{\infty}\! \frac{ds}{(e^2\!+\!s)^{3/2}(1\!+\!s)^{2}}, \\[6pt] && I_{22} = I_{23} = I_{33} &= 2\pi e \!\int_0^{\infty}\! \frac{ds}{(1\!+\!s)^{3}\sqrt{e^2\!+\!s}}. \end{aligned} \tag{6}\]Since $a_2 = a_3$, the integrands for $I_{22}$ and $I_{23}$ are identical, so $I_{23} \equiv I_{22}$. The first-order integrals satisfy $I_1 + 2I_2 = 4\pi$.
Closed-form evaluation
For a prolate spheroid ($e > 1$, fiber/needle):
\[\begin{aligned} I_1 &= \frac{4\pi}{(e^2-1)^{3/2}}\Bigl[e\,\cosh^{-1}\!e - \sqrt{e^2-1}\,\Bigr] \\[6pt] I_2 &= \frac{2\pi e}{(e^2-1)^{3/2}}\Bigl[e\sqrt{e^2-1} - \cosh^{-1}\!e\,\Bigr] \end{aligned} \tag{7a}\]For an oblate spheroid ($e < 1$, disk/platelet), replace $\cosh^{-1}$ with $\cos^{-1}$:
\[\begin{aligned} I_1 &= \frac{4\pi}{(1-e^2)^{3/2}}\Bigl[\sqrt{1-e^2} - e\,\cos^{-1}\!e\,\Bigr] \\[6pt] I_2 &= \frac{2\pi e}{(1-e^2)^{3/2}}\Bigl[\cos^{-1}\!e - e\sqrt{1-e^2}\,\Bigr] \end{aligned} \tag{7a′}\]Both forms satisfy $I_1 + 2I_2 = 4\pi$ and reduce to $I_1 = I_2 = 4\pi/3$ for the sphere ($e = 1$).
The second-order integrals follow from the recurrence relations $I_{ij} = (I_j - I_i)/(a_i^2 - a_j^2)$ and $3I_{ii} + I_{ij} + I_{ik} = 4\pi/a_i^2$:
\[I_{12} = \frac{I_2 - I_1}{e^2 - 1}, \qquad I_{11} = \frac{\dfrac{4\pi}{e^2} - 2\,I_{12}}{3}, \qquad I_{22} = \frac{4\pi - I_{12}}{4}. \tag{7b}\]In the interactive panel above, these are computed by numerical quadrature (midpoint rule, 256 panels) rather than the closed forms, which avoids catastrophic cancellation near $e = 1$.
Component formulas
\[\begin{aligned} S_{1111} &= \frac{3\,e^2\,I_{11} + (1-2\nu)\,I_1}{8\pi(1-\nu)} \\[6pt] S_{2222} = S_{3333} &= \frac{3\,I_{22} + (1-2\nu)\,I_2}{8\pi(1-\nu)} \\[6pt] S_{1122} = S_{1133} &= \frac{I_{12} - (1-2\nu)\,I_1}{8\pi(1-\nu)} \\[6pt] S_{2211} = S_{3311} &= \frac{e^2\,I_{12} - (1-2\nu)\,I_2}{8\pi(1-\nu)} \\[6pt] S_{2233} = S_{3322} &= \frac{I_{22} - (1-2\nu)\,I_2}{8\pi(1-\nu)} \\[6pt] S_{1212} = S_{1313} &= \frac{(e^2+1)\,I_{12} + (1-2\nu)(I_1+I_2)}{16\pi(1-\nu)} \end{aligned} \tag{8}\]Transverse isotropy demands $S_{2323} = (S_{2222} - S_{2233})/2$, which can be verified by substituting equations (8). All other components vanish.
Full 6×6 Voigt matrix
In Voigt notation ($1{=}11$, $2{=}22$, $3{=}33$, $4{=}23$, $5{=}13$, $6{=}12$) the Eshelby tensor is
\[\mathbb{S} = \begin{pmatrix} S_{1111} & S_{1122} & S_{1122} & 0 & 0 & 0 \\ S_{2211} & S_{2222} & S_{2233} & 0 & 0 & 0 \\ S_{2211} & S_{2233} & S_{2222} & 0 & 0 & 0 \\ 0 & 0 & 0 & S_{2323} & 0 & 0 \\ 0 & 0 & 0 & 0 & S_{1212} & 0 \\ 0 & 0 & 0 & 0 & 0 & S_{1212} \end{pmatrix} \tag{9}\]Note that $\mathbb{S}$ is not symmetric: the off-diagonal blocks satisfy $S_{1122} \neq S_{2211}$ whenever $e \neq 1$. Rows 2 and 3 of the upper-left $3\times3$ block are related by the 2–3 swap ($S_{2222} \leftrightarrow S_{2233}$), reflecting the transverse isotropy of the 2–3 plane.
Relation to the Hill tensor
A closely related tensor introduced by Hill [8] is
\[\mathbb{P} = \mathbb{S}\,\mathbb{C}_m^{-1}, \qquad \text{equivalently} \quad \mathbb{S} = \mathbb{P}\,\mathbb{C}_m.\]While the Eshelby tensor lacks major symmetry ($S_{ijkl} \neq S_{klij}$), the Hill tensor possesses it ($P_{ijkl} = P_{klij}$), making it the natural choice for variational bounds and self-consistent schemes [8].
Practical applications of eigenstrain
- Thermal residual stresses: $\boldsymbol{\varepsilon}^{\ast} = (\alpha_f - \alpha_m)\,\Delta T\,\mathbf{I}$ after cure.
- Phase transformations: martensitic transformation eigenstrain (dilatational + deviatoric) drives stress-induced plasticity.
- Precipitation hardening: misfitting precipitates in Al alloys interact with dislocations via eigenstrain fields.
- Residual stress reconstruction: invert measured strain fields to recover manufacturing residual stresses.
Limitations
- Infinite matrix: no boundary effects. Corrections needed for finite RVEs or thin films.
- Ellipsoidal inclusions: uniformity holds only for ellipsoids.
- Linear elasticity: nonlinear or viscoelastic phases require incremental extensions.
References
[1] 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.
[2] Mura, T. (1987). Micromechanics of Defects in Solids (2nd ed.). Martinus Nijhoff.
[3] Nemat-Nasser, S. & Hori, M. (1999). Micromechanics: Overall Properties of Heterogeneous Materials (2nd ed.). Elsevier.
[4] Clyne, T. W. & Withers, P. J. (1993). An Introduction to Metal Matrix Composites. Cambridge University Press.
[5] Mori, T. & Tanaka, K. (1973). Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica, 21(5), 571–574.
[6] Benveniste, Y. (1987). A new approach to the application of Mori-Tanaka’s theory in composite materials. Mechanics of Materials, 6(2), 147–157.
[7] Eshelby, J. D. (1959). The elastic field outside an ellipsoidal inclusion. Proceedings of the Royal Society A, 252(1271), 561–569.
[8] Parnell, W. J. (2016). The Hill and Eshelby tensors for ellipsoidal inhomogeneities in the Newtonian potential problem and linear elastostatics. Journal of Elasticity, 125(1), 87–142.