Axisymmetric necking versus Treloar–Kearsley instability in a hyperelastic sheet under equibiaxial stretching

Abstract

We consider bifurcations from the homogeneous solution of a circular or square hyperelastic sheet that is subjected to equibiaxial stretching under either force- or displacement-controlled edge conditions. We derive the condition for axisymmetric necking and show, for the class of strain-energy functions considered, that the critical stretch for necking is greater than the critical stretch for the Treloar–Kearsley (TK) instability and less than the critical stretch for the limiting-point instability. An amplitude equation for the bifurcated necking solution is derived through a weakly nonlinear analysis and is used to show that necking initiation is generally sub-critical. Abaqus simulations are conducted to verify the bifurcation conditions and the expectation that the TK instability should occur first under force control, but when the edge displacement is controlled, the TK instability is suppressed, and it is the necking instability that will be observed. It is also demonstrated that axisymmetric necking follows a growth/propagation process typical of all such localization problems.

1. Introduction

Necking is a phenomenon commonly observed in the tension test of ductile metals. It has also been observed in soft materials such as semi-crystalline and glassy polymers [1–6] and hydro-gels [7–9]. Numerical simulations have been conducted [10–14], and there exists a large body of literature devoted to the analysis of necking based on one-dimensional (1D) models (see, for instance, [15–27]).

Under the framework of nonlinear elasticity, the initiation condition for necking was traditionally derived as the condition for a periodic mode satisfying the end conditions to appear, in the same manner as for a buckling mode under compression (see, for instance, [28–33]). However, there also exists an alternative point of view based on the dynamical systems theory whereby localization is viewed to initiate from a “zero wavenumber” mode [34–36]. For a simple and purely algebraic demonstration of how a localized mode bifurcates from a uniform state, we refer to Fu et al. [37]. Although the critical loads associated with the two approaches are usually close to each other, the post-buckling behaviors predicted by the two approaches become qualitatively different. Whereas the “zero wavenumber” point of view predicts a post-bifurcation profile that is more and more localized away from the bifurcation point, exactly as observed in numerical simulations and experiments, the other approach predicts a periodic profile that has a vanishingly small range of validity [38]. The “zero wavenumber” approach has recently been applied to study elasto-capillary necking and has been shown to give predictions in excellent agreement with numerical simulations [39–41]. This approach has previously been applied to study localized bulging in inflated rubber tubes (see [42–44], and the references therein).

The main motivation for this study is to provide insight into necking of a dielectric membrane that is subjected to the combined action of an electric field and mechanical forces. This phenomenon has been addressed in the papers by Puglisi and Zurlo [45], Zurlo et al. [46], Fu et al. [47,48], but has otherwise received very little attention in the vast literature on electroelasticity. Instead, most literature on stability and bifurcation has focused on wrinkling [49–55] or marginal violation of the Hessian criterion [56–59].

To focus on the explanation of ideas without the extra complications arising from the incorporation of electric effect, we consider in this paper the purely mechanical case when a circular hyperelastic sheet is subjected to an all-round tension along its curved edge. The resulting trivial primary deformation is the same as in a square sheet that is subjected to equibiaxial tension. Such a square sheet is known to be susceptible to both Treloar–Kearsley (TK) and limiting-point instabilities. The former refers to the existence of deformation field in which the two in-plane stretches become unequal although the associated tensions are equal, whereas the latter refers to the fact that the tension as a function of the stretch on the trivial loading path may reach a maximum. The TK instability was first noted by Kearsley [60] based on the experimental results of Treloar [61], with its plane-strain counterpart analyzed slightly earlier by Ogden [62]. It has subsequently been studied by MacSithigh [63], Ogden [64], Chen [65], MacSithigh and Chen [66], Batra et al. [67], and Steigmann [68]. However, our focus in this study will be on a third instability, that associated with axisymmetric necking. Based on our previous studies on localized bulging of an inflated tube, it is tempting to conjecture that the bifurcation condition for necking might be that the Jacobian determinant of the two stretching forces, as functions of the two corresponding stretches, is zero on the trivial loading path (i.e., when the two in-plane stretches are equal). However, it turns out that the Jacobian determinant equal to zero may always be factorized into two separate equations, and each equation may have zero, one, or two roots for the stretch depending on the material model used. When each equation has at least one root (which is the case for the class of material models considered), the two smallest roots correspond to the TK and limiting-point instabilities, respectively. Guided by our previous studies on localized bulging of inflated rubber tubes, we postulate that localized necking may occur when a static axially symmetric extensional mode may exist in the limit of zero wavenumber (the extensional mode is the mode whose out-of-plane displacement is symmetric with respect to the mid-plane). The associated critical stretch then lies between the two smallest roots of the Jacobian determinant, and this is verified by Abaqus simulations.

The rest of this paper is organized as follows. After summarizing the governing equations in the next section, we present a linear analysis in section 3 where we derive the bifurcation conditions for the static and axially symmetric, extensional and flexural modes, respectively. Each bifurcation condition is an implicit equation relating the stretch to the wavenumber in the radial direction. As for localized bulging in inflated rubber tubes, the initiation condition for necking is postulated to be the condition for the extensional mode with zero wavenumber to exist. This condition is further discussed in relation to the TK and limiting-point instabilities. In section 4, we conduct a weakly nonlinear analysis and derive the amplitude equation in a small neighborhood of the bifurcation point and show that the necking instability is generally subcritical. In section 5, we conduct Abaqus-based numerical simulations of the TK and necking instabilities and show that for the class of constitutive models considered, TK instability appears first under force control, and it gives way to necking instability under displacement control. Section 6 is devoted to an analytical description of necking propagation. We determine the solutions for the two “phases” analytically and compare them with numerical simulations. The paper is then concluded in the final section with some additional remarks.

2. Governing equations

We consider a sufficiently large circular or square hyperelastic sheet that is subjected to an equibiaxial tension in its plane. We denote by $B_0$ and $B_e$ the undeformed and the finitely deformed configurations, respectively. The sheet thicknesses in $B_0$ and $B_e$ are denoted by H and h, respectively. We shall look for axisymmetric deformations, and so cylindrical polar coordinates will be employed, with coordinates $\theta$ , z, and r in $B_e$ corresponding to 1-, 2-, and 3-directions, respectively. The r and $\theta$ coordinates define positions in the plane, whereas z measures distance in the out-of-plane direction such that $z=0$ corresponds to the mid-plane and the surfaces to $z=\pm h/2$ . Thus, corresponding to the equibiaxial deformation, the three principal stretches are given by

$${\lambda }_1={\lambda }_3=\lambda ,{\lambda }_2={\lambda }^{-2},$$

(1)

in terms of a single parameter $\lambda$ due to the constraint of incompressibility.

To introduce the incremental equations, we denote the deformation gradient from $B_0$ to $B_e$ by $\overline{F}$ and the associated nominal stress by $\overline{S}$ . A small amplitude axially symmetric perturbation is now applied to $B_e$ , giving rise to the final configuration $B_t$ , and the associated incremental displacement field $\delta \mathit{x}$ is given by

$$\delta \mathit{x}=u(r,z){\mathit{e}}_r+v(r,z){\mathit{e}}_z,$$

(2)

where ${\mathit{e}}_r$ and ${\mathit{e}}_z$ are the basis vectors in the r- and z-directions, respectively, and u and v are the associated displacement components. The deformation gradient corresponding to the deformation $B_0\to B_t$ is denoted by F and the associated nominal stress by S. We write $F=(I+\eta )\overline{F}$ so that $\eta$ denotes the deformation gradient associated with the incremental deformation $B_e\to B_t$ . The divergence operator with respect to coordinates in $B_0$ and $B_e$ will be denoted by Div and div, respectively.

The incremental equilibrium equation can best be expressed in terms of the incremental stress tensor $\chi$ defined by

$${\chi }^{\mathrm{T}}={\overline{J}}^{-1}\overline{F}(S-\overline{S}),$$

(3)

where the superscript T stands for transpose, and $\overline{J}$ denotes the determinant of $\overline{F}$ (which is unity in the current case but is kept in the formula to maintain the generality of the formula). With the use of the identity $\mathrm{d}\mathrm{iv}{\overline{J}}^{-1}\overline{F}=0$ and the equilibrium equations $\mathrm{D}\mathrm{iv}\overline{S}=0$ and $\mathrm{D}\mathrm{iv}S=0$ , we obtain the incremental equilibrium equation

$$\mathrm{d}\mathrm{iv}{\chi }^{\mathrm{T}}=0.$$

(4)

For the current axisymmetric deformation, only the equations corresponding to $i=2,3$ are not satisfied automatically, and they are given by

$${\chi }_{3j,j}+\frac{1}{r}({\chi }_{33}-{\chi }_{11})=0,{\chi }_{2j,j}+\frac{1}{r}{\chi }_{23}=0.$$

(5)

For our linear and weakly nonlinear analyses, we need expansions of ${\chi }_{\mathrm{ij}}$ up to the quadratic order, and they are given by

$${\chi }_{\mathrm{ij}}=B_{\mathrm{jilk}}{\eta }_{\mathrm{kl}}+\overline{p}{\xi }_{\mathrm{ji}}-p^{*}({\delta }_{\mathrm{ji}}-{\xi }_{\mathrm{ji}})+\frac{1}{2}{B}_{\mathrm{jilknm}}^2{\eta }_{\mathrm{kl}}{\eta }_{\mathrm{mn}}+\cdots ,$$

(6)

where

$${\xi }_{\mathrm{ji}}={\delta }_{\mathrm{ji}}-F_{\mathrm{Ai}}^{-1}{\overline{F}}_{\mathrm{jA}}={\eta }_{\mathrm{ji}}-{\eta }_{\mathrm{jm}}{\eta }_{\mathrm{mi}}+\cdots .$$

(7)

See, for example, Fu et al. [69]. In equation (6), the $\overline{p}$ and $p^{*}$ are the Lagrangian multipliers associated with the deformations $B_0\to B_e$ and $B_e\to B_t$ , respectively. The ${\eta }_{\mathrm{kl}}$ are the components of $\eta$ given by

$$\eta =\left[\begin{array}{ccc}\frac{u}{r}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& v_z\hfill & v_r\hfill \\ \multicolumn{1}{c}{0}& u_z\hfill & u_r\hfill \end{array}\right],v_z\equiv \frac{\partial v}{\partial z},v_r\equiv \frac{\partial v}{\partial r},\mathrm{and}\mathrm{so}\mathrm{on},$$

(8)

and $B_{\mathrm{jilk}}$ and $B_{\mathrm{jilknm}}^2$ are the first- and second-order instantaneous elastic moduli, the expressions of which can be found in Ogden [70] or Fu and Ogden [69].

Due to the introduction of the extra variable $p^{*}$ , the equilibrium equations are augmented by the incompressibility condition which can be expanded as

$${\eta }_{\mathrm{ii}}-\frac{1}{2}{\eta }_{\mathrm{mn}}{\eta }_{\mathrm{nm}}+\cdots =0.$$

(9)

For the current problem, both the top and bottom surfaces are traction-free. Thus, the associated incremental boundary conditions are

$${\chi }_{22}=0,{\chi }_{32}=0,\mathrm{on}z=\pm h/2.$$

(10)

3. Linear analysis

With linearization, the incompressibility condition (9) for the current problem reduces to

$$\frac{\partial \left(\mathrm{ru}\right)}{\partial r}+\frac{\partial \left(\mathrm{rv}\right)}{\partial z}=0.$$

(11)

It is then convenient to satisfy this condition automatically by introducing a “stream function” $\varphi (r,z)$ such that

$$u=\frac{1}{r}{\varphi }_z,v=-\frac{1}{r}{\varphi }_r,$$

(12)

where as in equation (8), a subscript signifies differentiation (e.g., ${\varphi }_z=\partial \varphi /\partial z$ ). The non-zero stress components are given by

$${\chi }_{11}=(B_{1111}+\overline{p})\frac{u}{r}+B_{1122}{v}_z+B_{1133}{u}_r-p^{*},$$

(13)

$${\chi }_{22}=B_{1122}\frac{u}{r}+(B_{2222}+\overline{p})v_z+B_{2233}{u}_r-p^{*},$$

(14)

$${\chi }_{33}=B_{1133}\frac{u}{r}+B_{2233}{v}_z+(B_{3333}+\overline{p})u_r-p^{*},$$

(15)

$${\chi }_{23}=B_{3232}{v}_r+(B_{3223}+\overline{p})u_z,$$

(16)

$${\chi }_{32}=B_{2323}{u}_z+(B_{2332}+\overline{p})v_r.$$

(17)

On substituting these expressions together with equation (12) into the linearized forms of equation (5) and then eliminating $p^{*}$ by cross-differentiation, we obtain

$$\alpha ({\varphi }_{\mathrm{rrrr}}-\frac{2}{r}{\varphi }_{\mathrm{rrr}}+\frac{3}{r^2}{\varphi }_{\mathrm{rr}}-\frac{3}{r^3}{\varphi }_r)+2\beta ({\varphi }_{\mathrm{rrzz}}-\frac{1}{r}{\varphi }_{\mathrm{rzz}})+\gamma {\varphi }_{\mathrm{zzzz}}=0,$$

(18)

where

$$\alpha =B_{1212},2\beta =B_{1111}+B_{2222}-2B_{1122}-2B_{1221},\gamma =B_{2121}.$$

(19)

Note that these constants are the same as those in Dowaikh and Ogden [71], Fu and Rogerson [72] for plane-strain deformations. The $\alpha$ , for instance, is given by

$$\alpha =\frac{({\sigma }_1-{\sigma }_2){\lambda }_1^2}{{\lambda }_1^2-{\lambda }_2^2},$$

(20)

with ${\sigma }_1$ and ${\sigma }_2$ denoting the principal Cauchy stresses in the $x_1$ - and $x_2$ -directions. We make the assumption that both $\alpha$ and $\gamma$ are positive material constants, which is consistent with the Baker–Ericksen inequality $({\sigma }_1-{\sigma }_2)({\lambda }_1-{\lambda }_2)>0$ $({\lambda }_1\ne {\lambda }_2)$ [73].

Equation (18) admits a “normal mode” buckling solution of the form

$$\varphi (r,z)=r{J}_1\left(\zeta r\right)S\left(\zeta z\right),$$

(21)

where $\zeta$ is a constant playing the role of wavenumber, $J_1(x)$ is Bessel’s function of the first kind, and the other function $S(\zeta z)$ is to be determined. This is similar to looking for periodic bifurcation solutions proportional to ${\mathrm{e}}^{\mathrm{i}k{x}_1}$ in rectangular coordinates. In fact, $J_1(x)$ is oscillatory, although it also decays like $1/\sqrt{x}$ for large values of x due to geometric spreading.

On substituting equation (21) into equation (18) and simplifying by making use of the identity

$$J_{\nu }\left(x\right)=\frac{2(\nu -1)}{x}{J}_{\nu -1}\left(x\right)-J_{\nu -2}\left(x\right),$$

we obtain $J_1(\zeta r)\{\gamma S^{(4)}(\zeta z)-2\beta S″(\zeta z)+\alpha S(\zeta z)\}=0$ . It then follows that

$$\gamma S^{(4)}\left(\zeta z\right)-2\beta S″\left(\zeta z\right)+\alpha S\left(\zeta z\right)=0,$$

(22)

which has the same form as in the plane-strain case. Trying a solution of the form $S(x)={\mathrm{e}}^{x/\omega }$ , we find that $\omega$ must satisfy the bi-quadratic equation

$$\alpha {\omega }^4-2\beta {\omega }^2+\gamma =0,\Rightarrow {\omega }^2=\frac{1}{\alpha }(\beta \pm \sqrt{{\beta }^2-\alpha \gamma }).$$

(23)

Thus, a general solution is given by

$$S\left(\zeta z\right)=C_1\sinh \frac{\zeta }{\sqrt{k_1}}z+C_2\sinh \frac{\zeta }{\sqrt{k_2}}z+C_3\cosh \frac{\zeta }{\sqrt{k_1}}z+C_4\cosh \frac{\zeta }{\sqrt{k_2}}z,$$

(24)

where

$$k_1=\frac{1}{\alpha }(\beta -\sqrt{{\beta }^2-\alpha \gamma }),k_2=\frac{1}{\alpha }(\beta +\sqrt{{\beta }^2-\alpha \gamma }).$$

(25)

The boundary conditions (10) take the form

$$B_{2323}{u}_z+(B_{2332}+\overline{p})v_r=0,\mathrm{on}z=\pm h/2,$$

(26)

$$B_{1122}\frac{u}{r}+(B_{2222}+\overline{p})v_z+B_{2233}{u}_r-p^{*}=0,\mathrm{o}\mathrm{n}z=\pm h/2.$$

(27)

The $p^{*}$ in equation (27) can be eliminated by differentiating equation (27) with respect to r and then using equation (5)₁ to eliminate $p_r^{*}$ . We then obtain

$$\begin{array}{c}(B_{1122}-B_{1111}-\overline{p})(u_{\mathrm{rr}}+\frac{u_r}{r}-\frac{u}{r^2})\hfill \\ -B_{2121}{u}_{\mathrm{zz}}+(B_{2222}-B_{1221}-B_{1122})v_{\mathrm{rz}}=0,\mathrm{o}\mathrm{n}z=\pm \frac{h}{2}.\hfill \end{array}$$

(28)

On substituting equations (21) and (24) into the four boundary conditions (26) and (28), we obtain four algebraic equations. By simple addition and subtraction, the four equations can be decoupled into two equations for $C_1$ and $C_2$ and another two equations for $C_3$ and $C_4$ . Since the determinants of the two coefficient matrices cannot vanish simultaneously, when the equations for $C_1$ and $C_2$ have non-trivial solutions, $C_3$ and $C_4$ must vanish. In this case, $\varphi$ , and hence the vertical displacement v, are odd functions of z. The associated modes are thus extensional. However, when the equations for $C_3$ and $C_4$ have non-trivial solutions, $C_1$ and $C_2$ must vanish, and the associated modes are flexural. The bifurcation conditions for extensional and flexural modes are thus given by

$$\begin{array}{c}\sqrt{k_1}(1+k_1)(k_2(2\beta +\gamma )-\gamma )\tanh \left(\frac{\zeta h}{2\sqrt{k_1}}\right)\hfill \\ -\sqrt{k_2}(1+k_2)(k_1(2\beta +\gamma )-\gamma )\tanh \left(\frac{\zeta h}{2\sqrt{k_2}}\right)=0,\hfill \end{array}$$

(29)

$$\begin{array}{c}\sqrt{k_1}(1+k_1)(k_2(2\beta +\gamma )-\gamma )\tanh \left(\frac{\zeta h}{2\sqrt{k_2}}\right)\hfill \\ -\sqrt{k_2}(1+k_2)(k_1(2\beta +\gamma )-\gamma )\tanh \left(\frac{\zeta h}{2\sqrt{k_1}}\right)=0,\hfill \end{array}$$

(30)

respectively. Expanding equation (29) to order $(\zeta h{)}^2$ , we obtain

$$\gamma (\beta +\gamma )+\frac{1}{24}\{\alpha \gamma -{(2\beta +\gamma )}^2\}{\left(\zeta h\right)}^2+\cdots =0,$$

(31)

where we have used equation (25) to eliminate $k_1$ and $k_2$ .

From the discussion in Fu [36], we may postulate that the bifurcation condition for localized necking can be obtained by setting the leading order term in equation (31) to zero, that is $\beta +\gamma =0$ since $\gamma >0$ , or equivalently

$$B_{1111}+B_{2222}+2B_{2121}-2B_{1221}-2B_{1122}=0.$$

(32)

In terms of the strain-energy function $W({\lambda }_1,{\lambda }_2,{\lambda }_3)$ , this bifurcation condition takes the form

$$2{\lambda }^2{W}_2^{(0)}+W_{22}^{(0)}-2{\lambda }^3{W}_{12}^{(0)}+{\lambda }^6{W}_{11}^{(0)}=0,$$

(33)

where $W_2=\partial W/\partial {\lambda }_2$ , $W_{12}={\partial }^{2}W/\partial {\lambda }_1\partial {\lambda }_2$ , and so on, and a superscript “(0)” signifies evaluation at $({\lambda }_1,{\lambda }_2,{\lambda }_3)=(\lambda ,{\lambda }^{-2},\lambda )$ .

This bifurcation condition has the following interpretations. The nominal stresses in the $x_3$ - and $x_2$ -directions (i.e., r- and z-directions) are given by

$$s_3=W_3-{\lambda }_3^{-1}\overline{p},s_2=W_2-{\lambda }_2^{-1}\overline{p}.$$

On eliminating $\overline{p}$ using the condition $s_2=0$ , we obtain

$$s_3=W_3-\frac{{\lambda }_2}{{\lambda }_3}{W}_2=\frac{\partial }{\partial {\lambda }_3}W({\lambda }_1,{\lambda }_1^{-1}{\lambda }_3^{-1},{\lambda }_3),$$

where in the last term, we are viewing W as a function of ${\lambda }_1$ and ${\lambda }_3$ . Viewing $s_3$ as a function of ${\lambda }_1$ and ${\lambda }_3$ , we obtain

$${\frac{\partial s_3}{\partial {\lambda }_3}|}_{{\lambda }_1\mathrm{f}\mathrm{ixed}}=W_{33}-\frac{1}{{\lambda }_1{\lambda }_3^2}{W}_{32}+\frac{2}{{\lambda }_1{\lambda }_3^3}{W}_2-\frac{1}{{\lambda }_1{\lambda }_3^2}[W_{32}-\frac{1}{{\lambda }_1{\lambda }_3^2}{W}_{22}].$$

On evaluation at $({\lambda }_1,{\lambda }_2,{\lambda }_3)=(\lambda ,{\lambda }^{-2},\lambda )$ , we obtain

$${\lambda }^6{\frac{\partial s_3}{\partial {\lambda }_3}|}_{{\lambda }_1\mathrm{f}\mathrm{ixed}}={\lambda }^6{W}_{33}^{(0)}-2{\lambda }^3{W}_{32}^{(0)}+2{\lambda }^2{W}_2^{(0)}+W_{22}^{(0)}.$$

(34)

Since indexes 1 and 3 can be exchanged, this may be replaced by

$${\lambda }^6{\frac{\partial s_3}{\partial {\lambda }_3}|}_{{\lambda }_1\mathrm{f}\mathrm{ixed}}={\lambda }^6{W}_{11}^{\left(0\right)}-2{\lambda }^3{W}_{12}^{\left(0\right)}+2{\lambda }^2{W}_2^{\left(0\right)}+W_{22}^{\left(0\right)}.$$

(35)

It is seen that setting $\partial s_3/\partial {\lambda }_3=0$ gives the bifurcation condition (32).

In the case of localized bulging of an inflated rubber tube, the bifurcation condition was shown to be equivalent to the Jacobian determinant of the internal pressure P and resultant axial force N equal to zero when P and N are each viewed as functions of two deformation variables such as the internal volume and axial stretch [42]. In the current setting, it is natural to choose the nominal stresses $s_1$ and $s_3$ as the force variables and ${\lambda }_1$ and ${\lambda }_3$ as the deformation variables with ${\lambda }_2$ eliminated using the incompressibility condition ${\lambda }_2={\lambda }_1^{-1}{\lambda }_3^{-1}$ . Thus, we have

$$s_1=\frac{\partial }{\partial {\lambda }_1}W({\lambda }_1,{\lambda }_1^{-1}{\lambda }_3^{-1},{\lambda }_3),s_3=\frac{\partial }{\partial {\lambda }_3}W({\lambda }_1,{\lambda }_1^{-1}{\lambda }_3^{-1},{\lambda }_3).$$

(36)

The associated Jacobian determinant can be evaluated as follows

$$\frac{\partial s_1}{\partial {\lambda }_1}\frac{\partial s_3}{\partial {\lambda }_3}-\frac{\partial s_1}{\partial {\lambda }_3}\frac{\partial s_3}{\partial {\lambda }_1}={\left(\frac{\partial s_3}{\partial {\lambda }_3}\right)}^2-{\left(\frac{\partial s_3}{\partial {\lambda }_1}\right)}^2=(\frac{\partial s_3}{\partial {\lambda }_3}+\frac{\partial s_3}{\partial {\lambda }_1})(\frac{\partial s_3}{\partial {\lambda }_3}-\frac{\partial s_3}{\partial {\lambda }_1}),$$

(37)

where all the partial derivatives are evaluated at ${\lambda }_1={\lambda }_3=\lambda$ (and hence indexes 1 and 3 can be exchanged). Thus, the Jacobian determinant equal to zero implies

$$\frac{\partial s_3}{\partial {\lambda }_3}-\frac{\partial s_3}{\partial {\lambda }_1}=0,\mathrm{o}\mathrm{r}\frac{\partial s_3}{\partial {\lambda }_3}+\frac{\partial s_3}{\partial {\lambda }_1}=0,$$

(38)

which may be compared with the condition $\partial s_3/\partial {\lambda }_3=0$ for necking. The two equations in equation (38) are obviously equivalent to

$${\left\{\underset{{\lambda }_3\to {\lambda }_1}{lim}\frac{s_3-s_1}{{\lambda }_3-{\lambda }_1}\right\}}_{{\lambda }_1=\lambda }=0,\mathrm{o}\mathrm{r}\frac{d}{d\lambda }{s}_3(\lambda ,\lambda )=0,$$

(39)

which are simply the conditions for the TK and limiting-point instabilities, respectively. Thus, for the current problem, the Jacobian determinant equal to zero is not the bifurcation condition for necking.

For the purpose of illustration, we consider the following class of strain-energy functions

$$W=\frac{2{\mu }_1}{m_1^2}({\lambda }_1^{m_1}+{\lambda }_2^{m_1}+{\lambda }_3^{m_1}-3)+\frac{2{\mu }_2}{{m_2}^2}({\lambda }_1^{m_2}+{\lambda }_2^{m_2}+{\lambda }_3^{m_2}-3),$$

(40)

where ${\mu }_1,{\mu }_2$ , $m_1,m_2$ are the material constants. We then have

$${\left\{\frac{\partial s_1}{\partial {\lambda }_3}\right\}}_{{\lambda }_1={\lambda }_3=\lambda }={\left\{\frac{\partial s_3}{\partial {\lambda }_1}\right\}}_{{\lambda }_1={\lambda }_3=\lambda }=\frac{2}{{\lambda }^2}({\mu }_1{\lambda }^{-2m_1}+{\mu }_2{\lambda }^{-2m_2}),$$

(41)

which is positive provided both ${\mu }_1$ and ${\mu }_2$ are positive. Assuming that the latter is true, it then follows from equation (38) that necking would occur after the TK instability but before the limiting-point instability. For instance, for the case when ${\mu }_2=0,m_1=1/2$ , we have

$${\lambda }_{\mathrm{T}\mathrm{K}}=2^{2/3},{\lambda }_{\mathrm{n}\mathrm{ecking}}=3^{2/3},{\lambda }_{\mathrm{l}\mathrm{imiting}}=4^{2/3},$$

(42)

where the three subscripts in $\lambda$ signify association with the TK, necking, and limiting-point instabilities, respectively. However, when ${\mu }_2={\mu }_1/12,m_1=1/2,m_2=2$ with which the material exhibits moderate stiffening for large values of $\lambda$ , we have

$$\begin{array}{c}{\lambda }_{\mathrm{T}\mathrm{K}1}=1.698,{\lambda }_{\mathrm{T}\mathrm{K}2}=7.803,{\lambda }_{\mathrm{n}\mathrm{ecking}1}=2.318,\hfill \\ {\lambda }_{\mathrm{n}\mathrm{ecking}2}=7.487,{\lambda }_{\mathrm{l}\mathrm{imiting}1}=2.957,{\lambda }_{\mathrm{l}\mathrm{imiting}2}=7.102,\hfill \end{array}$$

where we have two bifurcation values for each of the three instabilities. Finally, when ${\mu }_2={\mu }_1/80,m_1=1/2,m_2=4$ which will be considered further in our numerical simulations later, we have

$${\lambda }_{\mathrm{T}\mathrm{K}1}=1.647,{\lambda }_{\mathrm{T}\mathrm{K}2}=3.403,{\lambda }_{\mathrm{n}\mathrm{ecking}1}=2.439,{\lambda }_{\mathrm{n}\mathrm{ecking}2}=2.919,{\lambda }_{\mathrm{l}\mathrm{imiting}}\to \infty ,$$

(43)

where we have two bifurcation values for both the TK instability and necking, but no limit points exist due to the stronger stiffening behavior of the material model. To be more precise, the second bifurcation value for necking is really the value for bulging when unloading takes place. We expect, however, that the TK instability can only occur under tension control even if the bifurcation condition (38)₁ is satisfied. Under displacement control, the two in-plane principal stretches are forced to be equal, and we expect that the TK instability should give way to necking if equation (33) is satisfied. These facts will indeed be confirmed by our numerical simulations to be conducted in section 5.

4. Weakly nonlinear analysis

Having determined the bifurcation condition for necking, we now denote the critical stretch by ${\lambda }_{\mathrm{c}r}$ and conduct a weakly nonlinear analysis to derive an amplitude equation for the necking solution. To simplify notation, we scale r, z, and all displacement components by h, and use the same notations for the scaled quantities. As a result, we have $h=1$ .

We use ${\lambda }_1={\lambda }_3\equiv \lambda$ as the control parameter in our near-critical analysis. Guided by the analysis in Fu [36], we may write

$$\lambda ={\lambda }_{\mathrm{cr}}+\epsilon {\lambda }_0,$$

(44)

and define a far distance variable s through

$$s=\sqrt{\epsilon }r,$$

(45)

where ${\lambda }_0$ is a constant and $\epsilon$ is a small positive parameter characterizing the order of deviation of $\lambda$ from its critical value ${\lambda }_{\mathrm{cr}}$ . Note that as a result of equation (44), the moduli $B_{\mathrm{jilk}}$ and $B_{\mathrm{jilknm}}^2$ must all be expanded in terms of $\epsilon$ as well, but these expansions are not written out for the sake of brevity. The pressure $\overline{p}$ is eliminated using the identity $\overline{p}=B_{2121}-B_{2112}$ .

By analyzing the linear results in section 3 in the limit $\lambda \to {\lambda }_{\mathrm{c}\mathrm{r}}$ , we may deduce that $v=O(\sqrt{\epsilon }u),p^{*}=O(\sqrt{\epsilon }u)$ . From the fact that the ${\lambda }_0$ in equation (44) must appear at the order where a solvability condition should be imposed, we find that u must be of order $\sqrt{\epsilon }$ . Thus, we look for an asymptotic solution of the form

$$\begin{array}{c}u=\sqrt{\epsilon }\{u^{(1)}(s,z)+\epsilon u^{(2)}(s,z)+{\epsilon }^2{u}^{(3)}(s,z)+\cdots \},\hfill \\ v=\epsilon \{v^{(1)}(s,z)+\epsilon v^{(2)}(s,z)+{\epsilon }^2{v}^{(3)}(s,z)+\cdots \},\hfill \\ p^{*}=\epsilon \{p^{(1)}(s,z)+\epsilon p^{(2)}(s,z)+{\epsilon }^2{p}^{(3)}(s,z)+\cdots \},\hfill \end{array}$$

(46)

where all the functions on the right-hand sides are to be determined from successive approximations.

On substituting equation (46) into the equilibrium equations (5), the incompressibility condition (9), the boundary conditions (10), and then equating the coefficients of like powers of $\epsilon$ , we obtain a hierarchy of boundary value problems. The solutions at different orders are obtained as follows. In the following description, the two equilibrium equations in equation (5) are referred to as the r- and z-equilibrium equations, respectively.

At leading order, the r-equilibrium equation subject to the second boundary condition in equation (10) yields

$$u^{(1)}(s,z)=A\left(s\right),$$

(47)

where $A(s)$ is to be determined. Integrating the incompressibility condition then gives

$$v^{(1)}(s,z)=-z\frac{1}{s}{\left(\mathrm{sA}\left(s\right)\right)}^{\prime }+B\left(s\right),$$

(48)

where $B(s)$ is another function to be determined. Integrating the z-equilibrium equation subject to the first boundary condition at $z=\pm h/2$ yields the following expression for $p^{(1)}(s,z)$

$$p^{(1)}(s,z)=-(B_{2222}-B_{1122}+B_{2121}-B_{2112})\frac{1}{s}{\left(\mathrm{sA}\left(s\right)\right)}^{\prime },$$

(49)

where here and hereafter in this section all the moduli are evaluated at $\lambda ={\lambda }_{\mathrm{cr}}$ .

At second order, solving the r-equilibrium equation yields an expression for $u^{(2)}(s,z)$ that contains two new functions $C(s)$ and $D(s)$ in the form $C(s)+\mathrm{zD}(s)$ . Subtracting and adding the second boundary condition at $z=\pm h/2$ , respectively, we obtain

$$B_{2222}-2B_{1122}+B_{1111}+2B_{2121}-2B_{2112}=0,$$

(50)

and

$$D\left(s\right)=-B\prime \left(s\right).$$

(51)

The first result (50) may be shown to be equivalent to the bifurcation condition (33). Integrating the incompressibility condition at second order gives an expression for $v^{(2)}(s,z)$ that contains a new function $F(s)$ , whereas integrating the z-equilibrium yields an expression for $p_2(s,z)$ containing a new function $E(s)$ . Finally, subtracting and adding the second boundary condition at $z=\pm h/2$ , respectively, we obtain $\mathrm{sB}″(s)+B\prime (s)=0$ , and an expression for E(s). Solving the equation for $B(s)$ , we obtain $B(s)=d_1\ln s+d_2$ . Since $v^{(1)}$ should be bounded at $s=0$ , and without loss of generality, we may also impose the condition $v^{(1)}(0,0)=0$ , we must set $d_1=d_2=0$ . Thus, we may set $B(s)=0$ . This then completes the solution at second order.

At third order, the nonlinear terms come into play and a solvability condition must be imposed. The r-equilibrium equation can be solved to find an expression for $u^{(3)}(s,z)$ that contains two new functions $G(s)$ and $H(s)$ in the form $G(s)+\mathrm{zH}(s)$ . Subtracting and adding the second boundary condition evaluated at $z=\pm h/2$ , respectively, we obtain the amplitude equation for $A(s)$ and an expression for $H(s)$ . After some simplification, it is found that the amplitude equation takes the form

$$c_0\frac{d}{\mathrm{ds}}\frac{1}{s}\frac{d}{\mathrm{ds}}\mathrm{sP}\prime \left(s\right)+c_1{\lambda }_{0}P\prime \left(s\right)+c_2\frac{d}{\mathrm{ds}}{P}^2\left(s\right)+c_{3}A″\left(s\right)(A\prime \left(s\right)-\frac{1}{s}A\left(s\right))=0,$$

(52)

where a prime signifies differentiation, $P(s)$ is defined by

$$P\left(s\right)=\frac{1}{s}{\left(\mathrm{sA}\left(s\right)\right)}^{\prime },$$

(53)

and the three coefficients are given by

$$\begin{array}{c}{c}_0=\frac{1}{12}(B_{1212}-B_{2121}),\hfill \\ c_1=2{B\prime }_{1122}+2{B\prime }_{1221}-{B\prime }_{1111}-2{B\prime }_{2121}-{B\prime }_{2222},\hfill \\ c_2=\frac{1}{4}(-4B_{1111}-2B_{1122}+6B_{2222}-B_{111111}^2+4B_{111122}^2-B_{111133}^2\hfill \\ -6B_{112222}^2+2B_{112233}^2+2B_{222222}^2),\hfill \\ c_3=B_{1122}-B_{1111}+B_{111122}^2-B_{112233}^2-\frac{1}{2}{B}_{111111}^2+\frac{1}{2}{B}_{111133}^2.\hfill \end{array}$$

In the above expressions, $B{\prime }_{1122}$ denotes $d{B}_{1122}/d\lambda$ and so on, and we have used the bifurcation condition (50) to eliminate $B_{2112}$ . When the strain-energy function is given by equation (40) with ${\mu }_2=0$ , these constants take the form

$$\begin{array}{c}{c}_0=\frac{{\lambda }^{-2m_1}(2m_1({\lambda }_{\mathrm{c}\mathrm{r}}^6+1)({\lambda }_{\mathrm{c}\mathrm{r}}^{6m_1}-1)-m_1^2({\lambda }_{\mathrm{c}\mathrm{r}}^6-1){({\lambda }_{\mathrm{c}\mathrm{r}}^{3m_1}+1)}^2-({\lambda }_{\mathrm{c}\mathrm{r}}^6-1){({\lambda }_{\mathrm{c}\mathrm{r}}^{3m_1}-1)}^2)}{24m_1({\lambda }_{\mathrm{c}\mathrm{r}}^{3m_1}-1)},\hfill \\ c_1=2{\lambda }_{\mathrm{c}\mathrm{r}}^{-2m_1-1}({\lambda }_{\mathrm{c}\mathrm{r}}^{3m_1}-m_1({\lambda }_{\mathrm{c}\mathrm{r}}^{3m_1}-2)+2),c_3=-(m_1-1){\lambda }_{\mathrm{c}\mathrm{r}}^{m_1},\hfill \\ c_2=-\frac{(m_1-1)(m_1({\lambda }_{\mathrm{c}r}^{3m_1}-2)+2({\lambda }_{\mathrm{c}r}^{3m_1}-1))}{2m_1{\lambda }_{\mathrm{c}r}^{2m_1}},\hfill \end{array}$$

where ${\lambda }_{\mathrm{c}r}=\{(1+m)/(1-m){\}}^{1/(3m)}$ . More specifically when $m_1=1/2$ , we have

$$\frac{c_1}{c_0}=\frac{9\sqrt[3]{3}}{2},\frac{c_2}{c_0}=\frac{27}{8},\frac{c_3}{c_0}=\frac{9}{4}.$$

As a consistency check, we may neglect the nonlinear terms in equation (52) to obtain

$$c_0\frac{d}{\mathrm{ds}}\frac{1}{s}\frac{d}{\mathrm{ds}}\mathrm{sP}\prime \left(s\right)+c_1{\lambda }_{0}P\prime \left(s\right)=0.$$

(54)

On substituting a solution of the form $A(s)=J_1(\mathrm{ks})$ into equation (53) and the resulting expression into equation (54) and simplifying, where k is a constant, we obtain

$$c_1{\lambda }_0-c_0{k}^2=0.$$

(55)

However, expanding equation (31) around $\lambda ={\lambda }_{\mathrm{c}\mathrm{r}}$ and taking $\zeta =k\sqrt{\epsilon },h=1$ , we obtain

$${\left\{\frac{d}{d\lambda }\gamma (\beta +\gamma )\right\}}_{\mathrm{c}\mathrm{r}}{\lambda }_0+\frac{1}{24}{\{\alpha \gamma -{(2\beta +\gamma )}^2\}}_{\mathrm{c}\mathrm{r}}{k}^2=0,$$

(56)

where the subscripts “cr” signify evaluation at $\lambda ={\lambda }_{\mathrm{c}\mathrm{r}}$ . We have verified that equation (55) is indeed consistent with equation (56).

As another consistency check, we may expand equation (53) fully out and omit all the terms that are divided by powers of s to obtain its planar counterpart

$$c_0{A}^{(4)}\left(s\right)+c_1{\lambda }_{0}A″\left(s\right)+c_2^{*}A\prime \left(s\right)A″\left(s\right)=0,$$

(57)

where

$$c_2^{*}=2c_2+c_3=3B_{2222}-3B_{1111}-B_{111111}^2+3B_{111122}^2-3B_{112222}^2+B_{222222}^2.$$

(58)

Equation (57) applies to a finitely deformed sheet that is subjected to plane-strain incremental deformations, and recovers the static counterpart of the amplitude equation (4.41) in Fu [36]. It has an exact solution given by

$$A\left(s\right)=\frac{6c_0}{c_2^{*}}\sqrt{\frac{-c_1{\lambda }_0}{c_0}}\tan h\left(\frac{1}{2}\sqrt{\frac{-c_1{\lambda }_0}{c_0}}s\right).$$

(59)

This solution has the property $A\prime (s)\to 0$ as $s\to \infty$ and is the localized necking solution in the two-dimensional (2D) case.

Although the solution equation (59) for the plane-strain case tends to a constant as $s\to \infty$ , we expect its counterpart for the axisymmetric case to decay algebraically due to geometrical spreading. Thus, for large s, we may neglect the quadratic part in equation (52) to obtain

$$\frac{d}{\mathrm{ds}}\frac{1}{s}\frac{d}{\mathrm{ds}}\mathrm{sP}\prime \left(s\right)-{\kappa }^{2}P\prime \left(s\right)=0,$$

(60)

where ${\kappa }^2=-c_1{\lambda }_0/c_0$ . Integrating once and setting the constant to zero, we obtain

$$s^{2}P″\left(s\right)+\mathrm{sP}\prime \left(s\right)-{\kappa }^2{s}^{2}P\left(s\right)=0.$$

This is a modified Bessel’s equation. It has a decaying solution given by

$$P\left(s\right)=a_1{K}_0\left(\kappa s\right),$$

(61)

where $K_0$ is the modified Bessel function of the second kind and $a_1$ is a constant. On substituting equation (61) into equation (53) and integrating, we obtain

$$A\left(s\right)=A_{\infty }\left(s\right)\equiv \frac{a_2}{s}-\frac{a_1}{\kappa }{K}_1\left(\kappa s\right),$$

(62)

where $a_1$ and $s_2$ are the constants, and $K_1$ is the modified Bessel function of the second kind that has the asymptotic behavior

$$K_1\left(x\right)~\sqrt{\frac{\pi }{2x}}[1+\frac{3}{8x}+\cdots ]{\mathrm{e}}^{-x},\mathrm{as}x\to \infty .$$

(63)

The asymptotic behavior (62) is consistent with our earlier assumption that $A(s)$ decays algebraically. In principle, the fourth-order differential equation (52) can be integrated subject to the conditions

$$A\left(s\right),A″\left(s\right)\to 0\mathrm{a}\mathrm{s}s\to 0,$$

(64)

and

$$A\left(s\right)\to A_{\infty }\left(s\right),A\prime \left(s\right)\to A_{\infty }\prime \left(s\right)\mathrm{a}\mathrm{s}s\to \infty .$$

(65)

The necessary numerical integration will not be carried out here. Instead, we shall proceed to a fully nonlinear numerical simulation of the necking solution in the next section. We conclude this section by observing that the decay behavior (62) requires $\kappa$ $(=\sqrt{-c_1{\lambda }_0/c_0})$ to be real. For all the material models considered, it is found that $c_1/c_0>0$ . Thus, bifurcation into a localized necking solution is generally subcritical (i.e., ${\lambda }_0<0$ ).

5. Fully nonlinear numerical simulations

To verify the various theoretical predictions presented in the previous sections, we now conduct nonlinear finite element (FE) simulations within Abaqus 2020/Standard [74] to investigate the effect of loading methods and material properties on the type of bifurcation and post-bifurcation behavior. In particular, we shall verify that when the loading is force controlled, TK instability will precede necking, whereas when the loading is displacement-controlled, TK instability will be suppressed and it is necking that will occur.

5.1. Equibiaxial stretching of a square plate

We first perform fully three-dimensional (3D) FE simulations for a square plate with $L=5H$ , where L is the length and width of the plate and H is the thickness of the plate (Figure 1(a)). The plate is meshed using eight-node fully integrated hybrid linear brick elements (Abaqus Element C3D8H). The following constraints are imposed on the opposite faces of the plate

$$u_{\alpha }^{L_i}=u_{\alpha }^{R_i}+u_{\alpha }^X,u_{\alpha }^{F_i}=u_{\alpha }^{B_i}+u_{\alpha }^Y,(\alpha =1,2,3,i=1,2,...,N),$$

(66)

where N is the number of nodes on each of the four faces (with L, R, F, and B denoting left, right, front, and back faces, respectively), and $u_{\alpha }^{L_i}$ denotes the displacement in the $\alpha -\mathrm{direction}$ at the i^th node of the left face and so on. The $u_{\alpha }^X$ and $u_{\alpha }^Y$ denote the displacements in the $\alpha -\mathrm{direction}$ at two arbitrarily chosen reference points X and Y that are used to apply the loading. In the above notations, the paired nodes $L_i$ and $R_i$ have the same coordinates in the 2–3 plane, whereas $F_i$ and $B_i$ have the same coordinates in the 1–2 plane. Moreover, to prevent rigid body motions, we fix the node at the center of the plate.

Figure 1.

Comparison between numerical results and theoretical predictions for the Treloar–Kearsley instability. (a) Geometrical parameters and boundary conditions for the square plate. (b) ${\lambda }_3$ versus ${\lambda }_1$ . (c) Nominal stress $S_1$ scaled by ${\mu }_1$ versus ${\lambda }_1$ .

The plate is loaded by applying the following concentrated forces and displacement constraints at the two reference points

$$\begin{array}{c}{f}_1^X=1.001*\overline{f},u_2^X=u_3^X=0,\mathrm{at}\mathrm{point}X,\hfill \\ f_3^Y=\overline{f},u_1^Y=u_2^Y=0,\mathrm{at}\mathrm{point}Y,\hfill \end{array}$$

(67)

where $f_1^X$ and $f_3^Y$ are the concentrated forces applied at the reference points X and Y in the 1- and 3-directions, respectively, and $\overline{f}$ represents the magnitude of the concentrated force. The factor 1.001 in $f_1^X$ can be viewed as an imperfection that will trigger the TK instability. Note that the near equibiaxial conditions can also be obtained by applying near equal pressures directly at the edges of the plate, but this will result in convergent issues in simulations. The imposed constraints (66) serve to circumvent convergent issues and ensure a homogeneous deformation.

It is known that the bifurcation associated with the TK instability can be either sub-critical or super-critical depending on the strain-energy function used. The corresponding test has been given by Ogden [64] in his equation (3.13). We thus consider two representative material models exhibiting these different behaviors. The first model is given by equation (40) with ${\mu }_1=1,{\mu }_2=0$ , and $m_1=1/2$ for which the bifurcation is sub-critical. We employ the modified RIKS algorithm [75] to capture the unloading behavior associated with the sub-critical bifurcation. We find that, initially, the common principal stretch in the 1- and 3-directions increases uniformly as the load increases. Then, when ${\lambda }_1$ reaches the critical value 1.58 given by equation (42), bifurcation occurs, resulting in unequal stretches and decline of the nominal stress in both directions. As Figure 1(b) and (c) shows, there is very good agreement with the theoretical predictions for both the stretches and nominal stress.

The second material model that we consider is given by equation (40) with $m_1=1/2$ , $m_2=4$ , and $0.012<{\mu }_2/{\mu }_1<0.035$ . The two limits of ${\mu }_2/{\mu }_1$ are such that if ${\mu }_2<0.012$ , the bifurcation will be sub-critical, whereas if the upper limit is exceeded, TK instability can no longer occur. The second term in equation (40) introduces straining-stiffening effect, and as a result, there are two bifurcation stretches associated with the TK instability. It is observed that the nonuniform stretches in the 1- and 3-directions initiating from the first bifurcation stretch would become uniform again at the second bifurcation stretch. This is demonstrated in Figure 2 where we compare the numerical and theoretical results for the principal stretch ${\lambda }_3$ as a function of ${\lambda }_1$ when ${\mu }_2/{\mu }_1=1/80$ . Again, we have excellent agreement between numerical results and theoretical predictions.

Figure 2.

Emergence and disappearance of the Treloar–Kearsley instability for the material model (40) with ${\mu }_2/{\mu }_1=1/80,m_1=1/2$ , $m_2=4$ .

Moreover, we remark that the TK instability can also be observed in an equibiaxially stretched square plate by introducing material imperfections. To demonstrate this, we introduce a small imperfection at the center of the plate by assuming that the shear modulus in a small rectangular area is 0.5% smaller than that for the surrounding material (see the schematic in Figure 3, where the blue region in the square plate has the modulus reduced). Then, by applying equal concentrated forces at the reference points X and Y, we can observe unequal stretches when the applied nominal stress $S_1$ is larger than 1.67 (see the solid black lines in Figure 3, which are obtained by RIKS simulation in Abaqus). However, if the reference points are stretched using displacement control (i.e., by applying equal displacements in both directions), TK instability will be suppressed, and instead, localized necking will appear at a larger stretch as indicated by the solid yellow line in Figure 3, which is obtained using STATIC simulation in Abaqus. The TK instability giving way to necking under displacement control will be discussed further in the next section. In addition, we note that a rectangular (anisometric) imperfection is essential for observing the TK instability in a square plate. If a square (isometric) defect was introduced instead, the TK instability would not be triggered even with force control.

Figure 3.

Numerical simulations of a square plate with material imperfections, modeled by the strain-energy function (40) with $m_1=1/2$ , $m_2=4$ , and ${\mu }_2/{\mu }_1=1/80$ . Solid yellow line: results obtained by applying equal displacements at the reference points using the STATIC module in Abaqus (in which case TK instability is suppressed and necking takes place at the kinked point in the black circle). Solid black line: results obtained by applying equal concentrated forces (force control) at the reference points using the RIKS module in Abaqus (in which case the TK instability occurs as a super-critical bifurcation).

5.2. All-round stretching of a circular plate

To demonstrate the occurrence of necking in a more convenient way than in the above square plate setting, we now consider a circular plate with $R_0=2H$ , where $R_0$ and H are respectively the radius and thickness of the plate in the undeformed configuration. We control the amount of stretch on the circular edge of the plate in order to suppress the occurrence of TK instability. In addition, since strain-stiffening is essential for observing the stable propagation stage of necking, we shall consider the two-terms material model (40) with ${\mu }_2\ne 0,m_2\ge 2$ . To reduce the cost of simulations, we only model a quarter of the plate and apply axisymmetric boundary conditions on the symmetric surfaces shown in Figure 4(a). Moreover, we note that the initiation location of necking is sensitive to geometric imperfections. In Figure 4(b) and (c), we showcase snapshots of our simulations with gradual thickness thinning toward the origin $(H_i=0.9999H_o)$ or the circular edge $(H_o=0.9999H_i)$ , respectively, where $H_i$ and $H_o$ denote the plate thicknesses at the center and outer edge. For the case with thinning toward the origin, necking initiates at the center of the plate and propagates outwards as the applied load increases (Figure 4(b)). In contrast, when thinning toward the outer edge is introduced, necking initiates at the outer edge and spreads toward the center as the applied load increases. At some stage, however, the necked region switches to the center of the plate and propagates outwards with increased load (Figure 4(c)). Thus, necking starting at the center seems to be the dominant deformation mode, and so from now on, we only consider plates with thinning toward the origin. In addition, the critical stretch for necking initiation is sensitive to the magnitude of geometric imperfections. For consistence, all the following numerical results correspond to the same geometric imperfection mentioned above. Due to the presence of near-critical snap-through behavior, we also add volume-proportional damping to the model to facilitate convergence (using the option STABILIZE in Abaqus), and set the dissipated energy fraction equal to ${\mathrm{e}}^{-10}$ and the maximum ratio of stabilization to strain energy equal to 0.01. Furthermore, we observe that although 2D axisymmetric simulations in Abaqus considering only a rectangular slice can reduce the cost of simulations significantly, the 2D model does not seem to be capable of characterizing the necking phenomenon in the sense that axisymmetric elements are prone to distortion when necking occurs, especially for materials with only a moderate strain-stiffening behavior.

Figure 4.

Simulations of necking in a circular plate subject to an all-round stretch, modeled by the strain-energy function (40) with $m_1=1/2$ , $m_2=4$ , and ${\mu }_2/{\mu }_1=1/120$ . (a) Geometric parameters for a quarter of the plate. (b and c) Numerical snapshots of the plate subjected to increased stretches with gradual thickness thinning toward the center (b) or the edge (c).

Figure 5 shows a typical plot of the necking amplitude $\hat{P}$ against the stretch $\lambda$ , where $\hat{P}$ denotes the absolute difference between the vertical displacements at points $(0,-H_i/2,0)$ and $(R_0,-H_i/2,0)$ normalized by $H_i$ . It is seen that it is appropriate to take the numerical value of critical stretch as the stretch at which an amplitude jump takes place (the snap-through behavior mentioned above). We observe, however, that the variation of amplitude against stretch should have a slight snap-back (indicated by the dashed line) that is not captured by our simulations because the snap-back section corresponds to an unstable solution (see Figure 3(b) in Fu and Xie [76] for a discussion of the same phenomenon in the context of inflation of a membrane tube). To verify our theoretical predictions, we fix $m_1=1/2$ , and compare the numerical and theoretical results in Figure 6 for the critical stretch as a function of ${\mu }_1/{\mu }_2$ for two typical values of $m_2$ . It is seen that for each given $m_2$ , the critical stretch is a monotonically decreasing function of the ratio ${\mu }_1/{\mu }_2$ , and the excellent agreement between the numerical and theoretical results confirms that the bifurcation condition for necking is indeed given by equation (32), not by the criterion based on the Jacobian determinant.

Figure 5.

Dependence of the necking amplitude on the applied stretch $\lambda$ . The solid curve represents numerical results for the amplitude evolution at the center of the plate. The dashed line is a sketch of how the actual solution should look like. The material is modeled by the strain-energy function (40) with $m_1=1/2$ , $m_2=4$ and ${\mu }_2/{\mu }_1=1/80$ .

Figure 6.

Dependence of the critical stretch on ${\mu }_1/{\mu }_2$ with (a) $m_1=1/2$ and $m_2=2$ or (b) $m_1=1/2$ and $m_2=4$ .

Furthermore, to quantify the post-bifurcation behavior of the necking deformation, we report evolution of the necking amplitude in Figure 7. The results shown in Figure 7(a) represent a plate with $m_1=1/2$ , $m_2=4$ , ${\mu }_2/{\mu }_1=1/80$ , and the six profiles correspond to $\lambda =2.40$ , 2.429, 2.430, 2.431, 2.710, and 3.10, respectively (the critical stretch being 2.439). The first two profiles correspond to tiny inhomogeneous deformations due to thickness thinning, whereas the third profile with $\lambda =2.430$ corresponds to the profile immediately after the jump in Figure 5 has taken place. Thus, the critical stretch found numerically is 2.430 which differs from the exact critical stretch 2.439 given by equation (43)₂ with a relative error less than 0.4%. It is observed that the necking amplitude first increases rapidly with respect to $\lambda$ , and then as a maximum amplitude is approached, the necked section begins to propagate outwards. This is reminiscent of the well-known growth/propagation behavior of localized bulges in an inflated rubber tube [44,77]. To demonstrate the effect of strain-stiffening on the necking amplitude, we report in Figure 7(b) the dependence on ${\mu }_2/{\mu }_1$ of the necking amplitude with $m_1=1/2$ , $m_2=4$ , and a fixed stretch of $\lambda =3.0$ . It can be seen that the necking amplitude increases as the value of ${\mu }_2/{\mu }_1$ decreases, which implies that the necking amplitude is significantly affected by the strain-stiffening response of the material.

Figure 7.

(a) Evolution of the necking profile with $m_1=1/2$ , $m_2=4$ , and ${\mu }_2/{\mu }_1=1/80$ . The six profiles correspond to $\lambda =2.40,2.429,2.430,2.431,2.710$ , and 3.10, respectively. (b) Evolution of the necking amplitude of the plate with $m_1=1/2$ , and $m_2=4$ when subjected to equal-biaxial stretch of $\lambda =3.0$ . The six curves represent ${\mu }_2/{\mu }_1=1/80,1/100,1/120,1/140,1/160$ , and $1/180$ , respectively.

Finally, we remark that although the radius/thickness ratio used in our simulations is only 2 initially ( $R_0/H=2$ ), this ratio has become as large as 7.62 when the stretch reaches its critical value 2.439. We have run simulations to verify that this choice of initial aspect ratio is sufficient for obtaining results that are almost free of edge effects. To demonstrate this, we have shown in Figure 8 two representative simulations corresponding to $R_0/H=5,7$ superimposed on the results for $R_0/H=2$ . Since $\lambda$ is the deformed radius divided by the undeformed radius, and is no longer a true characterization of local radial stretch after necking has taken place, we choose the three stretches shown in Figure 8(a) such that the maxima of V are the same (or closest numerically) and check whether the remaining profiles would then also coincide. It is seen that this is indeed the case for both V and U although the three profiles for each displacement terminate at different values of r. Notably, to demonstrate the negligible effect of radius/thickness ratio on the observed necking phenomenon with various imperfections, the results shown in Figure 8 correspond to a circular material imperfection whereby the shear modulus of the material within $R<H$ is 0.01% smaller than that of the surrounding material.

Figure 8.

Total displacement (i.e., displacement measured from the undeformed configuration) in the (a) z-direction and (b) radial direction, at the lower surface of the plate as a function of the radial coordinate r in the current configuration.

6. Characterization of necking propagation

To offer more insight into the propagation stage, we have shown in Figure 9 a typical set of simulation results for the evolution of ${\lambda }_z$ versus r as the outer radius is increased. It is seen that in the propagation stage (red curves), the ${\lambda }_z$ is a constant in the necked region, but seems to be slowly varying in the r-direction in the un-necked region. For easy reference, we refer to the necked region centered around the origin as the “− phase” and the other region with larger thickness as the “+ phase.” In this section, we give the two “phases” an analytical description by assuming that ${\lambda }_z$ is a constant in both “phases.” Although the two “phases” are connected by a thin transition region, we further assume that the transition region has zero thickness and simply refer to it as the interface. We use superscripts “+” and “−” to signify evaluations at the “+ phase” and “− phase” sides of the interface, respectively. Thus, for instance, ${\lambda }_1^{+}$ denotes the azimuthal stretch as the interface is approached from the “+ phase” side.

Figure 9.

Numerical simulation of the evolution of the principal stretch in the z-direction as the radius is increased. Thestrain-energy function is given by equation (40) with $m_1=1/2$ , $m_2=4$ , and ${\mu }_2/{\mu }_1=1/80$ , and the aspect ratio is $R_0/H=7$ .The blue and red curves correspond to the growth and propagation stages, respectively, and the six red curves correspond to ${\lambda }^{*}=2.412,2.456,2.536,2.616,2.696$ , and 2.776.

The “− phase” is homogeneous. The solution is the same as in section 3, and it is uniquely determined by a single parameter which we take to be the stretch ${\lambda }_z^{-}$ in the z-direction. We then have ${\lambda }_1={\lambda }_3=1/\sqrt{{\lambda }_z^{-}}$ and ${\lambda }_2={\lambda }_z^{-}$ .

For the “+ phase,” the solution is inhomogeneous although ${\lambda }_z^{+}$ is assumed to be a constant. In terms of cylindrical polar coordinates, the axisymmetric solution is given by

$$r=r\left(R\right),\theta =\mathrm{\Theta },z={\lambda }_z^{+}Z.$$

(68)

The corresponding deformation gradient is given by

$$F=r\prime \left(R\right){\mathit{e}}_r\otimes {\mathit{e}}_r+\frac{r}{R}{\mathit{e}}_{\theta }\otimes {\mathit{e}}_{\theta }+{\lambda }_z^{+}{\mathit{e}}_z\otimes {\mathit{e}}_z,$$

so that ${\lambda }_1=r/R$ , ${\lambda }_2={\lambda }_z^{+}$ , and ${\lambda }_3=r\prime (R)$ . To simplify notation, from now on we shall write ${\lambda }_z^{+}$ simply as ${\lambda }_z$ .

By integrating the incompressibility condition det $F=1$ , we obtain

$$r=\sqrt{{\lambda }_z^{-1}{R}^2+C},$$

(69)

where C is a constant of integration. This constant may be determined by evaluating this expression at the edge of the plate $R=R_0$

$$C=R_0^2({\lambda }^{*2}-{\lambda }_z^{-1}),$$

(70)

where ${\lambda }^{*}$ denotes the azimuthal stretch at $R=R_0$ (so that ${\lambda }^{*}=r^{*}/R_0$ with $r^{*}$ denoting the deformed radius). Once the hydrostatic pressure has been eliminated with the use of ${\sigma }_2=0$ , the principal Cauchy stresses ${\sigma }_1$ and ${\sigma }_3$ are given by

$${\sigma }_1={\lambda }_1{W}_1-{\lambda }_2{W}_2,{\sigma }_3={\lambda }_3{W}_3-{\lambda }_2{W}_2.$$

(71)

It then follows that

$${\sigma }_3-{\sigma }_1={\lambda }_1{w}_{{\lambda }_1},$$

(72)

where $w_{{\lambda }_1}=\mathrm{dw}/d{\lambda }_1$ and $w=w({\lambda }_1)$ is the reduced strain-energy function defined by

$$w\left({\lambda }_1\right)=W\left({\lambda }_1,{\lambda }_z,\frac{1}{{\lambda }_z{\lambda }_1}\right).$$

(73)

The only equilibrium equation that is not satisfied automatically is

$$\frac{d{\sigma }_3}{\mathrm{dr}}+\frac{{\sigma }_3-{\sigma }_1}{r}=0.$$

(74)

With use of the result (72) and the relation ${\lambda }_1=r/R$ , we may convert equation (74) to the form [78]

$$\frac{d{\sigma }_3}{d{\lambda }_1}=\frac{w_{{\lambda }_1}}{{\lambda }_z{\lambda }_1^2-1}.$$

(75)

On integrating this equation across the entire “+ phase,” we obtain

$${\sigma }_3^{*}-{\sigma }_3^{+}=\underset{{\lambda }_1^{+}}{\overset{{\lambda }_1^{*}}{\int }}\frac{w_{{\lambda }_1}}{{\lambda }_z{\lambda }_1^2-1}d{\lambda }_1,$$

(76)

where the superscript “*” signifies evaluation at the outer edge $r=r^{*}$ .

The jump conditions that need to be imposed across the interface are

$${\sigma }_3^{+}={\sigma }_3^{-},{\lambda }_1^{+}={\lambda }_1^{-},W^{+}-W^{-}={\sigma }_3^{-}(\frac{{\lambda }_z^{-}}{{\lambda }_z}-1).$$

(77)

The first two correspond to traction and displacement continuity, whereas the third corresponds to stationarity of the total energy with respect to perturbations of the interface position in the undeformed configuration [79–81].

With the use of equation (77)_1,2, equations (76) and (77)₃ may be rewritten as

$${\sigma }_3^{*}-{\sigma }_3^{-}=\underset{{\lambda }_1^{-}}{\overset{{\lambda }_1^{*}}{\int }}\frac{w_{{\lambda }_1}}{{\lambda }_z{\lambda }_1^2-1}d{\lambda }_1,W({\lambda }_1^{-},{\lambda }_z,\frac{1}{{\lambda }_1^{-}{\lambda }_z})-W(\frac{1}{\sqrt{{\lambda }_z^{-}}},{\lambda }_z^{-},\sqrt{{\lambda }_z^{-}})={\sigma }_3^{-}(\frac{{\lambda }_z^{-}}{{\lambda }_z}-1).$$

(78)

These two equations can be used to find ${\lambda }_z$ and ${\lambda }_z^{-}$ for each specified ${\lambda }_1^{*}$ (noting that ${\lambda }_1^{-}=1/\sqrt{{\lambda }_z^{-}}$ ). The location of the interface, $r=r_s$ say, may be expressed in terms of ${\lambda }_z^{-}$ and ${\lambda }_z$ by solving the equation ${\lambda }^{-}={\lambda }^{+}$ , i.e.

$$\frac{1}{\sqrt{{\lambda }_z^{-}}}=\frac{r_s}{\sqrt{{\lambda }_z{r}_s^2-{\lambda }_{z}C}},\Rightarrow {\left(\frac{r_s}{R_0}\right)}^2=\frac{{\lambda }_z{\lambda }_1^{*2}-1}{{\lambda }_z-{\lambda }_z^{-}}.$$

(79)

Since we must have $0<r_s<R_0{\lambda }_1^{-}$ , the ${\lambda }_1^{*}$ must satisfy

$$1/\sqrt{{\lambda }_z}<{\lambda }_1^{*}<1/\sqrt{{\lambda }_z^{-}}.$$

(80)

In Figure 10, we have shown the dependence of ${\lambda }_z$ $(={\lambda }_z^{+})$ and ${\lambda }_z^{-}$ on ${\lambda }_1^{*}$ , as predicted by equation (78), when the strain-energy function is given by equation (40) with $m_1=1/2$ , $m_2=4$ , and ${\mu }_2/{\mu }_1=1/80$ . On the same figure are shown six representative numerical simulation results (evaluated approximately at the center of the “+ phase”). It is seen that there is excellent agreement for the values of ${\lambda }_z^{-}$ with relative error less than 0.57%. The agreement for ${\lambda }_z^{+}$ is less satisfactory. We attribute this to the approximations that we have used, namely, that the thickness of the transition region is neglected and the ${\lambda }_z$ in the “+ phase” is independent of r.

Figure 10.

Comparison of theoretical (solid lines) and numerical simulation results (squares) for ${\lambda }_z^{+}$ and ${\lambda }_z^{-}$ as functions of ${\lambda }_1^{*}$ . The squares are corresponding to the six red lines in Figure 9.

7. Conclusion

For problems that have translational invariance in one spatial direction, e.g., the inflation of a homogeneous rubber tube or extension of a homogeneous plate in a state of plane strain, a rigorous analysis based on center manifold reduction can be conducted to establish the existence of a localized solution bifurcating from the trivial state [34,35,82]. A necessary condition for localization has the simple interpretation that a static extensional mode can exist in the small wavenumber limit. For the problem of localized bulging in an inflated rubber tube, this condition has been found to be equivalent to the Jacobian determinant of two force functions equal to zero. For problems defined over an infinite or semi-infinite domain that do not have translational invariance in the direction of localization, such as the current axially symmetric necking problem, no center manifold reductions seem to have been attempted. We postulated that localized necking is still governed by the condition that an extensional bifurcation mode exists in the small wavenumber limit. Although we have not yet been able to ground this condition on any rigorous theory, we have been able to verify this postulate with the aid of Abaqus simulations. We further showed that the bifurcation condition for necking does not correspond to the Jacobian determinant equal to zero. The latter condition can always be factorized and the two factors give the conditions for the TK instability and limiting-point instability, respectively. It is also demonstrated, first numerically and then analytically, that necking evolution follows a growth/propagation process that is typical of similar localization problems. It is hoped that the current analysis will guide future studies of elastic localizations under multiple fields such as axisymmetric necking of dielectric membranes that are subjected to both electrical and mechanical loading.