On the nonlinear Eshelby inclusion problem and its isomorphic growth limit

Abstract

In the late 1950s, Eshelby’s linear solutions for the deformation field inside an ellipsoidal inclusion and, subsequently, the infinite matrix in which it is embedded were published. The solutions’ ability to capture the behavior of an orthotropically symmetric shaped inclusion made it invaluable in efforts to understand the behavior of defects within, and the micromechanics of, metals and other stiff materials throughout the rest of the 20th century. Over half a century later, we wish to understand the analogous effects of microstructure on the behavior of soft materials, both organic and synthetic, but in order to do so, we must venture beyond the linear limit, far into the nonlinear regime. However, no solutions to these analogous problems currently exist for non-spherical inclusions. In this work, we present an accurate semi-inverse solution for the elastic field in an isotropically growing spheroidal inclusion embedded in an infinite matrix, both made of the same incompressible neo-Hookean material. We also investigate the behavior of such an inclusion as it grows infinitely large, demonstrating the existence of a non-spherical asymptotic shape and an associated asymptotic pressure. We call this the isomorphic limit, and the associated pressure the isomorphic pressure.

1. Introduction

Many of the most pressing problems in the physics of soft solids are inextricably linked to the physics of incompatibility [1, 2]. Growing biological objects such as tumors or biofilms are often constrained within, or against, other solid masses [3–7]. So much of the behavior of biological materials is dependent on the micromechanical incompatibilities and interactions between the cells and other microstructures that exist at a scale above the molecular but beneath our own [8–10]. Even in synthetic materials, differences in rates of curing can lead to incompatibility-driven stress concentrations within critical components and the resulting deformations can be large. Understanding these concentrations is especially important as we lean more heavily into the use of resins and other polymers as engineering materials [11–13].

The problem of an elastic inclusion embedded in an infinite matrix, as first formulated by John D. Eshelby in 1957, has served as a fundamental and instrumental example of incompatibility in the physics of solids [14, 15]. Eshelby set out to model any number of physical problems where “...the uniformity of an [infinite] elastic medium is disturbed by a region within it which has changed its form ...”. He called the changing region the inclusion and the rest of the infinite body the matrix. Ultimately he was able to solve this problem for an arbitrarily shaped inclusion and to find a particularly elegant solution for the case where this region took on the shape of an ellipsoid. “Fortunately”, remarks Eshelby in 1957, “the general ellipsoid is versatile enough to cover a wide variety of particular cases.”

The power and versatility of Eshelby’s solution is twofold: its ability to handle the case of an inclusion of an orthotropically symmetric shape and the fact that it is fully analytical. As it turns out, Eshelby was quite right. Over the years, a truly staggering number of examples were able to be well approximated as ellipsoidal inclusions (as they are convex, vaguely round, and longer in some directions than others). His exact solution has been instrumental in understanding the fields which develop around defects in stiff materials like metal, and through homogenization methods (which utilized the analytical nature of the solution) has laid the foundation for the development of micromechanical models for metals and stiff composite materials [16–18]. Extensions of his 1957 and 1959 works to other shapes and materials have been applied to better understand a wide variety of materials including metamaterials and piezoelectric materials [19–22].

One might imagine that in trying to understand the physics of incompatibility in soft solids, Eshelby’s solution to the inclusion problem may serve as a strong basis in doing so. However, although exact, Eshelby’s solution, and the thought experiment he utilized to attain it, relies on the geometric and constitutive linearity inherent to small deformations. Unfortunately, in the soft systems we wish to understand in the 21st century, we do not have the luxury of this linearity.

Generally speaking, exact analytical solutions to nonlinear problems, especially those lacking strong symmetries, are sparse at best. The inclusion problem is no exception. The extent to which exact solutions to the nonlinear inclusion problem have been found (largely by Yavari [23], Golgoon [24], and Goriely [25]) for commonly used constitutive relations is limited to specific spherical and circular cylindrical inclusions changing their shapes radially, axisymmetrically, torsionally, or azimuthally such that they remain spherical or circular cylindrical in their deformed shape. The problem of an elliptical cavity or rigid inclusion (both limits of the problems known as inhomogeneities) has also been solved in two dimensions using a specific two-dimensional constitutive relation called a harmonic material (whose physical applications are greatly limited) [26, 27]. However, no exact solutions exist for a general nonlinear ellipsoidal inclusion (or inhomogeneity) inside an infinite matrix.

This has led many to turn to numerical methods such as finite element analysis (FEA) as a first-line solution to approximately solve nonlinear problems. For instance, in 2000, Diani and Parks [28] examined the effect of geometric nonlinearity on the ellipsoidal inclusion problem using FEA to solve the problem with a logarithmic strain energy function. However, there are a few downsides to using these numerical methods, especially for a problem like the inclusion problem. Eshelby’s original statement of the problem imagines that the inclusion is embedded within an infinite matrix. When solving the problem numerically one cannot have a truly infinite matrix. As the matrix is made larger and larger it should converge to the solution of an infinite matrix, however, this also greatly increases the computational cost of the simulation, often taking hours of simulation time to reach moderate levels of stretch, and days to weeks of simulation time when performing parameter sweeps. Additionally, because a numerical solution is not expressed analytically, it is difficult to prove any sort of asymptotic behavior of the system and any homogenized solutions tend to be even more resource intensive.

Another set of methods that can be used to find approximate (or exact) solutions to nonlinear problems are inverse or semi-inverse methods. In an inverse solution, a fully determined, kinematically admissible deformation field is assumed and then the resulting stress field is checked to see whether it satisfies the necessary field and boundary equations. In a semi-inverse solution (as first developed by St. Venant), a form for the deformation field is assumed involving some set of unknowns, and the unknowns are then used to satisfy the relevant field equations and boundary conditions of the system. If these equations can be exactly satisfied, then you have found the exact solution to the problem. However, even if the equations are not exactly satisfied, the closest approximate solution which satisfies the initial kinematic assumptions can be found by minimizing the total elastic energy of the system [29]. Inverse and semi-inverse methods have been used readily in finding fairly accurate approximate solutions to problems involving spherical and spheroidal voids in nonlinear elastic and plastic materials [13, 30–34].

In this work, we utilize a semi-inverse method to formulate an approximate solution to the problem of a nonlinear elastic, incompressible, isotropically growing ellipsoidal inclusion embedded within an infinite matrix made of the same material. We begin by examining Eshelby’s solution in order to formulate a key assumption about the “shape” of the nonlinear deformation: that there exists a set of ellipsoids which fully cover the infinite body, and which remain ellipsoids after deformation. We then use these kinematic assumptions, along with incompressibility, to simplify our description of the field. We find that the entire field can be described by two functions which describe this set of ellipsoids. We then return to Eshelby’s solution to select this set of ellipsoids, generalizing such that we retain two “knobs” which we can tune to minimize the total elastic energy of the system, a la the semi-inverse method.

We apply this formulation to the problems of an isotropically growing spheroidal inclusion and demonstrate the striking accuracy of the resulting solutions through comparison with FEM simulations and Eshelby’s solution deep into the nonlinear range. Finally, we analyze the behavior of this solution as the inclusion grows infinitely large and show that for the case of spheroidal inclusions there exists a non-spherical shape which the inclusion tends toward. We call this limit the isomorphic growth limit and demonstrate its properties. We find that in the isomorphic growth limit the pressure inside the inclusion approaches a limit higher than that of spherical cavitation and call this pressure limit the isomorphic pressure. We end the work with a discussion of the application of this solution to homogenization and opine on the implications of the “shape” of the resulting fields for this approximate nonlinear solution to a canonical problem.

Figure 1.

The infinitely large body $\mathcal{B}$ in the reference (bottom left), transformed (top center), and deformed (bottom right) configurations. In each configuration, the region occupied by the inclusion (denoted by the subscript “I”) and matrix (denoted by the subscript “M”) are indicated. The inclusion has the three semi-axes $A_j$ in the reference configuration, ${\lambda }^{*}{A}_j$ in the transformed configuration, and $a_j$ in the deformed configuration aligned with the right-handed Cartesian coordinate system ( ${\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3$ ). The red region in the transformed configuration denotes the incompatibility present between the inclusion and matrix. The mappings ${\mathbf{\chi }}^e$ and $\mathbf{\chi }$ and their gradients ${\mathbf{F}}^{*}$ , ${\mathbf{F}}^e$ , and $\mathbf{F}$ represent the transformation, elastic part of the deformation, and total deformation, respectively.

2. The inclusion problem in finite elasticity

Consider an infinitely large body, $\mathcal{B}$ , described by a set of material points $\mathbf{X}\in {\mathbb{R}}^3$ . We divide this body into two material sub-regions, the inclusion - ${\mathcal{B}}_I$ , and the matrix within which it is enclosed - ${\mathcal{B}}_M$ , such that ${\mathcal{B}}_I\cup {\mathcal{B}}_M\equiv \mathcal{B}$ and ${\mathcal{B}}_I\cap {\mathcal{B}}_M$ forms a closed and simply connected surface we call the inclusion–matrix interface.

Initially, in their stress-free state, the inclusion and the matrix occupy the respective regions ${\mathcal{D}}_I^0$ and ${\mathcal{D}}_M^0$ (Figure 1), and the inclusion–matrix interface is denoted by $\partial {\mathcal{D}}_I^0$ . The inclusion then undergoes a permanent, inelastic transformation, such as growth or a phase transition, which without the constraint imposed on it by the surrounding matrix, would cause an arbitrary homogeneous transformation stretch. However, since the matrix and inclusion constrain each other, they both deform to accommodate the transformation. The corresponding regions occupied by the deformed inclusion and matrix are denoted by ${\mathcal{D}}_I$ and ${\mathcal{D}}_M$ , respectively, and the inclusion–matrix interface is denoted by $\partial {\mathcal{D}}_I$ . Material points in the deformed configuration are thus mapped to their new, deformed coordinates $\mathbf{x}=\mathbf{\chi }\left(\mathbf{X}\right)$ .

In the context of finite elasticity, we restrict our attention to homogeneous, isotropic, incompressible, and hyperelastic materials and assume that the matrix and inclusion are perfectly bonded at their interface, namely that the deformation is continuous ( $\mathbf{\chi }\left(\mathbf{X}\right)\in C^0$ ). Under these assumptions, we seek a solution to the deformation field for a given transformation stretch.

Or, in the words of Eshelby, we now task ourselves with solving for the “... elastic state of inclusion and matrix.”

2.1. Kinematic framework

We define the deformation gradient as follows:

$$\mathbf{F}=\frac{\partial \mathbf{\chi }\left(\mathbf{X}\right)}{\partial \mathbf{X}}.$$

(1)

However, because the transformation of the inclusion alters its stress-free configuration, the elastic free energy of the solid will no longer be dependent on $\mathbf{F}$ , but instead on only the part of the deformation which forces the two bodies together elastically.

In order to isolate this part of the deformation, we introduce a motion ${\mathbf{x}}^{*}={\mathbf{\chi }}^{*}\left(\mathbf{X}\right)$ we call the transformation which maps material point $\mathbf{X}$ to the point ${\mathbf{x}}^{*}$ which it occupies in the transformed body ${\mathcal{D}}^{*}$ . In this transformed configuration, the inclusion occupies the region ${\mathcal{D}}_I^{*}$ and the matrix occupies the region ${\mathcal{D}}_M^{*}$ (as shown at the top of Figure 1). We note that the transformed configuration is incompatible, meaning that gaps and/or overlaps exist between ${\mathcal{D}}_I^{*}$ and ${\mathcal{D}}_M^{*}$ (as denoted by the shaded red region in Figure 1). This implies that ${\mathbf{\chi }}^{*}$ is not a bijective map.

We now decompose $\mathbf{\chi }$ into the transformation ${\mathbf{\chi }}^{*}$ and the elastic deformation $\mathbf{x}={\mathbf{\chi }}^e\left({\mathbf{x}}^{*}\right)$ as $\mathbf{x}={\mathbf{\chi }}^e\left({\mathbf{\chi }}^{*}\left(\mathbf{X}\right)\right)$ , separating the motion into two steps. The first transforms the inclusion’s stress-free configuration, and the second deforms the bodies to fit together. Consequently, we can decompose the deformation gradient as

$$\mathbf{F}={\mathbf{F}}^e{\mathbf{F}}^{*},$$

(2)

where ${\mathbf{F}}^e=\partial {\mathbf{\chi }}^e/\partial {\mathbf{x}}^{*}$ is the elastic part of the deformation gradient and ${\mathbf{F}}^{*}=\partial {\mathbf{\chi }}^{*}/\partial \mathbf{X}$ is the transformation gradient. ¹

Because both the inclusion and matrix are made of an elastically incompressible material, we may restrict ourselves to isochoric deformations. This implies that the elastic part of the deformation must obey the constraint

$$\underset{\Omega }{\int }(J^e-1)\mathrm{dv}=0,$$

(3)

for any arbitrary sub-body in the matrix or inclusion, where $J^e=\left({\mathbf{F}}^e\right)$ . Since equation (3) must hold for any arbitrary sub-body Ω, the integrand of equation (3) vanishes everywhere, implying $J^e=1,\forall \mathbf{X}$ . Furthermore, from equation (2), we find that incompressibility implies

$$J=J^{*},\forall \mathbf{X},$$

(4)

where $J=\det \left(\mathbf{F}\right)$ and $J^{*}=\det \left({\mathbf{F}}^{*}\right)$ . We note that because we assume that the transformation is uniform and occurs only inside the inclusion, ${\mathbf{F}}^{*}$ and $J^{*}$ are piecewise constant with ${\mathbf{F}}^{*}=\mathbf{1}$ and $J^{*}=1$ for all $\mathbf{X}\in {\mathcal{D}}_M^0$ , where 1 denotes the identity tensor.

3. Ellipsoidal inclusions undergoing isotropic transformations

Having described the general finite elastic inclusion problem, we now restrict our attention to ellipsoidal inclusions which undergo isotropic transformations. The semi-axes of an inclusion are taken as $(A_1,A_2,A_3)$ , which align with the directions of a right handed Cartesian system $({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3)$ , respectively, and we set the origin of the coordinate system at the geometric center of the inclusion. We describe the body, $\mathcal{B}$ , as a continuous set of ellipsoidal surfaces where the dimensionless radial parameter, $\mathrm{\Lambda }=\overline{\mathrm{\Lambda }}\left(\mathbf{X}\right)$ , is constant on a given ellipsoid, and the functions, ${\mathrm{\Phi }}_2={\hat{\mathrm{\Phi }}}_2\left(\mathrm{\Lambda }\right)$ and ${\mathrm{\Phi }}_3={\hat{\mathrm{\Phi }}}_3\left(\mathrm{\Lambda }\right)$ , are the corresponding aspect ratios, to write the implicit relationship

$$X_1^2+\frac{X_2^2}{{\mathrm{\Phi }}_2^2}+\frac{X_3^2}{{\mathrm{\Phi }}_3^2}=A_1^2{\mathrm{\Lambda }}^2,$$

(5)

where $X_j\equiv \mathbf{X}·{\mathbf{e}}_j$ . According to this definition, Λ corresponds to the normalized distance from the origin to the intercept of the ellipsoid with ${\mathbf{e}}_1$ , such that $\mathrm{\Lambda }=X_1/A_1$ . The inclusion–matrix interface is thus at $\mathrm{\Lambda }=1$ and we enforce that this set of ellipsoidal surfaces includes the inclusion–matrix interface by setting ${\hat{\mathrm{\Phi }}}_2\left(1\right)={\mathrm{\Phi }}_{02}=A_2/A_1$ and ${\hat{\mathrm{\Phi }}}_3\left(1\right)={\mathrm{\Phi }}_{03}=A_3/A_1$ . Consequently,

$$\mathrm{for}\{\begin{array}{cc}0\le \mathrm{\Lambda }\le 1,\hfill & \mathbf{X}\in {\mathcal{D}}_I^0\hfill \\ \multicolumn{1}{c}{\mathrm{\Lambda }=1,}& \mathbf{X}\in \partial {\mathcal{D}}_I^0\hfill \\ \multicolumn{1}{c}{\mathrm{\Lambda }\ge 1,}& \mathbf{X}\in {\mathcal{D}}_M^0.\hfill \end{array}$$

(6)

Provided the above definitions, we write the volume of the inclusion as

$$V^0=\frac{4\pi }{3}{A}_1{A}_2{A}_3=\frac{4\pi }{3}{\mathrm{\Phi }}_{02}{\mathrm{\Phi }}_{03}{A}_1^3.$$

(7)

In the following sections, we use the summation convention of Mura [16], which operates by the following rules:

1. Repeated lowercase indices are summed from 1 to 3.

2. Uppercase indices take on the same value as the corresponding lowercase indices but are not summed.

Without loss of generality, we also define ${\mathrm{\Phi }}_1={\hat{\mathrm{\Phi }}}_1\left(\mathrm{\Lambda }\right)\equiv 1$ , and consequently, ${\hat{\mathrm{\Phi }}}_1\left(1\right)={\mathrm{\Phi }}_{01}=A_1/A_1=1$ , for use in calculations using this summation convention.

3.1. The linear limit and the confocal ellipsoids

Before attempting to solve the nonlinear inclusion problem, it is instructive to re-examine the small strain result. In his linear solution, Eshelby [14, 15] defines the inclusion problem in terms of a stress-free infinitesimal transformation strain (eigenstrain), related to ${\mathbf{F}}^{*}$ in the above-described nonlinear formulation as follows:

$${ϵ}^{*}=\frac{1}{2}[{\mathbf{F}}^{*}+{\left({\mathbf{F}}^{*}\right)}^T]-1.$$

(8)

He showed that in the linear limit, when ${ϵ}^{*}$ (and by extension ${\mathbf{F}}^{*}$ ) is constant for all $\mathbf{X}\in {\mathcal{D}}_I^0$ , the deformation is uniform inside the inclusion. This implies that in the linear limit, $\mathbf{F}$ should be constant inside the inclusion, and the inclusion should remain ellipsoidal.

From Eshelby’s solution, as formulated in Eshelby [14] and Eshelby [15], Mura [16], for the case of an isotropic transformation strain of the form

$${ϵ}^{*}={ϵ}^{*}\mathbf{1},\forall \mathbf{X}\in {\mathcal{D}}_I^0,$$

(9)

it can be trivially shown that if we choose the set of ellipses confocal with $\partial {\mathcal{D}}_I^0$

$${\mathrm{\Phi }}_J={\hat{\mathrm{\Phi }}}_J^C\left(\mathrm{\Lambda }\right)=\frac{\sqrt{{\mathrm{\Phi }}_{0J}^2+({\mathrm{\Lambda }}^2-1)}}{\mathrm{\Lambda }},\forall \mathbf{X}\in {\mathcal{D}}_M^0,$$

(10)

then the deformation takes the form

$$x_j={\hat{f}}_J\left(\mathrm{\Lambda }\right)X_j,\forall \mathbf{X}\in {\mathcal{D}}_M^0,$$

(11)

where $x_j\equiv \mathbf{x}·{\mathbf{e}}_j$ is the Cartesian coordinate of $\mathbf{x}$ in the ${\mathbf{e}}_j$ direction and the functions ${\hat{f}}_J$ are determined via the equations of linear elasticity. The above equation implies that for isotropic transformation strains, each of these confocal ellipsoidal surfaces is transformed into a new ellipsoidal surface in the deformed configuration. ²

3.2. Matching the linear limit: ellipsoids to ellipsoids

From equations (8) and (9), we find that an isotropic transformation strain corresponds to an isotropic transformation gradient of the form

$${\mathbf{F}}^{*}={\lambda }^{*}\mathbf{1},\forall \mathbf{X}\in {\mathcal{D}}_I^0,$$

(12)

where ${\lambda }^{*}=1+{ϵ}^{*}$ is the transformation stretch. Notice also that this implies $J^{*}=J_I^{*}\equiv {\left({\lambda }^{*}\right)}^3$ for all $\mathbf{X}\in {\mathcal{D}}_I^0$ . Since Eshelby’s solution is exact in the linear limit, the kinematics of any exact solution to the nonlinear problem of an ellipsoidal inclusion under an isotropic transformation must reduce to ellipsoids which are transformed to ellipsoids in the linear limit.

Therefore, in order to find an approximate solution to the nonlinear ellipsoidal inclusion problem under an isotropic transformation, we choose a kinematically admissable set of deformations in which ellipsoids of constant Λ, including $\partial {\mathcal{D}}_I^0$ , are transformed to a new set of ellipsoids via the mapping

$$\mathbf{x}=\mathbf{\chi }\left(\mathbf{X}\right)=\{\frac{\alpha }{\mathrm{\Lambda }}{X}_1,\frac{{\phi }_2\alpha }{{\mathrm{\Phi }}_2\mathrm{\Lambda }}{X}_2,\frac{{\phi }_3\alpha }{{\mathrm{\Phi }}_3\mathrm{\Lambda }}{X}_3\},$$

(13)

where we define the field variables $\alpha =\hat{\alpha }\left(\mathrm{\Lambda }\right)$ , ${\phi }_2={\hat{\phi }}_2\left(\mathrm{\Lambda }\right)$ , and ${\phi }_3={\hat{\phi }}_3\left(\mathrm{\Lambda }\right)$ as the deformed dimensionless radial parameter and deformed aspect ratios, respectively. ³ The semi-axes of the deformed inclusion are taken as $(a_1,a_2,a_3)$ , which align with $({\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3)$ , respectively, and the origin of the coordinate system remains at the geometric center of the inclusion.

As such, $\alpha \left(1\right)=a_1/A_1$ , ${\hat{\phi }}_2\left(1\right)={\phi }_{02}\equiv a_2/a_1$ , and ${\hat{\phi }}_3\left(1\right)={\phi }_{03}\equiv a_3/a_1$ , respectively, since this set of ellipsoidal surfaces includes the inclusion–matrix interface (with aspect ratios ${\phi }_{02}$ and ${\phi }_{03}$ ). Note that we do not yet enforce a specific form of ${\mathrm{\Phi }}_j$ , only the values of the undeformed aspect ratios at the matrix–inclusion interface. The volume of ${\mathcal{D}}_I$ is given as

$$V^{*}=\frac{4\pi }{3}{a}_1{a}_2{a}_3=J_I^{*}{V}^0,$$

(14)

since the inclusion has grown isotropically by a factor of ${\lambda }^{*}$ and the elastic part of the deformation is isochoric. We also define ${\phi }_1={\hat{\phi }}_1\left(\mathrm{\Lambda }\right)\equiv 1$ , and consequently, ${\hat{\phi }}_1\left(1\right)={\phi }_{01}\equiv a_1/a_1=1$ for use in calculations using summation convention. In this way,

$$x_j=\frac{{\phi }_J}{{\mathrm{\Phi }}_J}\frac{\alpha }{\mathrm{\Lambda }}{X}_j.$$

(15)

By inserting equation (15) into equation (1), and remembering that Λ is a function of $\mathbf{X}$ , we find that $F_{\mathrm{jk}}\equiv \mathbf{F}:({\mathbf{e}}_j\otimes {\mathbf{e}}_k)$ is

$$F_{\mathrm{jk}}=\frac{\partial x_j}{\partial X_k}=\frac{{\phi }_J}{{\mathrm{\Phi }}_J}(\overline{\lambda }{\delta }_{\mathrm{jk}}+\frac{(\lambda +{\beta }_J)-\overline{\lambda }(1+{\gamma }_J)}{A^2{\mathrm{\Lambda }}^2\mathrm{\Gamma }}·\frac{X_j{X}_k}{{\mathrm{\Phi }}_K^2}),$$

(16)

where ${\delta }_{\mathrm{ij}}$ is the Kronecker delta,

$$\mathrm{\Gamma }=1+{\gamma }_N\frac{X_n{X}_n}{{\mathrm{\Phi }}_N^2},$$

(17)

$\overline{\lambda }=\alpha /\mathrm{\Lambda }$ , $\lambda ={\alpha }^{\text{'}}$ , ${\beta }_J={\phi }_J^{\text{'}}\alpha /{\phi }_J$ , and ${\gamma }_J={\mathrm{\Phi }}_J^{\text{'}}\mathrm{\Lambda }/{\mathrm{\Phi }}_J$ with ${•}^{\text{'}}$ denoting differentiation with respect to Λ. We can now find J by taking the determinant of $\mathbf{F}$ and substitute ${\mathrm{\Phi }}_1={\phi }_1=1$ . Setting $J=J^{*}$ equal to $J_I^{*}$ inside the inclusion and 1 inside the matrix yields

$$J=\frac{{\phi }_2{\phi }_3}{{\mathrm{\Phi }}_2{\mathrm{\Phi }}_3}{\overline{\lambda }}^2\left(\frac{\lambda +{\beta }_2\frac{X_2^2}{{\mathrm{\Phi }}_2^2}+{\beta }_3\frac{X_3^2}{{\mathrm{\Phi }}_3^2}}{1+{\gamma }_2\frac{X_2^2}{{\mathrm{\Phi }}_2^2}+{\gamma }_3\frac{X_3^2}{{\mathrm{\Phi }}_3^2}}\right)=\{\begin{array}{cc}{J}_I^{*},\hfill & \forall \mathbf{X}\in {\mathcal{D}}_I^0\hfill \\ \multicolumn{1}{c}{1,}& \forall \mathbf{X}\in {\mathcal{D}}_M^0\hfill \end{array}.$$

(18)

Since equation (4) must hold for all $\mathbf{X}\in {\mathcal{D}}^0$ (arbitrary choice of $X_2$ and $X_3$ ) we find three constraints due to incompressibility

$$\frac{{\phi }_2{\phi }_3}{{\mathrm{\Phi }}_2{\mathrm{\Phi }}_3}{\overline{\lambda }}^2\lambda =\{\begin{array}{cc}{J}_I^{*},\hfill & \forall \mathbf{X}\in {\mathcal{D}}_I^0\hfill \\ \multicolumn{1}{c}{1,}& \forall \mathbf{X}\in {\mathcal{D}}_M^0\hfill \end{array},$$

(19a)

$${\beta }_2=\lambda {\gamma }_2,$$

(19b)

$${\beta }_3=\lambda {\gamma }_3.$$

(19c)

For ease of calculations, it is also useful to employ the conservation of volume of each ellipsoidal subregion of $\mathcal{B}$ as

$${\phi }_2{\phi }_3{\alpha }^3=\{\begin{array}{cc}{J}_I^{*}{\mathrm{\Phi }}_2{\mathrm{\Phi }}_3{\mathrm{\Lambda }}^3,\hfill & \forall \mathbf{X}\in {\mathcal{D}}_I^0\hfill \\ \multicolumn{1}{c}{({\mathrm{\Phi }}_2{\mathrm{\Phi }}_3{\mathrm{\Lambda }}^3+\mathrm{\Delta }),}& \forall \mathbf{X}\in {\mathcal{D}}_M^0\hfill \end{array}$$

(20)

where

$$\mathrm{\Delta }=\frac{3}{4\pi A_1^3}(V^{*}-V^0)={\mathrm{\Phi }}_{02}{\mathrm{\Phi }}_{03}(J_I^{*}-1).$$

(21)

Inserting equation (19a) into equation (20) for all $\mathbf{X}\in {\mathcal{D}}_I^0$ yields the integrable ordinary differential equation ${\alpha }^{\text{'}}=\alpha /\mathrm{\Lambda }$ , and thus, we find

$$\hat{\alpha }\left(\mathrm{\Lambda }\right)=C_1\mathrm{\Lambda },\forall \mathbf{X}\in {\mathcal{D}}_I^0,$$

(22)

where $C_1$ is an arbitrary integration constant. This implies that the deformation must be uniform inside the inclusion in order to satisfy incompressibility, a result which matches Eshelby’s solution in the linear limit. Because the deformation is uniform within the inclusion, any choice of ${\hat{\mathrm{\Phi }}}_J$ ( ${\hat{\mathrm{\Phi }}}_J\left(1\right)={\mathrm{\Phi }}_{0J}$ ) is kinematically admissible, and thus, we choose the convenient ${\mathrm{\Phi }}_J={\mathrm{\Phi }}_{0J}$ and ${\phi }_J={\phi }_{0J}$ for all $\mathbf{X}\in {\mathcal{D}}_I^0$ . In this way, from equation (20), for the inclusion, we find

$$\hat{\alpha }\left(\mathrm{\Lambda }\right)={\lambda }^{*}{\left(\frac{{\mathrm{\Phi }}_{02}{\mathrm{\Phi }}_{03}}{{\phi }_{02}{\phi }_{03}}\right)}^{\frac{1}{3}}\mathrm{\Lambda },\forall \mathbf{X}\in {\mathcal{D}}_I^0.$$

(23)

Inserting equation (23) into equation (16) and then applying equation (2) yields the elastic deformation gradient

$$F_{\mathrm{jk}}^e=\frac{{\phi }_{0J}}{{\mathrm{\Phi }}_{0J}}{\left(\frac{{\mathrm{\Phi }}_{02}{\mathrm{\Phi }}_{03}}{{\phi }_{02}{\phi }_{03}}\right)}^{\frac{1}{3}}{\delta }_{\mathrm{jk}},\forall \mathbf{X}\in {\mathcal{D}}_I^0.$$

(24)

For the matrix, inserting equation (19a) into equation (20) and integrating yields

$$\ln \left(\alpha \right)=\int \frac{{\mathrm{\Phi }}_2{\mathrm{\Phi }}_3{\mathrm{\Lambda }}^2}{\mathrm{\Delta }+{\mathrm{\Phi }}_2{\mathrm{\Phi }}_3{\mathrm{\Lambda }}^3}d\mathrm{\Lambda },\forall \mathbf{X}\in {\mathcal{D}}_M^0,$$

(25)

and inserting equations (19b) and (19c) into equation (16) yields

$$F_{\mathrm{jk}}^e=\frac{{\phi }_J}{{\mathrm{\Phi }}_J}(\overline{\lambda }{\delta }_{\mathrm{jk}}+\frac{(\lambda -\overline{\lambda })(1+{\gamma }_J)}{A^2{\mathrm{\Lambda }}^2\mathrm{\Gamma }}\frac{X_j{X}_k}{{\mathrm{\Phi }}_K^2}),\forall \mathbf{X}\in {\mathcal{D}}_M^0.$$

(26)

We have now reduced the problem to the determination of the two undeformed aspect ratio functions ${\hat{\mathrm{\Phi }}}_2\left(\mathrm{\Lambda }\right)$ and ${\hat{\mathrm{\Phi }}}_3\left(\mathrm{\Lambda }\right)$ . We now task ourselves with determining forms of these functions which accurately approximate the true solution.

3.3. A generalized set of ellipses

In the following, we will apply this general framework for ellipsoids to the specific case of spheroidal inclusions ( $A_1=A_2$ ) growing isotropically. We will begin choosing a set of confocal ellipsoids ( ${\mathrm{\Phi }}_J={\hat{\mathrm{\Phi }}}_J^C$ ). In this case, the deformation will be fully constrained by incompressibility and we can find a fully analytical set of expressions that define the deformation. Although a deformation of the form in equation (10) has been used in the past, we will show that it is a poor approximation of the true solution at even moderate levels of volume growth. Instead, we will show that generalizing ${\mathrm{\Phi }}_J$ to the form

$${\mathrm{\Phi }}_J={\hat{\mathrm{\Phi }}}_J^G\left(\mathrm{\Lambda }\right)=\frac{\sqrt{{\mathrm{\Phi }}_{0J}^2+{\mathrm{\Phi }}_{\infty J}^2({\mathrm{\Lambda }}^2-1)}}{\mathrm{\Lambda }},\forall \mathbf{X}\in {\mathcal{D}}_M^0,$$

(27)

leads to a more accurate solution, especially at large deformations. This generalized form describes a set non-confocal ellipsoids (formed by scaling the confocal ellipses uniformly, but anisotropically) which include the inclusion–matrix interface, where the parameter ${\mathrm{\Phi }}_{\infty J}$ describes the limiting value of their aspect ratios in the far-field ( $\mathrm{\Lambda }\to \infty$ ). We note that ${\hat{\mathrm{\Phi }}}_J^G\left(\mathrm{\Lambda }\right)$ is a monotonic function for all ${\mathrm{\Phi }}_{0J}$ , ${\mathrm{\Phi }}_{\infty J}$ . An illustration of the differences between the aspect ratio functions ${\hat{\mathrm{\Phi }}}^C\left(\mathrm{\Lambda }\right)$ and ${\hat{\mathrm{\Phi }}}^G\left(\mathrm{\Lambda }\right)$ is shown in Figure 2. In all cases, to remain consistent with equation (5), we will take ${\mathrm{\Phi }}_{\infty 1}=1$ . This set of ellipsoids reduces to the confocal set of ellipsoids in the case where ${\mathrm{\Phi }}_{\infty J}=1$ . Thus, we establish that this generalization can properly match the exact solution in the linear limit. For ease of calculations, it can also be shown from equation (27) that for these generalized ellipsoids

$${\gamma }_J=\frac{{\mathrm{\Phi }}_{\infty J}^2}{{\mathrm{\Phi }}_J^2}-1.$$

(28)

We have now further reduced the problem to the determination of the two far field aspect ratio parameters ${\mathrm{\Phi }}_{\infty 2}$ and ${\mathrm{\Phi }}_{\infty 3}$ . However, ${\mathrm{\Phi }}_{\infty 2}$ and ${\mathrm{\Phi }}_{\infty 3}$ cannot be determined by incompressibility, and so we turn to energy minimization to find the most accurate approximate solution to the problem which satisfies our kinematic assumptions.

Figure 2.

(a) A visualization of the difference between a confocal set of ellipses, characterized by the aspect ratio function ${\hat{\mathrm{\Phi }}}^C\left(\mathrm{\Lambda }\right)$ , and the generalized set of ellipses, characterized by the aspect ratio function ${\hat{\mathrm{\Phi }}}^G\left(\mathrm{\Lambda }\right)$ with ${\mathrm{\Phi }}_{\infty }=2$ . The gray region represents an undeformed inclusion with ${\mathrm{\Phi }}_0=5$ . (b) The aspect ratio functions of (a), ${\hat{\mathrm{\Phi }}}^C\left(\mathrm{\Lambda }\right)$ and ${\hat{\mathrm{\Phi }}}^G\left(\mathrm{\Lambda }\right)$ , versus Λ.

3.4. The elastic free energy and energy minimization

We begin by defining an elastic free energy per unit transformed volume ${\overline{\psi }}^{*}\left({\mathbf{F}}^e\right)$ . For an incompressible hyperelastic medium, the free energy may be expressed as functions of the first and second invariants of the elastic left Cauchy-Green deformation tensor $\mathbf{B}={\mathbf{F}}^e{\left({\mathbf{F}}^e\right)}^T$ , as ${\hat{\psi }}^{*}(I_1,I_2)$ , where $I_1=\mathrm{tr}\left(\mathbf{B}\right)=F_{\mathrm{jk}}^e{F}_{\mathrm{jk}}^e$ and $I_2=({\mathrm{tr}\left(\mathbf{B}\right)}^2-\mathrm{tr}\left({\mathbf{B}}^2\right))/2$ . For this problem, however, it is most convenient to use the free energy per unit reference volume (since it can be integrated over the continuous ${\mathcal{D}}^0$ ), which can be expressed as

$${\hat{\psi }}^0(I_1,I_2,J^{*})=J^{*}\left({\hat{\psi }}^{*}(I_1,I_2)\right),$$

(29)

since $J^{*}$ represents the ratio of the volume of an element in the transformed configuration to its volume in the reference configuration. In this work, we will limit our analysis to a neo-Hookean finite-elastic medium with an elastic free energy density function of the form

$${\hat{\psi }}^{*}\left(I_1\right)=\frac{\mu }{2}(I_1-3),$$

(30)

where μ is the ground state shear modulus of the material. ⁴ The total elastic free energy in the system for a neo-Hookean inclusion-matrix system can thus be given from equations (29) and (30), remembering that $J^{*}=1$ for all $\mathbf{X}\in {\mathcal{D}}_M^0$ , as

$$\mathrm{\Psi }=\underset{{\mathcal{D}}^0}{\int }{\psi }^0{dv}^0=\stackrel{{\mathrm{\Psi }}_M}{\overbrace{\frac{\mu }{2}\underset{{\mathcal{D}}_M^0}{\int }(I_1-3){dv}^0}}+\stackrel{{\mathrm{\Psi }}_I}{\overbrace{J_I^{*}\frac{\mu }{2}\underset{{\mathcal{D}}_I^0}{\int }(I_1-3){dv}^0}},$$

(31)

where ${\mathrm{\Psi }}_M$ and ${\mathrm{\Psi }}_I$ are the total elastic free energy in the inclusion and matrix, respectively.

Choosing ${\mathrm{\Phi }}_J={\hat{\mathrm{\Phi }}}_J^G\left(\mathrm{\Lambda }\right)$ and non-dimensionalizing equation (31) yields a non-dimensional elastic free energy of the form

$$\frac{2\mathrm{\Psi }}{\mu V^0}=\overline{\mathrm{\Psi }}({\lambda }^{*},{\mathrm{\Phi }}_{02},{\mathrm{\Phi }}_{\infty 2},{\mathrm{\Phi }}_{03},{\mathrm{\Phi }}_{\infty 3})={\overline{\mathrm{\Psi }}}_M+{\overline{\mathrm{\Psi }}}_I.$$

(32)

where ${\overline{\mathrm{\Psi }}}_M$ and ${\overline{\mathrm{\Psi }}}_I$ are the non-dimensional elastic-free energies of the matrix and inclusion, respectively. The closest approximate solution for an ellipsoidal inclusion of initial aspect ratios ${\mathrm{\Phi }}_{0J}$ with an isotropic eigenstretch ${\lambda }^{*}$ , satisfying our kinematic assumptions, can thus be found via the energy minimization conditions

$$\frac{\partial }{\partial {\mathrm{\Phi }}_{\infty 2}}\overline{\mathrm{\Psi }}({\lambda }^{*},{\mathrm{\Phi }}_{02},{\mathrm{\Phi }}_{\infty 2},{\mathrm{\Phi }}_{03},{\mathrm{\Phi }}_{\infty 3})=0\mathrm{and}\frac{\partial }{\partial {\mathrm{\Phi }}_{\infty 3}}\overline{\mathrm{\Psi }}({\lambda }^{*},{\mathrm{\Phi }}_{02},{\mathrm{\Phi }}_{\infty 2},{\mathrm{\Phi }}_{03},{\mathrm{\Phi }}_{\infty 3})=0.$$

(33)

3.5. Stress state inside the inclusion

For an incompressible neo-Hookean material, the Cauchy stress tensor $\mathbf{T}$ can be calculated as

$$\mathbf{T}=\mu {\mathbf{B}}^d-p\mathbf{1},$$

(34)

where ${·}^d$ represents the deviatoric part of a tensor, defined as ${\mathbf{A}}^d\equiv \mathbf{A}-(\mathrm{tr}\left(\mathbf{A}\right)/3)·\mathbf{1}$ , and p is an arbitrary hydrostatic pressure field associated with the incompressibility constraint. From equation (22) we can see that, given our kinematic assumptions, the deformation is uniform inside the inclusion. This, along with mechanical equilibrium ( $\mathrm{div}\left(\mathbf{T}\right)=\mathbf{0}$ ) also implies that the stress $\mathbf{T}$ is constant inside the inclusion and can be determined by solving for the spatially constant inclusion pressure $p=p_I$ .

By equating the change in internal energy of the matrix with the mechanical work done on it by the applied tractions at the inclusion–matrix interface, we can write the following configurational force balance:

$$\frac{\partial {\mathrm{\Psi }}_M}{\partial J_I^{*}}=-{\int }_{\partial {\mathcal{D}}_I}(\mathbf{Tn}·\frac{\partial \mathbf{x}}{\partial J_I^{*}})\mathrm{da}.$$

(35)

Remembering that $p=p_I$ is constant along $\partial {\mathcal{D}}_I$ and noticing that from mass conservation

$${\int }_{\partial {\mathcal{D}}_I}(\mathbf{n}·\frac{\partial \mathbf{x}}{\partial J_I^{*}})\mathrm{da}={\int }_{{\mathcal{D}}_I}\mathrm{div}\left(\frac{\partial \mathbf{x}}{\partial J_I^{*}}\right)\mathrm{dv}=\frac{\partial V^{*}}{\partial J_I^{*}}=V^0,$$

(36)

inserting equation (34) into equation (35) yields

$$p_I=\frac{1}{V^0}[\frac{\partial {\mathrm{\Psi }}_M}{\partial J_I^{*}}+\mu {\int }_{\partial {\mathcal{D}}_I}({\mathbf{B}}^d\mathbf{n}·\frac{\partial \mathbf{x}}{\partial J_I^{*}})\mathrm{da}].$$

(37)

3.6. Spheroids ( $A_1=A_2$ )

The following is a summary of the equations which can be used to solve for the field in a spheroidal inclusion and the surrounding matrix during isotropic transformations, given our kinematic assumptions and choosing the generalized version of the undeformed aspect ratio function (27). Detailed derivations of these equations can be found in Appendix 1.

For the case of a spheroid, from axisymmetry around ${\mathbf{e}}_3$ , we assume that ${\hat{\mathrm{\Phi }}}_2\left(\mathrm{\Lambda }\right)={\hat{\mathrm{\Phi }}}_1\left(\mathrm{\Lambda }\right)\equiv 1$ and ${\hat{\phi }}_2\left(\mathrm{\Lambda }\right)={\hat{\phi }}_1\left(\mathrm{\Lambda }\right)\equiv 1$ . As such, we are left with a single aspect ratio ${\mathrm{\Phi }}_3$ in the reference configuration and ${\phi }_3$ in the deformed configuration. For simplicity in the section on spheroids, we will drop the subscript “3” in the following section. That is to say that $({\mathrm{\Phi }}_3,{\mathrm{\Phi }}_{03},{\mathrm{\Phi }}_{\infty 3},{\phi }_3,{\phi }_{03},{\gamma }_3)\equiv (\mathrm{\Phi },{\mathrm{\Phi }}_0,{\mathrm{\Phi }}_{\infty },\phi ,{\phi }_0,\gamma )$ .

The deformation $\mathbf{\chi }\left(\mathbf{X}\right)$ can be described via the following functions for the undeformed aspect ratio

$$\hat{\mathrm{\Phi }}\left(\mathrm{\Lambda }\right)=\{\begin{array}{cc}{\mathrm{\Phi }}_0,\hfill & 0<\mathrm{\Lambda }<1\hfill \\ \multicolumn{1}{c}{\frac{\sqrt{({\mathrm{\Phi }}_0^2-{\mathrm{\Phi }}_{\infty }^2)+{\mathrm{\Phi }}_{\infty }^2{\mathrm{\Lambda }}^2}}{\mathrm{\Lambda }},}& \mathrm{\Lambda }>1,\hfill \end{array}$$

(38a)

and the corresponding deformed radial parameter α reads

$$\hat{\alpha }\left(\mathrm{\Lambda }\right)=\{\begin{array}{cc}{\lambda }^{*}{\left(\frac{{\mathrm{\Phi }}_0}{{\phi }_0}\right)}^{\frac{1}{3}}\mathrm{\Lambda },\hfill & 0<\mathrm{\Lambda }<1\hfill \\ \multicolumn{1}{c}{\frac{1}{{\mathrm{\Phi }}_{\infty }}\underset{R}{\Pi }{(\mathrm{\Phi }\mathrm{\Lambda }-R)}^{\frac{R^2}{3R^2-\delta }},}& \mathrm{\Lambda }>1,\hfill \end{array}$$

(38b)

where $\delta =({\mathrm{\Phi }}_0^2-{\mathrm{\Phi }}_{\infty }^2)$ and each R is a root of the depressed cubic polynomial $(R^3-\delta R+{\mathrm{\Phi }}_{\infty }^2\mathrm{\Delta }=0)$ . Note that we have not decided a priori whether or not the spheroidal inclusion is prolate ( $A_3>A_1=A_2$ ) or oblate ( $A_3<A_1=A_2$ ). Our formulation is agnostic to this choice, changing only the number of real R values. In either case, α will be entirely real. In the rest of this section, we find it convenient to define ${\alpha }_0\equiv \hat{\alpha }\left(1\right)$ and to use the relation

$$J_I^{*}={\left({\lambda }^{*}\right)}^3=\frac{{\phi }_0{\alpha }_0^3}{{\mathrm{\Phi }}_0},$$

(38c)

to calculate ${\phi }_0$ where necessary.

For an inclusion made of a neo-Hookean material, the total dimensionless elastic free energy of the inclusion–matrix system can be given as

$$\begin{array}{c}\overline{\mathrm{\Psi }}=\underset{1}{\overset{\infty }{\int }}\frac{\mathrm{\Phi }}{{\mathrm{\Phi }}_0}{\mathrm{\Lambda }}^2[2\left({\lambda }^2-{\overline{\lambda }}^2\right)\left(1-\frac{1}{{\lambda }^2{\overline{\lambda }}^4}\right)+\left(2+\frac{{\mathrm{\Phi }}_{\infty }^2}{{\mathrm{\Phi }}^2}\right)\left(2{\overline{\lambda }}^2+\frac{1}{{\overline{\lambda }}^4}-3\right)\\ \underset{{\overline{\mathrm{\Psi }}}_M}{\underbrace{+\frac{\left({\mathrm{\Phi }}_{\infty }^2-1\right)}{{\mathrm{\Phi }}^2}\left(\frac{{\mathrm{\Phi }}_{\infty }^2}{{\lambda }^2{\overline{\lambda }}^4}-1\right)\frac{{\left(\lambda -\overline{\lambda }\right)}^2}{{\gamma }^3}\left(2{\gamma }^2+3\gamma +\frac{3}{2}\frac{{\mathrm{\Phi }}_{\infty }^2}{{\mathrm{\Phi }}^2}\sqrt{-\gamma }\ln \left(\frac{1+\sqrt{-\gamma }}{1-\sqrt{-\gamma }}\right)\right)]d\mathrm{\Lambda }}}\\ +{\left({\lambda }^{\ast }\right)}^3\left[2{\left(\frac{{\mathrm{\Phi }}_0}{{\phi }_0}\right)}^{\frac{2}{3}}+{\left(\frac{{\phi }_0}{{\mathrm{\Phi }}_0}\right)}^{\frac{4}{3}}-3\right],\end{array}$$

(38d)

where ${\overline{\mathrm{\Psi }}}_M$ is the dimensionless elastic free energy of only the matrix. This expression is agnostic to whether the spheroidal inclusion is prolate or oblate (while $\sqrt{-\gamma }$ might be imaginary, $\overline{\mathrm{\Psi }}$ will always be real). ⁵

The parameter ${\mathrm{\Phi }}_{\infty }$ can now be determined for a given ${\lambda }^{*}$ and ${\mathrm{\Phi }}_0$ by satisfying equation (33) via the single condition

$$\frac{\partial }{\partial {\mathrm{\Phi }}_{\infty }}\overline{\mathrm{\Psi }}({\lambda }^{*},{\mathrm{\Phi }}_0,{\mathrm{\Phi }}_{\infty })=0.$$

(38e)

Finally, the pressure $p_I$ inside the inclusion can be determined by the following equation:

$$\frac{p_I}{\mu }=\frac{1}{2}\frac{\partial {\overline{\mathrm{\Psi }}}_M}{\partial J_I^{*}}+2{\left({\lambda }^{*}\right)}^2[{\left(\frac{{\mathrm{\Phi }}_0}{{\phi }_0}\right)}^{\frac{1}{3}}-{\left(\frac{{\phi }_0}{{\mathrm{\Phi }}_0}\right)}^{\frac{5}{3}}](\frac{\partial {\alpha }_0}{\partial J_I^{*}}-\frac{1}{3}\frac{{\alpha }_0}{J_I^{*}}).$$

(38f)

Thus, we have fully determined the stress state inside the inclusion.

Once ${\mathrm{\Phi }}_{\infty }$ is determined, the form of the deformation field is entirely analytic. However, because the authors were unable to analytically evaluate the integral necessary to calculate ${\overline{\mathrm{\Psi }}}_M$ , ${\mathrm{\Phi }}_{\infty }$ was determined numerically. Maybe the reader will have the luck or skill necessary to attain fully analytical results.

In order to validate the results attained for spheroids, we performed axisymmetric FEA of an isotropically growing inclusion using the open source software FEniCS, and the framework found in the work by Senthilnathan [35] (see section 5.5.1, Appendices C and Appendix E therein in his thesis) with the following, nearly incompressible neo-Hookean elastic free energy function

$${\psi }^{*}(I_1,J^e)=\frac{\mu }{2}(I_1-3-2\ln \left(J^e\right))+\frac{\kappa }{2}{\left(\ln \left(J^e\right)\right)}^2,$$

(39)

where we set $\kappa /\mu =1000$ to approximate incompressible behavior. The FEA results are compared with Eshelby’s linear solution [14, 15], a nonlinear solution using the confocal ellipses ${\hat{\mathrm{\Phi }}}^C$ (as used in the works by Gologanu et al. [31, 32], Avazmohammadi and Naghdabadi [33]), and the generalized solution formulated in this work using ${\hat{\mathrm{\Phi }}}^G$ .

Figure 3.

The deformed aspect ratio of the inclusion ${\phi }_0$ (a) and volume averaged dimensionless inclusion pressure ${\overline{p}}_I/\mu$ (b) versus the volume growth ratio ( $J_I^{*}-1$ ) of an incompressible and, neo-Hookean, isotropically growing prolate spheroidal inclusion–matrix system with undeformed aspect ratio ${\mathrm{\Phi }}_0=5$ . The dashed grey lines in (a) and (b) denote the isomorphic aspect ratio of the inclusion–matrix system ${\phi }_{\propto }=1.788$ and the dimensionless isomorphic inclusion pressure of the inclusion-matrix system $p_{\propto }/\mu =2.802$ , respectively (as derived in section 4). The numerical labels in (a) denote the volume growth ratios of the deformation magnitude plots in Figure 4.

Figure 4.

Full field plots of the normalized displacement magnitude $\left|\mathbf{u}\right|/A_1$ in the normalized reference coordinates $\mathbf{X}/A_1$ of an incompressible, neo-Hookean, isotropically growing prolate spheroidal inclusion–matrix system with undeformed aspect ratio ${\mathrm{\Phi }}_0=5$ . The volume growth ratios in each set of subplots are (a) $J_I^{*}=1.01$ , (b) $J_I^{*}=10$ , (c) $J_I^{*}=100$ , and (d) $J_I^{*}=26424$ . Note that the colorbar scale and zoom level vary between each set of plots.

Figure 5.

The deformed aspect ratio of the inclusion ${\phi }_0$ (a) and volume averaged dimensionless inclusion pressure ${\overline{p}}_I/\mu$ (b) versus the volume growth ratio ( $J_I^{*}-1$ ) of an incompressible and, neo-Hookean, isotropically growing prolate spheroidal inclusion–matrix system with undeformed aspect ratio ${\mathrm{\Phi }}_0=1/5$ . The dashed grey lines in (a) and (b) denote the isomorphic aspect ratio of the inclusion–matrix system ${\phi }_{\propto }=0.471$ and the dimensionless isomorphic inclusion pressure of the inclusion–matrix system $p_{\propto }/\mu =3.045$ , respectively (as derived in section 4). The numerical labels in (a) denote the volume growth ratios of the deformation magnitude plots in Figure 6.

Figure 6.

Full field plots of the normalized displacement magnitude $\left|\mathbf{u}\right|/A_1$ in the normalized reference coordinates $\mathbf{X}/A_1$ of an incompressible, neo-Hookean, isotropically growing oblate spheroidal inclusion-matrix system with undeformed aspect ratio ${\mathrm{\Phi }}_0=1/5$ . The volume growth ratios in each set of subplots are (a) $J_I^{*}=1.01$ , (b) $J_I^{*}=10$ , (c) $J_I^{*}=100$ , and (d) $J_I^{*}=45289$ . Note that the colorbar scale and zoom level vary between each set of plots.

The results for a prolate spheroid with ${\mathrm{\Phi }}_0=5$ , are shown in Figures 3 and 4, and the results for an oblate spheroid with ${\mathrm{\Phi }}_0=1/5$ are shown in Figures 5 and 6. The results for spheroids of aspect ratios ${\mathrm{\Phi }}_0=10$ , 2, $1/2$ , and $1/10$ are shown in Appendix 3, Figures 8–15. Figures 3 and 5 (a) present plots of the deformed aspect ratio of the inclusion ${\phi }_0$ as the inclusion grows (increasing $J_I^{*}$ ). Figures 3 and 5 (b) present plots of the volume averaged pressure inside the inclusion ${\overline{p}}_I$ nondimensionalized by the ground state shear modulus of the material μ as the inclusion grows. ⁶ In each plot, the solution presented in this work demonstrates striking accuracy and outperforms both the linear and confocal model both qualitatively and quantitatively.

Most importantly, note that both the FEA and the ${\hat{\mathrm{\Phi }}}^G$ models tend toward some non-spherical shape with an aspect ratio between ${\mathrm{\Phi }}_0$ and 1, and some pressure greater than the cavitation limit [36] of a neo-Hookean material of $5\mu /2$ . This behavior is not captured by any other existing model for the system. The Eshelby solution incorrectly inverts from a prolate to oblate (or vice versa) deformed inclusion shape and the inclusion pressure increases boundlessly with $J_I^{*}$ . The confocal solution, on the other hand, incorrectly approaches a spherical shape and the spherical cavitation pressure exactly, even for non spherical inclusions. Note that another useful generalization of the confocal set, which was introduced in the work by Zhang et al. [13], has the same limitations.

Furthermore, Figures 4 and 6, which plot the magnitude of the normalized displacement $\left|\mathbf{u}\right|/A_1=|\mathbf{x}-\mathbf{X}|/A_1$ in the normalized reference coordinates $\mathbf{X}/A_1$ at several orders of magnitude of $J_I^{*}$ , demonstrate the ability of our theory to capture the behavior of the FEA simulations in a way that the linear solution and confocal solution do not. For small $J_I^{*}$ , the models all agree. Energy minimization shows that in this linear limit, ${\mathrm{\Phi }}_{\infty }\to 1$ and thus the three fields are equivalent as they all must obey incompressibility and satisfy energy minimization. However, as $J_I^{*}$ increases, the linear solution has a much smaller field of influence on the body, as it can only scale with the transformation strain of the inclusion. The confocal solution, on the other hand, has a similarly sized field of influence of the deformation to that of the FEA but becomes spherical as the inclusion grows larger and larger. The shape of the field of influence of the deformation of FEA and solution presented in this work begin with a multi-lobed shape and approach a shape with aspect ratio ${\mathrm{\Phi }}_{\infty }\ne 1$ far from the inclusion. Similar observations were noted in the work by Zhang et al. [13]. This means that ${\mathrm{\Phi }}_{\infty }\ne 1$ dictates the “shape” of the deformation in the far field. Regardless of the value of ${\mathrm{\Phi }}_{\infty }$ , the deformation will vanish as $\mathrm{\Lambda }\to \infty$ , but the way in which it does so is dictated by this parameter. The reasons for this will be elucidated in the following section.

4. The isomorphic limit

In the following section, we investigate the behavior of an isotropically growing inclusion as its volume tends toward infinity, implying $J_I^{*}\to \infty$ and by extension $\mathrm{\Delta }\to \infty$ . We can equivalently imagine this as the limit where the inclusions initial volume goes to zero, i.e., $V^0/V^{*}\to 0$ . We find that for a neo-Hookean material an arbitrary ellipsoidal inclusion asymptotically approaches a non-spherical shape. Because in this limit the inclusion approaches this constant shape, we call this limit the isomorphic limit and the associated aspect ratios and maximum pressure (analogous to the well-established cavitation pressure) the isomorphic aspect ratios and isomorphic pressure, respectively.

4.1. The isomorphic aspect ratios

In this limit, because we assume ${\mathrm{\Phi }}_2$ and ${\mathrm{\Phi }}_3$ are finite, ${\mathrm{\Phi }}_2{\mathrm{\Phi }}_3{\mathrm{\Lambda }}^2$ is negligible unless $O\left({\mathrm{\Lambda }}^2\right)\ge O\left(\mathrm{\Delta }\right)$ . As $\mathrm{\Lambda }\to \infty$ , ${\mathrm{\Phi }}_J\to {\mathrm{\Phi }}_{\infty J}$ , thus, as $\mathrm{\Delta }\to \infty$ , equation (25) becomes

$$\underset{\mathrm{\Delta }\to \infty }{\lim }\ln \left(\alpha \right)=\underset{\mathrm{\Delta }\to \infty }{\lim }\int \frac{{\mathrm{\Phi }}_{\infty 2}{\mathrm{\Phi }}_{\infty 3}{\mathrm{\Lambda }}^2}{\mathrm{\Delta }+{\mathrm{\Phi }}_{\infty 2}{\mathrm{\Phi }}_{\infty 3}{\mathrm{\Lambda }}^3}d\mathrm{\Lambda },\forall \mathbf{X}\in {\mathcal{D}}_M^0.$$

(40)

Notice that this implies that the deformation field inside the matrix no longer depends on the shape of the region which the inclusion initially occupied. Thus, in this limit, our kinematic assumptions imply the deformation field $\underset{\mathrm{\Delta }\to \infty }{\lim }\mathbf{\chi }\left(\mathbf{X}\right)$ can be described via the following functions for the undeformed aspect ratios

$$\underset{\mathrm{\Delta }\to \infty }{\lim }{\hat{\mathrm{\Phi }}}_J\left(\mathrm{\Lambda }\right)=\{\begin{array}{cc}{\mathrm{\Phi }}_{0J},\hfill & 0<\mathrm{\Lambda }<1\hfill \\ \multicolumn{1}{c}{{\phi }_{\propto J},}& \mathrm{\Lambda }>1,\hfill \end{array}$$

(41a)

deformed aspect ratios

$$\underset{\mathrm{\Delta }\to \infty }{\lim }{\hat{\phi }}_J\left(\mathrm{\Lambda }\right)={\phi }_{\propto J},\forall \mathrm{\Lambda },$$

(41b)

and deformed radial parameter

$$\underset{\mathrm{\Delta }\to \infty }{\lim }\hat{\alpha }\left(\mathrm{\Lambda }\right)=\{\begin{array}{cc}{\lambda }^{*}{\left(\frac{{\mathrm{\Phi }}_{02}{\mathrm{\Phi }}_{03}}{{\phi }_{\propto 2}{\phi }_{\propto 3}}\right)}^{\frac{1}{3}}\mathrm{\Lambda },\hfill & 0<\mathrm{\Lambda }<1\hfill \\ \multicolumn{1}{c}{{({\mathrm{\Lambda }}^3+\frac{\mathrm{\Delta }}{{\phi }_{\propto 2}{\phi }_{\propto 3}})}^{\frac{1}{3}},}& \mathrm{\Lambda }>1.\hfill \end{array}$$

(41c)

Because the deformation field can be described by a constant deformed aspect ratio everywhere, we call this large growth limit the isomorphic limit and the corresponding aspect ratios ${\phi }_{\propto J}$ as the isomorphic aspect ratios. The value of each ${\phi }_{\propto J}$ can now be determined via energy minimization.

Here, we also find an explanation for why the confocal solution becomes spherical as $J_I^{*}\to \infty$ : it is constrained to do so by incompressibility since ${\phi }_{\propto }={\mathrm{\Phi }}_{\infty }=1$ . Therefore, the solution cannot capture the non-spherical asymptotic behavior seen in the simulations.

4.2. Energy minimization and the isomorphic pressure

Through the calculations in Appendix 2 it can be shown that for an arbitrary ellipsoidal inclusion, in the isomorphic limit, the total dimensionless free energy in the inclusion–matrix system is

$$\begin{array}{cc}\tilde{\mathrm{\Psi }}=\underset{\mathrm{\Delta }\to \infty }{\lim }\frac{2\mathrm{\Psi }}{\mu V^{*}}=\hfill & \underset{\underset{{\tilde{\mathrm{\Psi }}}_M}{︸}}{\frac{1}{5}[\frac{{({\phi }_{\propto 2}^2-1)}^2}{{\phi }_{\propto 2}^2}+\frac{{({\phi }_{\propto 3}^2-1)}^2}{{\phi }_{\propto 3}^2}+{(\frac{{\phi }_{\propto 2}}{{\phi }_{\propto 3}}-\frac{{\phi }_{\propto 3}}{{\phi }_{\propto 2}})}^2]+5}\hfill \\ \multicolumn{1}{c}{}& +[{\left(\frac{{\mathrm{\Phi }}_{02}{\mathrm{\Phi }}_{03}}{{\phi }_{\propto 2}{\phi }_{\propto 3}}\right)}^{\frac{2}{3}}+{\left(\frac{{\phi }_{\propto 2}^2{\mathrm{\Phi }}_{03}}{{\mathrm{\Phi }}_{02}^2{\phi }_{\propto 3}}\right)}^{\frac{2}{3}}+{\left(\frac{{\mathrm{\Phi }}_{02}{\phi }_{\propto 3}^2}{{\phi }_{\propto 2}{\mathrm{\Phi }}_{03}^2}\right)}^{\frac{2}{3}}-3],\hfill \end{array}$$

(42)

where ${\tilde{\mathrm{\Psi }}}_M$ is the total dimensionless free energy inside the matrix. Notice that by normalizing by the grown volume, the formulation for an arbitrary ellipsoidal inclusion becomes scale invariant and no longer depends on $J_I^{*}$ , ${\lambda }^{*}$ , Δ, or any other related quantity. In order to identify the values of ${\phi }_{\propto 2}$ and ${\phi }_{\propto 3}$ which minimize the elastic free energy, we find the fixed point where $\partial \tilde{\mathrm{\Psi }}/\partial {\phi }_{\propto 2}=\partial \tilde{\mathrm{\Psi }}/\partial {\phi }_{\propto 3}=0$ . In Appendix 3, it is further shown that as long as ${\mathrm{\Phi }}_{02}\ne 0$ or ${\mathrm{\Phi }}_{03}\ne 0$ , there exists a non-spherical shape which minimizes the total free energy of the system.

Any exact solution must have lower free energy than our approximate system, and since ${\phi }_{\propto 2}={\phi }_{\propto 3}=1$ is not the minimum energy solution for even our approximate form of deformation, the true solution as $J_I^{*}\to \infty$ must be different from that of isotropic spherical expansion.

As the inclusion grows, in a neo-Hookean material, it will also reach a maximum pressure (analogous to a cavitation pressure). We call this maximum inclusion pressure $p_{\propto }=\underset{\mathrm{\Delta }\to \infty }{\lim }{p}_I$ the isomorphic pressure and for an arbitrary ellipsoidal inclusion it can be found as follows:

$$\frac{p_{\propto }}{\mu }=\frac{{\tilde{\mathrm{\Psi }}}_M}{2}=\frac{5}{2}+\frac{1}{10}[\frac{{({\phi }_{\propto 2}^2-1)}^2}{{\phi }_{\propto 2}^2}+\frac{{({\phi }_{\propto 3}^2-1)}^2}{{\phi }_{\propto 3}^2}+{(\frac{{\phi }_{\propto 2}}{{\phi }_{\propto 3}}-\frac{{\phi }_{\propto 3}}{{\phi }_{\propto 2}})}^2].$$

(43)

In the case of a spherical inclusion ( ${\phi }_{\propto 2}={\phi }_{\propto 3}=1$ ), the isomorphic pressure reduces to the cavitation pressure $5\mu /2$ . And in any other case, the isomorphic pressure is higher than the cavitation pressure. This means that a growing ellipsoidal inclusion must be able to sustain a higher pressure than the cavitation pressure to be able to grow indefinitely large, even in a neo-Hookean material.

4.3. Isomorphic limit of a spheroidal inclusion

It can trivially be seen that for a spheroidal inclusion, described by only one isomorphic aspect ratio ${\phi }_{\propto }$ , equation (42) reduces to

$$\tilde{\mathrm{\Psi }}=\underset{\mathrm{\Delta }\to \infty }{\lim }\frac{2\mathrm{\Psi }}{\mu V^{*}}=\underset{{\tilde{\mathrm{\Psi }}}_M}{\frac{2}{5}\frac{{({\phi }_{\propto }^2-1)}^2}{{\phi }_{\propto }^2}+5}+2{\left(\frac{{\mathrm{\Phi }}_0}{{\phi }_{\propto }}\right)}^{\frac{2}{3}}+{\left(\frac{{\phi }_{\propto }}{{\mathrm{\Phi }}_0}\right)}^{\frac{4}{3}}-3,$$

(44)

and equation (43) reduces to

$$\frac{p_{\propto }}{\mu }=\frac{{\tilde{\mathrm{\Psi }}}_M}{2}=\frac{5}{2}+\frac{1}{5}\frac{{({\phi }_{\propto }^2-1)}^2}{{\phi }_{\propto }^2}.$$

(45)

By setting $\partial \tilde{\mathrm{\Psi }}/\partial {\phi }_{\propto }=0$ , we find that ${\phi }_{\propto }$ can be found as any real positive root of the polynomial

$$P\left({\phi }_{\propto }\right)=27{\mathrm{\Phi }}_0^4{({\phi }_{\propto }^4-1)}^3+125{\phi }_{\propto }^4{({\phi }_{\propto }^2-{\mathrm{\Phi }}_0^2)}^3.$$

(46)

From equation (46), it can be seen by inspection that for the spherical case ( ${\mathrm{\Phi }}_0=1$ ) that ${\phi }_{\propto }$ =1 is a real positive root of P, and that for any other positive value of ${\phi }_{\propto }$ , $P>0$ or $P<0$ . Thus, for a spherical inclusion, we find that the isomorphic limit is also always spherical.

From equation (46), it can be seen by inspection that for oblate spheroids ( ${\mathrm{\Phi }}_0<1$ ), if ${\phi }_{\propto }\le {\mathrm{\Phi }}_0$ then $P<0$ and if ${\phi }_{\propto }\ge 1$ then $P>0$ . Because P is a continuous function of ${\phi }_{\propto }$ , this means that P must have at least one positive real root in the range ${\phi }_{\propto }\in ({\mathrm{\Phi }}_0,1)$ . Additionally, the first derivative of P is given as follows:

$$P^{\text{'}}\left({\phi }_{\propto }\right)=324{\mathrm{\Phi }}_0^4{\phi }_{\propto }^3{({\phi }_{\propto }^4-1)}^2+500{\phi }_{\propto }^3{({\phi }_{\propto }^2-{\mathrm{\Phi }}_0^2)}^3+750{\phi }_{\propto }^5{({\phi }_{\propto }^2-{\mathrm{\Phi }}_0^2)}^2.$$

(47)

By inspection it is clear that for an oblate spheroid with ${\phi }_{\propto }\in ({\mathrm{\Phi }}_0,1)$ , P is monotonically increasing. Thus, for an oblate spheroid, we find that the isomorphic limit is also oblate with an aspect ratio between ${\mathrm{\Phi }}_0$ and 1.

Finally, from equation (46), it can be seen by in inspection that for prolate spheroids ( ${\mathrm{\Phi }}_0>1$ ), if ${\phi }_{\propto }\ge {\mathrm{\Phi }}_0$ , then $P>0$ and if ${\phi }_{\propto }\le 1$ then $P<0$ . Because P is a continuous function of ${\phi }_{\propto }$ , this means that P must have at least one positive real root in the range ${\phi }_{\propto }\in (1,{\mathrm{\Phi }}_0)$ . It can further be proven via Sturm’s theorem that there is only one real root within this range. Thus, for a prolate spheroid, we find that the isomorphic limit is also prolate with an aspect ratio between 1 and ${\mathrm{\Phi }}_0$ .

Here we make two important notes:

1. This analysis no longer depends on the assumption of an isotropically growing inclusion, only that the growth is homogeneous inside the inclusion. The isomorphic results depend on only the final, grown, shape of the inclusion.

2. In the limit where ${\mathrm{\Phi }}_0\to \infty$ , the problem becomes similar to that of a growing cylindrical cavity. For a cylindrical cavity, cavitation does not occur. We see similarly that there is no finite $p_{\propto }$ as ${\mathrm{\Phi }}_0\to \infty$ . This analysis gives some understanding of a connection between the two problems.

The analytical values of the isomorphic aspect ratios ${\phi }_{\propto }$ and associated isomorphic pressures $p_{\propto }$ plotted on Figures 3 and 5, agree with the limits of the results in the previous section and match quantitatively well with the FEA results. The analytical values of ${\phi }_{\propto }$ and $p_{\propto }$ versus initial aspect ratio ${\mathrm{\Phi }}_0$ are shown in Figure 7. The limiting values of these functions are also shown in Figure 7.

Figure 7.

(a) Isomorphic aspect ratio ${\phi }_{\propto }$ and (b) non-dimensional isomorphic inclusion pressure $p_{\propto }$ of a an incompressible, neo-Hookean, spheroidal inclusion–matrix system versus undeformed aspect ratio ${\mathrm{\Phi }}_0$ .

5. Conclusion

The accurate approximate solution in this work can be expressed in an analytical form (once ${\mathrm{\Phi }}_{\infty }$ is determined) and is able to capture the behavior of an isotropically growing neo-Hookean inclusion in a way that no previous theory has. Therefore, it may be useful in predicting the physics of incompatibility in synthetic and biological soft matter systems, such as the growth of bacterial biofilms and tumors, or the curing induced stress fields within light cured polymers. To attain this level of accuracy, however, the solution is limited to isotropic growth and relies on the ground state material properties of the matrix and inclusion being the same. In conjunction with homogenization methods, this and similarly derived solutions to related problems (fluid inclusions, inhomogeneities, stiffening inclusion/matrix materials, anisotropic growth, non-spheroidal geometries, external loads, etc.) may serve as a basis to better understand the micromechanics of soft materials in the future.

Furthermore, this work elucidates the asymptotic behavior of an ellipsoidal inclusion as it grows infinitely large. An analysis of this heretofore unexamined physics would be functionally impossible with only numerical simulations. This limit, associated with a higher asymptotic pressure than that of the commonly examined spherical cavitation limit, gives insight into how the shape of an inclusion affects the maximum stress it must be able to sustain in order to grow indefinitely.

The authors would like to emphasize a principal insight garnered from the methods used in this work: there is no reason to believe that the simplest way to express the solution to a nonlinear problem is the same at all levels of deformation. Expressing the solution presented in this work in a standard spheroidal coordinate system (convenient in the linear limit) would completely obscure the simplicity of the solution made possible by choosing a more generalized set of scaling spheroids. As the field looks toward understanding related problems on the incompatibility of soft systems, we reiterate that a carefully chosen arbitrary coordinate system may result in the discovery of more elegant and interpretable solutions.