Cavitation of a spherical body under mechanical and self-gravitational forces

Abstract

In this paper, we look for minimizers of the energy functional for isotropic compressible elasticity taking into consideration the effect of a gravitational field induced by the body itself. We consider two types of problems: the displacement problem in which the outer boundary of the body is subjected to a Dirichlet-type boundary condition, and the one with zero traction on the boundary but with an internal pressure function. For a spherically symmetric body occupying the unit ball $\mathcal {B}\in \mathbb {R}^3$, the minimization is done within the class of radially symmetric deformations. We give conditions for the existence of such minimizers, for satisfaction of the Euler–Lagrange equations, and show that for large displacements or large internal pressures, the minimizer must develop a cavity at the centre. We discuss a numerical scheme for approximating the minimizers for the displacement problem, together with some simulations that show the dependence of the cavity radius and minimum energy on the displacement and mass density of the body.

1. Introduction

The study of the shape of self-gravitating bodies is extensive and dates back to the time of Newton itself. It is well known that depending on the density of a dying star, there are several possibilities for the resulting object: white dwarf, neutron star, black hole, etc. The case of a black hole forming is also referred to as gravitational collapse. The literature on these phenomena is extensive and we refer to [6] and [8] for a historical account.

In this paper, we consider the problem of a self-gravitating spherical body. Apart from its apparent ‘simplicity’, this problem plays an important role on the study of the more complex phenomena described above. The proposed variational model combines both mechanical and gravitational responses, the mechanical part based on a model from non-linear elasticity which allows for the characterization of large deformations. Under certain mathematically physical conditions, the extrema of the corresponding energy functional can be characterized via the Euler–Lagrange equations. This combined model has been used by [5], [6] and [8] among others. In [5] the existence of solutions to the Euler–Lagrange equations, with a zero dead load boundary condition on the outer boundary of the body, is established via the implicit function theorem and is valid for ‘small bodies’ of arbitrary shape.

The work in [8] is for spherically symmetric deformations with a zero dead load boundary condition on the outer boundary as well, and combines asymptotic analysis with numerics to get results for varying densities and reference configuration body radius. They used a stored energy function of the form

(1.1)\begin{equation} W(\mathbf{F})=\frac{\mu}{2}\left(\left\Vert \mathbf{F}\right\Vert^2-3-2f_\alpha(\det\mathbf{F})\right) +\frac{\beta}{2}\left(\det\mathbf{F}-1\right)^2, \end{equation}

where $\alpha \ge 0$, $\beta$ and $\mu$ are positive constants, and $f_\alpha (d)=\ln (d)-\alpha d^{-\alpha }(d-1)^4$. This material corresponds to a ‘soft’ compressible material for $\alpha =0$ or small, and to a ‘strong’ compressible material otherwise. For constant reference configuration density $\rho _0$, the authors in [8] show numerically that for $\alpha$ small there exists a critical density $\rho _0^*$ such that if $\rho _0\le \rho _0^*$, then the Euler–Lagrange equations (cf. (3.8)) can have multiple solutions, most of them unstable, while if $\rho _0>\rho _0^*$, then there are no solutions, a result which could be interpreted as gravitational collapse. Moreover for $\alpha$ large, there are solutions for all densities $\rho _0$, which appear to be unique.

By adapting the techniques in [1] for polyconvex stored energy functions, the authors in [6] show the existence of minimizers for the resulting energy functional, now for large deformations and arbitrary bodies, and for both, zero dead load and displacement boundary conditions. The stored energy functions used in [6] could be classified as corresponding to ‘strong’ compressible materials (cf. eqns. (13) and (14) in [6]).

In this paper, we look for minimizers of the energy functional for isotropic compressible elasticity and taking into consideration the effect of a gravitational field induced by the body itself. We consider both, the displacement problem in which the outer boundary of the body is subjected to a Dirichlet-type boundary condition, and the one with zero traction on the boundary but with a specified internal pressure function. For a spherically symmetric body occupying the unit ball $\mathcal {B}\in \mathbb {R}^3$ and with radially symmetric mass density, the minimization is done within the class of radially symmetric deformations. Contrary to previous works, the deformations we consider belong to $W^{1,p}(\mathcal {B})$ with $p<3$, and thus may develop singularities. For the particular case of radially symmetric deformations, we study the occurrence or initiation of a cavitation at the centre of the ball and its dependence on the boundary displacement, internal pressure and gravitational-related constants.

In § 2 we introduce the basic model, energy functionals and admissible function spaces, for radial deformations (cf. (2.5)) of a spherically symmetric body. These deformations are characterized by a function $r\,:\,[0,\,1]\rightarrow [0,\,\infty )$. For the displacement problem, the boundary condition takes the form $r(1)=\lambda$. We now show in § 3 that under certain growth conditions on the stored energy function (cf. (3.1) with H1–H3) and for any (radial) reference configuration density function $\rho _0$ that is bounded, non-negative and bounded away from zero, a minimizer of the energy functional (2.7)–(2.9) exists over the admissible set (2.13).

For the internal pressure problem, we consider an internal pressure function $P:[0,\,\infty )\rightarrow [0,\,\infty )$ satisfying

(1.2)\begin{equation} P(A)=P_*\mbox{~for~} 0\le A\le\varepsilon, \quad\mbox{and}\quad A^3P(A)\rightarrow c\quad\mbox{as}\quad A \rightarrow \infty, \end{equation}

where $\varepsilon,\,c>0$ are constants. The motivation for this model comes from experiments performed by Gent and Tompkins [7], where samples of polymers were left in a gas under pressure, over a period of time (hours). The pressurizing gas then induces a negative hydrostatic pressure inside the polymer. It was observed that upon release, small holes formed inside the polymer if the initial confining pressure was sufficiently large. Under the same growth conditions as in the displacement problem, we show in § 3 that a minimizer of the energy functional¹ (2.14)–(2.15), which includes the model pressure function (1.2), exists over the admissible set (2.16).

Under the additional constitutive assumption (3.6), these minimizers satisfy the Euler–Lagrange equations (3.8) where either $r(0)=0$ or $r(0)>0$ (cavitation) with zero Cauchy stress at the origin in the displacement problem, and with zero modified Cauchy stress² at the origin for the internal pressure problem. In § 4 we show that for $\lambda$ or $P_*$ sufficiently large, these minimizers must satisfy $r(0)>0$. The result for the displacement problem is an adaptation to the problem with self-gravity of a similar result in [11] for compressible inhomogeneous materials, while that for the internal pressure problem is again an adaption to the self-gravity problem of the corresponding problem treated in [12].

In § 5 we collect several results for the displacement problem with $\lambda$ small where the minimizers must have the centre intact. In addition we show in Theorem 5.3 that any minimizer which leaves the centre intact must have strains at the origin less than the critical boundary displacement corresponding to a homogeneous isotropic material made of the material at the centre of the original body. Once again this result is an adaptation to the problem with self-gravity of a similar result in [11] for compressible inhomogeneous materials.

Finally in § 6 we present a numerical scheme for the computation of the minimizers of the energy functional for the displacement problem. This method is based on a combination of a gradient flow iteration which works as a predictor, together with a shooting method to solve the EL-equations, that works as a corrector. For constant reference configuration densities, we present several simulations that show the dependence of the cavity radius and minimum energy on the displacement $\lambda$ and density $\rho _0$. In addition, we include a simulation for a spherical body with composition resembling that of the planet Mercury.

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 [6]) 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 gravitational³ 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.1 The function $\rho _0(\cdot )$ is the so-called natural density which refers to the density of the material in a stress-free configuration without gravity. On the other hand, the current density is the mass density function due to the deformation $r(\cdot )$ with gravity, and is given by $\rho (r(R))=\rho _0(R)/[r^\prime (R)(r(R)/R)^2]$. Thus, the assumption that the current density is constant implies that $r^3(R)=KM_R+r^3(0)$ for some positive constant $K$, which essentially determines the deformation.

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)$.

Remark 2.3 Including a term like $P_*\,r^\prime (R)(r(R)/R)^2$ inside the functional $I$, specifying a Cauchy stress of $P_*$ at the centre of the ball, would lead to a functional unbounded from below (see [3]). Growth condition (1.2) not only yields a functional bounded from below as we will show, but it is consistent with the original experiments reported in Gent and Tompkins [7] for the so-called ideal gases.

The following result from [12, 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 [3] due to the presence of the potential energy functional (2.9). We do emphasize that they are not direct consequence of those in [6] 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 [6] (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:

Proposition 3.2 Under growth assumption (3.1) with H1 and H3, the functionals $I(\cdot )$ and $I_P(\cdot )$ are bounded below on $\mathcal {A}_\lambda$ and $\mathcal {A}$, respectively.

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. [2]), 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. [4]): 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 [3] 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 [11] to establish a similar fact for compressible inhomogeneous materials.

Theorem 4.1 Let $r$ be a minimizer of functional (2.7) over (2.13) and assume that (4.1) holds where $\tilde {\Phi }$ satisfies (3.1) with H1–H4. Then for $\lambda$ sufficiently large, we must have that $r(0)>0$.

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 [10, 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. [3]) 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 [12] 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 [11] to establish a similar fact for compressible homogeneous materials.

Theorem 4.2 Let $r_P$ be a minimizer of (2.15) over the set (2.16) with $P(\cdot )$ satisfying (1.2). Assume that (4.1) holds where $\tilde {\Phi }$ satisfies (3.1) with H1–H4. Then for $P_*$ sufficiently large we must have that $r_P(0)>0$.

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 [3] with stored energy function $\tilde {\Phi }$.

Theorem 5.1 Let $r$ be a minimizer of (5.1) over $\mathcal {A}_\lambda$. If $\lambda <\lambda _c^h$, then $r(0)=0$.

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 [3] 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$.

Theorem 5.2 Let $r$ be a minimizer of (2.7) over $\mathcal {A}_\lambda$ corresponding to the stored energy function (5.2) with $\beta \equiv 0$ and $\gamma \equiv 1$. Assume that $\alpha ^\prime (R)\le 0$ for all $R$ and that $\phi ^\prime (t)\ge 0$ for all $t$. Then if $\lambda ^3< d_0$ we must have that $r(0)=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 [3, 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 [11, 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 [11, 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 [11, 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. [9]) 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. [10, 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 [8], 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 [8] and is a special case of one of the problems treated in [6]. 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 [8], 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 [8] 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.