A Fourier-Chebyshev Spectral Method for Cavitation Computation in Nonlinear Elasticity

TL;DR

This paper introduces a Fourier-Chebyshev spectral method for accurately solving cavitation problems in nonlinear elasticity, with proven convergence and efficient algorithms demonstrated through numerical experiments.

Contribution

The paper develops a novel spectral method combining Fourier and Chebyshev techniques for cavitation in nonlinear elasticity, including error analysis and convergence proof.

Findings

The method achieves high accuracy in cavitation simulations.

Numerical experiments confirm the efficiency and precision of the approach.

The approach converges reliably for the nonlinear elasticity cavitation problem.

Abstract

A Fourier-Chebyshev spectral method is proposed in this paper for solving the cavitation problem in nonlinear elasticity. The interpolation error for the cavitation solution is analyzed, the elastic energy error estimate for the discrete cavitation solution is obtained, and the convergence of the method is proved. An algorithm combined a gradient type method with a damped quasi-Newton method is applied to solve the discretized nonlinear equilibrium equations. Numerical experiments show that the Fourier-Chebyshev spectral method is efficient and capable of producing accurate numerical cavitation solutions.

Tables

Table 1. Table 1: The minima of ( I N ​ M ​ P ) ρ subscript superscript 𝐼 𝑁 𝑀 𝑃 𝜌 (I^{NM}P)_{\rho} in Ω ′ superscript Ω ′ \Omega^{\prime} in various cases (independent of N 𝑁 N )

$\lambda \times \gamma$	$(\epsilon ,\gamma )$	$M$
		4	6	8	10	12	14
2	$({10}^{-2},1)$	2.66e-3	2.86e-3	2.86e-3	2.86e-3	2.86e-3	2.86e-3
	$({10}^{-3},1)$	8.57e-5	2.92e-4	2.88e-4	2.88e-4	2.88e-4	2.88e-4
	$({10}^{-4},1)$	-1.75e-4	3.27e-5	2.88e-5	2.89e-5	2.89e-5	2.89e-5
1.25	$({10}^{-3},{10}^{-1})$	3.97e-5	3.97e-5	3.97e-5	3.97e-5	3.97e-5	3.97e-5
	$({10}^{-4},{10}^{-2})$	3.96e-7	3.96e-7	3.96e-7	3.96e-7	3.96e-7	3.96e-7
	$({10}^{-5},{10}^{-3})$	3.96e-9	3.96e-9	3.96e-9	3.96e-9	3.96e-9	3.96e-9

Table 2. Table 2: The regressed quantities and parameters for the non-radially-symmetric case.

$q$	$c_1$	$c_2$	${\nu }_1$	${\nu }_2$	$q^{\infty }$
energy	1.45e+24	-3.10e+6	25	6.1	22.85959048
semi-major axis	-2.02e+14	-1.37e-1	17	2.1	1.67481624
semi-minor axis	-1.61e+14	-8.42e-2	16	2.1	1.42872097

Table 3. Table 3: The regressed quantities and parameters for the radially-symmetric case.

$\epsilon$	$q$	$c_2$	${\nu }_2$	$q^{\infty }$	ODE solution
${10}^{-3}$	energy	-1.38e+8	9.5	18.61960091	18.61960090
	cavity radius	-3.18e-2	2.0	1.26772534	1.26772534
${10}^{-4}$	energy	-2.22e+6	6.1	18.87582778	18.87582146
	cavity radius	-1.49e-1	2.1	1.25228561	1.25228643

Equations

$$E(\mathbf{u})={\int }_{{\Omega }_{\epsilon }}W(\nabla \mathbf{u}(\mathbf{x}))\mathrm{d}\mathbf{x},$$

$${\mathcal{A}}_{\epsilon }=\left\{\mathbf{u}\in W^{1,p}({\Omega }_{\epsilon }):\mathbf{u}\text{leavevmode}\text{is one-to-one}\text{leavevmode}a.e.,\text{leavevmode}\mathbf{u}|{}_{\partial \Omega }={\mathbf{u}}_0,\text{leavevmode}\det \nabla \mathbf{u}>0\text{leavevmode}a.e.\right\},$$

$$\{\begin{array}{ll}\mathrm{div}\frac{\partial W(\nabla \mathbf{u})}{\partial \nabla \mathbf{u}}=0,& \text{in}\text{leavevmode}{\Omega }_{\epsilon };\text{vspace}\ast 3mm\\ \frac{\partial W(\nabla \mathbf{u})}{\partial \nabla \mathbf{u}}\cdot \mathbf{n}=0,& \text{on}\text{leavevmode}\partial {\mathbb{B}}_{\epsilon }({\mathbf{x}}_0);\text{vspace}\ast 1mm\\ \mathbf{u}(\mathbf{x})={\mathbf{u}}_0(\mathbf{x}),& \text{on}\text{leavevmode}\partial \Omega ,\end{array}$$

$$W(\nabla \mathbf{u})=\kappa |\nabla \mathbf{u}{|}^p+h(\det \nabla \mathbf{u}),\nabla \mathbf{u}\in {\mathbb{M}}_{+}^{2\times 2},\text{leavevmode}1<p<2,$$

$$h(t)\to +\infty \text{leavevmode}\text{as}\text{leavevmode}t\to 0^{+},\text{leavevmode}\text{and}\text{leavevmode}\frac{h(t)}{t}\to +\infty \text{leavevmode}\text{as}\text{leavevmode}t\to +\infty .$$

$$\begin{array}{rl}{\mathcal{A}}_{(\epsilon ,\gamma )}({\mathbf{u}}_0)=\{\mathbf{u}\in W^{1,p}({\Omega }_{(\epsilon ,\gamma )}):\text{leavevmode}\mathbf{u}\text{leavevmode}\text{is one-to-one}\text{leavevmode}a.e.,\text{leavevmode}& \mathbf{u}|{}_{\partial {\mathbb{B}}_{\gamma }(0)}={\mathbf{u}}_0,\\ & \det \nabla \mathbf{u}>0\text{leavevmode}a.e.\}.\end{array}$$

$$D(\mathbf{u}):=\det \nabla \mathbf{u},F(\mathbf{u}):=\frac{1}{2}|\nabla \mathbf{u}{|}^2,g(t):=\kappa {\left(\sqrt{2t}\right)}^p,$$

$$E(\mathbf{u})={\int }_{{\Omega }_{(\epsilon ,\gamma )}}\left[g\left(F(\mathbf{u})\right)+h\left(D(\mathbf{u})\right)\right]\mathrm{d}\mathbf{x}.$$

$$\{\begin{array}{l}{x}_1=r\cos \theta ,\\ x_2=r\sin \theta ,\end{array}\text{and}\{\begin{array}{l}{u}_1=R(r,\theta )\cos \Theta (r,\theta ),\\ u_2=R(r,\theta )\sin \Theta (r,\theta ),\end{array}$$

$$\{\begin{array}{l}r=\frac{\gamma +\epsilon }{2}+\frac{\gamma -\epsilon }{2}\rho ,\\ \theta =\varphi ,\end{array}\text{and}\{\begin{array}{l}R(r,\theta )=P(\rho ,\varphi ),\\ \Theta (r,\theta )=Q(\rho ,\varphi )+\varphi ,\end{array}$$

$$D(P,Q)$$

$$F(P,Q)$$

$$E(P,Q)={\int }_{{\Omega }^{\prime }}\left[g\left(F(P,Q)\right)+h\left(D(P,Q)\right)\right]\frac{r}{{\rho }_r}\mathrm{d}\rho \mathrm{d}\varphi ,$$

$$\begin{split}\mathcal{A}_{\Omega^{\prime}}(\mathbf{u}_{0})=\{(P,Q):\leavevmode\nobreak\ \exists\leavevmode\nobreak\ \mathbf{u}\in\mathcal{A}_{(\varepsilon,\gamma)}(\mathbf{u}_{0}),\,&\text{{\it s.t.} }\mathbf{u}\leavevmode\nobreak\ \text{is mapped to}\leavevmode\nobreak\ (P,Q)\\ &\text{by the transformations \eqref{coordinate:polar} and \eqref{coordinate:rho_phi}}\}.\end{split}$$

$$(P,Q)=\arg \underset{(P,Q)\in {\mathcal{A}}_{{\Omega }^{\prime }}({\mathbf{u}}_0)}{\min }E(P,Q),$$

$$\{\begin{array}{r}{\int }_{{\Omega }^{\prime }}{f}_1(P,Q;\bar{P})\mathrm{d}\rho \mathrm{d}\varphi =0,\\ {\int }_{{\Omega }^{\prime }}{f}_2(P,Q;\bar{Q})\mathrm{d}\rho \mathrm{d}\varphi =0,\end{array}\forall \text{leavevmode}(\bar{P},\bar{Q})\in {\mathcal{A}}_{{\Omega }^{\prime }}(0),$$

$$f_1$$

$$+{\bar{P}}_{\rho }\left[g^{\prime }(F)r{\rho }_r{P}_{\rho }+h^{\prime }(D)P(Q_{\varphi }+1)\right]+{\bar{P}}_{\varphi }\left[g^{\prime }(F)\frac{1}{r{\rho }_r}{P}_{\varphi }-h^{\prime }(D)P{Q}_{\rho }\right],$$

$$f_2$$

$$P^{NM}(\rho ,\varphi )$$

$$Q^{NM}(\rho ,\varphi )$$

$$T_j(x)=\cos (j\arccos x),$$

$$T_{j+1}(x)=2x{T}_j(x)-T_{j-1}(x),j\ge 1.$$

$$\{\begin{array}{r}{\int }_{{\Omega }^{\prime }}{f}_1(P^{NM},Q^{NM};\bar{P})\mathrm{d}\rho \mathrm{d}\varphi =0,\\ {\int }_{{\Omega }^{\prime }}{f}_2(P^{NM},Q^{NM};\bar{Q})\mathrm{d}\rho \mathrm{d}\varphi =0,\end{array}\forall \text{leavevmode}(\bar{P},\bar{Q})\in {\mathcal{B}}_0^{NM},$$

$${\mathcal{B}}^{NM}$$

$$P^{NM}(1,{\varphi }_n)=P_0(1,{\varphi }_n),\text{leavevmode}{Q}^{NM}(1,{\varphi }_n)=Q_0(1,{\varphi }_n),\text{leavevmode}0\le n\le N-1\},$$

$${\mathcal{B}}_0^{NM}$$

$$P^{NM}(1,{\varphi }_n)=0,\text{leavevmode}{Q}^{NM}(1,{\varphi }_n)=0,\text{leavevmode}0\le n\le N-1\},$$

$$\displaystyle\begin{array}[]{lll}\rho_{m^{\prime}}=\cos\displaystyle\frac{(2m^{\prime}+1)\pi}{2M^{\prime}+2},&\omega^{C}_{m^{\prime}}=\displaystyle\frac{\pi}{M^{\prime}+1},&0\leq m^{\prime}\leq M^{\prime},\\ \phi_{n^{\prime}}=\displaystyle\frac{2\pi n^{\prime}}{N^{\prime}},&\omega^{F}_{n^{\prime}}=\displaystyle\frac{2\pi}{N^{\prime}},&0\leq n^{\prime}\leq N^{\prime}-1,\end{array}$$

$$\{\begin{array}{r}\sum _{n^{\prime }=0}^{N^{\prime }-1}\sum _{m^{\prime }=0}^{M^{\prime }}{f}_1(P^{NM}({\rho }_{m^{\prime }},{\varphi }_{n^{\prime }}),Q^{NM}({\rho }_{m^{\prime }},{\varphi }_{n^{\prime }});\bar{P}({\rho }_{m^{\prime }},{\varphi }_{n^{\prime }}))\sqrt{1-{\rho }_{m^{\prime }}^2}{\omega }_{m^{\prime }}^C{\omega }_{n^{\prime }}^F=0,\\ \sum _{n^{\prime }=0}^{N^{\prime }-1}\sum _{m^{\prime }=0}^{M^{\prime }}{f}_2(P^{NM}({\rho }_{m^{\prime }},{\varphi }_{n^{\prime }}),Q^{NM}({\rho }_{m^{\prime }},{\varphi }_{n^{\prime }});\bar{Q}({\rho }_{m^{\prime }},{\varphi }_{n^{\prime }}))\sqrt{1-{\rho }_{m^{\prime }}^2}{\omega }_{m^{\prime }}^C{\omega }_{n^{\prime }}^F=0.\end{array}$$

$$\displaystyle\begin{array}[]{cc}\begin{aligned} \alpha_{k,0}&=-\sum_{j=1}^{M}\alpha_{k,j}+a_{k},\\ \beta_{k,0}&=-\sum_{j=1}^{M}\beta_{k,j}+b_{k},\end{aligned}&\begin{aligned} \xi_{k,0}&=-\sum_{j=1}^{M}\xi_{k,j}+c_{k},\quad 0\leq k\leq N/2,\\ \eta_{k,0}&=-\sum_{j=1}^{M}\eta_{k,j}+d_{k},\quad 1\leq k\leq N/2-1.\end{aligned}\end{array}$$

$${\left\{\cos k\varphi \cdot (T_j(\rho )-1)\right\}}_{1\le j\le M}^{0\le k\le N/2},{\left\{\sin k\varphi \cdot (T_j(\rho )-1)\right\}}_{1\le j\le M}^{1\le k\le N/2-1},$$

$${\mathbf{B}}_{k+1}={\mathbf{B}}_k+\frac{\left({\mathbf{s}}_k-{\mathbf{B}}_k{\mathbf{z}}_k\right){\mathbf{s}}_k^T{\mathbf{B}}_k}{{\mathbf{s}}_k^T{\mathbf{B}}_k{\mathbf{z}}_k}.$$

$${\Vert \psi \Vert }_{H_{\omega }^{\sigma }(I)}={\left(\sum _{j=0}^{\sigma }{\int }_I{\left|\frac{d^j\psi }{d{\rho }^j}\right|}^2\omega \mathrm{d}\rho \right)}^{1/2},$$

Full text

Key words: Fourier-Chebyshev spectral method, cavitation, nonlinear elasticity, interpolation error analysis, energy error estimate, convergence.

1 Introduction

In 1958, Gent and Lindley [1] established the well known defective model for the cavitation in nonlinear elasticity characterizing the phenomenon as material instability associated to the dramatic growth of pre-existing micro voids under large hydrostatic tensions, which very well matched the experimental observation of sudden void formation in vulcanized rubber. Using the defective model, Gent et.al. [2], Lazzeri et.al. [3, 4], and many other researchers studied the cavitation phenomenon in elastomers containing rigid spherical inclusions as well as in the standard model problems. In 1982, Ball [5] established the famous perfect model, in which cavitations form in an originally intact body as an absolute energy minimizing bifurcation solution, and produced the same cavitation criterion. The profound relationship of the two models are studied by Sivaloganathan et.al. [6, 7] and Henao [8].

Since the perfect model is known to be seriously challenged by the Lavrentiev phenomenon [9], the defective model is chosen by most researchers in numerical studies of the cavitation phenomenon, using mainly a variety of the finite element methods (see Xu and Henao [10], Lian and Li [11, 12], Su and Li [13] among many others). A spectral collocation method [14], which approximates the cavitation solution with truncated Fourier series in the circumferential direction and finite differences in the radial direction, is also found some success.

In a typical 2-dimensional defective model with a prescribed displacement boundary condition, one considers to minimize the stored energy of the form

[TABLE]

in the set of admissible deformations

[TABLE]

where $\Omega_{\varepsilon}=\Omega\setminus\overline{\mathbb{B}_{\varepsilon}(\mathbf{x}_{0})}\subset\mathbb{R}^{2}$ is a domain occupied by the compressible hyperelastic material in its reference configuration, with $\Omega$ being a regular simply-connected domain and $\mathbb{B}_{\varepsilon}(\mathbf{x}_{0})=\{\mathbf{x}\in\mathbb{R}^{2}:|\mathbf{x}|<\varepsilon\}$ being a pre-existing circular defect of radius $\varepsilon\ll 1$ centered at $\mathbf{x}_{0}$, and where $W:\mathbb{M}_{+}^{2\times 2}\rightarrow\mathbb{R}^{+}$ is the stored energy density function of the hyperelastic material, and $\mathbb{M}_{+}^{2\times 2}$ denotes the set of $2\times 2$ matrices with positive determinant.

The Euler-Lagrange equation of the above minimization problem is the following displacement/traction boundary value problem:

[TABLE]

where $\mathbf{n}$ is the unit exterior normal with respect to $\Omega_{\varepsilon}$.

In the present paper, without loss of generality [5, 15], we consider the stored energy density function $W(\cdot)$ of the form

[TABLE]

where $\kappa$ is a positive material constant, $|\cdot|$ denotes the Frobenius norm of a matrix and $h\in\mathcal{C}^{3}((0,+\infty))$ is a strictly convex function satisfying

[TABLE]

Since the cavitation solution is generally considered to have high regularity except in a neighborhood of the defects, where the material experiences large expansion dominant deformations, we restrict ourselves to a simplified reference configuration $\Omega_{(\varepsilon,\gamma)}=\mathbb{B}_{\gamma}(\mathbf{0})\setminus\overline{\mathbb{B}_{\varepsilon}(\mathbf{0})}$ $(0<\varepsilon\ll\gamma\leq 1)$, and denote

[TABLE]

Taking the advantages of the smoothness of the cavitation solutions in the defective model when $\mathbf{u}_{0}$ is sufficiently smooth and the high efficiency and accuracy of spectral methods in approximating smooth solutions of partial differential equations (see Li and Guo [16], Shen [17, 18] etc.), we develop a Fourier-Chebyshev spectral method to solve the Euler-Lagrange equation (1.3), which approximates the cavitation solution with truncated Fourier series in the circumferential direction and truncated Chebyshev series in the radial direction. The interpolation error for the cavitation solution is analyzed, the elastic energy error estimate for the discrete cavitation solution is derived, and the convergence of the method is proved. An algorithm combined a gradient type method with a damped quasi-Newton method is applied to solve the discretized nonlinear equilibrium equations. Numerical experiments show that the Fourier-Chebyshev spectral method is efficient and capable of producing highly accurate numerical cavitation solutions. We would like to point out here, even though the reference domain is restricted to a circular ring $\Omega_{(\varepsilon,\gamma)}$, to further exploring its highly efficient feature in a neighborhood of a cavity surface, our method can be coupled with a domain decomposition method, especially in combining with some finite element methods to extend the application to more general situations with multiple pre-existing tiny voids.

The structure of the rest of the paper is as follows. In §2, we rewrite the Euler-Lagrange equation of the cavitation problem in a proper computing coordinates. In §3, the Fourier-Chebyshev spectral method is applied, the corresponding discrete equilibrium equation is derived, and an algorithm to solve the nonlinear equation is presented. §4 is devoted to the analysis of the interpolation error of the cavitation solution, the elastic energy error bound and the convergence of the discrete cavitation solution. In §5, numerical experiments and results are presented to show the efficiency and accuracy of our method.

2 The Euler-Lagrange Equation

In the Cartesian coordinate system, an admissible deformation $\mathbf{u}\in\mathcal{A}_{(\varepsilon,\gamma)}(\mathbf{u}_{0})$ is written as $\mathbf{u}(\mathbf{x})=\left[u_{1}(x_{1},x_{2}),u_{2}(x_{1},x_{2})\right]^{T}$. Denote

[TABLE]

and to futher simplify the notation, $D(\mathbf{u})$ and $F(\mathbf{u})$ will be denoted below as $D$, $F$ wherever no ambiguity is caused. For the elastic energy density function $W(\cdot)$ given by (1.4) and the elastic energy $E(\cdot)$ given by (1.1), we have

[TABLE]

For the convenience of the implementation of the Fourier-Chebyshev spectral method, we introduce a $(\rho,\phi)$-coordinate system defined on the computational domain $\Omega^{\prime}:=(-1,1)\times(0,2\pi)$, by coupling the Cartesian to polar coordinates transformation

[TABLE]

defined on the domain $(\varepsilon,\gamma)\times(0,2\pi)$, with a transformation defined by

[TABLE]

defined on the computational domain $\Omega^{\prime}=(-1,1)\times(0,2\pi)$.

In $(\rho,\phi)$-coordinates, $D(\mathbf{u})=\det\nabla\mathbf{u}$, $F(\mathbf{u})=|\nabla\mathbf{u}|^{2}/2$ defined in (2.1) can be rewritten as functions of $P(\rho,\phi)$, $Q(\rho,\phi)$:

[TABLE]

where $\rho_{r}=2/(\gamma-\varepsilon)$; the elastic energy $E(\mathbf{u})$ in (2.2) can be expressed as

[TABLE]

and the set of admissible deformation $\mathcal{A}_{(\varepsilon,\gamma)}(\mathbf{u}_{0})$ (see (1.6)) is reformulated as

[TABLE]

Thus, in $(\rho,\phi)$-coordinates, the cavitation solution $(P,Q)\in\mathcal{A}_{\Omega^{\prime}}(\mathbf{u}_{0})$ is characterized as the minimizer of $E(P,Q)$ in $\mathcal{A}_{\Omega^{\prime}}(\mathbf{u}_{0})$, i.e.

[TABLE]

or alternatively, as the solution to the Euler-Lagrange equation of (2.8):

[TABLE]

where, by the definition and direct calculations, we have

[TABLE]

3 The Fourier-Chebyshev Spectral Method

To discretize the Euler-Lagrange equation (2.9) defined on $\Omega^{\prime}=(-1,1)\times(0,2\pi)$ in $(\rho,\phi)$-coordinates, we first approximate the unknowns $(P(\rho,\phi),Q(\rho,\phi))$ by the finite Fourier-Chebyshev polynomials:

[TABLE]

where $T_{j}$ is the Chebyshev polynomial of the first kind of degree $j$, defined as

[TABLE]

with $T_{0}(x)=1,T_{1}(x)=x$ and satisfying the recurrence relation [18]

[TABLE]

Remark 1**.**

We use the trigonometric polynomials to approximate $Q=\Theta-\theta$ instead of $\Theta$ (see (2.4) and (3.1b)) so that the Gibbs phenomenon can be avoided (see e.g. [19]), since the periodic extension of $Q=\Theta-\theta$ from $[0,2\pi)$ to $\mathbb{R}^{1}$ is smooth, while that of $\Theta$ is a sawtooth function with jump discontinuities at $2k\pi$, $k=0,1,\dotsc$.

The discretized problem of solving the Euler-Lagrange equation (2.9) is then read as: find $(P^{NM},Q^{NM})\in\mathcal{B}^{NM}$ such that

[TABLE]

where $\mathcal{B}^{NM}$ and $\mathcal{B}^{NM}_{0}$ are the discrete trial and test function spaces defined as

[TABLE]

where, in (3.3), $\phi_{n}=2\pi n/N,\leavevmode\nobreak\ 0\leq n\leq N-1$, and the Dirichlet boundary condition $(P_{0}(1,\phi_{n}),Q_{0}(1,\phi_{n}))$ is defined by $\mathbf{u}_{0}$ via the coordinates transformations (2.3) and (2.4). To solve the equation (3.2) numerically, we need to replace the integrals in (3.2) by proper numerical quadratures. Let $\left\{\rho_{m^{\prime}},\omega^{C}_{m^{\prime}}\right\}_{m^{\prime}=0}^{M^{\prime}}$ and $\left\{\phi_{n^{\prime}},\omega^{F}_{n^{\prime}}\right\}_{n^{\prime}=0}^{N^{\prime}-1}$ be the sets of Gauss-Chebyshev and Fourier quadrature nodes and weights respectively, i.e. [18]

[TABLE]

then we are led to the following discretized Euler-Lagrange equation: find $(P^{NM},Q^{NM})\in\mathcal{B}^{NM}$ such that for all $(\bar{P},\bar{Q})\in\mathcal{B}^{NM}_{0}$

[TABLE]

Let $\{a_{k},b_{k}\}$ and $\{c_{k},d_{k}\}$ be the discrete Fourier coefficients of $P_{0}(1,\phi)$ and $Q_{0}(1,\phi)$ respectively, then the boundary condition in (3.3a) can be expressed as

[TABLE]

Noticing also that the following $N\times M$ functions

[TABLE]

form a set of bases for $\mathcal{B}^{NM}_{0}$, we conclude that the discrete Euler-Lagrange equation (3.7) consists of $2NM$ nonlinear algebraic equations, which, for the simplicity of the notations, will be denoted as $\mathbf{f}(\mathbf{y})=\mathbf{0}$, with $2NM$ unknowns $\mathbf{y}=$ $\left\{\alpha_{k,j},\xi_{k,j}\right\}_{1\leq j\leq M}^{0\leq k\leq N/2}\cup$ $\left\{\beta_{k,j},\eta_{k,j}\right\}_{1\leq j\leq M}^{1\leq k\leq N/2-1}$. Denote $E(\mathbf{y})$ as the discrete elastic energy defined by replacing the integral in $E(P^{NM},Q^{NM})$ (see (2.6)) with the numerical quadrature, then $\mathbf{f}(\mathbf{y})$ may be viewed as the gradient of the discrete elastic energy $E(\mathbf{y})$.

In our numerical experiments, the discrete equilibrium equations $\mathbf{f}(\mathbf{y})=\mathbf{0}$, i.e. (3.7), are solved by an algorithm combined a gradient type method with a damped quasi-Newton method [20]. More specifically, we use a gradient type method, which calculates a descent direction of the energy and conducts a incomplete line search in each iteration, to provide an appropriate initial cavity deformation for a damped quasi-Newton method with Broyden’s correction, which will then produce a reasonably accurate numerical cavity solution. The algorithm is summarized as follows, where the determinant of the deformation $(P^{NM},Q^{NM})$ corresponding to $\mathbf{y}$ is denoted as $D(\mathbf{y}):=D(P^{NM},Q^{NM})$ (see (2.5a)).

Algorithm:

Step 1

Given $\mathbf{y}^{G}_{0}$, set $TOL=10^{-1}$, compute $\mathbf{f}(\mathbf{y}^{G}_{0})$ and $E(\mathbf{y}^{G}_{0})$.

Step 2

If $TOL<10^{-10}$, then output $\mathbf{y}^{G}_{0}$ and stop; else, set $t^{G}_{-1}=1$ and $j:=0$.

Step 3

For $j\geq 0$, if $|\mathbf{f}(\mathbf{y}^{G}_{j})|<TOL$, then go to Step 6; else, set $t^{G}_{j}=4\cdot t^{G}_{j-1}$.

Step 4

Set $\mathbf{y}^{G}_{j+1}=\mathbf{y}^{G}_{j}-t^{G}_{j}\cdot\mathbf{f}(\mathbf{y}^{G}_{j})$, compute $\mathbf{f}(\mathbf{y}^{G}_{j+1})$, $E(\mathbf{y}^{G}_{j+1})$ and $D(\mathbf{y}^{G}_{j+1})$.

Step 5

If $t^{G}_{j}<10^{-16}$, then output $\mathbf{y}^{G}_{j}$ and stop; else if $E(\mathbf{y}^{G}_{j+1})<E(\mathbf{y}^{G}_{j})$ and $D(\mathbf{y}^{G}_{j+1})>0$, then set $j:=j+1$ and go to Step 3; else, set $t^{G}_{j}:=t^{G}_{j}/2$ and go to Step 4.

Step 6

Set $\mathbf{y}^{N}_{0}:=\mathbf{y}^{G}_{j}$, compute $\mathbf{f}(\mathbf{y}^{N}_{0})$ and $\mathbf{B}_{0}=[\nabla\mathbf{f}(\mathbf{y}^{N}_{0})]^{-1}$, set $t^{N}_{-1}=1$ and $k:=0$.

Step 7

For $k\geq 0$, if $|\mathbf{f}(\mathbf{y}^{N}_{k})|<10^{-10}$, then output $\mathbf{y}^{N}_{k}$ as the solution and stop; else, set $t^{N}_{k}=4\cdot t^{N}_{k-1}$.

Step 8

Set $\mathbf{y}^{N}_{k+1}=\mathbf{y}^{N}_{k}-t^{N}_{k}\ cdot\mathbf{B}_{k}\cdot\mathbf{f}(\mathbf{y}^{N}_{k})$, compute $\mathbf{f}(\mathbf{y}^{N}_{k+1})$ and $D(\mathbf{y}^{N}_{k+1})$.

Step 9

If $t^{N}_{k}<10^{-16}$, then go to Step 2 with $\mathbf{y}^{G}_{0}:=\mathbf{y}^{N}_{k}$ and $TOL:=TOL/10$; else if $|\mathbf{f}(\mathbf{y}^{N}_{k+1})|<|\mathbf{f}(\mathbf{y}^{N}_{k})|$ and $D(\mathbf{y}^{N}_{k+1})>0$, then go to Step 10; else, set $t^{N}_{k}:=t^{N}_{k}/2$ and go to Step 8.

Step 10

Compute $\mathbf{s}_{k}=\mathbf{y}^{N}_{k+1}-\mathbf{y}^{N}_{k}$, $\mathbf{z}_{k}=\mathbf{f}(\mathbf{y}^{N}_{k+1})-\mathbf{f}(\mathbf{y}^{N}_{k})$, and

[TABLE]

Set $k:=k+1$ and go to Step 7.

4 Error Analysis and the Convergence Theorem

In this section, we analyze the interpolation error of the discrete Fourier-Chebyshev spectral method for the cavitation solutions, which will enable us to derive the elastic energy error estimate for the discrete cavitation solution, and prove the convergence of the method.

Before analyzing the interpolation error of a cavitation solution, we first introduce some notations. Let $\mathcal{B}:=\{(P,Q)\in C^{1}(\overline{\Omega^{\prime}}):P(1,\phi)=P_{0}(1,\phi),\leavevmode\nobreak\ Q(1,\phi)=Q_{0}(1,\phi)\}$, let $\mathcal{B}_{+}:=\{(P,Q)\in\mathcal{B}:D(P,Q)>0\}$ (see (2.5a)), and denote $\mathcal{B}^{NM}_{+}=\mathcal{B}^{NM}\cap\mathcal{B}_{+}$ (see (3.3a)). Let $\omega(\rho):=(1-\rho^{2})^{-1/2}$, $\Lambda:=(0,2\pi)$, $I:=(-1,1)$, and recall that $(\rho,\phi)\in\Omega^{\prime}=I\times\Lambda$. For given integers $\sigma\geq 0$ and $\mu\geq 0$, denote $H^{\sigma}_{\omega}(I)=\{\psi:\left\lVert\psi\right\rVert_{H^{\sigma}_{\omega}(I)}<\infty\}$ the weighted Hilbert space with the norm defined as

[TABLE]

and denote $H^{\mu}(\Lambda;H^{\sigma}_{\omega}(I))=\{v:\left\lVert v\right\rVert_{H^{\mu}(\Lambda;H^{\sigma}_{\omega}(I))}<\infty\}$ the Hilbert space equipped with the norm defined as

[TABLE]

Definition 1**.**

Define the interpolation operator $I^{NM}:\mathcal{B}\rightarrow\mathcal{B}^{NM}$ as

[TABLE]

where $\rho_{m}=\cos(m\pi/M)$, $\phi_{n}=2n\pi/N$.

The interpolation operator $I^{NM}$ is shown to have the following error estimates (see Lemma 5 in [16]).

Lemma 1**.**

[16*]**

If $v\in H^{\beta}(\Lambda;H^{\sigma}_{\omega}(I))\cap H^{\mu}(\Lambda;H^{\alpha}_{\omega}(I))\cap H^{\mu^{\prime}}(\Lambda;H^{\sigma^{\prime}}_{\omega}(I)),\leavevmode\nobreak\ 0\leq\alpha\leq\sigma,\sigma^{\prime},\leavevmode\nobreak\ 0\leq\beta\leq\mu,\mu^{\prime},\leavevmode\nobreak\ \sigma,\sigma^{\prime}>\frac{1}{2}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \mu,\mu^{\prime}>1$, then, there exists a constant $c>0$ independent of $P$, $Q$, $M$ and $N$, such that*

[TABLE]

where $q(\beta)=0$ for $\beta>1$ and $q(\beta)=1$ for $\beta\leq 1$.

Theorem 1**.**

Let $(P,Q)\in\mathcal{B}\cap H^{s}(\Lambda;H^{l}_{\omega}(I))$, $l>2$, $s>1$. Then there exists a constant $c>0$ independent of $P$, $Q$, $M$ and $N$, such that

[TABLE]

where $\alpha,\beta=0,1$ and $\|\cdot\|_{*}$ is a norm defined by:

[TABLE]

Furthermore, if $(P,Q)\in\mathcal{B}\cap H^{l}(\Lambda;H^{l}_{\omega}(I))$ with $l>6$, then

[TABLE]

Proof.

The first half of the theorem is a direct consequence of Lemma 1 by setting $\sigma=\sigma^{\prime}=l$, $\mu=\mu^{\prime}=s$, and taking $\alpha=0$ or $1$ and $\beta=0$ or $1$ respectively.

By (2.5a), the error estimate (LABEL:det_error) follows from Lemma 1 with $\sigma=\sigma^{\prime}=\mu=\mu^{\prime}=l$, $\alpha=\beta=3$, and the fact that $H^{3}_{\omega}(\Omega^{\prime})\hookrightarrow H^{3}(\Omega^{\prime})\hookrightarrow C^{1}(\overline{\Omega^{\prime}})$. ∎

In what follows below, we always assume that, for a cavitation solution $(P,Q)$, the following hypotheses hold:

(H1)

$(P,Q)\in\mathcal{B}_{+}\cap H^{l}(\Lambda;H^{l}_{\omega}(I))$ with $l>6$ and $(P,Q)$ is the energy minimizer in $\mathcal{B}_{+}$.

(H2)

there exists constants $c_{F}>1$ and $c_{D}>1$ such that $F$ and $D$ (see (2.5)) satisfies

[TABLE]

Remark 2**.**

Notice that, by (2.5)

[TABLE]

and for a cavitation solution $P\geq P_{0}>0$, and in the radially symmetric case $Q_{\phi}=0$. the hypothesis $c_{F}^{-1}\leq 2r^{2}F$ is not too harsh a requirement on a general solution. While the other bounds are the direct consequences of $(P,Q)\in C^{1}(\overline{\Omega^{\prime}})$.

To estimate the error on the elastic energy of the interpolation function of a cavitation solution, we will making use of an auxiliary grid in radial direction on which the elastic energy of the cavitation solution is radially quasi-equi-distributed in the sense given in Lemma 3. The properties of such grids are given by the following two lemmas.

Lemma 2**.**

Let $(P,Q)\in\mathcal{B}_{+}$. Let $D$ and $F$ be defined by (2.5). Then, there exists a constant $c\geq 1$ such that, for all grid $\varepsilon=r_{0}<r_{1}<\cdots<r_{K}=\gamma$, the elastic energy

[TABLE]

satisfies

[TABLE]

Proof.

It follows from the convexity of $h(\cdot)$ and the hypothesis (H1) that

[TABLE]

Since $r_{i-1}<\gamma\leq 1$ and $1<p<2$, the hypothesis (H2) implies

[TABLE]

Hence, the conclusion (4.3) follows by taking

[TABLE]

∎

Lemma 3**.**

Let $K\gg\varepsilon^{-1}$ be a sufficiently large integer. Let $\varepsilon=r_{0}<r_{1}<\cdots<r_{K}=\gamma$ be a given grid satisfying

[TABLE]

where $\tau_{i}:=r_{i}-r_{i-1}$, $1\leq i\leq K$. Then, we have

[TABLE]

and the energy is radially quasi-equi-distributed on the grid, i.e.

[TABLE]

where $c$ is the same constant in Lemma 2.

Proof.

By (4.4), we have

[TABLE]

where denote $\Upsilon_{i}:=\frac{\tau_{i}}{r_{i-1}}>0$, $1\leq i\leq K-1$. Notice that, by definition,

[TABLE]

i.e. $\Upsilon_{i}$ is a strictly deceasing function of $i$. On the other hand, $\frac{\tau_{i+1}}{\tau_{i}}=(1+\Upsilon_{i})^{p-1}>1$ implies that $\tau_{i}$ is a strictly increasing function of $i$, and as a consequence, we have $K\Upsilon_{1}\varepsilon=K\tau_{1}<\sum_{i=1}^{K}\tau_{i}=\gamma-\varepsilon<K\tau_{K}$, which yields, since $K\gg\varepsilon^{-1}$,

[TABLE]

By (4.7) and (4.8), we also have, for all $1\leq i\leq K-1$,

[TABLE]

Now, express the left-hand side of (4.5) as

[TABLE]

For the first term on the right hand side of (4.10), by (4.7), (4.8) and (4.9), we have

[TABLE]

Since $r_{K}/r_{K-1}<2$ (see (4.4)), by (4.9) and (4.10), this yields (4.5).

Next, it follows from (4.3), (4.4) and $\Upsilon_{i}\ll 1$, $\forall\leavevmode\nobreak\ i$ that

[TABLE]

This proves the second inequality of (4.6).

Notice that, by (4.8) and (4.9), we have

[TABLE]

Thus, by (4.3) and (4.4),

[TABLE]

Denote $E_{max}=\max_{1\leq i\leq K}E_{(r_{i-1},r_{i})}$, then $2^{p-1}c^{2}KE_{(r_{i-1},r_{i})}>KE_{max}\geq E_{(\varepsilon,\gamma)}$. This proves the first inequality of (4.6). ∎

Remark 3**.**

For $K$ sufficiently large, it is not difficult to show that there exists an auxiliary grid $\varepsilon=r_{0}<r_{1}<\cdots<r_{K}=\gamma$ such that (4.4) holds.

The theorem below gives the relative and absolute errors of the elastic energy $E(P,Q)=E_{(\varepsilon,\gamma)}$ when a cavitation solution $(P,Q)$ is replaced by its interpolation functions $(I^{NM}P,I^{NM}Q)$.

Theorem 2**.**

Let $(P,Q)$ be a cavitation solution satisfying the hypotheses (H1) and (H2). Then, for $M$, $N$ sufficiently large, there exists a constant $C$ such that

[TABLE]

Proof.

To simplify the notation, we denote (see (2.5))

[TABLE]

By (2.1) and (2.6), the energy error can be bounded as follows

[TABLE]

By the hypotheses (H1), (H2), and as a consequence of Theorem 1, we have $(P,Q)$, $(\widetilde{P},\widetilde{Q})$ and their first order derivatives are all bounded, and

[TABLE]

where $\|\cdot\|_{\omega,\Omega^{\prime}}$ is the weighted $L^{2}$-norm on $\Omega^{\prime}$.

Thus, by (4.2) and recalling $\rho_{r}=2/(\gamma-\varepsilon)$, we have

[TABLE]

and as a consequence, it follows from (4.12) that

[TABLE]

By hypothesis (H1), (H2) and Theorem 1, both $D>0$ and $\widetilde{D}>0$ are bounded away from [math] and $+\infty$, hence, by (4.13a), we have

[TABLE]

where $\vartheta_{1}$ is between $\widetilde{D}$ and $D$, and thus $h^{\prime}(\vartheta_{1})$ is bounded.

On the other hand, let $\varepsilon=r_{0}<r_{1}<\cdots<r_{K}=\gamma$ be an auxiliary grid in radial direction satisfying the conditions of Lemma 3, then, by (4.3) and (4.6), we have

[TABLE]

Since hypotheses (H1), (H2) and Lemma 1 implies that both $r^{2}F$ and $r^{2}\widetilde{F}$ are bounded, by the Hölder inequality, we have

[TABLE]

where $\vartheta_{2}$ is between $2r^{2}\widetilde{F}$ and $2r^{2}F$, and thus $|\vartheta_{2}^{p/2-1}|$ is also bounded. Substituting this into (4) and applying the Hölder inequality, we have, by (4.5) and (4.13b),

[TABLE]

The proof is completed by combining(4), (4) and (4). ∎

Notice that $\mathcal{B}_{+}^{NM}\subset\mathcal{B}_{+}$, the result of Theorem 2 allows us to obtain the following elastic energy error estimate for the discrete cavitation solution.

Theorem 3**.**

Let a cavitation solution $(P,Q)$ satisfy the hypotheses (H1), (H2) and be a global energy minimizer of $E(\cdot,\cdot)$ in $\mathcal{B}_{+}$. Let $(P^{NM},Q^{NM})$ be a global energy minimizer of $E(\cdot,\cdot)$ in $\mathcal{B}_{+}^{NM}$. Then, for $M$, $N$ sufficiently large, there exists a constant $C>0$ such that

[TABLE]

Proof.

The first inequality is a direct consequence of $\mathcal{B}_{+}^{NM}\subset\mathcal{B}_{+}$. By (LABEL:det_error) and for $M$, $N$ sufficiently large, we have $(I^{NM}P,I^{NM}Q)\in\mathcal{B}_{+}^{NM}$. Then, the second inequality follows from Theorem 2 and $E(P^{NM},Q^{NM})\leq E(I^{NM}P,I^{NM}Q)$. ∎

Let $(P,Q)\in\mathcal{B}_{+}$ be a global energy minimizer of $E(\cdot,\cdot)$ in $\mathcal{B}_{+}$, let $(I^{NM}P,I^{NM}Q)\in\mathcal{B}_{+}^{NM}$ be its interpolation functions. Let $(P^{NM},Q^{NM})\in\mathcal{B}_{+}^{NM}$ be global energy minimizers of $E(\cdot,\cdot)$ in $\mathcal{B}_{+}^{NM}$. Denote $\mathbf{\bar{u}}$, $\mathbf{\bar{u}}^{NM}$ and $\mathbf{U}^{NM}$ as the corresponding functions on $\Omega_{(\varepsilon,\gamma)}$ defined by $(P,Q)$, $(I^{NM}P,I^{NM}Q)$ and $(P^{NM},Q^{NM})$ respectively via the coordinates transformations (2.3) and (2.4). Then, (4.17) can be rewritten as

[TABLE]

According to Theorem 4.9 in [21] and its proof, for a conforming discrete approximation method of the cavitation problem, the inequality (4.18) implies the convergence of the discrete cavitation solutions. Thus, we have the following convergence theorem. For the convenience of the readers, we sketch its proof below.

Theorem 4**.**

Let a cavitation solution $(P,Q)$ satisfy the hypotheses (H1), (H2) and be a global energy minimizer of $E(\cdot,\cdot)$ in $\mathcal{B}_{+}$. Let $(P^{NM},Q^{NM})$ be global energy minimizers of $E(\cdot,\cdot)$ in $\mathcal{B}_{+}^{NM}$. Let $\mathbf{U}^{NM}$ correspond to $(P^{NM},Q^{NM})$ under the transformations (2.3) and (2.4). Then, there exists a subsequence, still denoted as $\{\mathbf{U}^{NM}\}$, and a function $\mathbf{u}\in\mathcal{A}_{(\varepsilon,\gamma)}(\mathbf{u}_{0})$, such that $\mathbf{U}^{NM}\rightarrow\mathbf{u}$ in $W^{1,p}(\Omega_{(\varepsilon,\gamma)})$ and $\mathbf{u}$ is a global energy minimizer of $E(\cdot)$ in $\mathcal{A}_{(\varepsilon,\gamma)}(\mathbf{u}_{0})$.

Proof.

Since $h>0$ is a convex function satisfying the growth conditions (1.5), $1<p<2$ and $\mathbf{U}^{NM}$ satisfies the Direchlet boundary condition, by (4.18) and the De La Vallée Poussin theorem [23], we conclude that $\{|\nabla\mathbf{U}^{NM}|^{p}\}_{N,M\rightarrow\infty}$ and $\{\det\nabla\mathbf{U}^{NM}\}_{N,M\rightarrow\infty}$ are equi-integrable. As a consequence, there exists a subsequence, still denoted as $\{\mathbf{U}^{NM}\}$, a function $\mathbf{u}\in W^{1,p}(\Omega_{(\varepsilon,\gamma)})$ and a function $\zeta\in L^{1}(\Omega_{(\varepsilon,\gamma)})$ such that

[TABLE]

Hence, by $\det\nabla\mathbf{U}^{NM}>0$, $a.e.$, we have $\zeta\geq 0$, $a.e.$. We conclude that $\zeta>0$, $a.e.$. Suppose otherwise, i.e. $\zeta=0$ on a set $S$ with positive measure, then there exists a subsequence, still denoted as $\{\mathbf{U}^{NM}\}$, such that $\int_{S}|\det\nabla\mathbf{U}^{NM}|\mathrm{d}\mathbf{x}\rightarrow 0$ and $\det\nabla\mathbf{U}^{NM}\rightarrow 0$, $a.e.$ on the set $S$, which, by (1.5), implies $h(\det\nabla\mathbf{U}^{NM})\rightarrow\infty$, $a.e.$ on the set $S$. Thus, by the Fatou lemma, we have and $E(\mathbf{U}^{NM})\rightarrow\infty$, which contradicts to $\displaystyle\varlimsup_{N,M\rightarrow\infty}E(\mathbf{U}^{NM})<\infty$.

Thanks to Theorem 3 in [25] and Theorem 3 in [24], as a consequence of (4.19), $\zeta>0$, $a.e.$ and the continuity of $\mathbf{U}^{NM}$, we have $\zeta=\det\nabla\mathbf{u}$, $a.e.$ and $\mathbf{u}$ is one-to-one a.e.. In addition, it is easily verified that $\mathbf{u}|_{\partial\Omega}=\mathbf{u}_{0}$. Hence $\mathbf{u}\in\mathcal{A}_{(\varepsilon,\gamma)}(\mathbf{u}_{0})$. On the other hand, since $\displaystyle E(\mathbf{u})\leq\varliminf_{N,M\rightarrow\infty}E(\mathbf{U}^{NM})$, due to the weakly lower semi-continuity of $E(\cdot)$ on $W^{1,p}(\Omega_{(\varepsilon,\gamma)})$ (see the theorem 5.4 in [26]), we conclude from $\mathbf{u}\in\mathcal{A}_{(\varepsilon,\gamma)}(\mathbf{u}_{0})$ and (4.18) that $\mathbf{u}$ is a global minimizer of $E(\cdot)$ in $\mathcal{A}_{(\varepsilon,\gamma)}(\mathbf{u}_{0})$ and $\displaystyle E(\mathbf{u})=\lim_{N,M\rightarrow\infty}E(\mathbf{U}^{NM})$.

Since $h$ is a convex function, it follows from $\det\nabla\mathbf{U}^{NM}\rightharpoonup\det\nabla\mathbf{u}$ in $L^{1}(\Omega_{(\varepsilon,\gamma)})$ that

[TABLE]

which leads to $\displaystyle\varlimsup_{N,M\rightarrow\infty}\left\lVert\mathbf{U}^{NM}\right\rVert_{W^{1,p}(\Omega_{(\varepsilon,\gamma)})}\leq\left\lVert\mathbf{u}\right\rVert_{W^{1,p}(\Omega_{(\varepsilon,\gamma)})}$. Thus, by the uniform convexity of $W^{1,p}(\Omega_{(\varepsilon,\gamma)})$ (see [27]) and $\mathbf{U}^{NM}\rightharpoonup\mathbf{u}$ in $W^{1,p}(\Omega_{(\varepsilon,\gamma)})$, it follows from proposition 3.30 in [28] that $\mathbf{U}^{NM}\rightarrow\mathbf{u}$ in $W^{1,p}(\Omega_{(\varepsilon,\gamma)})$. This completes the proof. ∎

5 Numerical Experiments and Results

In our numerical experiments, the stored energy density function $W(\cdot)$ is taken of the form (1.4) with

[TABLE]

The reference configuration is $\Omega_{(\varepsilon,\gamma)}=\mathbb{B}_{\gamma}(\mathbf{0})\setminus\overline{\mathbb{B}_{\varepsilon}(\mathbf{0})}$, $(0<\varepsilon\ll\gamma\leq 1)$. We consider $\mathbf{u}_{0}(\mathbf{x})=\lambda\mathbf{x},\leavevmode\nobreak\ \mathbf{x}\in\partial\mathbb{B}_{\gamma}(\mathbf{0})$, $\lambda>1$ in the radially-symmetric case and $\mathbf{u}_{0}(\mathbf{x})=\left[\lambda_{1}x_{1},\lambda_{2}x_{2}\right]^{T}$, $\mathbf{x}\in\partial\mathbb{B}_{\gamma}(\mathbf{0})$, $\lambda_{1},\lambda_{2}>1$ in the non-radially-symmetric case.

By (LABEL:det_error) and the hypotheses (H1), (H2), we expect to have $(I^{NM}P,I^{NM}Q)\in\mathcal{B}_{+}^{NM}$ for sufficiently large $N$ and $M$. Before proceeding to the numerical experiments, we first check in Table 1 the orientation preservation condition $D(I^{NM}P,I^{NM}Q)>0$ for the exact cavitation solution $(P,Q)$ in the radially-symmetric case for incompressible elastic materials, since in such a case the cavity solutions have a simple explicit form $R(r)=\sqrt{\lambda^{2}+r^{2}-\gamma^{2}}$ in the polar coordinate systems. By $D(I^{NM}P,I^{NM}Q)=\frac{\rho_{r}}{r}\cdot I^{NM}P\cdot(I^{NM}P)_{\rho}$, we only need to check whether $(I^{NM}P)_{\rho}>0$ is satisfied. In fact, whenever $(I^{NM}P)_{\rho}>0$ is satisfied, we have $D(I^{NM}P,I^{NM}Q)\approx 1$. The data shown in Table 1 suggest that the orientation preservation condition $D>0$ should not impose much real additional restrictions on the choice of $N$ and $M$ in practical computations.

Next, we investigate the effect of the number of quadrature points $N^{\prime},M^{\prime}$ used in (3.7). Figure 1 shows the convergence behavior of the errors on the cavity radius for various $N^{\prime}/N$ (with $M=32$ and $M^{\prime}=8M$ fixed) and $M^{\prime}/M$ (with $N=16$ and $N^{\prime}=2N$ fixed), where for the non-radially-symmetric case $\Omega_{(\varepsilon,\gamma)}=\Omega_{(10^{-4},1)}$ with $\lambda_{1}=2.4$ and $\lambda_{2}=2$, and for the radially-symmetric case $\Omega_{(\varepsilon,\gamma)}=\Omega_{(10^{-2},1)}$ with $\lambda=2$. To balance the accuracy and computational cost, we set in our numerical experiments below $N^{\prime}=2N$ and $M^{\prime}=8M$.

5.1 Radially-Symmetric Case

In the radially-symmetric case with $\mathbf{u}_{0}(\mathbf{x})=\lambda\mathbf{x}$, $\lambda>1$, the cavitation solution $\mathbf{u}$ can be written in polar coordinates systems as

[TABLE]

where $s(r)$ satisfies $s(\varepsilon)>0$ and $s(\gamma)=\lambda\cdot\gamma$. Since in theory the numerical solution is independent of the circumferential DOF (degree of freedom) $N$, we fix $N=16$ and examine the effect of the radial DOF $M$ on the numerical performance of our method. In comparison, high precision numerical solutions to the equivalent 1-dimensional ODE boundary value problems [29], obtained by the ode15s routine in MATLAB with the tolerance $10^{-16}$, are taken as the exact solutions.

For the standard case of $\gamma=1$, $\lambda=2$, and $\varepsilon=10^{-3}$, $10^{-4}$, the convergence behavior of our numerical cavitation solutions $\mathbf{U}^{NM}$ is shown in Figure 2(a), 2(b), where $L^{2}_{\omega}$ and $W^{1,p}$ represent $L^{2}_{\omega}(\Omega^{\prime})$-norm and $W^{1,p}(\Omega_{(\varepsilon,\gamma)})$-semi-norm respectively.

For $\gamma=1$ and $\varepsilon=10^{-2},10^{-3},10^{-4}$, we show in Figure 4 the numerical results obtained with $M=32$ on the $\lambda$-$R^{NM}(\varepsilon)$ (i.e. the expansion rate on the outer boundary against the cavity radius on the inner boundary) graph. Figure 4 shows, for $\varepsilon=10^{-4}$, the convergence of our numerical results to that of the 1-dimensional ODE solution.

To explore the potential of the method in coupling with a domain decomposition method, especially when combining with a finite element method in a multi-defects problem, we examine the convergence behavior of our algorithm on a small neighbourhood of the defect. Taking $\Omega_{(10^{-3},10^{-1})}$, $\Omega_{(10^{-4},10^{-2})}$ and $\Omega_{(10^{-5},10^{-3})}$ as the reference configurations and setting $\lambda\cdot\gamma=1.25$, we show in Figure 5 the energy error and cavity radius error as a function of $M$ (with $N=16$ fixed), where it is clearly seen that high precision numerical results can be obtained with rather small $M$.

5.2 Non-radially Symmetric Case

For the non-radially symmetric case, we consider the circular ring reference configuration $\Omega_{(\varepsilon,\gamma)}$ with oval boundary stretch $\mathbf{u}_{0}(\mathbf{x})=\left[\lambda_{1}x_{1},\lambda_{2}x_{2}\right]^{T}$, $\lambda_{1}$, $\lambda_{2}>1$. Assuming that the error of the numerical solution $\mathbf{U}^{NM}$ satisfies

[TABLE]

where $q^{\infty}$ and $q^{NM}$ represent the exact and numerical results of a specific quantity, such as the elastic energy, semi-major axis and semi-minor axis etc., and $c_{1}$, $c_{2}$, $v_{1}$, $v_{2}$ are the corresponding parameters to be determined by the least squares data fitting.

For $\Omega_{(\varepsilon,\gamma)}=\Omega_{(10^{-4},1)}$, $\lambda_{1}=2.4$ and $\lambda_{2}=2$, we show in Figure 6(a) the errors between $\mathbf{U}^{N}$ and $\mathbf{U}^{1.5N}$ with $M=32$ fixed, and in Figure 6(b) the errors between $\mathbf{U}^{M}$ and $\mathbf{U}^{1.25M}$ with $N=16$ fixed, where $L^{2}_{\omega}$ and $W^{1,p}$ represent $L^{2}_{\omega}(\Omega^{\prime})$-norm and $W^{1,p}(\Omega_{(\varepsilon,\gamma)})$-semi-norm respectively. The regressed quantities and parameters are shown in Table 2.

As a comparison, we show in Figure 7 the corresponding errors obtained in the same way for the radially-symmetric case with $\lambda=2$, and show in Table 3 the regressed quantities and parameters. It is clearly seen that the regressed formula (5.2) produces quite sharp numerical results in the radially-symmetric case.

To see how well the regressed formula (5.2) fits the data, we show in Figure 8(a) the errors on the cavity dimensions, i.e. the cavity radius in the radially-symmetric case and the cavity major and minor axes in the non-radially-symmetric case, and in 8(b) the errors on the elastic energy, between the corresponding quantities produced by the numerical solution $\mathbf{U}^{NM}$ and the regressed formula (5.2) respectively. In particular, compare also to Figure 2(b), it is clearly seen that the regressed formula (5.2) is highly accurate and reliable.

To examine how the axial ratio of the oval stretch affect the critical displacement, we show in Figure 9 the semi-major and semi-minor axes of the numerical cavity formed as functions of $\lambda_{1}$ for $\lambda_{1}/\lambda_{2}=1.2,1.3,1.4$, with $\Omega_{(\varepsilon,\gamma)}=\Omega_{(10^{-4},1)}$, $N=16$ and $M=32$, where it is obviously seen that both are monotonously increasing functions.

Taking $\Omega_{(\varepsilon,\gamma)}=\Omega_{(10^{-3},10^{-1})}$, $\Omega_{(10^{-4},10^{-2})}$, $\Omega_{(10^{-5},10^{-3})}$ as the reference configuration respectively, setting the oval boundary data $\lambda_{1}\gamma=1.67$ and $\lambda_{2}\gamma=1.42$ (see Table 2), for $N=16$ fixed, the errors of the corresponding non-radially-symmetric numerical solutions $\mathbf{U}^{NM}$ are shown in with Figure 10, where the numerical solution with $M=32$ is taken as the exact solution. It is clearly seen that high precision numerical results can also be obtained in a neighborhood of a defect with rather small $M$ in non-radially-symmetric case.