2. Problem formulation

Consider a body which in its reference configuration occupies the region

(2.1)\begin{equation} \mathcal{B}=\{ \mathbf{x}\in \mathbb{R} ^3\,:\, \left\Vert \mathbf{x}\right\Vert<1\}, \end{equation}

where $\left \Vert \cdot \right \Vert$ denotes the Euclidean norm. Let $\mathbf {u}:\mathcal {B}\rightarrow \mathbb {R}^3$ denote a deformation of the body and let its deformation gradient be

(2.2)\begin{equation} \nabla\mathbf{u}(\mathbf{x})=\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}\mathbf{x}}(\mathbf{x}). \end{equation}

For smooth deformations, the requirement that $\mathbf {u}(\mathbf {x})$ is locally invertible and preserves orientation takes the form

(2.3)\begin{equation} \det \nabla\mathbf{u}(\mathbf{x})>0,\quad\mathbf{x}\in\mathcal{B}. \end{equation}

Let $W:\mathcal {B}\times M_+^{3\times 3}\rightarrow \mathbb {R}$ be the stored energy function of the material of the body where $M_+^{3\times 3}=\{ \mathbf {F}\in M^{3\times 3}\,|\,\det \mathbf {F}>0\}$ and $M^{3\times 3}$ denotes the space of real 3 by 3 matrices. Since we are interested in modelling large deformations, we assume for fixed $\mathbf {x}\in \mathcal {B}$ that the stored energy function $W$ satisfies that $W(\mathbf {x},\,\mathbf {F})\rightarrow \infty$ as either $\det \mathbf {F}\rightarrow 0^+$ or $\left \Vert \mathbf {F}\right \Vert \rightarrow \infty$.

We assume that the stored energy function, in units of energy per unit volume, describing the mechanical response of the body is given by

(2.4)\begin{equation} W(\mathbf{x},\mathbf{F})=\Phi(\mathbf{x},v_1,v_2,v_3),\quad\mathbf{F}\in M_+^{3\times 3},\quad \mathbf{x}\in\mathcal{B}, \end{equation}

for some function $\Phi :\mathcal {B}\times \mathbb {R}_+^3\rightarrow \mathbb {R}_+$ symmetric in its last three arguments, and where $v_1,\,v_2,\,v_3$ are the eigenvalues of $(\mathbf {F}^t\mathbf {F})^{1/2}$, known as the principal stretches. Note that for any fixed $\mathbf {x}$, the material response $W(\mathbf {x},\,\cdot )$ corresponds to an isotropic and frame-indifferent material.

We now restrict attention to the special case in which the deformation $\mathbf {u}(\cdot )$ is radially symmetric, so that

(2.5)\begin{equation} \mathbf{u}(\mathbf{x})=r(R)\,\frac{\mathbf{x}}{R},\quad\mathbf{x}\in \mathcal{B}, \end{equation}

for some scalar function $r(\cdot )$, where $R=\left \Vert \mathbf {x}\right \Vert$. In this case, one can easily check that

(2.6)\begin{equation} v_1=r^\prime(R),\quad v_2=v_3=\frac{r(R)}{R}. \end{equation}

In this case, we assume that the dependence of $\Phi$ on $\mathbf {x}$ in (2.4) is only on $R=\left \Vert \mathbf {x}\right \Vert$. We employ throughout this paper the following notation:

\[ \Phi(R,r(R))=\Phi\left(R,r^\prime(R),\frac{r(R)}{R}, \frac{r(R)}{R}\right), \]

with similar notation for any partial derivative $\Phi _{,i}=\frac {\partial \Phi }{\partial v _i}$ of $\Phi$. We shall drop the explicit dependence on $R$ above, whenever $\Phi$ does not depend explicitly on $R$.

The total stored energy functional due to internal mechanical and gravitational forces (see [Reference Calogero and Leonori6]) is given, up to a multiplicative constant of $4\pi$, by:

(2.7)\begin{equation} I(r)=I_{\mbox{mec}}(r)-I_{\mbox{pot}}(r), \end{equation}

where

(2.8)\begin{align} I_{\mbox{mec}}(r)& =\int_0^1 \Phi(R,r(R))\,R^2\,\mathrm{d} R, \end{align}

(2.9)\begin{align} I_{\mbox{pot}}(r)& =\int_0^1 \rho_0(R)\dfrac{M_R}{r(R)}\,R^2\,\mathrm{d} R, \end{align}

are the mechanical and gravitationalFootnote ³ potential energy functionals respectively. Here $\rho _0$ is the mass density of the body (mass per unit volume) in the reference configuration, and

\[ M_R=4\pi\int_0^R\rho_0(u)u^2\,\mathrm{d} u, \]

is the mass of the ball of radius $R$ in the reference configuration and centred at the origin. We assume that

(2.10)\begin{equation} k_0\le\rho_0(R)\le k_1,\quad 0\le R\le1, \end{equation}

for some positive constants $k_0$ and $k_1$.

Remark 2.2 We note that the stored energy function in (2.8) includes as a special case that in which

\[ \Phi\left(R,r^\prime(R),\frac{r(R)}{R},\frac{r(R)}{R}\right)= \rho_0(R)\phi\left(R,r^\prime(R),\frac{r(R)}{R},\frac{r(R)}{R}\right), \]

where now $\phi$ is the energy per unit mass, also called the specific energy density.

In accord with (2.3), we have the inequalities

(2.11)\begin{equation} r^\prime(R),\,\frac{r(R)}{R}>0,\quad 0< R<1. \end{equation}

In the boundary displacement problem, we require that the deformation $\mathbf {u}$ satisfies the boundary condition:

\[ \mathbf{u}(\mathbf{x})=\lambda\mathbf{x},\quad\mathbf{x}\in\partial\mathcal{B}, \]

where $\lambda >0$ is given. For (2.5), this reduces to:

(2.12)\begin{equation} r(1)=\lambda. \end{equation}

The (radial) boundary displacement problem now is to minimize the functional $I(\cdot )$ over the set

(2.13)\begin{equation} \mathcal{A}_\lambda=\big\{r\in W^{1,1}(0,1)\,|\,r(0)\ge0,\,r(1)=\lambda,\,r^\prime(R)>0\ \mbox{a.e.}, \, I_{\mbox{mec}}(r)<\infty\big\}. \end{equation}

Note that $\mathcal {A}_\lambda \ne \emptyset$ as $r_\lambda \in \mathcal {A}_\lambda$ where $r_\lambda (R)=\lambda R$.

To state the internal pressure problem, we let $P:[0,\,\infty )\rightarrow [0,\,\infty )$ be a continuous function satisfying (1.2), and let

(2.14)\begin{equation} H_P(s)=\int_0^{s}P(t)t^2\,\mathrm{d} t,\quad s\ge0. \end{equation}

The internal pressure problem is now to minimize the functional given by

(2.15)\begin{equation} I_P(r)=I(r)-H_P(r(0)), \end{equation}

over the set

(2.16)\begin{equation} \mathcal{A}=\left\{r\in W^{1,1}(0,1)\,|\,r(0)\ge0,\,r^\prime(R)>0\mbox{ a.e.} ,\, I_{\mbox{mec}}(r)<\infty\right\}. \end{equation}

Note that $\mathcal {A}\ne \emptyset$ as $\hat {r}(R)\equiv R$ belongs to $\mathcal {A}$. The term $H_P(r(0))$ in (2.15) represents the work done by the internal pressure to open up a cavity of size $r(0)$.

The following result from [Reference Sivaloganathan12, Lemma 4.1] will be useful for the analysis of this problem in the next section.

Lemma 2.4 let $P:[0,\,\infty )\rightarrow [0,\,\infty )$ be a continuous function satisfying (1.2). Then there exist constants $\alpha >0$ and $\beta$ such that

\[ \frac{1}{3}\,P_*\,\min\left\{s^3,\varepsilon^3\right\}\le H_P(s)\le\alpha\log(s)+\beta,\quad s\ge0. \]

3. Existence of minimizers

In this section, we show that the functionals $I(\cdot )$ in (2.7) and $I_P$ in (2.15) have minimizers over the sets $\mathcal {A}_\lambda$ in (2.13) and $\mathcal {A}$ in (2.16), respectively. We also give physically reasonable conditions for these minimizers to satisfy the Euler–Lagrange equations for the corresponding functionals. The proofs of the results in this section are adaptations of the corresponding ones in [Reference Ball3] due to the presence of the potential energy functional (2.9). We do emphasize that they are not direct consequence of those in [Reference Calogero and Leonori6] as these are for maps in Sobolev spaces $W^{1,p}$ with $p>3$ and thus they represent continuous deformations.

Throughout this section and the rest of the paper, we assume that the stored energy function $\Phi$ in (2.8) satisfies that

(3.1)\begin{equation} \Phi(R,v_1,v_2,v_3)\ge \phi(v_1)+\phi(v_2)+\phi(v_3)+h(v_1v_2v_3),\quad R\in[0,1], \end{equation}

where $\phi,\,h\,:\,(0,\,\infty )\rightarrow (0,\,\infty )$ are strictly convex and such that

1. H1: $\phi (v)\ge C v^\gamma$ for some positive constant $C$ and $1<\gamma <3$;

2. H2: $\dfrac {h(d)}{d}\rightarrow \infty,\,\quad \mbox {as }d\rightarrow \infty ;$

3. H3: $h(d)\ge Kd^{-s},\,\quad d>0,\,$ for some positive constant $K$ and $s\ge \gamma ^*={\gamma }/{\gamma -1}$.

If we let

\[ \delta_r(R)=r^\prime(R)\left(\frac{r(R)}{R}\right)^2, \]

then the specialization of equation (31) in [Reference Calogero and Leonori6] (which follows from the Hardy–Littlewood–Sobolev inequality) to radial map (2.5), together with (2.10) gives that

(3.2)\begin{equation} \left\vert \int_0^1 \rho_0(R)\dfrac{M_R}{r(R)}\,R^2\,\mathrm{d} R\right\vert\le C\left(\int_0^1\delta_r(R)^{{-}s}R^2\,\mathrm{d} R\right)^{{1}/{3s}}, \end{equation}

for some positive constant $C$ independent of $r\in \mathcal {A}_\lambda$ or $r\in \mathcal {A}$.

The following lemma will be needed to prove the next proposition.

Lemma 3.1 Let $u:[0,\,1]\rightarrow \mathbb {R}$ be a non-negative function. Then if $\phi$ satisfies H1, there exist constants $k_1>0$ and $k_2$ such that

\[ \int_0^1\phi(u(R))R^2\,\mathrm{d} R\ge k_1\int_0^1u(R)R^2\,\mathrm{d} R+k_2. \]

Moreover, $k_1$ can be chosen arbitrarily large by properly adjusting $k_2$.

Proof. It follows from H1 that $\phi (t)/t\rightarrow \infty$ as $t\rightarrow \infty$. Thus, for any given $k_1>0$, there exists $M>0$ such that $\phi (t)\ge k_1t$ whenever $t\ge M$. Hence,

\begin{align*} \int_0^1\phi(u(R))R^2\,\mathrm{d} R& =\int_{\left\{u< M\right\}}\phi(u(R))R^2\,\mathrm{d} R+ \int_{\left\{u\ge M\right\}}\phi(u(R))R^2\,\mathrm{d} R,\\ & \ge k_1\int_{\left\{u\ge M\right\}}u(R)R^2\,\mathrm{d} R,\\ & = k_1\int_0^1u(R)R^2\,\mathrm{d} R-k_1\int_{\left\{u< M\right\}}u(R)R^2\,\mathrm{d} R. \end{align*}

But

\[ \int_{\left\{u< M\right\}}u(R)R^2\,\mathrm{d} R\le M\int_{\left\{u< M\right\}}R^2\,\mathrm{d} R \le M\int_0^1R^2\,\mathrm{d} R=M/3. \]

Combining this with the previous inequality, we get that

\[ \int_0^1\phi(u(R))R^2\,\mathrm{d} R\ge k_1\int_0^1u(R)R^2\,\mathrm{d} R- \frac{k_1M}{3}. \]

The result follows upon setting $k_2=-({k_1M}/{3})$.

We now have the following:

Proof. Combining (2.10), (3.1) with H3, and (3.2) we get for some positive constants $K_1,\,K_2$ that

\[ I(r)\ge K_1\int_0^1\delta_r(R)^{{-}s}R^2\,\mathrm{d} R-K_2 \left(\int_0^1\delta_r(R)^{{-}s}R^2\,\mathrm{d} R\right)^{{1}/{3s}}, \]

for all $r\in \mathcal {A}_\lambda$. Since $s\ge \frac {3}{2}$, the function $g(x)=K_1x-K_2x^{\frac {1}{3s}}$ is bounded below for $x\ge 0$, and the result for $I(\cdot )$ follows.

For the functional $I_P(\cdot )$, in reference to (2.7) and (2.15), we can write that

\[ I_P(r)=\left[\frac{1}{2}I_{\mbox{mec}}(r)-I_{\mbox{pot}}(r)\right] +\left[\frac{1}{2}I_{\mbox{mec}}(r)-H_P(r(0))\right]. \]

The first bracketed term can be bounded below as in the first part of this proof. For the second bracket, we begin by noting that using (3.1) we get that

\[ I_{\mbox{mec}}(r)\ge \int_0^1(\phi(r'(R))+\phi(r(R)/R))R^2\,\mathrm{d} R. \]

Using lemma 3.1 twice, we can write

\begin{align*} \int_0^1\phi(r'(R))R^2\,\mathrm{d} R& \ge k_1\int_0^1r'(R)R^2\,\mathrm{d} R+k_2,\\ \int_0^1\phi(r(R)/R)R^2\,\mathrm{d} R& \ge 2k_1\int_0^1(r(R)/R)R^2\,\mathrm{d} R +\hat{k}_2, \end{align*}

where $k_1>0$ is chosen such that $k_1/2>\alpha$ with $\alpha$ is as in lemma 2.4. Combining these inequalities with lemma 2.4, we get that

\begin{align*} \frac{1}{2}I_{\mbox{mec}}(r)-H_P(r(0))& \ge \frac{1}{2} \left[k_1\int_0^1r'(R)R^2\,\mathrm{d} R+k_2\right]\\ & \quad +\frac{1}{2}\left[ 2k_1\int_0^1(r(R)/R)R^2\,\mathrm{d} R +\hat{k}_2\right] -\alpha\log(r(0))-\beta,\\ & \ge \frac{1}{2} \left[k_1\left(r(1)-2\int_0^1r(R)R\,\mathrm{d} R\right)+k_2\right]\\ & \quad +\frac{1}{2}\left[ 2k_1\int_0^1r(R)R\,\mathrm{d} R +\hat{k}_2\right] -\alpha r(1)-\beta,\\ & =(k_1/2-\alpha)r(1)+(k_2+\hat{k}_2)/2-\beta, \end{align*}

where we used, for the second inequality, that $\log (t)< t$ and that $r(0)\le r(1)$. The result for $I_P$ now follows.

Using this we can now establish the existence of minimizers for $I$ over $\mathcal {A}_\lambda$.

Theorem 3.3 Let the stored energy function $\Phi$ in (2.8) satisfy (3.1) with H1–H3. Then there exists $r_\lambda \in \mathcal {A}_\lambda$ such that

\[ I(r_\lambda)=\inf_{r\in\mathcal{A}_\lambda}I(r). \]

Proof. Since $\mathcal {A}_\lambda \ne \emptyset$, it follows from proposition 3.2 that $\inf _{r\in \mathcal {A}_\lambda }I(r)\in \mathbb {R}$. Let $(r_j)$ with $r_j\in \mathcal {A}_\lambda$ for all $j$ be an infimizing sequence, i.e.

\[ \inf_{r\in\mathcal{A}_\lambda}I(r)=\lim_{j\rightarrow{\infty}} I(r_j). \]

Since $(I(r_j))$ is bounded, it follows from the proof of proposition 3.2 that the sequence

(3.3)\begin{equation} \left(K_1\int_0^1\delta_{r_j}(R)^{{-}s}R^2\,\mathrm{d} R-K_2 \left(\int_0^1\delta_{r_j}(R)^{{-}s}R^2\,\mathrm{d} R\right)^{{1}/{3s}}\right), \end{equation}

is bounded. Hence, the sequence

\[ \left(\int_0^1\delta_{r_j}(R)^{{-}s}R^2\,\mathrm{d} R\right), \]

must be bounded as well, and thus from (3.2) that $(I_{\mbox {pot}}(r_j))$ is bounded. From this and the boundedness of $(I(r_j))$, we get that $(I_{\mbox {mec}}(r_j))$ is bounded.

From the boundedness of $(I_{\mbox {mec}}(r_j))$ and (3.1), we get that

\[ \left(\int_0^1h(\delta_{r_j}(R))R^2\,\mathrm{d} R\right), \]

is bounded. Let $\rho =R^3$ and $u_j(\rho )=r_j^3(\rho ^{1/3})$. It follows now that

\[ \dot{u}_j(\rho)=\frac{\mathrm{d}u_j}{\mathrm{d}\rho}(\rho)=\delta_{r_j}(\rho^{1/3}), \]

and that the sequence

\[ \left(\int_0^1h(\dot{u}_j(\rho))\,\mathrm{d} \rho\right), \]

is bounded. It follows now from H2 and De La Vallée–Poussin Criterion that for some subsequence $(\dot {u}_k)$ of $(\dot {u}_j)$, we have $\dot {u}_k\rightharpoonup w$ in $L^1(0,\,1)$ for some $w\in L^1(0,\,1)$, and that $(\dot {u}_j)$ is equi-integrable. Using H3 is easy to show that $w>0$ a.e. Letting

\[ u(\rho)=\lambda^3-\int_\rho^1w(s)\,\mathrm{d} s, \]

we get from the equi-integrability of $(\dot {u}_j)$ that $u_k\rightarrow u$ in $C[0,\,1]$. Thus, $r_k\rightarrow r_\lambda$ in $C[0,\,1]$ where $r_\lambda (R)=u(R^3)^{1/3}$. From these we can conclude $r_k\rightharpoonup r_\lambda$ in $W^{1,1}(\varepsilon,\,1)$ and that $\delta _{r_j}\rightharpoonup \delta _{r_\lambda }$ in $L^1(\varepsilon,\,1)$ for any $\varepsilon \in (0,\,1)$. By the weak lower semi-continuity properties of $I_{\mbox {mec}}(\cdot )$ (cf. [Reference Ball, Currie and Olver2]), we get that

\begin{align*} \int_\varepsilon^1\!\! \Phi\left(R,r_\lambda(R) \right)R^2\,\mathrm{d} R& \le\liminf_k \int_\varepsilon^1\!\! \Phi\left(R,r_k(R)\right]R^2\,\mathrm{d} R,\\ & \le \liminf_k I_{\mbox{mec}}(r_k)<\infty. \end{align*}

We get now from the Monotone Convergence Theorem and the arbitrariness of $\varepsilon$ that

(3.4)\begin{equation} I_{\mbox{mec}}(r_\lambda)\le \liminf_k I_{\mbox{mec}}(r_k). \end{equation}

This together with the facts that $r_\lambda (0)\ge 0$, $r_\lambda ^\prime (R)\ge 0$ a.e., and $r_\lambda (1)=\lambda$, shows that $r_\lambda \in \mathcal {A}_\lambda$.

To get that $r_\lambda$ is a minimizer of $I$ over $\mathcal {A}_\lambda$, we still have to deal with the potentials $(I_{\mbox {pot}}(r_k))$. First note that

(3.5)\begin{align} & \left\vert I_{\mbox{pot}}(r_k)-I_{\mbox{pot}}(r_\lambda)\right\vert\nonumber\\ & =\left\vert \int_0^1\dfrac{\rho_0(R)M_RR^2}{r_k(R)r_\lambda(R)}(r_\lambda(R)-r_k(R))\,\mathrm{d} R\right\vert,\nonumber\\ & \le\left\Vert r_\lambda-r_k\right\Vert_{C[0,1]} \left[\int_0^1\dfrac{\rho_0(R)M_RR^2}{r_k^2(R)}\mathrm{d} R\right]^{{1}/{2}} \left[\int_0^1\dfrac{\rho_0(R)M_RR^2}{r_\lambda^2(R)}\mathrm{d} R\right]^{{1}/{2}}, \end{align}

where in the last step we used a weighted Holder’s inequality with weight $\rho _0(R)M_RR^2$. We now show that each of the two integrals on the right-hand side of this inequality are bounded. Upon recalling (2.10), we can take $M_R\le CR^2$ for some constant $C$. Hence,

\begin{align*} \int_0^1\dfrac{\rho_0(R)M_RR^2}{r_k^2(R)}\mathrm{d} R& \le \mbox{(const)} \int_0^1\dfrac{R^4}{r_k^2(R)}\,\mathrm{d} R= \mbox{(const)} \int_0^1\dfrac{R^2}{\delta_{r_k}(R)}\,r_k^\prime(R)\mathrm{d} R,\\ & \le \mbox{(const)}\left[\int_0^1\dfrac{R^2}{\delta_{r_k}(R)^{\gamma^*}}\mathrm{d} R \right]^{{1}/{\gamma^*}}\left[\int_0^1R^2(r_k^\prime(R))^\gamma\mathrm{d} R \right]^{{1}/{\gamma}}. \end{align*}

However, by the weighted Holder’s inequality (with weight $R^2$),

\[ \left[\int_0^1\dfrac{R^2}{\delta_{r_k}(R)^{\gamma^*}}\mathrm{d} R\right]^{\frac{1}{\gamma^*}} \le \mbox{(const)} \left[\int_0^1\dfrac{R^2}{\delta_{r_k}(R)^{s}}\mathrm{d} R\right]^{\frac{1}{s}}, \]

with the sequence on the right-hand side bounded. From (3.1) and H1, it follows that

\[ C\int_0^1R^2(r_k^\prime(R))^\gamma\mathrm{d} R\le \int_0^1R^2\phi(r_k^\prime(R))\,\mathrm{d} R\le I_{\mbox{mec}}(r_k), \]

with the sequence on the right-hand side bounded. Combining these results, we can conclude that the sequence

\[ \left(\int_0^1\dfrac{\rho_0(R)M_RR^2}{r_k^2(R)}\mathrm{d} R\right), \]

is bounded. Using this in (3.5) we get that

\[ \left\vert I_{\mbox{pot}}(r_k)-I_{\mbox{pot}(r_\lambda)}\right\vert\le\mbox{(const)} \left\Vert r_\lambda-r_k\right\Vert_{C[0,1]} \rightarrow 0, \]

as $k\rightarrow \infty$, which together with (3.4) implies that

\[ I(r_\lambda)\le \liminf_k I(r_k)=\inf_{r\in\mathcal{A}_\lambda}I(r), \]

i.e. that $r_\lambda$ is a minimizer.

We now give the corresponding result for the existence of minimizers for $I_P$ over $\mathcal {A}$.

Theorem 3.4 Let the stored energy function $\Phi$ in (2.8) satisfy (3.1) with H1–H3 and $P(\cdot )$ be a continuous function satisfying (1.2). Then there exists $r_P\in \mathcal {A}$ such that

\[ I_P(r_P)=\inf_{r\in\mathcal{A}}I_P(r). \]

Proof. Since $\mathcal {A}\ne \emptyset$, it follows from proposition 3.2 that $\inf _{r\in \mathcal {A}}I_P(r)\in \mathbb {R}$. Let $(r_j)$ with $r_j\in \mathcal {A}$ for all $j$, be an infimizing sequence, i.e.

\[ \inf_{r\in\mathcal{A}}I_P(r)=\lim_{j\rightarrow{\infty}} I_P(r_j). \]

Since $(I_P(r_j))$ is bounded, it follows from the proof of proposition 3.2 that the sequence $(r_j(1))$ is bounded, and that

\[ \left(K_1\int_0^1\delta_{r_j}(R)^{{-}s}R^2\,\mathrm{d} R-K_2 \left(\int_0^1\delta_{r_j}(R)^{{-}s}R^2\,\mathrm{d} R\right)^{{1}/{3s}}\right), \]

is bounded as well. Hence, the sequence $(H(r_j(0)))$ must be bounded, and

\[ \left(\int_0^1\delta_{r_j}(R)^{{-}s}R^2\,\mathrm{d} R\right), \]

must be bounded as well, and thus from (3.2) that $(I_{\mbox {pot}}(r_j))$ is bounded. From this and the boundedness of $(I_P(r_j))$, we get that $(I_{\mbox {mec}}(r_j))$ is bounded. The rest of the proof is similar to that of theorem 3.3 upon setting

\[ u(\rho)=u(1)-\int_\rho^1w(s)\,\mathrm{d} s, \]

and $r_P(R)=u(R^3)^{1/3}$, and by the continuity of the functional $H(r(0))$ over $C[0,\,1]$.

For our next result, we shall need the following assumption on the stored energy function $\Phi$ (cf. [Reference Ball4]): there exist constants $M,\, \varepsilon _0\in (0,\,\infty )$ such that

(3.6)\begin{equation} \left\vert \frac{\partial\Phi}{\partial v_k}(R,\alpha_1v_1,\alpha_2v_2,\alpha_3v_3)v_k\right\vert\le M\left[\Phi(R,v_1,v_2,v_3)+1\right], \end{equation}

for all $R\in [0,\,1]$, $k=1,\,2,\,3$, and $\left \vert \alpha _i-1\right \vert <\varepsilon _0$ for $i=1,\,2,\,3$. Henceforth, we use the notation

(3.7)\begin{equation} T(R,r(R))=\dfrac{R^2}{r^2(R)}\,\Phi_{,1}(R,r(R)), \end{equation}

for the radial component of the Cauchy stress. The techniques in [Reference Ball3] can now be adapted to show the following result.

Theorem 3.5 Assume that the function $\Phi$ satisfies (3.6) and (3.1) with H1–H3. Let $P(\cdot )$ be a continuous function satisfying (1.2). If $r$ is any minimizer of $I$ over $\mathcal {A}_\lambda$ or $I_P$ over $\mathcal {A}$, we have that $r\in C^1(0,\,1]$, $r^\prime (R)>0$ for all $R\in (0,\,1]$, $R^{2}\Phi _{,1}(R,\,r(R))$ is $C^1(0,\,1]$, and

(3.8)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}R}\left[R^{2}\Phi_{,1}(R,r(R))\right]= 2R\Phi_{,2}(R,r(R))+ R^2\,\dfrac{\rho_0(R)M_R}{r^2(R)}, \quad0< R<1. \end{equation}

For the functional $I$, the minimizer $r$ satisfies (2.12) while for $I_P$ it must satisfy that

(3.9)\begin{equation} T(1,r(1))=0. \end{equation}

Moreover $r(0)\ge 0$ and if $r(0)>0$, then

(3.10a)\begin{align} \lim_{R\rightarrow{0^+}}T(R,r(R))& =0,\mbox{ for }I, \end{align}

(3.10b)\begin{align} \lim_{R\rightarrow{0^+}}[T(R,r(R))+P(r(R))]& =0,\mbox{ for }I_P. \end{align}

We note that using (3.8) together with definition (3.7), we get that

(3.11)\begin{align} \frac{\mathrm{d}T}{\mathrm{d}R}(R,r(R))& =2\dfrac{R^2}{r^3(R)}\left[\dfrac{r(R)}{R} \Phi_{,2}(R,r(R))-r^\prime(R)\Phi_{,1}(R,r(R))\right]\nonumber\\ & \quad +R^2\,\dfrac{\rho_0(R)M_R}{r^4(R)}. \end{align}

If $\rho (r)$ and $p(r)$ denote the density and pressure (due to gravity), respectively, on the deformed configuration, then

\[ \mathrm{d} p(r)=\dfrac{\rho(r)(4\pi r^2\mathrm{d} r)}{4\pi r^2}\,g_r=\rho(r)g_r\mathrm{d} r, \]

where $g_r$ is the gravity at radius $r$ in the deformed configuration. Since $g_r=\frac {M_r}{r^2}$ (setting the gravitational constant $G$ to one), we get

\[ \mathrm{d} p(r)=\rho(r)\,\frac{M_r}{r^2}\,\mathrm{d} r. \]

Here $M_r$ represents the mass within radius $r$ in the deformed configuration. In terms of the reference configuration, since $M_{r(R)}=M_R$ by conservation of mass, we have

\[ \mathrm{d} p_0(R)= \frac{\rho_0(R)}{r^\prime(R)(r(R)/R)^2}\,\frac{M_R}{r^2(R)}\,r^\prime(R)\, \mathrm{d} R =\rho_0(R)\frac{M_RR^2}{r^4(R)}\,\mathrm{d} R, \]

where $p_0(R)=p(r(R))$. Integrating this expression from $R$ to $1$, and using that $p_0(1)=0$, we get that

\[ p_0(R)={-}\int_R^1\rho_0(U)\frac{M_UU^2}{r^4(U)}\,\mathrm{d} U. \]

Using this we can write (3.11) now as

(3.12)\begin{align} \frac{\mathrm{d}}{\mathrm{d}R}\left[T(R,r(R))-p_0(R)\right]=2\dfrac{R^2}{r^3(R)}\left[\dfrac{r(R)}{R} \Phi_{,2}(R,r(R))-r^\prime(R)\Phi_{,1}(R,r(R))\right]. \end{align}

By the Baker–Ericksen inequality, the right-hand side of this equation is positive whenever $r^\prime (R)<{r(R)}/{R}$.

The material of the body $\mathcal {B}$ is homogeneous if $\rho _0(R)$ is constant, still denoted by $\rho _0$, and

(3.13)\begin{equation} \Phi(R,v_1,v_2,v_3)=\tilde{\Phi}(v_1,v_2,v_3). \end{equation}

In this case (3.8) reduces to

(3.14)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}R}\left[R^{2}\tilde{\Phi}_{,1}(r(R))\right]= 2R\tilde{\Phi}_{,2}(r(R))+ \frac{4\pi}{3}\, \rho_0^2\,\dfrac{R^5}{r^2(R)}, \quad0< R<1, \end{equation}

and (3.11) reduces to

(3.15)\begin{equation} \frac{\mathrm{d}\tilde{T}}{\mathrm{d}R}(r(R))=2\dfrac{R^2}{r^3(R)}\left[\dfrac{r(R)}{R} \tilde{\Phi}_{,2}(r(R))-r^\prime(R)\tilde{\Phi}_{,1}(r(R))\right] +\frac{4\pi}{3}\, \rho_0^2\,\dfrac{R^5}{r^4(R)}. \end{equation}

4. Cavitation

In this section, we show that for large values of $\lambda$ in the displacement problem, and for large values of $P_*$ in the internal pressure problem, the corresponding minimizer $r$ of $I(\cdot )$ over $\mathcal {A}_{\lambda }$ or $I_P(\cdot )$ over $\mathcal {A}$ has to have $r(0)>0$.

We assume for some positive constants $c_0$ and $c_1$,

(4.1)\begin{equation} c_0\tilde{\Phi}(v_1,v_2,v_3)\le\Phi(R,v_1,v_2,v_3)\le c_1\tilde{\Phi}(v_1,v_2,v_3), \end{equation}

for all $R\in [0,\,1]$ and where $\tilde {\Phi }$ corresponds to an isotropic and frame-indifferent material that satisfies (3.1) with H1–H3. We further assume that:

1. H4: For any $\eta >1$,

   \[ \dfrac{v^2}{(v^3-1)^2}\,\tilde{\Phi}\left(\frac{1}{v^2},v,v\right)\in L^1(\eta,\infty). \]

We now state and prove the above mentioned result for the displacement problem. The proof given here of this fact is an adaptation (to the problem with self-gravity) of the technique used in [Reference Sivaloganathan11] to establish a similar fact for compressible inhomogeneous materials.

Proof. We consider an incompressible deformation given by

\[ r_{inc}(R)=\sqrt[3]{R^3+\lambda^3-1}. \]

It follows that $r_{inc}\in \mathcal {A}_\lambda$ for $\lambda \ge 1$. If $r_\lambda$ is any minimizer of $I(\cdot )$ over $\mathcal {A}_\lambda$, then

\begin{align*} \Delta I& =I(r_{inc})-I(r_\lambda)\\ & \le c_1\int_0^1\tilde{\Phi}(r_{inc}(R))R^2\,\mathrm{d} R-c_0 \int_0^1\tilde{\Phi}(r_{\lambda}(R))R^2\,\mathrm{d} R\\ & \quad -\int_0^1 \rho_0(R)\dfrac{M_R}{r_{inc}(R)}\,R^2\,\mathrm{d} R+ \int_0^1 \rho_0(R)\dfrac{M_R}{r_\lambda(R)}\,R^2\,\mathrm{d} R. \end{align*}

By [Reference Sivaloganathan10, Proposition 4.10] we have that for $\lambda _1\le \lambda _2$

\[ r_{\lambda_1}(R)\le r_{\lambda_2}(R),\quad0\le R\le1. \]

Hence, for $\lambda _0$ fixed, we get that for $\lambda \ge \lambda _0$,

\[ \int_0^1 \rho_0(R)\dfrac{M_R}{r_\lambda(R)}\,R^2\,\mathrm{d} R\le \int_0^1 \rho_0(R)\dfrac{M_R}{r_{\lambda_0}(R)}\,R^2\,\mathrm{d} R. \]

It follows now that

\begin{align*} \Delta I& \le c_1\int_0^1\tilde{\Phi}(r_{inc}(R))R^2\,\mathrm{d} R-c_0 \int_0^1\tilde{\Phi}(r_{\lambda}(R))R^2\,\mathrm{d} R\\ & \quad +\int_0^1 \rho_0(R)\dfrac{M_R}{r_{\lambda_0}(R)}\,R^2\,\mathrm{d} R. \end{align*}

If $r_\lambda (0)=0$, it follows (cf. [Reference Ball3]) that

\[ \int_0^1\tilde{\Phi}(r_{\lambda}(R))R^2\,\mathrm{d} R\ge \int_0^1\tilde{\Phi}(r_h(R))R^2\,\mathrm{d} R, \]

where $r_h(R)=\lambda R$. Hence,

\begin{align*} \Delta I& \le c_1\int_0^1\tilde{\Phi}(r_{inc}(R))R^2\,\mathrm{d} R-c_0 \int_0^1\tilde{\Phi}(r_h(R))R^2\,\mathrm{d} R\\ & \quad +\int_0^1 \rho_0(R)\dfrac{M_R}{r_{\lambda_0}(R)}\,R^2\,\mathrm{d} R. \end{align*}

Since the third integral on the right-hand side of this inequality is fixed, the result now follows as in [Reference Sivaloganathan12] using H2 and H4.

We now state and prove the corresponding result for the internal pressure problem. The proof given here of this fact is an adaptation to the problem for inhomogeneous materials and with self-gravity, of the technique used in [Reference Sivaloganathan11] to establish a similar fact for compressible homogeneous materials.

Proof. We consider an incompressible deformation given by

\[ r_{\varepsilon}(R)=\sqrt[3]{R^3+\varepsilon^3}, \]

whwere $\varepsilon$ is as in (1.2). It follows that $r_{\varepsilon }\in \mathcal {A}$. If $r_P$ is any minimizer of $I_P(\cdot )$ over $\mathcal {A}$, then

\begin{align*} \Delta I_P& =I_P(r_{\varepsilon})-I_P(r_P)\\ & \le c_1\int_0^1\tilde{\Phi}(r_{\varepsilon}(R))R^2\,\mathrm{d} R-c_0 \int_0^1\tilde{\Phi}(r_P(R))R^2\,\mathrm{d} R\\ & \quad -\int_0^1 \rho_0(R)\dfrac{M_R}{r_{\varepsilon}(R)}\,R^2\,\mathrm{d} R+ \int_0^1 \rho_0(R)\dfrac{M_R}{r_P(R)}\,R^2\,\mathrm{d} R\\ & \quad -H_P(r_\varepsilon(0))+H_P(r_P(0)). \end{align*}

Since $r_\varepsilon (0)=\varepsilon$, we have that $H_P(r_\varepsilon (0))\ge {1}/{3}P_*\varepsilon ^3$. Moreover if $r_P(0)=0$, then $H_P(r_p(0))=0$. Using these in the above inequality and the non-negativity of $\tilde {\Phi }$, we have that

\[ \Delta I_P \le c_1\int_0^1\tilde{\Phi}(r_{\varepsilon}(R))R^2\,\mathrm{d} R +\int_0^1 \rho_0(R)\dfrac{M_R}{r_P(R)}\,R^2\,\mathrm{d} R -\frac{1}{3}P_*\varepsilon^3. \]

Thus, provided that the potential integral for $r_P$ (second term on the right) is bounded by some constant independent of $P_*$, we get that the right-hand side of the inequality above becomes negative for $P_*$ sufficiently large, contradicting the minimality of $r_P$. Thus, we must have that $r_P(0)>0$.

To show the uniform boundedness of the potential integral for $r_P$, we argue by contradiction. That is, assume that for a sequence of functions $\left \{P_k\right \}$ with the stated properties, we have that

(4.2)\begin{equation} \int_0^1 \rho_0(R)\dfrac{M_R}{r_{P_k}(R)}\,R^2\,\mathrm{d} R\rightarrow\infty,\quad k\rightarrow\infty, \end{equation}

where $r_{P_k}$ is a minimizer of $I_{P_k}$ with $r_{P_k}(0)=0$. From (3.2) (which holds for $r\in \mathcal {A}$ as well) and (H3), we get

(4.3)\begin{equation} \left\vert \int_0^1 \rho_0(R)\dfrac{M_R}{r_{P_k}(R)}\,R^2\,\mathrm{d} R\right\vert\le C\left(\int_0^1\tilde{\Phi}(r_{P_k}(R))R^2\,\mathrm{d} R\right)^{\frac{1}{3s}}. \end{equation}

Thus, using (4.2) we get that

(4.4)\begin{equation} \int_0^1\tilde{\Phi}(r_{P_k}(R))R^2\,\mathrm{d} R\rightarrow\infty,\quad k\rightarrow\infty. \end{equation}

It follows from (4.1) that

\[ I_{P_k}(r_{P_k})\ge c_0 \int_0^1\tilde{\Phi}(r_{P_k}(R))R^2\,\mathrm{d} R- \int_0^1 \rho_0(R)\dfrac{M_R}{r_{P_k}(R)}\,R^2\,\mathrm{d} R, \]

where we used that $H_{P_k}(r_{P_k}(0))=0$ because $r_{P_k}(0)=0$. Combining this with (4.3) gives

(4.5)\begin{equation} I_{P_k}(r_{P_k})\ge c_0 \int_0^1\tilde{\Phi}(r_{P_k}(R))R^2\,\mathrm{d} R- C\left(\int_0^1\tilde{\Phi}(r_{P_k}(R))R^2\,\mathrm{d} R\right)^{{1}/{3s}}. \end{equation}

If we let $\hat {r}(R)\equiv R$, then as $\hat {r}\in \mathcal {A}$, we get

(4.6)\begin{equation} I_{P_k}(r_{P_k})\le I_{P_k}(\hat{r})=I(\hat{r}). \end{equation}

But by (H3) we have that ${1}/{3s}<1$, and thus by (4.4) and (4.5) that $I_{P_k}(r_{P_k})\rightarrow \infty$, which contradicts (4.6). Thus, the potential integrals must be bounded uniformly among minimizers with $r_P(0)=0$.

5. No cavitation results

We now collect some results for the displacement problem under which the minimizer of $I(\cdot )$ over $\mathcal {A}_\lambda$ for $\lambda$ sufficiently small must satisfy that $r(0)=0$.

We first consider the case of a homogeneous material for which (3.13) holds. The functional $I$ is now given by

(5.1)\begin{equation} I(r)=\int_0^1 \left[\tilde{\Phi}\left(r(R) \right)-\frac{4\pi}{3}\,\rho_0^2\dfrac{R^3}{r(R)}\right]\,R^2\, \mathrm{d} R. \end{equation}

We denote by $\lambda _c^h$ the critical boundary displacement for the cavitation problem considered in [Reference Ball3] with stored energy function $\tilde {\Phi }$.

Proof. To argue by contradiction, assume that $r(0)>0$. For $\hat {\lambda }=(\lambda +\lambda _c^h)/2$, we let

\[ R_0=\inf\left\{R\,|\,r(R)=\hat{\lambda}R\right\}. \]

Since $\lambda <\lambda _c^h$ and $\frac {r(R)}{R}\rightarrow \infty$ as $R \rightarrow 0^+$, we get that $R_0$ is well defined and $R_0>0$. We define

\[ \hat{r}(R)=\left\{\begin{array}{@{}rcl}\hat{\lambda}R & , & R\le R_0,\\ r(R) & , & R>R_0.\end{array}\right. \]

If follows that $\hat {r}\in \mathcal {A}_\lambda$. Now

\begin{align*} \Delta I& =I(r)-I(\hat{r})= \int_0^{R_0} \left[\tilde{\Phi}\left(r(R)\right)- \tilde{\Phi}\left(\hat{\lambda},\hat{\lambda}, \hat{\lambda}\right)\right]R^2\,\mathrm{d} R\\ & \quad +\frac{4\pi}{3}\,\rho_0^2\int_0^{R_0}\left[\dfrac{1}{\hat{r}(R)} -\dfrac{1}{r(R)}\right]\,R^5\,\mathrm{d} R,\\ & \ge\int_0^{R_0} \left[\tilde{\Phi}\left(r(R)\right)- \tilde{\Phi}\left(\hat{\lambda},\hat{\lambda}, \hat{\lambda}\right)\right]R^2\,\mathrm{d} R, \end{align*}

as $\hat {r}(R)\le r(R)$ for $R\le R_0$. Since $\hat {\lambda }<\lambda _c^h$ and $r(0)>0$, the results in [Reference Ball3] imply that

\[ \int_0^{R_0}\tilde{\Phi}\left(r(R)\right)R^2\,\mathrm{d} R> \int_0^{R_0} \tilde{\Phi}\left(\hat{\lambda},\hat{\lambda}, \hat{\lambda}\right)R^2\,\mathrm{d} R, \]

and thus that $\Delta I>0$ which contradicts the minimality of $r$.

We now consider inhomogeneous materials of the form

(5.2)\begin{equation} \Phi(R,v_1,v_2,v_3)=\alpha(R)\sum_i\phi(v_i)+\beta(R)\sum_{i< j}\psi(v_iv_j)+ \gamma(R)h(v_1v_2v_3), \end{equation}

where $\alpha$, $\beta$ and $\gamma$ are smooth positive functions over $[0,\,1]$, and $\phi$, $\psi$, $h$ are non-negative convex functions, with $h$ strictly convex satisfying H2 and H3. We let $d_0$ be the value of the argument at which $h$ assumes its global minimum value. Note that $h^\prime (d)<0$ for $d< d_0$.

Proof. As in the proof of theorem 5.1, we argue by contradiction. Thus, we assume that $r(0)>0$ and let $\hat {\lambda }=(\lambda +d_0^{\frac {1}{3}})/2$. Define $R_0$ and $\hat {r}$ as in the proof of theorem 5.1. Then

\begin{align*} \Delta I& =I(r)-I(\hat{r})= \int_0^{R_0} \left[\Phi\left(R,r(R)\right)- \Phi\left(R,\hat{\lambda}R\right)\right]R^2\,\mathrm{d} R\\ & \quad +\int_0^{R_0}\rho_0(R)M_R\left[\dfrac{1}{\hat{r}(R)} -\dfrac{1}{r(R)}\right]\,R^2\,\mathrm{d} R,\\ & \ge\int_0^{R_0} \left[\Phi\left(R,r(R)\right)- \Phi\left(R,\hat{\lambda}R\right)\right]R^2\,\mathrm{d} R, \end{align*}

as $\hat {r}(R)\le r(R)$ for $R\le R_0$. Using the convexity of $\phi$, and $h$ in (5.2), we have now that (see [Reference Ball3, p. 589]):

\begin{align*} \left[\Phi\left(R,r(R)\right)- \Phi\left(R,\hat{\lambda}R\right)\right]R^2& \ge \alpha(R)\phi^\prime(\hat{\lambda})(r(R)R^2-\hat{\lambda}R^3)^\prime \\ & \quad +\frac{1}{3}h^\prime(\hat{\lambda}^3)(r^3(R)-\hat{\lambda}^3R^3)^\prime. \end{align*}

It follows now, after integrating by parts, that

\[ \int_0^{R_0} \left[\Phi\left(R,r(R)\right)- \Phi\left(R,\hat{\lambda}R\right)\right]R^2\,\mathrm{d} R\ge -\frac{1}{3}h^\prime(\hat{\lambda}^3)r^3(0)>0. \]

Thus, $\Delta I>0$ which contradicts the minimality of $r$.

For the stored energy function (5.2), it is easy to check that for some constant $C>0$,

(5.3)\begin{equation} \left\vert \frac{\partial\Phi}{\partial R}(R,v_1,v_2,v_3)\right\vert\le C\,\tilde{\Phi}(v_1,v_2,v_3), \end{equation}

where

(5.4)\begin{equation} \tilde{\Phi}(v_1,v_2,v_3)=\sum_i\phi(v_i)+\sum_{i< j}\psi(v_iv_j)+h(v_1v_2v_3). \end{equation}

We now denote by $\lambda _c^h$ the critical boundary displacement corresponding to the stored energy function $\Phi (0,\,v_1,\,v_2,\,v_3)$. Note that any deformation with finite $\Phi (0,\,v_1,\,v_2,\,v_3)$ energy has finite $\tilde {\Phi }$ energy as well. The following result is reminiscent to [Reference Sivaloganathan11, Proposition 12]. It shows that any minimizer which leaves the centre intact must have strains at the origin less than $\lambda _c^h$.

Theorem 5.3 Let $r\in \mathcal {A}_\lambda$ satisfy $r(0)=0$ and that $\ell \in (0,\,\infty ]$ where

\[ \lim_{R\rightarrow{0^+}}r^\prime(R)=\ell. \]

If $\ell >\lambda _c^h$, then $r$ can not be a minimizer of $I(\cdot )$ over $\mathcal {A}_\lambda$.

Proof. The following proof is similar to the one of [Reference Sivaloganathan11, Proposition 12] except for the treatment of the gravitational potential and the specific dependence on $R$ of the stored energy function.

For $\varepsilon \in (0,\,1)$, we let

\[ \lambda(\varepsilon)=\dfrac{r(\varepsilon)}{\varepsilon}. \]

From the given hypotheses, it follows that

\[ \lambda(\varepsilon)\rightarrow\ell,\quad\mbox{as }\varepsilon \rightarrow 0^+. \]

Assume for the moment that $\ell$ is finite, and let $r_c$ be a cavitating extrema corresponding to $\Phi (0,\,v_1,\,v_2,\,v_3)$. We define

(5.5)\begin{equation} r_\varepsilon(R)=\left\{\begin{array}{@{}rcl} \alpha_\varepsilon r_c\left(\dfrac{R}{\alpha_\varepsilon}\right), & R\in[0,\varepsilon],\\ r(R), & R\in(\varepsilon,1], \end{array}\right. \end{equation}

where $\alpha _\varepsilon$ is such that $\alpha _\varepsilon r_c(\varepsilon /\alpha _\varepsilon )=\lambda (\varepsilon )\varepsilon$. That $\alpha _\varepsilon$ exists follows from the fact that $\lambda (\varepsilon )>\lambda _c^h$ for $\varepsilon$ sufficiently small. Now

\begin{align*} \Delta I& =I(r_\varepsilon)-I(r)\\ & =\int_0^\varepsilon R^2\bigg[\Phi(R,r_\varepsilon(R))) -\Phi(R,r(R))\bigg]\,\mathrm{d} R\\ & \quad -\int_0^\varepsilon\rho_0(R)\dfrac{M_R}{r_\varepsilon(R)} \,R^2\,\mathrm{d} R +\int_0^\varepsilon\rho_0(R)\dfrac{M_R}{r(R)}\,R^2\,\mathrm{d} R \end{align*}

With the change of variables $R=\varepsilon U$ with $U\in [0,\,1]$, we can rewrite the first term above to get that

\begin{align*} \Delta I& =\varepsilon^3\int_0^1 U^2\bigg[\Phi\left(\varepsilon U,r_\varepsilon(\varepsilon U)\right) -\Phi\left(\varepsilon U,r(\varepsilon U)\right)\bigg]\,\mathrm{d} U\\ & \quad -\int_0^\varepsilon\rho_0(R)\dfrac{M_R}{r_\varepsilon(R)} \,R^2\,\mathrm{d} R +\int_0^\varepsilon\rho_0(R)\dfrac{M_R}{r(R)}\,R^2\,\mathrm{d} R \end{align*}

We first examine the gravitational integrals. For this, we use that $\rho _0(\cdot )$ is non-negative and bounded above, and that $M_R$ is bounded by a constant times $R^3$. Since

\[ \dfrac{r_\varepsilon(\varepsilon)}{\varepsilon}=\dfrac{r_c(\varepsilon/\alpha_\varepsilon)}{\varepsilon/\alpha_\varepsilon} =\lambda(\varepsilon) \rightarrow\ell\quad\mbox{as }\varepsilon \rightarrow 0^+, \]

we get that

\[ \dfrac{\varepsilon}{\alpha_\varepsilon}\rightarrow\mu, \]

where $\mu >0$ and $r_c(\mu )/\mu =\ell$. Upon recalling that ${r_c(S)}/{S}$ is a decreasing function of $S$, we have that:

\begin{align*} \int_0^\varepsilon\rho_0(R)\dfrac{M_R}{r_\varepsilon(R)} \,R^2\,\mathrm{d} R& = \varepsilon^3\int_0^1\rho_0(\varepsilon U)\dfrac{M_{\varepsilon U}}{r_\varepsilon(\varepsilon U)} \,U^2\,\mathrm{d} U\\ & \le K_1\varepsilon^5\int_0^1\dfrac{U^4}{r_\varepsilon(\varepsilon U)/\varepsilon U} \,\mathrm{d} U\le \dfrac{K\varepsilon^5}{5\ell}. \end{align*}

As $r(0)=0$ and $\ell >0$, the function $r(R)/R$ is positive and continuous in $[0,\,1]$. Thus, if $v_0$ is its minimum value, we have that for some positive constant $L$:

\[ \int_0^\varepsilon\rho_0(R)\dfrac{M_R}{r(R)}\,R^2\,\mathrm{d} R= \varepsilon^3\int_0^1 \rho_0(\varepsilon U)\dfrac{M_{\varepsilon U}}{r(\varepsilon U)}\,U^2\,\mathrm{d} U\le \frac{L\varepsilon^5}{5v_0}. \]

Thus, both gravitational potential terms go to zero faster than $\varepsilon ^3$.

We now examine the mechanical potential terms in $\Delta I$. For this, we note that

(5.6)\begin{align} & \int_0^1 U^2\bigg[\Phi\left(\varepsilon U,r_\varepsilon(\varepsilon U)\right) -\Phi\left(\varepsilon U,r(\varepsilon U)\right)\bigg]\,\mathrm{d} U\nonumber\\ & \quad =\int_0^1 U^2\bigg\{\bigg[\Phi\left(\varepsilon U,r_\varepsilon(\varepsilon U)\right)-\Phi\left(0,r_\varepsilon(\varepsilon U)\right)\bigg]\nonumber\\ & \qquad +\bigg[\Phi\left(0,r_\varepsilon(\varepsilon U)\right)-\Phi(0,\lambda(\varepsilon))\bigg] +\bigg[\Phi(0,\lambda(\varepsilon)) -\Phi\left(\varepsilon U,r(\varepsilon U)\right)\bigg]\bigg\}\,\mathrm{d} U \end{align}

where for simplicity, we have written

\begin{align*} \Phi\left(0,r_\varepsilon(\varepsilon U)\right)& \equiv \Phi\left(0,r_\varepsilon^\prime(\varepsilon U),\frac{r_\varepsilon(\varepsilon U)}{\varepsilon U},\frac{r_\varepsilon(\varepsilon U)}{\varepsilon U}\right),\\ \Phi(0,\lambda(\varepsilon))& \equiv\Phi(0,\lambda(\varepsilon),\lambda(\varepsilon),\lambda(\varepsilon)). \end{align*}

From (5.3) and Taylor’s Theorem, we have that

\[ \int_0^1U^2[\Phi\left(\varepsilon U,r_\varepsilon(\varepsilon U)\right)-\Phi\left(0,r_\varepsilon(\varepsilon U)\right)]\mathrm{d} U \le C_1\varepsilon \displaystyle\int_0^1 U^3\tilde{\Phi}\left(r_\varepsilon(\varepsilon U)\right)\mathrm{d} U\le C_2\varepsilon, \]

where we used that $r_c$ has finite $\tilde {\Phi }$ energy. Thus, the first bracketed term in (5.6) goes to zero with $\varepsilon$. For the third term, we note that the functions $r^\prime (S)$ and ${r(S)}/{S}$ are $C[0,\,1]$ and positive. Hence, for some $M>0$,

\[ \left\vert \Phi\left(\varepsilon U,r(\varepsilon U\right)\right\vert\le M,\quad\forall\ U, \]

and since

\[ \Phi\left(\varepsilon U,r(\varepsilon U)\right)\rightarrow\Phi(0,\ell,\ell,\ell), \]

pointwise, we get by the Lebesgue-dominated convergence theorem that

\[ \int_0^1U^2\Phi\left(\varepsilon U,r(\varepsilon U\right)\,\mathrm{d} U\rightarrow \int_0^1U^2\Phi(0,\ell,\ell,\ell)\,\mathrm{d} U, \]

as $\varepsilon \rightarrow 0^+$. This together with $\lambda (\varepsilon )\rightarrow \ell$ yields that

\[ \int_0^1U^2\bigg[\Phi(0,\lambda(\varepsilon)) -\Phi\left(\varepsilon U,r(\varepsilon U\right)\bigg]\,\mathrm{d} U \rightarrow 0, \]

as $\varepsilon \rightarrow 0^+$. Thus, the third term in (5.6) goes to zero with $\varepsilon$ as well.

For the second term in (5.6), note that with the change of variables $Z=(\varepsilon /\alpha _\varepsilon )U$, we get

\begin{align*} \int_0^1U^2\,\Phi\left(0,r_\varepsilon(\varepsilon U)\right)\,\mathrm{d} U& = \left(\frac{\alpha_\varepsilon}{\varepsilon}\right)^3 \int_0^{\varepsilon/\alpha_\varepsilon}Z^2\,\Phi\left(0,r_c(Z)\right)\, \mathrm{d} Z\\ & \rightarrow\frac{1}{\mu^3} \int_0^{\mu}Z^2\,\Phi\left(0,r_c(Z)\right)\,\mathrm{d} Z =\int_0^1U^2\,\Phi\left(0,r_c(\mu U)\right)\,\mathrm{d} U. \end{align*}

It follows now that

\begin{align*} & \int_0^1U^2\bigg[\Phi\left(0,r_\varepsilon(\varepsilon U)\right)-\Phi(0,\lambda(\varepsilon))\bigg] \,\mathrm{d} R\\ & \quad \rightarrow\int_0^1U^2\bigg[\Phi\left(0,r_c(\mu U)\right) -\Phi(0,\ell,\ell,\ell)\bigg]\,\mathrm{d} R<0, \end{align*}

where the last inequality follows since $\ell >\lambda _c^h$ and $\tilde {r}(U)=\mu ^{-1}r_c(\mu U)$ is the minimizer for the functional with stored energy $\Phi (0,\,v_1,\,v_2,\,v_3)$ and boundary condition $\tilde {r}(1)=\ell$. Collecting all of the intermediate results so far, we get that $\varepsilon ^{-3}\Delta I<0$ for $\varepsilon$ sufficiently small, which contradicts the minimality of $r$.

The case $\ell =\infty$ can be handled in a similar fashion using a suitable incompressible deformation on $[0,\,\varepsilon ]$ in (5.5). See [Reference Sivaloganathan11, Proposition 12] for details.

6. Numerical results

We now give some numerical examples illustrating several of the results from previous sections for the displacement problem. For the calculations, we employ a combination of two numerical schemes: (the predictor) a descent method for the minimization of (2.7) based on a gradient flow iteration; and (the corrector) a shooting method that solves directly the Euler–Lagrange boundary value problem (2.12), (3.8) and (3.10a). The use of adaptive ODE solvers in the shooting method allows for a more precise computation of the equilibrium states, especially near $R=0$ where both strains in our problem tend to develop boundary layers.

A gradient flow iteration (cf. [Reference Neuberger9]) assumes that $r$ depends on a flow parameter $t$, and that $r(R,\,t)$ satisfies

(6.1)\begin{align} \displaystyle\frac{\mathrm{d}^2}{\mathrm{d}R^2}(r_t(R,t))& ={-}\frac{\mathrm{d}}{\mathrm{d}R}\left[R^{2}\Phi_{,1} (R,r(R,t))\right] +2R\Phi_{,2}(R,r(R,t))\\ & \quad +R^2\,\dfrac{\rho_0(R)M_R}{r^2(R,t)},\quad0< R<1,~t>0,\nonumber \end{align}

(6.2)\begin{align} r(1,t)& =\lambda,\quad \lim_{R\rightarrow{0^+}}\left[\frac{\mathrm{d}}{\mathrm{d}R}(r_t(R,t))+R^2\Phi_{,1}(R,r(R,t))\right]=0, \, t\ge0. \end{align}

(Here $r_t=\frac {\partial r }{\partial t }$.) The gradient flow equation leads to a descent method for the minimization of (2.7) over (2.13). After discretization of the partial derivative with respect to ‘$t$’, one can use a finite element method to solve the resulting flow equation. In particular, if we let $\Delta t>0$ be given, and set $t_{i+1}=t_i+\Delta t$ where $t_0=0$, we can approximate $r_t(R,\,t_i)$ with:

\[ z_i(R)= \dfrac{r_{i+1}(R)-r_i(R)}{\Delta t}, \]

where $r_i(R)=r(R,\,t_i)$, etc. (We take $r_0(R)$ to be some initial deformation satisfying the boundary condition at $R=1$ , e.g. $\lambda R$.) Inserting this approximation into the weak form of (6.1), (6.2), and evaluating the right-hand side of (6.1) at $r=r_i$, we arrive at the following iterative formula:

(6.3)\begin{align} & \int_0^1z^\prime_i(R)v^\prime(R)\,\mathrm{d} R+ \int_0^1\Big[R^{2}\Phi_{,1}(R,r_i(R))v^\prime(R) \nonumber\\ & \quad +\left(2R\Phi_{,2}(R,r_i(R))+ R^2\,\dfrac{\rho_0(R)M_R}{r_i^2(R)}\right)v(R)\Big]\,\mathrm{d} R=0, \end{align}

for all functions $v$ with $v(1)=0$ and sufficiently smooth so that the integrals above are well defined. Given $r_i$, one can solve the above equation for $z_i$ via some finite element scheme, and then set $r_{i+1}=r_i+\Delta t\,z_i$. This process is repeated for $i=0,\,1,\,\ldots$, until $r_{i+1}-r_i$ is “small’ enough, or some maximum value of ‘$t$’ is reached, declaring the last $r_i$ as an approximate minimizer of (2.7) over (2.13).

In the shooting method technique, for given $\nu >0$, we solve the initial value problem

(6.4a)\begin{align} & \frac{\mathrm{d}}{\mathrm{d}R}\left[R^{2}\Phi_{,1}(R,r(R))\right]= 2R\Phi_{,2}(R,r(R))+ R^2\,\dfrac{\rho_0(R)M_R}{r^2(R)}, \quad0< R<1, \end{align}

(6.4b)\begin{align} & r(1)=\lambda,\quad r^\prime(1)=\nu, \end{align}

from $R=1$ to $R=0$. The value of $\nu$ is adjusted so that

(6.5)\begin{equation} \lim_{R\rightarrow{0^+}}\dfrac{R^2}{r^2(R)}\Phi_{,1}\left(R,\nu,\dfrac{r(R)}{R}, \dfrac{r(R)}{R}\right)=0. \end{equation}

In actual calculations, we solve (6.4) from $R=1$ to $R=\varepsilon$, where $\varepsilon >0$ is small, and replace (6.5) with

(6.6)\begin{equation} \dfrac{\varepsilon^2}{r^2(\varepsilon)}\Phi_{,1}\left(\varepsilon,\nu,\dfrac{r(\varepsilon)}{\varepsilon}, \dfrac{r(\varepsilon)}{\varepsilon}\right)=0. \end{equation}

This equation is solved for $\nu$ via a secant-type iteration, where the initial guess for $\nu$ is computed using the approximate solution from the predictor step. The solution of (6.6) requires repeated solutions of the initial value problem (6.4) from $R=1$ to $R=\varepsilon$. These intermediate initial value problems are solved with the routine ode45 of the MATLAB^TM ode suite.

For the simulations, we consider the homogeneous case (5.1) for which the stored energy function $\tilde {\Phi }$ is given by

\[ \tilde{\Phi}(v_1,v_2, v_3)=\frac{\kappa}{p}\,(v_1^p+v_2^p+v_3^p)+h(v_1v_2v_3), \]

with

\[ h(d)=C\,d^{\gamma}+D\,d^{-\delta}, \]

where $p<3$, $C\ge 0$, $D\ge 0$ and $\gamma,\,\delta >0$. The reference configuration is mechanically stress-free provided:

\[ D=\dfrac{\kappa+C\gamma}{\delta}. \]

The mass density function $\rho _0$ is taken to be constant. In this case, the differential equation in (6.4a) reduces to (3.14). For the first set of simulations, we used the following values of the mechanical parameters in $\tilde {\Phi }$:

(6.7)\begin{equation} p=2,\quad\kappa=1,\quad C=1,\quad\gamma=\delta=2, \end{equation}

with a value of $\varepsilon =0.001$ in (6.6) and in the shooting method. The gradient flow iteration was used as a predictor for the shooting method, with the integrals in (6.3) computed over $(\varepsilon,\,1)$ as well.

In our first simulation, we show the approximate cavity radius as a function of $\rho _0$ and $\lambda$, for values of $\rho _0\in [0.5,\,1.5]$ and $\lambda \in [0.9,\,1.2]$. The resulting surface is shown in figure 1. We note that for fixed values of $\rho _0$, the cavity size is essentially zero up to some certain value of $\lambda$ (the critical boundary displacement corresponding to $\rho _0$), after which the graph becomes concave. This critical boundary displacement appears to be an increasing function of $\rho _0$. In figure 2, we show the corresponding surface for the energies of the approximate minimizers. For fixed values of $\lambda$, the energy is an increasing function of $\rho _0$, while for $\rho _0$ constant the energy is non-monotone with respect to $\lambda$ with a convex shape.

Figure 1.

Cavity surface.

Figure 2.

Energy surface.

Our next simulations are for fixed values of $\lambda$ and $\rho _0$. In the first case $\lambda =1$ and $\rho _0=1$. In figure 3, we show the computed minimizer $r$ compared to the affine deformation $\lambda R$. The value of $r(0.001)$ is $7.7078\times 10^{-4}$ with an energy of $0.49074$. Figure 4 shows plots of the strains $r^\prime (R)$ and ${r(R)}/{R}$, and the determinant $r^\prime (R)(r(R)/R)^2$ in this case. Also in figure 5, we show the graph of the corresponding Cauchy stress. We note the boundary layer close to $R=\varepsilon$ in these plots. This boundary layer is a numerical artefact since the numerical scheme tries to make the Cauchy stress zero, while for this value of $\lambda$ the value of $r(0)$ should be zero. Still in this case the numerical scheme converges to the solution with $r(0)=0$ as $\varepsilon \rightarrow 0^+$ (cf. [Reference Sivaloganathan10, Section 5]).

Figure 3.

Radial displacement for $\lambda =1$ and $\rho _0=1$.

Figure 4.

Strains and determinant for $\lambda =1$ and $\rho _0=1$.

Figure 5.

Cauchy stress for $\lambda =1$ and $\rho _0=1$.

In our next simulation for the values (6.7), we take $\lambda =1.15$ and $\rho _0=1$. In figure 6, we show the computed minimizer $r$ which corresponds to a cavitating solution. The value of $r(0.001)$ is $0.48346$ with an energy of $0.91034$. Figure 7 shows plots of the strains $r^\prime (R)$ and ${r(R)}/{R}$, and the determinant $r^\prime (R)(r(R)/R)^2$, and in figure 8, we show the graph of the corresponding Cauchy stress. The boundary layer in $r(R)/R$ is now a ‘true’ one associated with the computed cavitating solution. Finally, in figure 9, we show a three-dimensional ‘cut’ of the deformed body with the colouring given by the density in the deformed configuration which is $\rho _0/[r^\prime (R)(r(R)/R)^2]$. Note that the density in the deformed configuration is higher near the centre, decreasing towards the outer radius.

Figure 6.

Radial displacement for $\lambda =1.15$ and $\rho _0=1$.

Figure 7.

Strains and determinant for $\lambda =1.15$ and $\rho _0=1$.

Figure 8.

Cauchy stress for $\lambda =1.15$ and $\rho _0=1$.

Figure 9.

Deformed configuration with colouring given by the density in the deformed configuration for $\lambda =1.15$ and $\rho _0=1$.

Our next numerical example uses values characteristic of the composition of the planet Mercury. Differential equation (3.14) now is over the interval $[0,\,R_0]$, where $R_0$ is the planet average radius, and the coefficient $\rho _0^2$ in the last term must be replaced with $\rho _0^2G$, where $G$ is the gravitational constant. In the normalized variables $R\leftrightarrow R/R_0$ and $r\leftrightarrow r/R_0$, differential equation (3.14) becomes:

(6.8)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}R}\left[R^{2}\tilde{\Phi}_{,1}(r(R))\right]= 2R\tilde{\Phi}_{,2}(r(R))+ \frac{4\pi}{3}\, (\rho_0^2GR_0^2)\,\dfrac{R^5}{r^2(R)}, \quad 0< R<1. \end{equation}

Now using that $G\,{=}\,6.674\times 10^{-11}\,\mbox {m}^3\,\mbox {kg}^{-1}\,\mbox {s}^{-2}$, $R_0\,{=}\,2440\,\mbox {km}$ and $\rho _0\,{=}\,7860\,\mbox {kg}\,\mbox {m}^{-3}$ (for an iron planet), we get that

\[ \rho_0^2GR_0^2=24.5\,\mbox{GPa}. \]

For the mechanical parameters resembling an iron composition as used in [Reference Jia, Kodio, Chapman and Goriely8], we take

(6.9)\begin{equation} p=2,\quad\kappa=80\,\mbox{GPa},\quad C=60\,\mbox{GPa},\quad\gamma=\delta=2. \end{equation}

In the simulation, we compared the energies of the solution of (6.8) with boundary conditions (6.5) and $r(1)=\lambda$, to the solution of (6.8) with boundary conditions $r(0)=0$ and $r(1)=\lambda$. We denote these solutions as $r_Z$ (zero Cauchy stress at origin) and $r_I$ (centre intact), respectively. For $\lambda$ small, the solutions of these two boundary value problems are equal. However for the simulations, the first boundary value problem is solved on $[\varepsilon,\,1]$ with $\varepsilon$ small, with (6.5) replaced by (6.6). This causes an ‘artificial’ boundary layer close to $R=\varepsilon$ for the strains corresponding to $r_Z$, which disappears as $\varepsilon$ tends to zero and with $r_Z(\varepsilon )>0$ but small. The energies of the $r_Z$ for small $\lambda$’s are consequently slightly larger than those corresponding to $r_I$. However, for values of $\lambda$ larger, the energy of $r_Z$ is indeed lower than that of $r_I$, an indication that cavitation is energetically favourable. In figure 10, we show a plot of the differences in these energies for a range of values of $\lambda$. We recall that the horizontal axis in this graph is in normalized values. Thus, for instance, a value of $\lambda =1$ would correspond to $2440\,\mbox {km}$. The numerical results suggest a critical boundary displacement (for the initiation of cavitation) at about $\lambda =1.145$ or $2790\,\mbox {km}$ approximately.

Figure 10.

Cavity size vs boundary displacement (in normalized variables) for parameter values resembling the composition of the planet Mercury.

7. Final comments

The usual self-gravitating problem is that in which the centre of the body remains intact ($r(0)=0$) and no condition is explicitly prescribed on the outer boundary. This is the problem considered in [Reference Jia, Kodio, Chapman and Goriely8] and is a special case of one of the problems treated in [Reference Calogero and Leonori6]. It is straightforward to check that our results hold in this case as well where the admissible set is now given by

\[ \mathcal{A}=\big\{r\in W^{1,1}(0,1)\,|\,r(0)=0,\,r^\prime(R)>0\mbox{~a.e. for~}R\in(0,1),\,I_{\mbox{mec}}(r)<\infty\big\}. \]

In particular, minimizers exist for all densities $\rho _0$ and satisfy the Euler–Lagrange equation (3.8) with natural boundary condition at $R=1$ given by $T(1,\,r(1))=0$ (cf. (3.7)). In reference to the stored energy function (1.1) considered in [Reference Jia, Kodio, Chapman and Goriely8], growth condition (3.1) with H3 places our stored energy function under the ‘strong’ compressibility category. Thus, our result on the existence of minimizers with the centre intact for all densities $\rho _0$ is consistent with the results in [Reference Jia, Kodio, Chapman and Goriely8] which in turn suggest that there might be uniqueness of solutions of (3.8) in our problem as well.

Another interesting problem is that of a spinning body. In the three-dimensional case, radial solutions to this problem are not possible and a full three-dimensional variational problem must be considered. The two-dimensional problem of a spinning disk can be studied with radial symmetry. However, for cavitation, one would have condition H1 with $0<\gamma <2$, and it is easy to check that in this case functional (2.7) is in general unbounded below.