{'fraction_non_alphanumeric': 0.04421979455405185, 'fraction_numerical': 0.021685961193543126, 'mean_word_length': 4.509522098454905, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Growth induced instabilities in a circular hyperelastic plate

In this work, we have explored growth-induced mechanical instability in an isotropic circular hyperelastic plate. Consistent two-dimensional governing equations for a plate under a general finite strain are derived using a variational approach. The derived plate equations are solved using the compound matrix method for two cases of axisymmetric growth conditions - purely radial, and combined radial and circumferential growth. The effect of growth on the buckling behaviour of the plate (in particular, the critical growth factor and the associated buckling mode shapes) is investigated for different thickness values. These results are applicable to model growth induced deformation in planar soft tissues such as skin.


Introduction
Mechanical instabilities are ubiquitous in nature and often result in pattern formation in thin elastic structures. Classical plate theories like Kirchhoff-Love theory, Föppl-von Kármán theory, and Mindlin-Reissner theory have been widely used to study the instability behaviour of thin elastic structures (Coman and Haughton, 2006;Coman et al., 2015;Li et al., 2010). These theories are based on apriori kinematic assumptions which are suitable for solving small strain problems. Also, these theories when applied to plates under general loading conditions give inconsistent results due to the underlying assumptions of displacement variation along thickness of the plate. To overcome the inconsistencies in the classical plate theories, Kienzler (2002) developed the consistent plate theory using uniform approximation of unknown variables based on linear elasticity. The consistent plate theory does not incorporate apriori kinematic assumptions and all the coefficients are treated as independent unknown variables. However, small strain theories based on linear elasticity principles are not suitable for finite strain problems. To alleviate these problems, a consistent finite-strain plate theory was proposed by Dai and Song (2014). This approach was based on the principle of minimisation of potential energy under general three-dimensional loading conditions. They derived the two-dimensional plate vector equation by employing variational principle and series expansion of the independent variables about the bottom surface of a hyperelastic plate. Wang et al. (2016) extended this approach to incompressible hyperelastic materials with extra unknown variables to accommodate the incompressibility constraint. Mechanical instability is also a common phenomenon in morphoelastic structures (Ben Amar et al., 2011) and soft biological tissues , which exhibit nonlinear mechanical response due to growth.
Growth not only changes the mass and geometry of structures but can also alter their mechanical properties and stress state (Goriely, 2017). Growth can induce residual stresses inside the body which result in large deformations leading to instabilities such as wrinkling, folding, creasing, and crumpling . Residual stresses are self equilibrating (Hoger, 1986;Ben Amar and Goriely, 2005) generally arising due to the incompatibility of growth. Goriely and Ben Amar (2005) and Vandiver and Goriely (2009) studied the buckling of cylindrical elastic structures subjected to differential growth and residual stresses. Dervaux et al. (2009) discussed the nonlinear behaviour of thin elastic structures subjected to growth by considering Föppl-von Kármán plate theory. Moulton and Goriely (2011) investigated the circumferential instability in differentially growing cylindrical elastic tube subjected to uniform pressure and evaluated the critical pressure for buckling using an incremental theory. Papastavrou et al. (2013) investigated wrinkling in growing surface by considering membrane with zero thickness in the derivation of potential energy. Wu and Ben Amar (2015) performed the bifurcation analysis of uni-directional growing disk reinforced with fibres and determined the effect of fibre anisotropy on critical growth factor. Recently, Wang et al. (2018) derived a consistent finite-strain plate theory for growth-induced large deformations and investigated the buckling and post-buckling behaviour of a thin rectangular hyperelastic plate under axial growth.
Consistent plate theories with finite-strain have a wide range of applications as they incorporate bending as well as stretching. They are suitable to approximate the behaviour of soft biological tissue such as skin (Tepole et al., 2011) and their bifurcation properties under residual stress Gupta, 2015, 2016). Skin undulation near the edges of wound has been observed in healing experiments on mice (Nassar et al., 2012;Wang et al., 2013). Beyond biomedical applications, the mechanics of instability with large deformations have important implications on wrinkling analysis of gossamer space structures (Wang et al., 2009;Deng et al., 2019) when exposed to temperature gradients. Since the last decade, research in understanding buckling/wrinkling instabilities during micro-fabrication in the field of flexible/stretchable electronics such as sensory skins used in robotics and wearable communication devices (Rogers et al., 2010;Wei and Zhao, 2014) has increased.
In this manuscript, we have used the consistent finite-strain plate theory given by Wang et al. (2018) to derive the governing differential equations for circular isotropic hyperelastic plates. The symmetry of circular geometry allows us to transform the resulting partial differential equations (PDEs) to ordinary differential equations (ODEs) while still retaining the key aspects of mechanics. Our current focus is to comprehend the underlying mechanics of such systems. The formulation is general and can be applied to other geometries by solving the resulting PDEs using numerical techniques such as finite element method. Following the multiplicative decomposition approach proposed by Rodriguez et al. (1994), the deformation gradient tensor is decomposed into growth tensor, describing change of mass or growth laws and elastic deformation tensor that ensures compatibility (no overlaps) and integrity (no cavitation) (Goriely and Ben Amar, 2007). The principle of minimum total potential energy is applied to derive the 3-D governing partial differential equations in the polar coordinate system. These equations are reduced to two dimensions using a series expansion of unknown functions in terms of the thickness variable. After establishing the governing equations, we study the instability behaviour of neo-Hookean circular plates growing in radial as well as combined (radial and circumferential) direction. The traction conditions at the bottom surface of plate is applied through a Winkler support in both the cases. The resulting system of ODEs is stiff and standard shooting methods are not able to evaluate the bifurcation results accurately. To resolve this problem, we have used the compound matrix method Reid, 1979, 1985) to determine the critical point of buckling.
The remainder of this paper is organised as follows. In Section 2, a general formulation for three-dimensional circular plate is established using a variational formulation. The two-dimensional plate governing equations are derived by eliminating the dependence on Z variable using Taylor's expansion. In Section 3, we discuss two examples of growthinduced instability in incompressible neo-Hookean circular plate. A detailed discussion of numerical results is provided along with the comparison of existing analytical results for rectangular plate. Finally, we present our conclusions in Section 4.

Notation used in this manuscript
Brackets: Three types of brackets are used. Round brackets ( ) are used to define the functions applied on parameters or variables. Square brackets [ ] are used to clarify the order of operations in an algebraic expression. Curly brackets { } are used to define a set. Square brackets are also used for matrices and tensors. At some places we use the square bracket to define the functional. Symbols: A variable typeset in a normal weight font represents a scalar. A lower-case bold weight fonts denotes a vector and bold weight upper-case denotes the tensor or matrices. Tensor product of two vectors a and b is defined as Tensor product of two second order tensors A and B is defined as either Higher order tensors are written in bold calligraphic font with a superscript as A A A (i) , where superscript 'i' tells that the function is differentiated i + 1 times. For example, ∂A∂A is a fourth order tensor. Operation of a fourth order tensor on a second order tensor is denoted as The symbol ∇ denotes the two-dimensional gradient operator. We use the word 'Div' to denote divergence in three dimensions. Functions: det(A) denote the determinant of the tensor A. tr(A) denote the trace of a tensor A. diag(a, b, c) denotes a second order tensor with only diagonal entries a, b and c. configuration are given by R, Θ and Z and in the deformed configuration by r, θ and z. Radius of the plate in the reference configuration is R 0 . Position vectors in the configurations B 0 and B t are denoted as X(R, Θ, Z) and x(r, θ, z), respectively. The deformation gradient for a circular plate is defined as (Dai and Song, 2014) where ζ = Re R + Θe Θ and k is the unit normal to the surface Ω in B 0 . The deformation gradient tensor can be decomposed as (Rodriguez et al., 1994) where G represents growth tensor and A represents the elastic deformation tensor that ensures compatibility. We also assume the plate to be incompressible and the incompressibility constraint is given by The energy density per unit volume φ of the material is where J G = det(G) = det(F) describes the local change in volume due to growth and φ 0 (FG −1 ) is the elastic strain energy density. The external work done by the traction is given as where q + (respectively q − ) represents the applied traction on top (respectively bottom) surface of region andq represents the traction on the lateral surface ∂Ω q × [0, 2h], ∂Ω q being the boundary along the lateral surface, and the symbol (·) denotes the inner product. If we neglect the body force, the total potential energy functional (ψ) for incompressible plate structure is where p(X) is the Lagrange multiplier, associated with the incompressibility constraint. We apply the principle of minimum potential energy and vanishing of the first variation with respect to x and p of the above functional results in the governing equations (2.7c) and the incompressibility constraint (2.3). Here, n is the unit outward normal to the boundary ∂Ω q and P = J G ∂φ 0 ∂A -p ∂L0 ∂A G −T is recognised as the first Piola Kirchhoff stress tensor. While deriving these equations, we have used the assumption that the rate of deformation of the growth process is very small compared to the elastic deformation (Ben Wang et al., 2018), therefore the growth tensor G is assumed to be constant in time. Auxiliary calculations are presented in Appendix A.

Specialisation to two dimensions
To obtain the 2-D formulation of circular plate, we perform the series expansion of x and p in terms of Z about the bottom surface, Z = 0 following the approach by Wang et al. (2018Wang et al. ( , 2019b x where (·) (n) = ∂ n (·) ∂Z n Z=0 (n = 0, 1, 2, 3, 4). Likewise, we write the Taylor's expansion of the deformation gradient F, the elastic deformation tensor A and the inverse transpose of the growth tensor G as The recursion relation for the expansion coefficients in (2.10a) is expressed as (2.11) Using (2.2) and (2.10a) -(2.10c) we obtain the expressions for higher derivatives of A as Similarly, the expansion of the first Piola Kirchhoff stress tensor P is αj , (2.14) . Further mathematical details of above calculations are provided in Appendix B.
The two-dimensional governing system has x (·) and p (·) as unknown functions. In order to form a closed system of equations, we write the boundary conditions (2.7b) at the bottom surface (Z = 0) and at top surface (Z = 2h) The stress equilibrium equation neglecting the body force and external traction is given by (2.7a) where ∇ is the two-dimensional differentiation operator. Upon use of (2.13), we obtain a recursion relation for the first Piola-Kirchhoff stress as The series expansion of unknown functions x(X) and p(X) have 19 unknowns with x n (n = 0, 1, 2, 3, 4) comprising 15 unknowns and p n (n = 0, 1, 2, 3) comprising 4 unknowns. Thus, a closed system of 19 equations for the solution of x and p is derived from equilibrium equation, traction (bottom and top surface) boundary conditions (2.7a -2.7c) and incompressibility condition (2.3) as . On subtracting the top (2.18) and bottom traction (2.17) condition we obtain the equilibrium equation P is the average stress over the thickness,q is the effective body force due to traction at top and bottom surface (see Appendix A). Using Taylor's expansion, the equilibrium equation (2.22) can be rewritten as (Wang et al., 2019a) where the subscript 't' represents the in-plane (or tangential) component of a vector or tensor, . Physically, Eq. (2.23) represents the balance of forces. The explicit expressions for x (2) , p (1) and x (3) , p (2) in terms of x (0) , x (1) and p (0) are given as (2.28) (2.31)

Growth induced deformation
In this section, we discuss two cases of circular hyperelastic plate growing in only radial direction (Wu and Ben Amar, 2015) and combined radial and circumferential directions. We assume a constant growth function in both the cases.

Radial growth
Consider a thin isotropic circular plate which undergoes axisymmetric deformation (θ = Θ) with the position vector x(r, θ, z) in deformed configuration. The series expansion of unknown functions r(R, Z), z(R, Z) and p(R, Z) in terms of Z is written as where we have used the notation r (n) = ∂ n r ∂Z n , z (n) = ∂ n z ∂Z n and p (n) = ∂ n p ∂Z n . We assume the plate to follow neo-Hookean elastic constitutive law given by where C 0 is the ground state shear modulus. On substitution of (2.2) in (3.4), the first Piola Kirchhoff stress P is obtained as (see Appendix C ) We can rewrite (2.23) in component form as We can substitute the expression of stress from (2.22) and omit the higher order terms in h to get Consider constant radial growth given by the growth tensor G = diag(λ, 1, 1) that results in where a superposed prime denotes partial derivative with respect to R. Substituting (3.10) in (3.9) we obtain the governing equation Note that for the sake of brevity, we do not explicitly represent the O(h) and O(h 2 ) terms here but they are utilized for calculations later in Section 3.1.2.

Plate supported by the Winkler foundation: Traction condition
We concern ourselves with a circular plate resting on a Winkler foundation as shown in Figure 2. Winkler foundation models the elastic support provided to the growing plate. The top surface (Z = 0) of the plate is assumed to be tractionfree and the bottom surface, is supported by the Winkler foundation which provides a transverse load q + 3 = −K 0 λW 0 , where K 0 is the elastic constant of the foundation and W 0 is the transverse (Z) component of the displacement. λ is the growth multiplier which represents the fact that the traction is applied in the current grown configuration, as mass and stress state of the plate changes due to growth.
For this case, the components of the effective body force (q) are given bȳ The governing equations are obtained by substituting (3.13) in (3.11) Assuming simply supported condition along the edge of circular plate, the boundary conditions at the center R = 0 and the edge R = R 0 are given as where the unknown variables in terms of r (0) and z (0) are given as ( 3.17) The explicit expressions for r (1) , z (1) , p (0) are derived in Appendix C.

Stability analysis
The principal solution for plate deformation is given by where the second equation ensures that the Z displacement of the lower surface vanishes, that is z (0) + 2hz (1) − 2h = 0 and in the homogeneous deformation we get the non zero quantity z (1) = λ (because r (0) ′ = 1 and z (0) ′ = 0). A small perturbation in the homogeneous equilibrium state by a parameter ǫ results in (3.20) To simplify our calculations, we eliminate the terms with P (2) in equation (2.23) by subtracting the divergence of (3.20) from (2.23) 2 (Wang et al., 2019a). We keep the terms that correspond to bending and obtain where m t is the tangential component of m. Physically, Eq. (3.20) represents the balance of moments. We define the dimensionless quantities where ρ ∈ [0, 1] and R 0 is radius of circular plate in the reference configuration. On substituting the ansatz (3.19) in (3.14) 1 and (3.21b), simplifying using (3.17), (3.22) and collecting only O(ǫ) terms The higher order derivative of in-plane displacement term is omitted for simplification. The plate boundary conditions (3.16) is given as Equations (3.23) and (3.24) are coupled ODEs and can be written as a system of first order ODEs by defining (3.26) The system of first order ordinary differential equations is written in the form of where Y = [∆U, ∆U ′ , ∆U ′′ , ∆W, ∆W ′ , ∆W ′′ ] = [y 1 y 2 y 3 y 4 y 5 y 6 ] and A A A is given by ( 3.28) where We first determine the critical value of growth factor (λ cr ) that results in the onset of a bifurcation and then we discuss the associated buckling modes. The system of first order ODEs (3.27) is treated as two-point boundary value problem. This stiff eigenvalue problem is solved using the compound matrix method Reid, 1979, 1985;Lindsay and Rooney, 1992;Haughton and Orr, 1997). Following the compound matrix approach, the system (3.27) is converted into 20 first order equations of the form Φ ′ = A A A * (ρ; λ,h, β)Φ (A detailed description of the solution procedure is given in Appendix D) The initial condition for the system of equations (3.29) is The target condition is achieved by having det(CM) = 0 in order to obtain the non-trivial solution, where matrix C corresponds to the boundary condition at the edge of circular plate (3.25) and M is the solution matrix For the current case, the corresponding initial conditions using (3.25) are given by Φ(12) = 1 and rest all are zero. Then, we integrate the system numerically in the interval of 0 < ρ ≤ 1 until we achieve the target condition on the other boundary which is given by det(CM)=(1,3,5)=Φ(6) = 0 (see Appendix D). The main objective of this optimization problem is to determine the critical value of growth factor λ cr for which the target value Φ(6) is zero.

Results and discussion
In this section, the buckling behaviour of a circular hyperelastic plate due to radial growth rested on Winkler foundation is presented. In order to validate our numerical scheme based on the compound matrix method, we first evaluate the onset of buckling of a rectangular plate under uni-directional growth for which an analytical solution has been provided by Wang et al. (2018).

a) Comparison of numerical results and analytical results for rectangular plate
Consider a rectangular plate of thickness 2h clamped at the ends X = ±1 and supported by a Winkler foundation as shown in Figure 3. The plate grows along the X-axis with a growth stretch λ. The compression effects due to the clamped boundary condition lead to buckling. The governing plate differential equation is given by (Wang et al., 2018) ψ 0 ∆W + ψ 2 ∆W ′′ + ψ 4 ∆W ′′′′ = 0, (3.32)  (Wang et al., 2018). b) Variation of the critical value of growth factor (λ cr ) against thickness to length ratio (h) evaluated using the compound matrix method. c) Direct comparison of (a) and (b).
where ψ 0 = α 2h , ψ 2 = 1 λ 2 + 3λ 6 λ 4 − 1 2 + [6 + hα]λ 4 , ψ 4 = 4h 2 3 + 9λ 4 3 + hα + [2 + 3hα]λ 4 + [3 + 2hα]λ 8 , The analytical results of bifurcation curves for various buckling modes is shown in Figure 4a. For thick plates (h > 0.15), the first fundamental mode dominates. As we lower the value of h, the higher modes (n = 2, 3, . . . ) are more stable. Figure 4b represents the variation of critical value of growth factor (λ cr ) with respect to thickness to length ratio (h) obtained by solving (3.32)-(3.33) using compound matrix method. Also, the numerical compound matrix method solution of (3.21a) and (3.21b) specialized to a rectangular plate geometry provides the same result. The λ cr variation is non-monotonous and has discontinuous derivatives at certain thickness values which suggests the phenomenon of mode jumping from high energy state to a lower one associated with changes in mode shape. We conclude that the peaks in Figure 4b that represent the transition of modes i.e., for thin plates, the higher modes are more stable than the lower modes. In Figure 4c, we superpose the numerical results on the analytical curves for direct comparison. It can be seen that the analytical and numerical results are in a near perfect agreement which shows the efficacy of the compound matrix method.

b) Numerical results for circular plate with radial growth
The compound matrix method is used to compute the buckling parameter (λ cr ) of the circular hyperelastic plate under radial growth condition. The variation of critical growth factor (λ cr ) with respect to normalized thickness is shown in Figure 5. The results show the non-monotonous nature of λ cr variation and the points at which the mode jump phenomenon happens at certain plate thicknesses. Based on the above results, the graph can be divided into 3 regions of thickness in which a particular mode is dominant over the other modes. In the region 3, which corresponds to lower plate thickness, a higher order mode and for higher plate thickness, the lower order mode dominates the buckling behaviour. Next, we determine the mode shapes corresponding to the 3 regions defined in Figure 5.
The mode shape results are obtained for a thickness value in the defined regions and the results are shown in Figure 6. In the region 3, a thicknessh = 0.024 is chosen for which λ cr = 1.0218 is used to evaluate the mode shape of the circular hyperelastic plate. At this value ofh, a higher order mode shape (mode 3) exists and variation of out of plane displacement amplitude in radial direction is shown in Figure 6. Similarly, the modes in the region 2 and 1 represented by mode 2 and mode 1 are evaluated by choosingh = 0.05, λ cr = 1.0350 andh = 0.2, λ cr = 1.0611, respectively. Thus, we conclude that the buckling parameter (λ cr ) of circular plate under radial growth is non-monotonic with higher order modes in the thin plate regime and lower order modes in the thicker plate regime.

Combined radial and circumferential growth
In this second case, we assume a constant isotropic growth function λ (i.e., λ rr = λ θθ = λ) in both radial and circumferential direction (Wu and Ben Amar, 2015). The isotropic growth tensor G = diag(λ, λ, 1) results in the following kinematics On substituting (3.35) in the governing equation (3.9), we obtain where the unknown variables in this case are calculated as The explicit expressions for r (1) , θ (1) , z (1) , p (0) are given in the Appendix E. Here, again we use the Winkler support on the bottom surface of the plate with traction given by (3.13).

Linear buckling analysis for combined growth model
We conduct the linear bifurcation analysis by perturbing the principal solution with small parameter (ǫ) to determine the onset of instability by assuming the form where m = 1, 2, 3... represents the mode number in the circumferential direction. On substituting (3.37) in (3.36), making use of (3.21b) along with (3.13) which have the updated traction component as q + 3 = −K 0 λ 2 W 0 and collecting only O(ǫ) terms, we obtain the coupled differential equations for dimensionless displacement functions U and W as  Figure 7: Variation of the critical growth parameter λ cr with normalised plate thicknessh and β = 0.2 for different circumferential mode number (m). (3.39)

Results and discussion
The compound matrix method is applied to solve the governing equations to determine variation of the critical growth factor (λ cr ) with respect to normalized thickness (h). The variation of critical buckling parameter (λ cr ) for different circumferential mode numbers are given in Figure 7. The results shows that the higher modes arise for lower value of λ cr for all thicknesses. Also, the phenomenon of mode jump is not observed here as there is no intersection of the bifurcation curves corresponding to different mode numbers. Furthermore, we note that the mode numbers in the case of combined radial/circumferential growth appear explicitly in the governing equations as opposed to interpretation of mode numbers in the purely radial growth case. The normalized out of plane and in-plane displacement amplitude variation in the radial direction corresponding to each circumferential mode number is shown in Figure 8 for a plate thickness ofh = 0.1. For this plate thickness, the bifurcation curve corresponding to m = 10 mode is more stable as we conclude by analysing Figure 7. Other mode shapes may not be achievable but are plotted for completeness. Analysing the mechanics of a circular plate under radial growth and combined growth shows that the buckling configuration of the plate changes with the increase of thickness in both the cases. The higher modes are more stable in thin regime of the plate in the former case and the bifurcation solution corresponding to higher modes are more stable for latter combined growth case regardless of the thickness. However, the critical value of bifurcation parameter (λ cr ) is low in combined growth case for higher modes as compared to the purely radial growth case.

Concluding remarks
Mechanical instabilities are often observed in thin films, elastic structures and in soft biological tissues. In this paper, we have used a consistent finite-strain plate theory to investigate the buckling behaviour of incompressible circular hyperelastic plate subjected to growth. A 3-D governing system of PDEs is converted into 2-D plate system using series expansion in terms of the thickness variable. We consider two examples of growth-induced instability in neo-Hookean circular plate under a) radial growth and b) combined radial/circumferential growth conditions. In both the cases, the circular plate is resting on a Winkler support and subjected to simply supported boundary conditions. The compound matrix method is used to evaluate the critical buckling growth factor of the hyperelastic plate. The numerical performance of compound matrix method is validated by comparing the results with the existing analytical solution for rectangular plate under uniform growth. The variation of critical growth factor (λ cr ) with normalised thickness (h) for both the cases is evaluated numerically and the results show that the plate buckles in different modes as we increase the thickness. In thin regime of the plate, higher modes are more stable and in thick regime, lower modes are more stable for purely radial growth case. However, in combined growth case the higher modes are more stable independent of the thickness of the plate.
The current work can be applied in modelling of thin soft biological tissues such as skin wrinkling during wound healing and ageing, modelling of wrinkling patterns in thin deployable space structures and stretchable electronics. Constitutive model considered in this paper is isotropic and biological materials are generally anisotropic heterogeneous materials. An extension of the current analysis to account for more generally applicable material models will be taken up as a future study. We have also restricted ourselves to determine the critical buckling load, but a post-buckling analysis may provide insights on the magnitude of the out of plane deformation and the associated mode-switching response. These avenues are currently under investigation and we will report our results in a suitable forum at a later stage.

A Appendix: Expression for Piola stress
The average stress is given as Subtracting the top and bottom surface traction condition using (2.17) and (2.18) Using equilibrium equation Now, the Piola stress (P) is described in terms of strain energy function (φ(F, G) = J G φ 0 (A)) as The series expansion for Piola Kirchhoff stress (P) is given as Using (2.10c) and (A.5) we have Substituting expansion of φ 0 (A) and L 0 (A) about A (0) we obtain Substituting expression for [A − A (0) ] using (2.10b) and compare the terms A=A (0) and L L L i (A (0) ) = ∂R 0 (A) ∂A i+1 where, and T ijkl = δ il δ jk . Note: In this current work, we assume incompressible material thus, det(A) = 1 C Appendix: Expression for x (1) and p (0) For an incompressible neo-Hookean material elastic strain energy function is given as where I 1 = tr(C) = tr(A T A) as elastic strain energy depends only on elastic deformation tensor (A) and tr(·) is trace of a tensor (·). Using (A.4), the Piola stress (P) for incompressible neo-Hookean material is Now, the first term P (0) in (A.5) for neo-Hookean material is rewritten as By using bottom traction condition, P (0) k = −q − and substituting expression for A (0) we have where J (0) = J G | Z=0 . Traction at bottom surface is given as where Cofac(F (0) )k = Cofac(F (0) ) ij k j e i . In this work, we use ∇x (0) * in place of Cofac(F (0) )k and for deformation gradient in polar coordinate system it is given by As the definition of determinant det(F (0) ) = r (1) e 1 + r (0) θ (1) e 2 + z (1) e 3 · Cofac(F (0) )k = x (1) · ∇x (0) * . (C.8) By combining (C.7) and (C.8) we have Using (C.5) we obtain the explicit expression for x (1) (C.10) To obtain the explicit expression for p (0) we substitute (C.10) into (C.9) which yields (C.11)

D Appendix: General Description on Compound Matrix Method
Consider a two-point boundary value problem expressed in first order ordinary differential equations subjected to boundary condition where λ is the eigenvalue or critical buckling parameter, Y is 1 × 2n vector, A A A is 2n × 2n matrix and B and C both are n × 2n full rank matrices i.e., n boundary conditions are given at x = a, b. Assume Using (D.1), the system of first order differential equation in terms of compound variable is (D.8) The system is now converted into ( 2n C n ) ordinary differential equations which is in the form of Subjected to initial condition at x = 0, Now, if we consider sixth order ODE system (n = 3) then the solution matrix is defined as On substituting the p (0) in (C.10) we obtain r (1) = p (0) ∇x 11 2C 0 λ 2 , θ (1) = p (0) ∇x 22 2C 0 λ 2 r (0) , z (1) = p (0) ∇x 33 2C 0 λ 2 . (E.2)