# Growth induced instabilities in a circular hyperelastic plate

Sumit Mehta  Department of Mechanical and Aerospace Engineering
Indian Institute of Technology Hyderabad, India  Gangadharan Raju  Department
of Mechanical and Aerospace Engineering
Indian Institute of Technology Hyderabad, India  Prashant Saxena Corresponding
author email<EMAIL_ADDRESS>Glasgow Computational Engineering
Centre, James Watt School of Engineering,
University of Glasgow, Glasgow G12 8LT, UK

###### Abstract

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.

Keywords: Growth, nonlinear elasticity, stability analysis, compound matrix
method.

Note: This is the author-generated version of paper to be published in the
International Journal of Solids and Structures (2021).

##  1 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 (Cao et al., 2012; Wu and Amar,
2015), which exhibit non-linear 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
(Li et al., 2012). 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 (Swain and
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
(Ng and 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 growth-induced 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.

###  1.1 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 $\mathbf{a}$ and
$\mathbf{b}$ is defined as
$[\mathbf{a}\otimes\mathbf{b}]_{ij}=[\mathbf{a}]_{i}[\mathbf{b}]_{j}$. Tensor
product of two second order tensors $\mathbf{A}$ and $\mathbf{B}$ is defined
as either
$[\mathbf{A}\otimes\mathbf{B}]_{ijkl}=[\mathbf{A}]_{ij}[\mathbf{B}]_{kl}$ or
$[\mathbf{A}\boxtimes\mathbf{B}]_{ijkl}=[\mathbf{A}]_{ik}[\mathbf{B}]_{jl}$.
Higher order tensors are written in bold calligraphic font with a superscript
as $\boldsymbol{\mathcal{A}}^{(i)}$, where superscript ‘$i$’ tells that the
function is differentiated $i+1$ times. For example,
$\boldsymbol{\mathcal{A}^{(1)}}=\frac{\partial
f(\mathbf{A})}{\partial\mathbf{A}\partial\mathbf{A}}$ is a fourth order
tensor. Operation of a fourth order tensor on a second order tensor is denoted
as
$[\boldsymbol{\mathcal{A}^{(1)}}:\mathbf{A}]_{ij}=[\boldsymbol{\mathcal{A}^{(1)}}]_{ijkl}[\mathbf{A}]_{kl}$
Inner product is defined as
$\mathbf{a}\cdot\mathbf{b}=[\mathbf{a}]_{i}[\mathbf{b}]_{i}$ and
$\mathbf{A}\cdot\mathbf{B}=[\mathbf{A}]_{ij}[\mathbf{B}]_{ij}$. The symbol
$\nabla$ denotes the two-dimensional gradient operator. We use the word ‘Div’
to denote divergence in three dimensions.

Functions: $\det(\mathbf{A})$ denote the determinant of the tensor
$\mathbf{A}$. $\text{tr}(\mathbf{A})$ denote the trace of a tensor
$\mathbf{A}$. $\text{diag}(a,b,c)$ denotes a second order tensor with only
diagonal entries $a,b~{}\text{and}~{}c$.

##  2 Governing equations for growing circular plate

Consider a thin circular plate with constant thickness $2h$ that occupies the
region $\Omega\times[0,2h]$ in the reference configuration
$\mathcal{B}_{0}\in\mathscr{R}^{3}$ and deforms to the current configuration
$\mathcal{B}_{t}\in\mathscr{R}^{3}$ as shown in Figure 1. Coordinates of a
point in the reference configuration are given by $R,\Theta$ and $Z$ and in
the deformed configuration by $r$, $\theta$ and $z$. Radius of the plate in
the reference configuration is $R_{0}$.

Figure 1: Schematic of a circular plate under growth (a) Finite deformation
due to growth (perspective) (b) Front view of the plate in reference
configuration with applied traction on the top, bottom, and side surfaces

.

Position vectors in the configurations $\mathcal{B}_{0}$ and $\mathcal{B}_{t}$
are denoted as $\mathbf{X}({R},~{}\Theta,~{}Z)$ and
$\mathbf{x}({r},~{}\theta,~{}z)$, respectively. The deformation gradient for a
circular plate is defined as (Dai and Song, 2014)

$\displaystyle\mathbf{F}$
$\displaystyle=\frac{\partial\mathbf{x}}{\partial\mathbf{X}}=\frac{\partial\mathbf{x}}{\partial{R}}\otimes\mathbf{e}_{R}+\frac{\partial\mathbf{x}}{\partial\Theta}\otimes\mathbf{e}_{\Theta}+\frac{\partial\mathbf{x}}{\partial
Z}\otimes\mathbf{k},$ (2.1a)
$\displaystyle=\frac{\partial\mathbf{x}}{\partial\mathbf{\zeta}}+\frac{\partial\mathbf{x}}{\partial
Z}\otimes\mathbf{k},$ (2.1b)

where $\mathbf{\zeta}=R\mathbf{e}_{R}+\Theta\mathbf{e}_{\Theta}$ and
$\mathbf{k}$ is the unit normal to the surface $\Omega$ in $\mathcal{B}_{0}$.
The deformation gradient tensor can be decomposed as (Rodriguez et al., 1994)

$\displaystyle\mathbf{F}=\mathbf{AG},$ (2.2)

where $\mathbf{G}$ represents growth tensor and $\mathbf{A}$ represents the
elastic deformation tensor that ensures compatibility. We also assume the
plate to be incompressible and the incompressibility constraint is given by

$\displaystyle
L({\mathbf{F,G}})=L_{0}(\mathbf{FG}^{-1})=L_{0}(\mathbf{A})=\text{det}(\mathbf{A})-1=0.$
(2.3)

The energy density per unit volume $\phi$ of the material is

$\displaystyle\phi(\mathbf{F},\mathbf{G})=J_{G}\phi_{0}(\mathbf{FG}^{-1}),$
(2.4)

where $J_{G}=\det(\mathbf{G})=\det(\mathbf{F})$ describes the local change in
volume due to growth and $\phi_{0}(\mathbf{FG}^{-1})$ is the elastic strain
energy density. The external work done by the traction is given as

$\displaystyle
V=\int_{\Omega}\mathbf{q}^{-}\cdot\mathbf{x}({\zeta},0)dA+\int_{\Omega}\mathbf{q}^{+}\cdot\mathbf{x}({\zeta},2h)dA+\int_{\partial\Omega_{q}}\int_{0}^{2h}{\mathbf{\tilde{q}}}\cdot\mathbf{x}(s,Z)dsdZ,$
(2.5)

where $\mathbf{q}^{+}(\text{respectively}~{}\mathbf{q}^{-})$ represents the
applied traction on top (respectively bottom) surface of region and
$\mathbf{\tilde{q}}$ represents the traction on the lateral surface
$\partial\Omega_{q}\times[0,2h]$, $\partial\Omega_{q}$ being the boundary
along the lateral surface, and the symbol ($\cdot$) denotes the inner product.
If we neglect the body force, the total potential energy functional ($\psi$)
for incompressible plate structure is

$\displaystyle\psi[\mathbf{x}(\mathbf{X}),p(\mathbf{X})]=\int_{\Omega}\int_{0}^{2h}J_{G}\phi_{0}(\mathbf{FG}^{-1})dV-\int_{\Omega}\int_{0}^{2h}\big{[}J_{G}~{}p(\mathbf{X})L_{0}(\mathbf{FG}^{-1})\big{]}dV-V,$
(2.6)

where $p(\mathbf{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 $\mathbf{x}$ and
$p$ of the above functional results in the governing equations

$\displaystyle\text{Div}\mathbf{P}=\mathbf{0},\hskip
198.02019pt\text{in}~{}\Omega\times[0,2h],$ (2.7a)
$\displaystyle\mathbf{Pk}\big{|}_{Z=0}=-\mathbf{q}^{-}(\mathbf{\zeta}),\hskip
43.36243pt\left.\mathbf{Pk}\right|_{Z=2h}=\mathbf{q}^{+}(\zeta),\hskip
36.135pt\text{on}~{}\Omega,$ (2.7b)
$\displaystyle\mathbf{Pn}=\mathbf{\tilde{q}}(s,Z),\hskip
184.28871pt\text{on}~{}\partial\Omega_{q}\times[0,2h].$ (2.7c)

and the incompressibility constraint (2.3). Here, $\mathbf{n}$ is the unit
outward normal to the boundary $\partial\Omega_{q}$ and
$\mathbf{P}=J_{G}\big{[}\displaystyle\frac{\partial\phi_{0}}{\partial\mathbf{A}}$
-$p\frac{\partial L_{0}}{\partial\mathbf{A}}\big{]}\mathbf{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 Amar and
Goriely, 2005; Wang et al., 2018), therefore the growth tensor $\mathbf{G}$ is
assumed to be constant in time. Auxiliary calculations are presented in
Appendix A.

###  2.1 Specialisation to two dimensions

To obtain the 2-D formulation of circular plate, we perform the series
expansion of $\mathbf{x}~{}\text{and}~{}p$ in terms of $Z$ about the bottom
surface, $Z=0$ following the approach by Wang et al. (2018, 2019b)

$\displaystyle\mathbf{x}(\mathbf{X})$
$\displaystyle=\mathbf{x}^{(0)}(\mathbf{\mathbf{\zeta}})+Z\mathbf{x}^{(1)}(\mathbf{\mathbf{\zeta}})+\frac{Z^{2}}{2}\mathbf{x}^{(2)}(\mathbf{\mathbf{\zeta}})+\frac{Z^{3}}{3!}\mathbf{x}^{(3)}(\mathbf{\mathbf{\zeta}})+\frac{Z^{4}}{4!}\mathbf{x}^{(4)}(\mathbf{\mathbf{\zeta}})+O(Z^{5}),$
(2.8) $\displaystyle p(\mathbf{X})$
$\displaystyle=p^{(0)}(\mathbf{\zeta})+Zp^{(1)}(\mathbf{\zeta})+\frac{Z^{2}}{2}p^{(2)}(\mathbf{\zeta})+\frac{Z^{3}}{3!}p^{(3)}(\mathbf{\zeta})+\frac{Z^{4}}{4!}p^{(4)}(\mathbf{\zeta})+O(Z^{5}),$
(2.9)

where $(\cdot)^{(n)}=\displaystyle\left.\frac{\partial^{n}{(\cdot)}}{\partial
Z^{n}}\right|_{Z=0}~{}(n=0,1,2,3,4)$. Likewise, we write the Taylor’s
expansion of the deformation gradient $\mathbf{F}$, the elastic deformation
tensor $\mathbf{A}$ and the inverse transpose of the growth tensor
$\mathbf{G}$ as

$\displaystyle\mathbf{F}$
$\displaystyle=\mathbf{F}^{(0)}(\zeta)+Z\mathbf{F}^{(1)}(\zeta)+\frac{Z^{2}}{2}\mathbf{F}^{(2)}(\zeta)+\frac{Z^{3}}{3!}\mathbf{F}^{(3)}(\zeta)+O(Z^{4}),$
(2.10a) $\displaystyle\mathbf{A}$
$\displaystyle=\mathbf{A}^{(0)}(\zeta)+Z\mathbf{A}^{(1)}(\zeta)+\frac{Z^{2}}{2}\mathbf{A}^{(2)}(\zeta)+\frac{Z^{3}}{3!}\mathbf{A}^{(3)}(\zeta)+O(Z^{4}),$
(2.10b) $\displaystyle\mathbf{G}^{-T}$
$\displaystyle=\bar{\mathbf{G}}^{(0)}(\zeta)+Z\bar{\mathbf{G}}^{(1)}(\zeta)+\frac{Z^{2}}{2}\bar{\mathbf{G}}^{(2)}(\zeta)+\frac{Z^{3}}{3!}\bar{\mathbf{G}}^{(3)}(\zeta)+O(Z^{4}).$
(2.10c)

The recursion relation for the expansion coefficients in (2.10a) is expressed
as

$\displaystyle\mathbf{F}^{(n)}=\nabla{\mathbf{x}^{(n)}}+\mathbf{x}^{(n+1)}\otimes\mathbf{k}.~{}~{}~{}({n=0,1,2,3}).$
(2.11)

Using (2.2) and (2.10a) - (2.10c) we obtain the expressions for higher
derivatives of $\mathbf{A}$ as

$\displaystyle\mathbf{A}^{(0)}$
$\displaystyle=\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}},$ (2.12a)
$\displaystyle\mathbf{A}^{(1)}$
$\displaystyle=\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(1)}^{T}}+\mathbf{F}^{(1)}\bar{\mathbf{G}}^{{(0)}^{T}},$
(2.12b) $\displaystyle\mathbf{A}^{(2)}$
$\displaystyle=\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(2)}^{T}}+2\mathbf{F}^{(1)}\bar{\mathbf{G}}^{{(1)}^{T}}+\mathbf{F}^{(2)}\bar{\mathbf{G}}^{{(0)}^{T}},$
(2.12c) $\displaystyle\mathbf{A}^{(3)}$
$\displaystyle=\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(3)}^{T}}+3\mathbf{F}^{(1)}\bar{\mathbf{G}}^{{(2)}^{T}}+3\mathbf{F}^{(2)}\bar{\mathbf{G}}^{{(1)}^{T}}+\mathbf{F}^{(3)}\bar{\mathbf{G}}^{{(0)}^{T}}.$
(2.12d)

Similarly, the expansion of the first Piola Kirchhoff stress tensor
$\mathbf{P}$ is

$\displaystyle\mathbf{P}(\mathbf{x},p)$
$\displaystyle=\mathbf{P}^{(0)}(\mathbf{x},p)+Z\mathbf{P}^{(1)}(\mathbf{x},p)+\frac{Z^{2}}{2!}\mathbf{P}^{(2)}(\mathbf{x},p)+\frac{Z^{3}}{3!}\mathbf{P}^{(3)}(\mathbf{x},p)+O(Z^{4}),$
(2.13)

which can be written in component form ($[\mathbf{P}]_{ij}=P_{ij}$) as

$\displaystyle P_{ij}^{(0)}$
$\displaystyle=J_{G}\bigg{[}\mathcal{A}^{(0)}-p^{(0)}\mathcal{L}^{(0)}\bigg{]}_{i\alpha}\bar{G}^{(0)}_{\alpha
j},$ (2.14) $\displaystyle P_{ij}^{(1)}$
$\displaystyle=J_{G}\left[\bigg{[}\left[\mathcal{A}^{(1)}-p^{(0)}\mathcal{L}^{(1)}\right]_{i\alpha
k\beta}A^{(1)}_{k\beta}-p^{(1)}\mathcal{L}^{(0)}_{i\alpha}\bigg{]}\bar{G}^{(0)}_{\alpha
j}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+}\left[\mathcal{A}^{(0)}-p^{(0)}\mathcal{L}^{(0)}\right]_{i\alpha}\bar{G}^{(1)}_{\alpha
j}\right],$ (2.15) $\displaystyle P^{(2)}_{ij}$
$\displaystyle=J_{G}\bigg{[}\bigg{[}\mathcal{A}^{(1)}_{ik\alpha\beta}A^{(2)}_{\alpha\beta}+\big{[}\mathcal{A}^{(2)}-p^{(0)}\mathcal{L}^{(2)}\big{]}_{ik\alpha\beta
mn}A^{(1)}_{\alpha\beta}A^{(1)}_{mn}-2p^{(1)}\mathcal{L}^{(1)}_{ik\alpha\beta}A^{(1)}_{\alpha\beta}-p^{(0)}\mathcal{L}^{(1)}_{ik\alpha\beta}A^{(2)}_{\alpha\beta}$
$\displaystyle-p^{(2)}\mathcal{L}^{(0)}_{ik}\bigg{]}\bar{G}^{(0)}_{kj}+\bigg{[}2\mathcal{A}^{(1)}_{ik\alpha\beta}A^{(1)}_{\alpha\beta}-2p^{(0)}\mathcal{L}^{(1)}_{ik\alpha\beta}A^{(1)}_{\alpha\beta}-2p^{(1)}\mathcal{L}^{(0)}_{ik}\bigg{]}\bar{G}^{(1)}_{kj}$
$\displaystyle+\bigg{[}\mathcal{A}^{(0)}-p^{(0)}\mathcal{L}^{(0)}\bigg{]}_{ik}\bar{G}^{(2)}_{kj}\bigg{]},$
(2.16)

with
$\displaystyle\boldsymbol{\mathcal{A}^{{i}}}(\mathbf{A}^{(0)})=\left.\frac{\partial^{i+1}\phi_{0}(\mathbf{A})}{\partial\mathbf{A}^{i+1}}\right|_{\mathbf{A}=\mathbf{A}^{(0)}}$
and
$\displaystyle\boldsymbol{\mathcal{L}}^{{i}}(\mathbf{A}^{(0)})=\left.\frac{\partial
L_{0}(\mathbf{A})}{\partial\mathbf{A}^{i+1}}\right|_{\mathbf{A}=\mathbf{A}^{(0)}}$.
Further mathematical details of above calculations are provided in Appendix B.

The two-dimensional governing system has $\mathbf{x}^{(\cdot)}$ and
${p}^{(\cdot)}$ 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$)

$\displaystyle\left.\mathbf{P}^{(0)}\mathbf{k}\right|_{Z=0}=\mathbf{P}^{(0)}(\mathbf{A})\mathbf{k}=-\mathbf{q}^{-},$
(2.17)

and at top surface ($Z=2h$)

$\displaystyle\left.\mathbf{P}\mathbf{k}\right|_{Z=2h}=\mathbf{P}^{(0)}\mathbf{k}+2h\mathbf{P}^{(1)}\mathbf{k}+2h^{2}\mathbf{P}^{(2)}\mathbf{k}+\frac{4}{3}h^{3}\mathbf{P}^{(3)}\mathbf{k}+O(h^{4})=\mathbf{q}^{+}.$
(2.18)

The stress equilibrium equation neglecting the body force and external
traction is given by (2.7a)

$\displaystyle\text{Div}\mathbf{P}=\mathbf{0},$ (2.19)
$\displaystyle\nabla\cdot{\mathbf{P}}+\frac{\partial}{\partial
Z}[{\mathbf{P}}\mathbf{k}]=\mathbf{0},$

where $\nabla$ is the two-dimensional differentiation operator. Upon use of
(2.13), we obtain a recursion relation for the first Piola–Kirchhoff stress as

$\displaystyle\nabla\cdot{\mathbf{P}}^{(n)}+{\mathbf{P}}^{(n+1)}\mathbf{k}=\mathbf{0}.$
(2.20)

The series expansion of unknown functions $\mathbf{x(X)}$ and $p(\mathbf{X})$
have 19 unknowns with $\mathbf{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 $\mathbf{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

$\displaystyle L_{0}(\mathbf{A}^{(0)})=0,$ (2.21a)
$\displaystyle\boldsymbol{\mathcal{L}}^{(0)}[\mathbf{A}^{(1)}]=0,$ (2.21b)
$\displaystyle\boldsymbol{\mathcal{L}}^{(0)}[\mathbf{A}^{(2)}]+\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(1)},~{}\mathbf{A}^{(1)}]=0,$
(2.21c)
$\displaystyle\boldsymbol{\mathcal{L}}^{(0)}[\mathbf{A}^{(3)}]+3\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(1)},~{}\mathbf{A}^{(2)}]+\boldsymbol{\mathcal{L}}^{(2)}[\mathbf{A}^{(1)},~{}\mathbf{A}^{(1)},~{}\mathbf{A}^{(1)}]=0,$
(2.21d)

where
$\boldsymbol{\mathcal{L}}^{(0)}[\mathbf{A}^{(1)}]=\boldsymbol{\mathcal{L}}^{(0)}:\mathbf{A}^{(1)}=\displaystyle\text{det}(\mathbf{A})\mathbf{A}^{-T}:\mathbf{A}^{(1)}$.
On subtracting the top (2.18) and bottom traction (2.17) condition we obtain
the equilibrium equation

$\displaystyle\nabla\cdot\bar{\mathbf{P}}=-\bar{\mathbf{q}},$ (2.22)

where

$\displaystyle\bar{\mathbf{P}}$
$\displaystyle=\frac{1}{2h}\int_{0}^{2h}\mathbf{P}dZ=\mathbf{P}^{(0)}+h\mathbf{P}^{(1)}+\frac{2}{3}h^{2}\mathbf{P}^{(2)}+O(h^{3}),$
$\displaystyle\bar{\mathbf{q}}$
$\displaystyle=\frac{\mathbf{q}^{+}+\mathbf{q}^{-}}{2h}.$

$\bar{\mathbf{P}}$ is the average stress over the thickness,
$\bar{\mathbf{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)

$\left.\begin{aligned}
\nabla\cdot\mathbf{P}^{(0)}_{t}+h\nabla\cdot\mathbf{P}^{(1)}_{t}+\frac{2}{3}h^{2}\nabla\cdot\mathbf{P}^{(2)}_{t}+O(h^{3})=-\bar{\mathbf{q}}_{t},\\\
\big{[}\nabla\cdot\bar{\mathbf{P}}\big{]}\cdot\mathbf{k}=\nabla\cdot\left[{\mathbf{P}^{(0)}}^{T}\mathbf{k}\right]+h\nabla\cdot\left[{\mathbf{P}^{(1)}}^{T}\mathbf{k}\right]+\frac{2}{3}h^{2}\nabla\cdot\left[{\mathbf{P}^{(2)}}^{T}\mathbf{k}\right]+O(h^{3})=-\bar{{q}}_{3},\end{aligned}\right\\}$
(2.23)

where the subscript ‘$t$’ represents the in-plane (or tangential) component of
a vector or tensor,
${\mathbf{q}}_{t}=q_{1}\mathbf{e}_{R}+q_{2}\mathbf{e}_{\Theta},~{}\bar{\mathbf{q}}_{t}=\displaystyle\frac{\mathbf{q}_{t}^{+}+\mathbf{q}_{t}^{-}}{2h}$
and $\bar{q}_{3}=\displaystyle\frac{q_{3}^{+}+q_{3}^{-}}{2h}$. Physically, Eq.
(2.23) represents the balance of forces. The explicit expressions for
$\mathbf{x}^{(2)}$, $p^{(1)}~{}\text{and}~{}\mathbf{x}^{(3)},~{}p^{(2)}$ in
terms of $\mathbf{x}^{(0)},~{}\mathbf{x}^{(1)}~{}\text{and}~{}p^{(0)}$ are
given as

$\displaystyle\mathbf{x}^{(2)}$
$\displaystyle=-\mathbf{B}^{-1}\mathbf{f}^{(2)}-D^{{(0)}^{-1}}\boldsymbol{\mathcal{L}}^{(0)}\left[\nabla\mathbf{x}^{(1)}\bar{\mathbf{G}}^{{(0)}^{T}}-\mathbf{B}^{-1}\mathbf{f}^{(2)}\otimes\bar{\mathbf{G}}^{(0)}\mathbf{k}+\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(1)}^{T}}\right]\mathbf{B}^{-1}\boldsymbol{\mathcal{L}}^{(0)}\widehat{\mathbf{G}}^{(0)}\mathbf{k},$
(2.24) $\displaystyle p^{(1)}$
$\displaystyle=-D^{{(0)}^{-1}}\boldsymbol{\mathcal{L}}^{(0)}\left[\nabla\mathbf{x}^{(1)}\bar{\mathbf{G}}^{{(0)}^{T}}-\mathbf{B}^{-1}\mathbf{f}^{(2)}\otimes\bar{\mathbf{G}}^{(0)}\mathbf{k}+\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(1)}^{T}}\right],$
(2.25) $\displaystyle\mathbf{x}^{(3)}$
$\displaystyle={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\mathbf{B}^{-1}\mathbf{f}^{(3)}+p^{(2)}\mathbf{B}^{-1}\boldsymbol{\mathcal{L}}^{(0)}\widehat{\mathbf{G}}^{(0)}\mathbf{k}+2p^{(1)}\mathbf{B}^{-1}\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(1)}]\widehat{\mathbf{G}}^{(0)}\mathbf{k}+2p^{(1)}\mathbf{B}^{-1}\boldsymbol{\mathcal{L}}^{(0)}\widehat{\mathbf{G}}^{(1)}\mathbf{k}},$
(2.26) $\displaystyle p^{(2)}$
$\displaystyle=-D^{{(0)}^{-1}}\bigg{[}\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(1)},\mathbf{A}^{(1)}]+\boldsymbol{\mathcal{L}}^{(0)}\left[-\mathbf{B}^{-1}\mathbf{f}^{(3)}\otimes\bar{\mathbf{G}}^{(0)}\mathbf{k}\right]$
$\displaystyle\quad+2p^{(1)}\boldsymbol{\mathcal{L}}^{(0)}\left[\mathbf{B}^{-1}\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(1)}]\widehat{\mathbf{G}}^{{(0)}}\mathbf{k}\otimes\bar{\mathbf{G}}^{{(0)}}\mathbf{k}\right]+2p^{(1)}\boldsymbol{\mathcal{L}}^{(0)}\left[\mathbf{B}^{-1}\boldsymbol{\mathcal{L}}^{(0)}\widehat{\mathbf{G}}^{(1)}\mathbf{k}\otimes\bar{\mathbf{G}}^{(0)}\mathbf{k}\right]$
$\displaystyle\quad+\boldsymbol{\mathcal{L}}^{(0)}\left[\nabla\mathbf{x}^{(2)}\bar{\mathbf{G}}^{{(0)}^{T}}+2\mathbf{F}^{(1)}\bar{\mathbf{G}}^{{(1)}^{T}}+\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(2)}^{T}}\right]\bigg{]},$
(2.27)

where $\widehat{\mathbf{G}}^{{(0)}}=J_{G}\bar{\mathbf{G}}^{{(0)}}$ and

$\displaystyle{B}_{\alpha\beta}=J_{G}\left[\boldsymbol{\mathcal{A}}^{(1)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(1)}\right]_{\alpha
i\beta
j}\left[\bar{\mathbf{G}}^{(0)}\mathbf{k}\right]_{i}\left[\bar{\mathbf{G}}^{(0)}\mathbf{k}\right]_{j},$
(2.28)
$\displaystyle{D}^{(0)}=\boldsymbol{\mathcal{L}}^{(0)}\left[J_{G}\mathbf{B}^{-1}\boldsymbol{\mathcal{L}}^{(0)}\bar{\mathbf{G}}^{(0)}\mathbf{k}\otimes\bar{\mathbf{G}}^{(0)}\mathbf{k}\right],$
(2.29)
$\displaystyle\mathbf{f}^{(2)}=\nabla~{}.~{}\mathbf{P}^{(0)}+\bigg{[}\left[\boldsymbol{\mathcal{A}}^{(1)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(1)}\right]\left[\nabla\mathbf{x}^{(1)}\bar{\mathbf{G}}^{{(0)}^{T}}+\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(1)}^{T}}\right]\bigg{]}\widehat{\mathbf{G}}^{(0)}\mathbf{k}+\left[\boldsymbol{\mathcal{A}}^{(0)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(0)}\right]\widehat{\mathbf{G}}^{(1)}\mathbf{k},$
(2.30)
$\displaystyle\mathbf{f}^{(3)}=\nabla~{}.~{}\mathbf{P}^{(1)}+\bigg{[}\left[\boldsymbol{\mathcal{A}}^{(1)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(1)}\right]\left[\nabla\mathbf{x}^{(2)}\bar{\mathbf{G}}^{{(0)}^{T}}+2\mathbf{F}^{(1)}\bar{\mathbf{G}}^{{(1)}^{T}}+\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(2)}^{T}}\right]\bigg{]}\widehat{\mathbf{G}}^{(0)}\mathbf{k}$
$\displaystyle\hskip
50.58878pt+\bigg{[}\left[\boldsymbol{\mathcal{A}}^{(2)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(2)}\right]\left[\mathbf{A}^{(1)},\mathbf{A}^{(1)}\right]\bigg{]}\widehat{\mathbf{G}}^{(0)}\mathbf{k}+\bigg{[}\left[\boldsymbol{\mathcal{A}}^{(1)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(1)}\right]\left[\mathbf{A}^{(1)}\right]\bigg{]}\widehat{\mathbf{G}}^{(1)}\mathbf{k}$
$\displaystyle\hskip
234.87749pt+\left[\boldsymbol{\mathcal{A}}^{(0)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(0)}\right]\widehat{\mathbf{G}}^{(2)}\mathbf{k}.$
(2.31)

##  3 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.

###  3.1  Radial growth

Consider a thin isotropic circular plate which undergoes axisymmetric
deformation ($\theta=\Theta$) with the position vector
$\mathbf{x}$($r,\theta,z$) in deformed configuration. The series expansion of
unknown functions $r(R,Z)$, $z(R,Z)~{}\text{and}~{}p(R,Z)$ in terms of $Z$ is
written as

$\displaystyle r(R,Z)$
$\displaystyle=r^{(0)}(R)+Zr^{(1)}(R)+\frac{1}{2!}Z^{2}r^{(2)}(R)+\frac{1}{3!}Z^{3}r^{(3)}(R)+\frac{1}{4!}Z^{4}r^{(4)}(R)+O(Z^{5}),$
(3.1) $\displaystyle z(R,Z)$
$\displaystyle=z^{(0)}(R)+Zz^{(1)}(R)+\frac{1}{2!}Z^{2}z^{(2)}(R)+\frac{1}{3!}Z^{3}z^{(3)}(R)+\frac{1}{4!}Z^{4}z^{(4)}(R)+O(Z^{5}),$
(3.2) $\displaystyle p(R,Z)$
$\displaystyle=p^{(0)}+Zp^{(1)}(R)+\frac{1}{2!}Z^{2}p^{(2)}(R)+\frac{1}{3!}Z^{3}p^{(3)}(R)+\frac{1}{4!}Z^{4}p^{(4)}(R)+O(Z^{5}),$
(3.3)

where we have used the notation
$r^{(n)}=\displaystyle\frac{\partial^{n}r}{\partial Z^{n}}$,
$z^{(n)}=\displaystyle\frac{\partial^{n}z}{\partial Z^{n}}$ and
$p^{(n)}=\displaystyle\frac{\partial^{n}p}{\partial Z^{n}}$.

We assume the plate to follow neo-Hookean elastic constitutive law given by

$\displaystyle\phi_{0}(\mathbf{A})=C_{0}\left[\text{tr}(\mathbf{A}^{T}\mathbf{A})-3\right],$
(3.4)

where $C_{0}$ is the ground state shear modulus. On substitution of (2.2) in
(3.4), the first Piola Kirchhoff stress $\mathbf{P}$ is obtained as (see
Appendix C )

$\displaystyle\mathbf{P}=J_{G}\bigg{[}2C_{0}\mathbf{A}-p\mathbf{A}^{-T}\bigg{]}\mathbf{G}^{-T}.$
(3.5)

We can rewrite (2.23) in component form as

$\displaystyle
P_{rR,R}^{(0)}+\frac{1}{R}P_{r\Theta,\Theta}^{(0)}+\frac{1}{R}\left[P_{rR}^{(0)}\right]+h\bigg{[}P_{rR,R}^{(1)}+\frac{1}{R}P_{r\Theta,\Theta}^{(1)}+\frac{1}{R}\left[P_{rR}^{(1)}\right]\bigg{]}+O(h^{2})=-\bar{{q}}_{1},$
(3.6) $\displaystyle P_{\theta
R,R}^{(0)}+\frac{1}{R}P_{\theta\Theta,\Theta}^{(0)}+\frac{1}{R}\left[P_{\theta
R}^{(0)}\right]+h\bigg{[}P_{\theta
R,R}^{(1)}+\frac{1}{R}P_{\theta\Theta,\Theta}^{(1)}+\frac{1}{R}\left[P_{\theta
R}^{(1)}\right]\bigg{]}+O(h^{2})=-\bar{{q}}_{2},$ (3.7) $\displaystyle
P_{zR,R}^{(0)}+\frac{1}{R}P_{z\theta,\Theta}^{(0)}+\frac{1}{R}\left[P_{zR}^{(0)}\right]+h\bigg{[}P_{zR,R}^{(1)}+\frac{1}{R}P_{z\Theta,\Theta}^{(1)}+\frac{1}{R}\left[P_{zR}^{(1)}\right]\bigg{]}+O(h^{2})=-\bar{{q}}_{3}.$
(3.8)

We can substitute the expression of stress from (2.22) and omit the higher
order terms in $h$ to get

$\displaystyle
2C_{0}\nabla\cdot\left[J_{G}\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}}\bar{\mathbf{G}}^{(0)}\right]-\nabla\cdot\left[p^{(0)}J_{G}\mathbf{F}^{{(0)}^{-T}}\bar{\mathbf{G}}^{(0)}\right]+O(h)=-\bar{\mathbf{q}}.$
(3.9)

Consider constant radial growth given by the growth tensor
$\mathbf{G}=\text{diag}(\lambda,1,1)$ that results in

$\displaystyle\big{[}\bar{\mathbf{G}}^{(0)}\big{]}=\begin{bmatrix}\displaystyle\frac{1}{\lambda}&0&0\\\
0&1&0\\\
0&0&1\end{bmatrix},\quad\big{[}\mathbf{F}^{(0)}\big{]}=\begin{bmatrix}{r^{(0)}}^{\prime}&0&r^{(1)}\\\
0&\displaystyle\frac{r^{(0)}}{R}&0\\\
{z^{(0)}}^{\prime}&0&z^{(1)}\end{bmatrix},\quad
J_{G}=\text{det}({\mathbf{G}})=\lambda.$ (3.10)

where a superposed prime denotes partial derivative with respect to $R$.
Substituting (3.10) in (3.9) we obtain the governing equation

$\displaystyle
2C_{0}\nabla\cdot\lambda\begin{bmatrix}{r^{(0)}}^{\prime}&0&r^{(1)}\\\
0&\displaystyle\frac{r^{(0)}}{R}&0\\\
{z^{(0)}}^{\prime}&0&z^{(1)}\end{bmatrix}\begin{bmatrix}\displaystyle\frac{1}{\lambda^{2}}&0&0\\\
0&1&0\\\ 0&0&1\end{bmatrix}~{}-~{}\nabla\cdot
p^{(0)}\lambda\frac{\text{cofac}(\mathbf{F}^{(0)})}{\det\mathbf{F}^{(0)}}+O(h)=-\bar{\mathbf{q}}.$
(3.11)

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.

####  3.1.1 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
traction-free and the bottom surface, is supported by the Winkler foundation
which provides a transverse load ${q}_{3}^{+}=-K_{0}\lambda W_{0}$, where
$K_{0}$ is the elastic constant of the foundation and $W_{0}$ is the
transverse ($Z$) component of the displacement. $\lambda$ 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.

Figure 2: Circular plate in reference configuration growing with growth factor
$\lambda$ resting on a Winkler foundation. (a) Perspective view (b) Front
view.

For this case, the components of the effective body force $(\bar{\mathbf{q}}$)
are given by

$\displaystyle\bar{q}_{1}$ $\displaystyle=0,$ (3.12)
$\displaystyle\bar{q}_{3}$
$\displaystyle=-\frac{1}{2h}K_{0}\lambda\bigg{[}z^{(0)}+2hz^{(1)}+\frac{1}{2}~{}[2h]^{2}z^{(2)}+\frac{1}{6}~{}[2h]^{3}z^{(3)}-2h\bigg{]}.$
(3.13)

The governing equations are obtained by substituting (3.13) in (3.11)

$\displaystyle\frac{\partial}{\partial
R}\left[2C_{0}\frac{{r^{(0)}}^{\prime}}{\lambda}-p^{(0)}\frac{r^{(0)}z^{(1)}}{R}\right]+\frac{1}{R}\bigg{[}P_{rR}^{(0)}-p^{(0)}\left[\frac{r^{(0)}z^{(1)}}{R}\right]\bigg{]}+O(h)=0,$
(3.14) $\displaystyle\frac{\partial}{\partial
R}\bigg{[}2C_{0}\frac{{z^{(0)}}^{\prime}}{\lambda}-p^{(0)}\left[\frac{-r^{(0)}r^{(1)}}{R}\right]\bigg{]}+\frac{1}{R}\bigg{[}P_{zR}^{(0)}-p^{(0)}\left[\frac{-r^{(0)}r^{(1)}}{R}\right]\bigg{]}+O(h)$
$\displaystyle\hskip
144.54pt=-K_{0}\lambda\bigg{[}\frac{z^{(0)}}{2h}+z^{(1)}+hz^{(2)}+\frac{2}{3}h^{2}z^{(3)}-1\bigg{]}.$
(3.15)

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

$r^{(0)}(0)=0,\quad r^{(0)}(R_{0})=R_{0},\quad
r^{(0)}(0)+[2h]r^{(1)}(0)+\frac{1}{2}[2h^{2}]r^{(2)}(0)=0,$ (3.16)

where the unknown variables in terms of $r^{(0)}$ and $z^{(0)}$ are given as

$\displaystyle
r^{(1)}=\displaystyle-\frac{\lambda{z^{(0)}}^{\prime}}{\displaystyle\frac{r^{(0)}}{R}\left[{z^{(0)}}^{\prime
2}+{r^{(0)}}^{\prime 2}\right]},\quad z^{(1)}$
$\displaystyle=\frac{\lambda{r^{(0)}}^{\prime}}{\displaystyle\frac{r^{(0)}}{R}\left[{z^{(0)}}^{\prime
2}+{r^{(0)}}^{\prime 2}\right]},\quad
p^{(0)}=\frac{2C_{0}\lambda^{2}}{{\left[\displaystyle\frac{r^{(0)}}{R}\right]^{2}\left[{z^{(0)}}^{\prime
2}+{r^{(0)}}^{\prime 2}\right]}}.$ (3.17)

The explicit expressions for $r^{(1)},~{}z^{(1)},~{}p^{(0)}$ are derived in
Appendix C.

####  3.1.2 Stability analysis

The principal solution for plate deformation is given by

$\displaystyle r^{(0)}(R)=R,\quad z^{(0)}(R)=-2h[\lambda-1],$ (3.18)

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)}=\lambda$ (because
${r^{(0)}}^{\prime}=1~{}\text{and}~{}{z^{(0)}}^{\prime}=0$). A small
perturbation in the homogeneous equilibrium state by a parameter $\epsilon$
results in

$\displaystyle r^{(0)}(R)=R+\epsilon\Delta U(R)\quad\text{and}\quad
z^{(0)}(R)=-2h[\lambda-1]+\epsilon\Delta W(R).$ (3.19)

We define

$\displaystyle\mathbf{m}=\frac{\mathbf{q}^{+}-\mathbf{q}^{-}}{2}=\mathbf{P}^{(0)}\mathbf{k}+h\mathbf{P}^{(1)}\mathbf{k}+h^{2}\mathbf{P}^{(2)}\mathbf{k}+O(h^{3}).$
(3.20)

To simplify our calculations, we eliminate the terms with $\mathbf{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

$\displaystyle\nabla\cdot\mathbf{P}^{(0)}_{t}+h\nabla\cdot\mathbf{P}^{(1)}_{t}=-\bar{\mathbf{q}}_{t},$
(3.21a)
$\displaystyle\nabla\cdot\left[({\mathbf{P}^{(0)}}^{T}\mathbf{k})-({\mathbf{P}^{(0)}}\mathbf{k})\right]+h\bigg{[}\nabla\cdot\left[({\mathbf{P}^{(1)}}^{T}\mathbf{k})-({\mathbf{P}^{(1)}}\mathbf{k})\right]\bigg{]}$
$\displaystyle\hskip
144.54pt+\frac{1}{3}h^{2}\nabla\cdot[\nabla\cdot{\mathbf{P}_{t}}^{(1)}]=-\bar{q}_{3}-\nabla\cdot\mathbf{m}_{t},$
(3.21b)

where $\mathbf{m}_{t}$ is the tangential component of $\mathbf{m}$.
Physically, Eq. (3.20) represents the balance of moments.

We define the dimensionless quantities

$\displaystyle\rho=\frac{R}{R_{0}},\hskip
7.22743pt\bar{h}=\frac{h}{R_{0}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}\hskip
7.22743ptU=\frac{\Delta U}{R_{0}},\hskip 7.22743ptW=\frac{\Delta W}{R_{0}},$
(3.22)

where $\rho\in[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(\epsilon)$ terms

$\displaystyle\frac{2}{\lambda}[1+3\lambda^{4}]{U^{\prime\prime}}+\frac{2}{\rho\lambda}[1+5\lambda^{4}]U^{\prime}-\bar{h}\bigg{[}4[1+\lambda^{4}]W^{\prime\prime\prime}+\frac{2}{\rho}\big{[}1+5\lambda^{4}\big{]}\bigg{]}=0.$
(3.23)
$\displaystyle\frac{2}{\lambda}[1-\lambda^{4}]{W^{\prime\prime}}+\frac{2}{\rho\lambda}[1-\lambda^{4}]W^{\prime}$
$\displaystyle-\bar{h}\bigg{[}\frac{4}{\rho\lambda}\big{[}4\lambda^{5}-4\lambda^{4}+\lambda-1\big{]}U^{\prime\prime}+\frac{4}{\rho^{2}}\big{[}\lambda^{4}-1\big{]}U^{\prime}-\frac{4}{\rho^{3}}\big{[}\lambda^{4}-1\big{]}U\bigg{]}$
$\displaystyle-\frac{2}{3}\bar{h}^{2}\bigg{[}2[1+\lambda^{4}]W^{iv}+\frac{1}{\rho}\big{[}3+7\lambda^{4}\big{]}W^{\prime\prime\prime}\bigg{]}$
$\displaystyle-\frac{\beta\lambda}{2\bar{h}}{W}+\beta\lambda^{2}U^{\prime}+\frac{\beta\lambda^{2}}{\rho}U-\frac{\bar{h}\beta[2\lambda^{4}-1]}{\rho\lambda}W^{\prime}-\bar{h}\beta\lambda^{3}W^{\prime\prime}$
$\displaystyle-\frac{2}{3}\bar{h}^{2}\beta\bigg{[}\frac{4}{\rho}\big{[}1+\lambda^{4}\big{]}U^{\prime\prime}-\frac{1}{\rho^{2}}\big{[}2\lambda^{4}-3\big{]}U^{\prime}+\frac{1}{\rho^{3}}\big{[}2\lambda^{4}-3\big{]}U\bigg{]}=0.$
(3.24)

where $\beta=\displaystyle\frac{K_{0}}{C_{0}}R_{0}$ is a non-dimensional
constant. The higher order derivative of in-plane displacement term is omitted
for simplification. The plate boundary conditions (3.16) is given as

$\displaystyle U(0)=W^{\prime}(0)=W^{\prime\prime\prime}(0)=0,$ (3.25)
$\displaystyle U(1)=W(1)=W^{\prime\prime}(1)=0.$

Equations (3.23) and (3.24) are coupled ODEs and can be written as a system of
first order ODEs by defining

$\displaystyle U=y_{1},\qquad U^{\prime}=y_{2},\qquad W=y_{3},\qquad
W^{\prime}=y_{4},\qquad W^{\prime\prime}=y_{5},\qquad
W^{\prime\prime\prime}=y_{6}.$ (3.26)

The system of first order ordinary differential equations is written in the
form of

$\displaystyle\mathbf{Y}^{\prime}=\boldsymbol{\mathscr{A}}(\rho;\lambda,\bar{h},\beta)\mathbf{Y},$
(3.27)

where $\mathbf{Y}=[\Delta U,~{}\Delta U^{\prime},~{}\Delta
U^{\prime\prime},~{}\Delta W,~{}\Delta W^{\prime},~{}\Delta W^{\prime\prime}]$
= $[y_{1}~{}y_{2}~{}y_{3}~{}y_{4}~{}y_{5}~{}y_{6}]$ and
$\boldsymbol{\mathscr{A}}$ is given by

$\displaystyle\boldsymbol{\mathscr{A}}=\begin{bmatrix}0&1&0&0&0&0\\\
\mathscr{A}_{21}&\mathscr{A}_{22}&\mathscr{A}_{23}&\mathscr{A}_{24}&\mathscr{A}_{25}&\mathscr{A}_{26}\\\
0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\
\mathscr{A}_{61}&\mathscr{A}_{62}&\mathscr{A}_{63}&\mathscr{A}_{64}&\mathscr{A}_{65}&\mathscr{A}_{66}\end{bmatrix},$
(3.28)

where

$\displaystyle\mathscr{A}_{21}=0,\quad\mathscr{A}_{22}=\frac{\lambda}{2(1+3\lambda^{4})}\bigg{[}-\frac{2}{\rho\lambda}[1+5\lambda^{4}]\bigg{]},\quad\mathscr{A}_{23}=0,\quad\mathscr{A}_{24}=0,$
$\displaystyle\mathscr{A}_{25}=\frac{\lambda}{2(1+5\lambda^{4})}\bigg{[}\frac{2\bar{h}}{\rho}\big{[}1+5\lambda^{4}\big{]}\bigg{]},\qquad\mathscr{A}_{26}=\frac{\lambda}{2(1+3\lambda^{4})}\bigg{[}4\bar{h}[1+\lambda^{4}]\bigg{]}$
$\displaystyle\mathscr{A}_{61}=\frac{3}{4\bar{h^{2}}[1+\lambda^{4}]}\Bigg{[}\frac{4\bar{h}}{\rho^{3}}[\lambda^{4}-1]+\frac{\beta\lambda^{2}}{\rho}-\frac{2}{3}\frac{\bar{h}^{2}\beta}{\rho^{3}}[2\lambda^{4}-3]\Bigg{]},$
$\displaystyle\mathscr{A}_{62}=\frac{3}{4\bar{h^{2}}[1+\lambda^{4}]}\Bigg{[}-\frac{4h}{\rho\lambda}\left[\frac{\lambda[4\lambda^{5}-4\lambda^{4}+\lambda-1]}{2[1+3\lambda^{4}]}\right]\bigg{[}-\frac{2}{\rho\lambda}[1+5\lambda^{4}]\bigg{]}-\frac{4h}{\rho^{2}}[\lambda^{4}-1]$
$\displaystyle\qquad\qquad+\beta\lambda^{2}-\frac{2}{3}\bar{h}^{2}\beta\bigg{[}\frac{4\lambda[1+\lambda^{4}]}{2\rho[1+3\lambda^{4}]}\bigg{[}-\frac{2}{\rho\lambda}[1+5\lambda^{4}]\bigg{]}\bigg{]}+\frac{2}{3}\bar{h}^{2}\frac{\beta}{\rho^{2}}[2\lambda^{4}-3]\Bigg{]},$
$\displaystyle\mathscr{A}_{63}=\frac{3}{4\bar{h^{2}}[1+\lambda^{4}]}\bigg{[}-\frac{\beta\lambda}{2\bar{h}}\bigg{]},\qquad\mathscr{A}_{64}=\frac{3}{4\bar{h^{2}}[1+\lambda^{4}]}\bigg{[}\frac{2}{\rho\lambda}[1-\lambda^{4}]-\frac{h\beta}{\rho\lambda}[2\lambda^{4}-1]\bigg{]},$
$\displaystyle\mathscr{A}_{65}=\frac{3}{4\bar{h^{2}}[1+\lambda^{4}]}\Bigg{[}\frac{2}{\lambda}[1-\lambda^{4}]-\frac{4h}{\rho\lambda}\left[\frac{\lambda[4\lambda^{5}-4\lambda^{4}+\lambda-1]}{2[1+3\lambda^{4}]}\right]\bigg{[}\frac{2\bar{h}}{\rho}[1+5\lambda^{4}]\bigg{]}-h\beta\lambda^{3}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\frac{2}{3}\bar{h}^{2}\beta\bigg{[}\frac{4\lambda[1+\lambda^{4}]}{2\rho[1+3\lambda^{4}]}\bigg{[}\frac{2\bar{h}}{\rho}[1+5\lambda^{4}]\bigg{]}\bigg{]}\Bigg{]},$
$\displaystyle\mathscr{A}_{66}=\frac{3}{4\bar{h^{2}}[1+\lambda^{4}]}\Bigg{[}-\frac{4h}{\rho\lambda}\left[\frac{\lambda[4\lambda^{5}-4\lambda^{4}+\lambda-1]}{2[1+3\lambda^{4}]}\right]\bigg{[}4\bar{h}[1+\lambda^{4}]\bigg{]}-\frac{2}{3}\frac{\bar{h}^{2}}{\rho}[3+7\lambda^{4}]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\frac{2}{3}\bar{h}^{2}\beta\bigg{[}\frac{4\lambda[1+\lambda^{4}]}{2\rho[1+3\lambda^{4}]}\bigg{]}\bigg{[}4\bar{h}[1+\lambda^{4}]\bigg{]}\Bigg{]}.$

We first determine the critical value of growth factor ($\lambda_{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 (Ng and 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
$\mathbf{\Phi}^{\prime}=\boldsymbol{\mathscr{A}}^{*}(\rho;\lambda,\bar{h},\beta)\mathbf{\Phi}$
(A detailed description of the solution procedure is given in Appendix D)

$\displaystyle\Phi_{1}^{\prime}$
$\displaystyle=\Phi_{2}+\mathscr{A}_{22}\Phi_{1}-\mathscr{A}_{24}\Phi_{5}-\mathscr{A}_{25}\Phi_{6}-\mathscr{A}_{26}\Phi_{7},$
$\displaystyle\Phi_{2}^{\prime}$
$\displaystyle=\Phi_{3}+\mathscr{A}_{22}\Phi_{2}+\mathscr{A}_{23}\Phi_{5}-\mathscr{A}_{25}\Phi_{8}-\mathscr{A}_{26}\Phi_{9},$
$\displaystyle\Phi_{3}^{\prime}$
$\displaystyle=\Phi_{4}+\mathscr{A}_{22}\Phi_{3}+\mathscr{A}_{23}\Phi_{6}+\mathscr{A}_{24}\Phi_{8}-\mathscr{A}_{26}\Phi_{10},$
$\displaystyle\Phi_{4}^{\prime}$
$\displaystyle=\mathscr{A}_{22}\Phi_{4}+\mathscr{A}_{23}\Phi_{7}+\mathscr{A}_{24}\Phi_{9}+\mathscr{A}_{25}\Phi_{10}+\mathscr{A}_{63}\Phi_{1}+\mathscr{A}_{64}\Phi_{2}+\mathscr{A}_{65}\Phi_{3}+\mathscr{A}_{66}\Phi_{4},$
$\displaystyle\Phi_{5}^{\prime}$ $\displaystyle=\Phi_{11}+\Phi_{6},$
$\displaystyle\Phi_{6}^{\prime}$ $\displaystyle=\Phi_{12}+\Phi_{8}+\Phi_{7},$
$\displaystyle\Phi_{7}^{\prime}$
$\displaystyle=\Phi_{13}+\Phi_{9}-\mathscr{A}_{62}\Phi_{1}+\mathscr{A}_{64}\Phi_{5}+\mathscr{A}_{65}\Phi_{6}+\mathscr{A}_{66}\Phi_{7},$
$\displaystyle\Phi_{8}^{\prime}$ $\displaystyle=\Phi_{14}+\Phi_{9},$
$\displaystyle\Phi_{9}^{\prime}$
$\displaystyle=\Phi_{15}+\Phi_{10}-\mathscr{A}_{62}\Phi_{2}-\mathscr{A}_{63}\Phi_{5}+\mathscr{A}_{65}\Phi_{8}+\mathscr{A}_{66}\Phi_{9},$
$\displaystyle\Phi_{10}^{\prime}$
$\displaystyle=\Phi_{16}-\mathscr{A}_{62}\Phi_{3}-\mathscr{A}_{63}\Phi_{6}-\mathscr{A}_{64}\Phi_{8}+\mathscr{A}_{66}\Phi_{10},$
(3.29) $\displaystyle\Phi_{11}^{\prime}$
$\displaystyle=\Phi_{12}+\mathscr{A}_{21}\Phi_{5}+\mathscr{A}_{22}\Phi_{11}+\mathscr{A}_{25}\Phi_{17}+\mathscr{A}_{26}\Phi_{18},$
$\displaystyle\Phi_{12}^{\prime}$
$\displaystyle=\Phi_{14}+\Phi_{13}+\mathscr{A}_{21}\Phi_{6}+\mathscr{A}_{22}\Phi_{12}-\mathscr{A}_{24}\Phi_{17}+\mathscr{A}_{26}\Phi_{19},$
$\displaystyle\Phi_{13}^{\prime}$
$\displaystyle=\Phi_{15}+\mathscr{A}_{21}\Phi_{7}+\mathscr{A}_{22}\Phi_{13}-\mathscr{A}_{24}\Phi_{18}-\mathscr{A}_{25}\Phi_{19}+\mathscr{A}_{61}\Phi_{1}+\mathscr{A}_{64}\Phi_{11}+\mathscr{A}_{65}\Phi_{12}+\mathscr{A}_{66}\Phi_{13},$
$\displaystyle\Phi_{14}^{\prime}$
$\displaystyle=\Phi_{15}+\mathscr{A}_{21}\Phi_{8}+\mathscr{A}_{22}\Phi_{14}+\mathscr{A}_{23}\Phi_{17}+\mathscr{A}_{26}\Phi_{20},$
$\displaystyle\Phi_{15}^{\prime}$
$\displaystyle=\Phi_{16}+\mathscr{A}_{21}\Phi_{9}+\mathscr{A}_{22}\Phi_{15}+\mathscr{A}_{23}\Phi_{18}-\mathscr{A}_{25}\Phi_{20}+\mathscr{A}_{61}\Phi_{2}-\mathscr{A}_{63}\Phi_{11}+\mathscr{A}_{65}\Phi_{14}+\mathscr{A}_{66}\Phi_{15},$
$\displaystyle\Phi_{16}^{\prime}$
$\displaystyle=\mathscr{A}_{21}\Phi_{10}+\mathscr{A}_{22}\Phi_{16}+\mathscr{A}_{23}\Phi_{19}+\mathscr{A}_{24}\Phi_{20}+\mathscr{A}_{61}\Phi_{3}-\mathscr{A}_{63}\Phi_{12}-\mathscr{A}_{64}\Phi_{14}+\mathscr{A}_{66}\Phi_{16},$
$\displaystyle\Phi_{17}^{\prime}$ $\displaystyle=\Phi_{18},$
$\displaystyle\Phi_{18}^{\prime}$
$\displaystyle=\Phi_{19}+\mathscr{A}_{61}\Phi_{5}+\mathscr{A}_{62}\Phi_{11}+\mathscr{A}_{65}\Phi_{17}+\mathscr{A}_{66}\Phi_{18},$
$\displaystyle\Phi_{19}^{\prime}$
$\displaystyle=\Phi_{20}+\mathscr{A}_{61}\Phi_{6}+\mathscr{A}_{62}\Phi_{12}-\mathscr{A}_{64}\Phi_{17}+\mathscr{A}_{66}\Phi_{19},$
$\displaystyle\Phi_{20}^{\prime}$
$\displaystyle=\mathscr{A}_{61}\Phi_{8}+\mathscr{A}_{62}\Phi_{14}+\mathscr{A}_{63}\Phi_{17}+\mathscr{A}_{66}\Phi_{20}.$

The initial condition for the system of equations (3.29) is

$\displaystyle\Phi(0)=\begin{bmatrix}\Phi_{1},\Phi_{2},\Phi_{3},\Phi_{4},\Phi_{5},\Phi_{6},\Phi_{7},\Phi_{8},\Phi_{9},\Phi_{10},\Phi_{11},\Phi_{12},\Phi_{13},\Phi_{14},\Phi_{15},\Phi_{16},\Phi_{17},\Phi_{18},\Phi_{19},\Phi_{20}\end{bmatrix}.$
(3.30)

The target condition is achieved by having $\det(\mathbf{CM})=0$ in order to
obtain the non-trivial solution, where matrix $\mathbf{C}$ corresponds to the
boundary condition at the edge of circular plate (3.25) and $\mathbf{M}$ is
the solution matrix

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbf{C}=\begin{bmatrix}1&0&0&0&0&0\\\
0&0&1&0&0&0\\\
0&0&0&0&1&0\end{bmatrix}}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}\mathbf{M}=\begin{bmatrix}y_{1}^{(1)}&y_{1}^{(2)}&y_{1}^{(3)}\\\
y_{2}^{(1)}&y_{2}^{(2)}&y_{2}^{(3)}\\\ y_{3}^{(1)}&y_{3}^{(2)}&y_{3}^{(3)}\\\
y_{4}^{(1)}&y_{4}^{(2)}&y_{4}^{(3)}\\\ y_{5}^{(1)}&y_{5}^{(2)}&y_{5}^{(3)}\\\
y_{6}^{(1)}&y_{6}^{(2)}&y_{6}^{(3)}\end{bmatrix}.$ (3.31)

For the current case, the corresponding initial conditions using (3.25) are
given by $\Phi(12)=1$ and rest all are zero. Then, we integrate the system
numerically in the interval of $0<\rho\leq 1$ until we achieve the target
condition on the other boundary which is given by
$\det(\mathbf{CM})$=(1,3,5)=$\Phi(6)=0$ (see Appendix D). The main objective
of this optimization problem is to determine the critical value of growth
factor $\lambda_{cr}$ for which the target value $\Phi(6)$ is zero.

####  3.1.3 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=\pm 1$
and supported by a Winkler foundation as shown in Figure 3. The plate grows
along the X-axis with a growth stretch $\lambda$. The compression effects due
to the clamped boundary condition lead to buckling.

Figure 3: Schematic of a rectangular plate resting on Winkler foundation and
clamped at $X=\pm 1$.

The governing plate differential equation is given by (Wang et al., 2018)

$\displaystyle\psi_{0}\Delta W+\psi_{2}\Delta W^{\prime\prime}+\psi_{4}\Delta
W^{\prime\prime\prime\prime}=0,$ (3.32)

where

$\displaystyle\psi_{0}$ $\displaystyle=\frac{\alpha}{2h},$
$\displaystyle\psi_{2}$
$\displaystyle=\frac{1}{\lambda^{2}+3\lambda^{6}}\bigg{[}\lambda^{4}-1\left[2+[6+h\alpha]\lambda^{4}\right]\bigg{]},$
$\displaystyle\psi_{4}$
$\displaystyle=\frac{4h^{2}}{3+9\lambda^{4}}\bigg{[}3+h\alpha+[2+3h\alpha]\lambda^{4}+[3+2h\alpha]\lambda^{8}\bigg{]},$

subjected to the boundary conditions

$\displaystyle\Delta W^{\prime}(-1)=\Delta
W^{\prime}(1)=0,~{}~{}~{}~{}~{}~{}~{}\Delta W^{\prime\prime\prime}(-1)=\Delta
W^{\prime\prime\prime}(1)=0.$ (3.33)

---
Figure 4:  a) Bifurcation curves for various modes of rectangular plate (Wang
et al., 2018). b) Variation of the critical value of growth factor
($\lambda_{cr}$) against thickness to length ratio ($h$) evaluated using the
compound matrix method. c) Direct comparison of (a) and (b).

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,\dotsc$) are
more stable. Figure 4b represents the variation of critical value of growth
factor ($\lambda_{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 $\lambda_{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

Figure 5:  a) Variation of critical growth factor $\lambda_{cr}$ with
dimensionless thickness $\bar{h}$ for $\beta=0.2$. b) Enlarged part of (a).
The thickness values $\bar{h}_{1}$ and $\bar{h}_{2}$ correspond to the
critical points of mode switching.

The compound matrix method is used to compute the buckling parameter
($\lambda_{cr}$) of the circular hyperelastic plate under radial growth
condition. The variation of critical growth factor ($\lambda_{cr}$) with
respect to normalized thickness is shown in Figure 5. The results show the
non-monotonous nature of $\lambda_{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.

To obtain the mode shapes, the Matlab ODE package bvp4c is used for solving
the equations (3.23) and (3.24) subjected to boundary condition (3.25) at the
critical point. Equations (3.23) and (3.24) are rewritten into a first order
form of $\mathbf{Y}^{\prime}=\boldsymbol{\mathscr{A}}\mathbf{Y}$, where
$\mathbf{Y}=[y1,~{}y2,~{}y3,~{}y4,~{}y5,~{}y6]=[U,~{}U^{\prime},~{}W,~{}W^{\prime},~{}W^{\prime\prime},~{}W^{\prime\prime\prime}]$
where $\boldsymbol{\mathscr{A}}$ is given by (3.28) and the boundary condition
(3.25) is rewritten as

$\displaystyle
y_{a}(1)=0,~{}y_{b}(1)=0,~{}y_{a}(4)=0,~{}y_{b}(3)=0,~{}y_{a}(6)=0,~{}y_{b}(5)=0.$
(3.34)

where $y_{a}$ (respectively, $y_{b}$) define the centre of the plate
(respectively, edge of the plate).

Figure 6:  Normalised bifurcation mode shapes for circular plate under radial
growth (a) Normalised out-of-plane displacement $W$. b) Normalised in-plane
displacement $U$.

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 thickness
$\bar{h}=0.024$ is chosen for which $\lambda_{cr}=1.0218$ is used to evaluate
the mode shape of the circular hyperelastic plate. At this value of $\bar{h}$,
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 choosing  $\bar{h}=0.05,\lambda_{cr}=1.0350$ and
$\bar{h}=0.2,\lambda_{cr}=1.0611,$ respectively. Thus, we conclude that the
buckling parameter ($\lambda_{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.

###  3.2 Combined radial and circumferential growth

In this second case,  we assume a constant isotropic growth function $\lambda$
(i.e., $\lambda_{rr}=\lambda_{\theta\theta}=\lambda$) in both radial and
circumferential direction (Wu and Ben Amar, 2015). The isotropic growth tensor
$\mathbf{G}=\text{diag}(\lambda,\lambda,1)$ results in the following
kinematics

$\displaystyle\big{[}\bar{\mathbf{G}}^{(0)}\big{]}=\begin{bmatrix}\displaystyle\frac{1}{\lambda}&0&0\\\
0&\displaystyle\frac{1}{\lambda}&0\\\
0&0&1\end{bmatrix},\quad\big{[}\mathbf{F}^{(0)}\big{]}=\begin{bmatrix}\displaystyle\frac{\partial
r^{(0)}}{\partial R}&\displaystyle\frac{1}{R}\frac{\partial
r^{(0)}}{\partial\Theta}&r^{(1)}\vspace{5pt}\\\ \displaystyle
r^{(0)}\frac{\partial\theta}{\partial
R}&\displaystyle\frac{r^{(0)}}{R}\frac{\partial\theta}{\partial\Theta}&r^{(0)}\theta^{(1)}\vspace{5pt}\\\
\displaystyle\frac{\partial z^{(0)}}{\partial
R}&\displaystyle\frac{1}{R}\frac{\partial
z^{(0)}}{\partial\Theta}&z^{(1)}\end{bmatrix},\quad
J_{G}=\text{det}({\mathbf{G}})=\lambda^{2}.$ (3.35)

On substituting (3.35) in the governing equation (3.9), we obtain

$\displaystyle
2C_{0}\nabla\cdot\lambda^{2}\begin{bmatrix}\displaystyle\frac{\partial
r^{(0)}}{\partial R}&\displaystyle\frac{1}{R}\frac{\partial
r^{(0)}}{\partial\Theta}&r^{(1)}\vspace{5pt}\\\ \displaystyle
r^{(0)}\frac{\partial\theta}{\partial
R}&\displaystyle\frac{r^{(0)}}{R}\frac{\partial\theta}{\partial\Theta}&r^{(0)}\theta^{(1)}\vspace{5pt}\\\
\displaystyle\frac{\partial z^{(0)}}{\partial
R}&\displaystyle\frac{1}{R}\frac{\partial
z^{(0)}}{\partial\Theta}&z^{(1)}\end{bmatrix}\begin{bmatrix}\displaystyle\frac{1}{\lambda^{2}}&0&0\\\
0&\displaystyle\frac{1}{\lambda^{2}}&0\\\
0&0&1\end{bmatrix}-\nabla\cdot\left[p^{(0)}\lambda^{2}\frac{\text{cofac}(\mathbf{F}^{(0)})}{\det\mathbf{F}^{(0)}}\right]$
$\displaystyle+~{}O(h)=-\bar{\mathbf{q}},$ (3.36)

where the unknown variables in this case are calculated as

$\displaystyle
p^{(0)}=\frac{2C_{0}\lambda^{4}}{\left|\nabla{\mathbf{x}^{(0)}}^{**}\right|^{2}},\quad
r^{(1)}=\displaystyle\frac{p^{(0)}\nabla
x_{11}}{2C_{0}\lambda^{2}},\quad\theta^{(1)}=\frac{p^{(0)}\nabla
x_{22}}{2C_{0}\lambda^{2}r^{(0)}},\quad z^{(1)}=\frac{p^{(0)}\nabla
x_{33}}{2C_{0}\lambda^{2}}.$

The explicit expressions for $r^{(1)},\theta^{(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).

####  3.2.1 Linear buckling analysis for combined growth model

We conduct the linear bifurcation analysis by perturbing the principal
solution with small parameter ($\epsilon$) to determine the onset of
instability by assuming the form

$\displaystyle{r}^{(0)}({R},\Theta)$ $\displaystyle={R}+\epsilon\Delta
U(R)\cos(m\Theta),$ $\displaystyle\theta$ $\displaystyle=\Theta,$ (3.37)
$\displaystyle z^{(0)}(R,\Theta)$
$\displaystyle=-2h[\lambda^{2}-1]+\epsilon\Delta W(R)\cos(m\Theta).$

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}\lambda^{2}W_{0}$ and collecting only $O(\epsilon)$ terms, we
obtain the coupled differential equations for dimensionless displacement
functions $U$ and $W$ as

$\displaystyle
2[1+3\lambda^{6}]U^{\prime\prime}+\frac{2}{\rho}[1+5\lambda^{6}]U^{\prime}-\frac{2}{\rho^{2}}[m^{2}]U$
$\displaystyle-\bar{h}\bigg{[}4\lambda^{2}[1+\lambda^{6}]W^{\prime\prime\prime}+\frac{2}{\rho}\lambda^{2}\big{[}1+5\lambda^{6}\big{]}W^{\prime\prime}+\frac{2}{\rho^{2}}m^{2}\lambda^{2}[\lambda^{6}-2]W^{\prime}$
$\displaystyle-\frac{2}{\rho^{3}}\lambda^{2}m^{2}[\lambda^{6}-1]W\bigg{]}=0,$
(3.38) $\displaystyle
2[1-\lambda^{6}]W^{\prime\prime}+\frac{2}{\rho}\left[1-\lambda^{6}\right]W^{\prime}-\frac{2}{\rho^{2}}m^{2}[1-\lambda^{6}]W-\bar{h}\Bigg{[}\frac{4}{\rho}\bigg{[}4\lambda^{8}-4\lambda^{6}+\lambda^{2}-1\bigg{]}U^{\prime\prime}$
$\displaystyle-\frac{2}{\rho^{2}}\bigg{[}2m^{2}\lambda^{8}-3m^{2}\lambda^{6}-2\lambda^{8}+2m^{2}\lambda^{2}-m^{2}+2\lambda^{2}\bigg{]}U^{\prime}$
$\displaystyle-\frac{2}{\rho^{3}}\bigg{[}\lambda^{2}\big{[}2m^{2}\lambda^{6}-3m^{2}\lambda^{4}+2\lambda^{6}+m^{2}-2\big{]}\bigg{]}U\Bigg{]}$
$\displaystyle+\frac{2}{3}\bar{h}^{2}\lambda^{2}\Bigg{[}-2[1+\lambda^{6}]W^{iv}-\frac{1}{\rho}[3+7\lambda^{6}]W^{\prime\prime\prime}+\frac{m^{2}}{\rho^{2}}\bigg{[}3+\lambda^{6}\bigg{]}W^{\prime\prime}$
$\displaystyle+\frac{3m^{2}}{\rho^{3}}\bigg{[}2\lambda^{6}-1\bigg{]}W^{\prime}+\frac{m^{2}}{\rho^{4}}\bigg{[}\lambda^{6}m^{2}-3\lambda^{6}-m^{2}+3\bigg{]}W\Bigg{]}$
$\displaystyle-\frac{\beta\lambda^{2}}{2\bar{h}}W+\beta\lambda^{4}+\frac{\beta\lambda^{4}}{\rho}U-\frac{\bar{h}}{\rho}\beta[2\lambda^{6}-1]-\bar{h}\beta\lambda^{6}W^{\prime\prime}$
$\displaystyle-\frac{1}{3}h^{2}\lambda^{2}\beta\bigg{[}\frac{1}{\rho}[7+6\lambda^{6}]U^{\prime\prime}-\frac{1}{\rho^{2}}[9\lambda^{6}+2m^{2}-6]U^{\prime}+\frac{1}{\rho^{3}}[\lambda^{6}+3m^{2}-4]U\bigg{]}=0.$
(3.39)

####  3.2.2  Results and discussion

The compound matrix method is applied to solve the governing equations to
determine variation of the critical growth factor ($\lambda_{cr}$) with
respect to normalized thickness ($\bar{h}$).

Figure 7: Variation of the critical growth parameter $\lambda_{cr}$ with
normalised plate thickness $\bar{h}$ and $\beta=0.2$ for different
circumferential mode number ($m$).

The variation of critical buckling parameter ($\lambda_{cr}$) for different
circumferential mode numbers are given in Figure 7. The results shows that the
higher modes arise for lower value of $\lambda_{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.

Figure 8: Normalised bifurcation mode shapes for circular plate under combined
radial and circumferential growth (a) Normalised out-of-plane displacement $W$
for mode numbers $m=2,4,6,10$ at $\bar{h}=0.1$. (b) Normalised in-plane
displacement $U$ for mode number $m=2,4,6,10$ at $\bar{h}=0.1$.

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 of $\bar{h}=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 ($\lambda_{cr}$) is low in combined growth case for higher modes as
compared to the purely radial growth case.

##  4 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 $(\lambda_{cr})$ with normalised
thickness $(\bar{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.

## Acknowledgements

Prashant Saxena acknowledges the support of startup funds from the James Watt
School of Engineering at the University of Glasgow.

## References

  * Ben Amar and Goriely (2005) Ben Amar M. and Goriely A.  “Growth and instability in elastic tissues”.  _Journal of the Mechanics and Physics of Solids_ , 53(10):2284–2319 (2005)
  * Ben Amar et al. (2011) Ben Amar M., Goriely A., Müller M.M., and Cugliandolo L.  _New Trends in the Physics and Mechanics of Biological Systems: Lecture Notes of the Les Houches Summer School: Volume 92, July 2009_ , volume 92.  Oxford University Press (2011)
  * Cao et al. (2012) Cao Y., Jiang Y., Li B., and Feng X.  “Biomechanical modeling of surface wrinkling of soft tissues with growth-dependent mechanical properties”.  _Acta Mechanica Solida Sinica_ , 25(5):483–492 (2012)
  * Coman and Haughton (2006) Coman C. and Haughton D.  “Localized wrinkling instabilities in radially stretched annular thin films”.  _Acta Mechanica_ , 185(3-4):179–200 (2006)
  * Coman et al. (2015) Coman C.D., Matthews M.T., and Bassom A.P.  “Asymptotic phenomena in pressurized thin films”.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 471(2182):20150471 (2015)
  * Dai and Song (2014) Dai H.H. and Song Z.  “On a consistent finite-strain plate theory based on three-dimensional energy principle”.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 470(2171):20140494 (2014)
  * Deng et al. (2019) Deng X., Xu Y., and Clarke C.  “Wrinkling modelling of space membranes subject to solar radiation pressure”.  _Composites Part B: Engineering_ , 157:266–275 (2019)
  * Dervaux et al. (2009) Dervaux J., Ciarletta P., and Ben Amar M.  “Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the föppl–von kármán limit”.  _Journal of the Mechanics and Physics of Solids_ , 57(3):458–471 (2009)
  * Goriely (2017) Goriely A.  _The mathematics and mechanics of biological growth_ , volume 45.  Springer (2017)
  * Goriely and Ben Amar (2005) Goriely A. and Ben Amar M.  “Differential growth and instability in elastic shells”.  _Physical review letters_ , 94(19):198103 (2005)
  * Goriely and Ben Amar (2007) Goriely A. and Ben Amar M.  “On the definition and modeling of incremental, cumulative, and continuous growth laws in morphoelasticity”.  _Biomechanics and Modeling in Mechanobiology_ , 6(5):289–296 (2007)
  * Haughton and Orr (1997) Haughton D. and Orr A.  “On the eversion of compressible elastic cylinders”.  _International journal of solids and structures_ , 34(15):1893–1914 (1997)
  * Hoger (1986) Hoger A.  “On the determination of residual stress in an elastic body”.  _Journal of Elasticity_ , 16(3):303–324 (1986)
  * Kienzler (2002) Kienzler R.  “On consistent plate theories”.  _Archive of Applied Mechanics_ , 72(4-5):229–247 (2002)
  * Li et al. (2012) Li B., Cao Y.P., Feng X.Q., and Gao H.  “Mechanics of morphological instabilities and surface wrinkling in soft materials: a review”.  _Soft Matter_ , 8(21):5728–5745 (2012)
  * Li et al. (2010) Li B., Huang S.Q., and Feng X.Q.  “Buckling and postbuckling of a compressed thin film bonded on a soft elastic layer: a three-dimensional analysis”.  _Archive of Applied Mechanics_ , 80(2):175 (2010)
  * Lindsay and Rooney (1992) Lindsay K. and Rooney C.  “A note on compound matrices”.  _Journal of Computational Physics_ , 103(2):472–477 (1992)
  * Moulton and Goriely (2011) Moulton D. and Goriely A.  “Circumferential buckling instability of a growing cylindrical tube”.  _Journal of the Mechanics and Physics of Solids_ , 59(3):525–537 (2011)
  * Nassar et al. (2012) Nassar D., Letavernier E., Baud L., Aractingi S., and Khosrotehrani K.  “Calpain activity is essential in skin wound healing and contributes to scar formation”.  _PloS one_ , 7(5):e37084 (2012)
  * Ng and Reid (1979) Ng B. and Reid W.  “A numerical method for linear two-point boundary-value problems using compound matrices”.  _Journal of Computational Physics_ , 33(1):70–85 (1979)
  * Ng and Reid (1985) Ng B. and Reid W.  “The compound matrix method for ordinary differential systems”.  _Journal of Computational Physics_ , 58(2):209–228 (1985)
  * Papastavrou et al. (2013) Papastavrou A., Steinmann P., and Kuhl E.  “On the mechanics of continua with boundary energies and growing surfaces”.  _Journal of the Mechanics and Physics of Solids_ , 61(6):1446–1463 (2013)
  * Rodriguez et al. (1994) Rodriguez E.K., Hoger A., and McCulloch A.D.  “Stress-dependent finite growth in soft elastic tissues”.  _Journal of biomechanics_ , 27(4):455–467 (1994)
  * Rogers et al. (2010) Rogers J.A., Someya T., and Huang Y.  “Materials and mechanics for stretchable electronics”.  _science_ , 327(5973):1603–1607 (2010)
  * Swain and Gupta (2015) Swain D. and Gupta A.  “Interfacial growth during closure of a cutaneous wound: stress generation and wrinkle formation”.  _Soft matter_ , 11(32):6499–6508 (2015)
  * Swain and Gupta (2016) Swain D. and Gupta A.  “Mechanics of cutaneous wound rupture”.  _Journal of biomechanics_ , 49(15):3722–3730 (2016)
  * Tepole et al. (2011) Tepole A.B., Ploch C.J., Wong J., Gosain A.K., and Kuhl E.  “Growing skin: A computational model for skin expansion in reconstructive surgery”.  _Journal of the Mechanics and Physics of Solids_ , 59(10):2177–2190 (2011)
  * Vandiver and Goriely (2009) Vandiver R. and Goriely A.  “Differential growth and residual stress in cylindrical elastic structures”.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 367(1902):3607–3630 (2009)
  * Wang et al. (2009) Wang C., Du X., Tan H., and He X.  “A new computational method for wrinkling analysis of gossamer space structures”.  _International Journal of Solids and Structures_ , 46(6):1516–1526 (2009)
  * Wang et al. (2019a) Wang F.F., Steigmann D.J., and Dai H.H.  “On a uniformly-valid asymptotic plate theory”.  _International Journal of Non-Linear Mechanics_ , 112:117–125 (2019a)
  * Wang et al. (2016) Wang J., Song Z., and Dai H.H.  “On a consistent finite-strain plate theory for incompressible hyperelastic materials”.  _International Journal of Solids and Structures_ , 78:101–109 (2016)
  * Wang et al. (2018) Wang J., Steigmann D., Wang F.F., and Dai H.H.  “On a consistent finite-strain plate theory of growth”.  _Journal of the Mechanics and Physics of Solids_ , 111:184–214 (2018)
  * Wang et al. (2019b) Wang J., Wang Q., Dai H.H., Du P., and Chen D.  “Shape-programming of hyperelastic plates through differential growth: an analytical approach”.  _Soft matter_ , 15(11):2391–2399 (2019b)
  * Wang et al. (2013) Wang X., Ge J., Tredget E.E., and Wu Y.  “The mouse excisional wound splinting model, including applications for stem cell transplantation”.  _Nature protocols_ , 8(2):302–309 (2013)
  * Wei and Zhao (2014) Wei K. and Zhao Y.  “Fabrication of anisotropic and hierarchical undulations by benchtop surface wrinkling”.  _In_ “2014 IEEE 27th International Conference on Micro Electro Mechanical Systems (MEMS)”, pages 474–477. IEEE (2014)
  * Wu and Amar (2015) Wu M. and Amar M.B.  “Growth and remodelling for profound circular wounds in skin”.  _Biomechanics and modeling in mechanobiology_ , 14(2):357–370 (2015)
  * Wu and Ben Amar (2015) Wu M. and Ben Amar M.  “Modelling fibers in growing disks of soft tissues”.  _Mathematics and Mechanics of Solids_ , 20(6):663–679 (2015)

##  Appendix A Appendix: Expression for Piola stress

The average stress is given as

$\displaystyle\bar{\mathbf{P}}$
$\displaystyle=\frac{1}{2h}\int_{0}^{2h}\mathbf{P}dZ=\frac{1}{2h}\int_{0}^{2h}\left[\mathbf{P}^{(0)}+Z\mathbf{P}^{(1)}+\frac{Z^{2}}{2!}\mathbf{P}^{(2)}+\frac{Z^{3}}{3!}\mathbf{P}^{(3)}+O(Z^{4})\right]dZ,$
$\displaystyle\bar{\mathbf{P}}$
$\displaystyle=\mathbf{P}^{(0)}+h\mathbf{P}^{(1)}+\frac{2}{3}h^{2}\mathbf{P}^{(3)}+O(h^{3}).$
(A.1)

Subtracting the top and bottom surface traction condition using (2.17) and
(2.18)

$\displaystyle\left.\mathbf{P}\mathbf{k}\right|_{Z=2h}-\left.\mathbf{P}\mathbf{k}\right|_{Z=0}=$
$\displaystyle\mathbf{P}^{(0)}\mathbf{k}+2h\mathbf{P}^{(1)}\mathbf{k}+2h^{2}\mathbf{P}^{(2)}\mathbf{k}+\frac{4}{3}h^{3}\mathbf{P}^{(3)}\mathbf{k}-\mathbf{P}^{(0)}\mathbf{k}+O(h^{3})=\mathbf{q}^{+}+\mathbf{q}^{-},$
$\displaystyle=$
$\displaystyle\mathbf{P}^{(1)}\mathbf{k}+h\mathbf{P}^{(2)}\mathbf{k}+\frac{2}{3}h^{2}\mathbf{P}^{(3)}\mathbf{k}=\frac{\mathbf{q}^{+}+\mathbf{q}^{-}}{2h}+O(h^{3})=\bar{\mathbf{q}}.$
(A.2)

Using equilibrium equation

$\displaystyle\nabla\cdot\bar{\mathbf{P}}+\frac{\partial}{\partial
Z}\left[\bar{\mathbf{P}}\mathbf{k}\right]=\mathbf{0},$ (A.3)

where

$\displaystyle\frac{\partial}{\partial
Z}\left[\bar{\mathbf{P}}\mathbf{k}\right]$
$\displaystyle=\frac{\partial}{\partial
Z}\left[\mathbf{P}^{(0)}\mathbf{k}+h\mathbf{P}^{(1)}\mathbf{k}+\frac{2}{3}h^{2}\mathbf{P}^{(3)}\mathbf{k}+O(h^{3})\right]=\mathbf{P}^{(1)}\mathbf{k}+h\mathbf{P}^{(2)}\mathbf{k}+\frac{2}{3}h^{2}\mathbf{P}^{(3)}\mathbf{k}+O(h^{3}).$

Now, the Piola stress $(\mathbf{P})$ is described in terms of strain energy
function $\left(\phi(\mathbf{F,G})=J_{G}\phi_{0}(\mathbf{A})\right)$ as

$\displaystyle\mathbf{P}$
$\displaystyle=\frac{\partial\phi(\mathbf{F,G})}{\partial\mathbf{F}}-p\frac{\partial{L}(\mathbf{F,G})}{\partial\mathbf{F}}=J_{G}\bigg{[}\frac{\partial\phi_{0}(\mathbf{A})}{\partial\mathbf{F}}-p\frac{\partial{L_{0}}(\mathbf{A})}{\partial\mathbf{F}}\bigg{]}$
(A.4)
$\displaystyle=J_{G}\bigg{[}\frac{\partial\phi_{0}}{\partial\mathbf{A}}~{}\frac{\partial\mathbf{A}}{\partial\mathbf{F}}-p\frac{\partial
L(\mathbf{A})}{\partial\mathbf{A}}\frac{\partial\mathbf{A}}{\partial\mathbf{F}}\bigg{]}=J_{G}\bigg{[}\frac{\partial\phi_{0}}{\partial{A}_{ij}}\frac{\partial[{FG}^{-1}]_{ij}}{\partial{F}_{lm}}-\frac{\partial\det(\mathbf{A})}{\partial{A}_{ij}}\frac{\partial[{FG}^{-1}]_{ij}}{\partial{F}_{lm}}\bigg{]}$
$\displaystyle=J_{G}\bigg{[}\frac{\partial\phi_{0}}{\partial{A}_{ij}}~{}\frac{\partial{F}_{ik}}{\partial{F}_{lm}}{G}^{-1}_{kj}-p\det(\mathbf{A})A_{ij}^{-T}~{}\frac{\partial{F}_{ik}}{\partial{F}_{lm}}{G}^{-1}_{kj}\bigg{]},$
$\displaystyle=J_{G}\bigg{[}\frac{\partial\phi_{0}}{\partial{A}_{ij}}~{}\delta_{il}\delta_{km}\mathbf{G}^{-1}_{kj}-pA_{ij}^{-T}~{}\frac{\partial{F}_{ik}}{\partial{F}_{lm}}{G}^{-1}_{kj}\bigg{]}=J_{G}\bigg{[}\frac{\partial\phi_{0}}{\partial{A}_{ij}}~{}\delta_{il}{G}^{-1}_{mj}-pA_{ij}^{-T}~{}\delta_{il}{G}^{-1}_{mj}\bigg{]}$
$\displaystyle=\bigg{[}J_{G}\left[\frac{\partial\phi_{0}(\mathbf{A})}{\partial{A}}\right]_{lj}~{}{G}^{-1}_{mj}-pA_{lj}^{-T}{G}^{-1}_{mj}\bigg{]}=J_{G}\bigg{[}\frac{\partial\phi_{0}(\mathbf{A})}{\partial\mathbf{A}}~{}\mathbf{G}^{-T}-p\mathbf{A}^{-T}\mathbf{G}^{-T}\bigg{]},$
$\displaystyle=J_{G}\left[\frac{\partial\phi_{0}(\mathbf{A})}{\partial\mathbf{A}}-p\mathbf{A}^{-T}\right]\mathbf{G}^{-T}.$

The series expansion for Piola Kirchhoff stress ($\mathbf{P}$) is given as

$\displaystyle\mathbf{P}(\mathbf{x},p)$
$\displaystyle=\mathbf{P}^{(0)}(\mathbf{x},p)+Z\mathbf{P}^{(1)}(\mathbf{x},p)+\frac{Z^{2}}{2}\mathbf{P}^{(2)}(\mathbf{x},p)+\frac{Z^{3}}{3!}\mathbf{P}^{(3)}(\mathbf{x},p)+O(Z^{4}).$
(A.5)

Using (2.10c) and (A.5) we have

$\displaystyle\mathbf{P}(\mathbf{x},p)=J_{G}\left[\frac{\partial\phi_{0}(\mathbf{A})}{\partial\mathbf{A}}-p\frac{\partial
L_{0}(\mathbf{A})}{\partial\mathbf{A}}\right]\left[\bar{\mathbf{G}}^{(0)}(\mathbf{\zeta})+Z\bar{\mathbf{G}}^{(1)}(\mathbf{\zeta})+\frac{Z^{2}}{2}\bar{\mathbf{G}}^{(2)}(\mathbf{\zeta})+\frac{Z^{3}}{3!}\bar{\mathbf{G}}^{(3)}(\mathbf{\zeta})+O(Z^{4})\right].$
(A.6)

Substituting expansion of
$\phi_{0}(\mathbf{A})~{}\text{and}~{}L_{0}(\mathbf{A})$ about
$\mathbf{A}^{(0)}$ we obtain

$\displaystyle\mathbf{P}(\mathbf{x},p)$
$\displaystyle=J_{G}\bigg{[}\frac{\partial}{\partial\mathbf{A}}\left[\phi_{0}(\mathbf{A})+\frac{\partial\phi_{0}}{\partial\mathbf{A}}[\mathbf{A}-\mathbf{A^{(0)}}]+\frac{1}{2!}\frac{\partial^{2}\phi_{0}}{\partial\mathbf{A}\partial\mathbf{A}}[\mathbf{A}-\mathbf{A^{(0)}},~{}\mathbf{A}-\mathbf{A^{(0)}}]\right.$
(A.7)
$\displaystyle\left.+\frac{1}{3!}\frac{\partial^{3}\phi_{0}}{\partial\mathbf{A}\partial\mathbf{A}\partial\mathbf{A}}[\mathbf{A}-\mathbf{A^{(0)}},~{}\mathbf{A}-\mathbf{A^{(0)}},~{}\mathbf{A}-\mathbf{A^{(0)}}\right]$
$\displaystyle-\left[p^{(0)}+Zp^{(1)}+\frac{Z^{2}}{2}p^{(2)}+\frac{Z^{3}}{3!}p^{(3)}+\frac{Z^{4}}{4!}p^{(4)}+O(Z^{4})\right]\times$
$\displaystyle\hskip
10.84006pt\bigg{[}\frac{\partial}{\partial\mathbf{A}}\left[L_{0}(\mathbf{A})+\frac{\partial
L_{0}}{\partial\mathbf{A}}[\mathbf{A}-\mathbf{A^{(0)}}]+\frac{1}{2!}\frac{\partial^{2}L_{0}}{\partial\mathbf{A}\partial\mathbf{A}}[\mathbf{A}-\mathbf{A^{(0)}},~{}\mathbf{A}-\mathbf{A^{(0)}}]\right.$
$\displaystyle\left.+\frac{1}{3!}\frac{\partial^{3}L_{0}}{\partial\mathbf{A}\partial\mathbf{A}\partial\mathbf{A}}[\mathbf{A}-\mathbf{A^{(0)}},~{}\mathbf{A}-\mathbf{A^{(0)}},~{}\mathbf{A}-\mathbf{A^{(0)}}]\right]\bigg{]}$
$\displaystyle\hskip
7.22743pt\left[\bar{\mathbf{G}}^{(0)}(\mathbf{\zeta})+Z\bar{\mathbf{G}}^{(1)}(\mathbf{\zeta})+\frac{Z^{2}}{2}\bar{\mathbf{G}}^{(2)}(\mathbf{\zeta})+\frac{Z^{3}}{3!}\bar{\mathbf{G}}^{(3)}(\mathbf{\zeta})+O(Z^{4})\right].$

Substituting expression for $[\mathbf{A}-\mathbf{A^{(0)}}]$ using (2.10b) and
compare the terms

$\displaystyle\mathbf{P}^{(0)}(\mathbf{x},p)$
$\displaystyle=J_{G}\left[\boldsymbol{\mathcal{A}}^{(0)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(0)}\right]\bar{\mathbf{G}}^{(0)},$
(A.8) $\displaystyle\mathbf{P}^{(1)}(\mathbf{x},p)$
$\displaystyle=J_{G}\left[\left[\boldsymbol{\mathcal{A}}^{(1)}[\mathbf{A}^{(1)}]-p^{(0)}\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(1)}]-p^{(1)}\boldsymbol{\mathcal{L}}^{(0)}\right]\bar{\mathbf{G}}^{(0)}+\left[\boldsymbol{\mathcal{A}}^{(0)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(0)}\right]\bar{\mathbf{G}}^{(1)}\right],$
$\displaystyle\mathbf{P}^{(2)}(\mathbf{x},p)$
$\displaystyle=J_{G}\bigg{[}\bigg{[}\boldsymbol{\mathcal{A}}^{(1)}[\mathbf{A}^{(2)}]+\boldsymbol{\mathcal{A}}^{(2)}[\mathbf{A}^{(1)},~{}\mathbf{A}^{(1)}]-2p^{(1)}\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(1)}]-p^{(0)}\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(2)}]$
$\displaystyle\hskip
14.45377pt-p^{(0)}\boldsymbol{\mathcal{L}}^{(2)}[\mathbf{A}^{(1)},\mathbf{A}^{(1)}]-p^{(2)}\boldsymbol{\mathcal{L}}^{(0)}\bigg{]}\bar{\mathbf{G}}^{(0)}$
$\displaystyle\hskip
17.34518pt+\bigg{[}2\boldsymbol{\mathcal{A}}^{(1)}[\mathbf{A}^{(1)}]-2p^{(0)}\boldsymbol{\mathcal{L}}^{(1)}[\mathbf{A}^{(1)}]-2p^{(1)}\boldsymbol{\mathcal{L}}^{(0)}\bigg{]}\bar{\mathbf{G}}^{(1)}+\left[\boldsymbol{\mathcal{A}}^{(0)}-p^{(0)}\boldsymbol{\mathcal{L}}^{(0)}\right]\bar{\mathbf{G}}^{(2)}\bigg{]}.$

where
$\displaystyle\boldsymbol{\mathcal{A}}^{{i}}(\mathbf{A}^{(0)})=\left.\frac{\partial^{i+1}\phi_{0}(\mathbf{A})}{\partial\mathbf{A}^{i+1}}\right|_{\mathbf{A}=\mathbf{A}^{(0)}}$
and
$\displaystyle\boldsymbol{\mathcal{L}}^{{i}}(\mathbf{A}^{(0)})=\left.\frac{\partial
R_{0}(\mathbf{A})}{\partial\mathbf{A}^{i+1}}\right|_{\mathbf{A}=\mathbf{A}^{(0)}}$

##  Appendix B Appendix: Expression for $\mathcal{L}^{i}$

$\displaystyle\boldsymbol{\mathcal{L}}^{(0)}(\mathbf{A})$
$\displaystyle=\frac{\partial\det(\mathbf{A})}{\partial\mathbf{A}}=\det(\mathbf{A})\mathbf{A}^{-T}=\det(\mathbf{A})A_{ab}^{-T}\mathbf{e}_{a}\otimes\mathbf{e}_{b},$
$\displaystyle\boldsymbol{\mathcal{L}}^{(1)}(\mathbf{A})$
$\displaystyle=\frac{\partial^{2}\det(\mathbf{A})}{\partial\mathbf{A}~{}\partial\mathbf{A}}=\text{det}(\mathbf{A})\bigg{[}\mathbf{A}^{-T}\otimes\mathbf{A}^{-T}\bigg{]}+\text{det}(\mathbf{A})\bigg{[}\mathbb{T}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}[-\mathbf{A}^{-1}\boxtimes\mathbf{A}^{-T}]}\bigg{]},$
$\displaystyle=\det(\mathbf{A})\bigg{[}A_{ab}^{-T}A_{cd}^{-T}-A_{bc}^{-1}A_{ad}^{-T}\bigg{]}\mathbf{e}_{a}\otimes\mathbf{e}_{b}\otimes\mathbf{e}_{c}\otimes\mathbf{e}_{d},$
$\displaystyle\boldsymbol{\mathcal{L}}^{(2)}(\mathbf{A})$
$\displaystyle=\frac{\partial^{3}\det(\mathbf{A})}{\partial\mathbf{A}~{}\partial\mathbf{A}~{}\partial\mathbf{A}}$
$\displaystyle=\text{det}(\mathbf{A})\bigg{[}A_{ab}^{-T}A_{cd}^{-T}A_{ef}^{-T}-A_{cb}^{-T}A_{ad}^{-T}A_{ef}^{-T}-A^{-T}_{eb}A^{-T}_{af}A^{-T}_{cd}$
$\displaystyle\hskip
14.45377pt-A^{-T}_{ab}A^{-T}_{ed}A^{-T}_{cf}+A^{-T}_{eb}A^{-T}_{cf}A^{-T}_{ad}+A^{-T}_{cb}A^{-T}_{ed}A^{-T}_{af}\bigg{]}\mathbf{e}_{a}\otimes\mathbf{e}_{b}\otimes\mathbf{e}_{c}\otimes\mathbf{e}_{d}\otimes\mathbf{e}_{e}\otimes\mathbf{e}_{f}.$

where,

$\displaystyle[\mathbf{A}\otimes\mathbf{B}]_{abcd}$
$\displaystyle=\big{[}A_{ab}A_{cd}\big{]}\mathbf{e}_{a}\otimes\mathbf{e}_{b}\otimes\mathbf{e}_{c}\otimes\mathbf{e}_{d},$
$\displaystyle[\mathbf{A}\boxtimes\mathbf{B}]_{abcd}$
$\displaystyle=\big{[}A_{ac}A_{bd}\big{]}\mathbf{e}_{a}\otimes\mathbf{e}_{b}\otimes\mathbf{e}_{c}\otimes\mathbf{e}_{d},$
$\displaystyle[{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{T}[\mathbf{A}\boxtimes\mathbf{B}]}]_{abcd}$
$\displaystyle=[A_{bc}A_{ad}]\mathbf{e}_{a}\otimes\mathbf{e}_{b}\otimes\mathbf{e}_{c}\otimes\mathbf{e}_{d}.$

and $\mathbb{T}_{ijkl}=\delta_{il}\delta_{jk}$.
Note: In this current work, we assume incompressible material thus,
$\det(\mathbf{A})=1$

##  Appendix C Appendix: Expression for $\mathbf{x}^{(1)}$ and $p^{(0)}$

For an incompressible neo-Hookean material elastic strain energy function is
given as

$\displaystyle\phi_{0}(\mathbf{A})$ $\displaystyle=C_{0}[I_{1}-3],$
$\displaystyle\phi_{0}(\mathbf{A})$
$\displaystyle=C_{0}[\text{tr}(\mathbf{A}^{T}\mathbf{A})-3],$ (C.1)

where $I_{1}=\text{tr}(\mathbf{C})=\text{tr}(\mathbf{A}^{T}\mathbf{A})$ as
elastic strain energy depends only on elastic deformation tensor
$(\mathbf{A})$ and $tr(\cdot)$ is trace of a tensor ($\cdot$). Using (A.4),
the Piola stress $(\mathbf{P})$ for incompressible neo-Hookean material is

$\displaystyle\mathbf{P}=J_{G}\bigg{[}2C_{0}[\mathbf{A}]-p\mathbf{A}^{-T}\bigg{]}\mathbf{G}^{-T},$
(C.2)

where

$\displaystyle\frac{\partial\phi_{0}(\mathbf{A})}{\partial\mathbf{A}}$
$\displaystyle=C_{0}\left[\frac{\partial}{\partial\mathbf{A}}\bigg{[}\text{tr}(\mathbf{A}^{T}\mathbf{A})-3\bigg{]}\right]=C_{0}\left[\frac{\partial}{\partial{A}_{lm}}({A}^{T}{A})_{ii})\right]=C_{0}\left[\frac{\partial}{\partial{A}_{lm}}({A}_{ij}A_{ij})\right],$
$\displaystyle=C_{0}\bigg{[}\delta_{il}\delta_{jm}{A}_{ij}+{A}_{ij}\delta_{il}\delta_{jm}\bigg{]}=2C_{0}\bigg{[}{A}_{ij}\delta_{il}\delta_{jm}\bigg{]}=2C_{0}~{}{A}_{lm}=2C_{0}~{}\mathbf{A}.$

Now, the first term $\mathbf{P}^{(0)}$ in (A.5) for neo-Hookean material is
rewritten as

$\displaystyle\mathbf{P}^{(0)}$
$\displaystyle=J_{G}\left[2C_{0}\mathbf{A}^{(0)}-p{\mathbf{A}^{(0)}}^{-T}\right]\bar{\mathbf{G}}^{(0)}.$
(C.3)

By using bottom traction condition,
$\mathbf{P}^{(0)}\mathbf{k}=-\mathbf{q}^{-}$ and substituting expression for
$\mathbf{A}^{(0)}$ we have

$\displaystyle\mathbf{P}^{(0)}\mathbf{k}$
$\displaystyle=2C_{0}J_{G}\mathbf{F}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}}\bar{\mathbf{G}}^{(0)}\mathbf{k}-J_{G}~{}p^{(0)}\mathbf{F}^{{(0)}^{-T}}\bar{\mathbf{G}}^{{(0)}^{-1}}\bar{\mathbf{G}}^{(0)}\mathbf{k},$
$\displaystyle-\mathbf{q}^{-}$
$\displaystyle=2C_{0}\bigg{[}\nabla\mathbf{x}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}}+\mathbf{x}^{(1)}\otimes\bar{\mathbf{G}}^{(0)}\mathbf{k}\widehat{\mathbf{G}}^{(0)}\mathbf{k}\bigg{]}-J^{(0)}p^{(0)}\mathbf{F}^{{(0)}^{-T}}\mathbf{k},$
$\displaystyle=2C_{0}\bigg{[}\nabla\mathbf{x}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}}(\widehat{\mathbf{G}}^{(0)})\mathbf{k}+\left[\bar{\mathbf{G}}^{(0)}\mathbf{k}\cdot\widehat{\mathbf{G}}^{(0)}\mathbf{k}\right]\mathbf{x}^{(1)}\bigg{]}-\left[J^{(0)}p^{(0)}\frac{\text{Cofac}(\mathbf{F}^{(0)})}{\text{det}(\mathbf{A}^{(0)}){\text{det}(\mathbf{G}^{(0)})}}\right]\mathbf{k},$
$\displaystyle=2C_{0}\nabla\mathbf{x}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}}\widehat{\mathbf{G}}^{(0)}\mathbf{k}+2C_{0}J^{(0)}\left|\bar{\mathbf{G}}^{(0)}\mathbf{k}\right|^{2}\mathbf{x}^{(1)}-\left[J^{(0)}p^{(0)}\frac{\text{Cofac}(\mathbf{F}^{(0)})}{(1)J^{(0)}}\right]\mathbf{k}.$
(C.4)

where $J^{(0)}=\left.J_{G}\right|_{Z=0}$. Traction at bottom surface is given
as

$\displaystyle-\mathbf{q}^{-}$
$\displaystyle=2C_{0}\nabla\mathbf{x}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}}\widehat{\mathbf{G}}^{(0)}\mathbf{k}+2C_{0}J_{G}\left|\bar{\mathbf{G}}^{(0)}\mathbf{k}\right|^{2}\mathbf{x}^{(1)}-p^{(0)}\text{Cofac}(\mathbf{F}^{(0)})\mathbf{k},$
(C.5)

where
$\text{Cofac}(\mathbf{F}^{(0)})\mathbf{k}=\text{Cofac}(F^{(0)})_{ij}k_{j}\mathbf{e}_{i}$.
In this work, we use $\nabla\mathbf{x}^{{(0)}^{*}}$ in place of
$\text{Cofac}(\mathbf{F}^{(0)})\mathbf{k}$ and for deformation gradient in
polar coordinate system it is given by

$\displaystyle\nabla\mathbf{x}^{{(0)}^{*}}=\displaystyle\frac{r^{(0)}}{R}\left[\frac{\partial\theta}{\partial
R}\displaystyle\frac{\partial
z^{(0)}}{\partial\Theta}-\displaystyle\frac{\partial\theta}{\partial\Theta}\displaystyle\frac{\partial
z^{(0)}}{\partial
R}\right]\mathbf{e}_{1}+\displaystyle\frac{1}{R}\left[\displaystyle\frac{\partial
r^{(0)}}{\partial\Theta}\displaystyle\frac{\partial z^{(0)}}{\partial
R}-\displaystyle\frac{\partial r^{(0)}}{\partial R}\displaystyle\frac{\partial
z^{(0)}}{\partial\Theta}\right]\mathbf{e}_{2}$
$\displaystyle+\displaystyle\frac{r^{(0)}}{R}\left[\displaystyle\frac{\partial
r^{(0)}}{\partial
R}\displaystyle\frac{\partial\theta}{\partial\Theta}-\displaystyle\frac{\partial\theta}{\partial
R}\displaystyle\frac{\partial r^{(0)}}{\partial\Theta}\right]\mathbf{k}.$
(C.6)

We know from (2.3)

$\displaystyle\text{det}(\mathbf{A}^{(0)})=1=\text{det}(\mathbf{F}^{(0)})\text{det}(\bar{\mathbf{G}}^{{(0)}^{T}}),$
$\displaystyle\text{det}(\mathbf{F}^{(0)})=\text{det}(\bar{\mathbf{G}}^{{(0)}^{-T}}).$
(C.7)

As the definition of determinant

$\displaystyle\text{det}({\mathbf{F}}^{(0)})$
$\displaystyle=\bigg{[}r^{(1)}\mathbf{e}_{1}+r^{(0)}\theta^{(1)}\mathbf{e}_{2}+z^{(1)}\mathbf{e}_{3}\bigg{]}\cdot\text{Cofac}(\mathbf{F}^{(0)})\mathbf{k}=\mathbf{x}^{(1)}\cdot\nabla\mathbf{x}^{{(0)}^{*}}.$
(C.8)

By combining (C.7) and (C.8) we have

$\displaystyle\mathbf{x}^{(1)}\cdot\nabla\mathbf{x}^{{(0)}^{*}}=\det\left(\bar{\mathbf{G}}^{{(0)}^{-T}}\right).$
(C.9)

Using (C.5) we obtain the explicit expression for $\mathbf{x}^{(1)}$

$\displaystyle\mathbf{x}^{(1)}$
$\displaystyle=\frac{-\mathbf{q}^{-}-2C_{0}\nabla\mathbf{x}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}}\widehat{\mathbf{G}}^{(0)}\mathbf{k}+p^{(0)}\nabla\mathbf{x}^{{(0)}^{*}}}{2C_{0}J_{G}\left|\bar{\mathbf{G}}^{(0)}\mathbf{k}\right|^{2}}.$
(C.10)

To obtain the explicit expression for $p^{(0)}$ we substitute (C.10) into
(C.9) which yields

$\displaystyle
p^{(0)}=\frac{2C_{0}J_{G}\left|\bar{\mathbf{G}}^{(0)}\mathbf{k}\right|^{2}}{\text{det}\bar{\mathbf{G}}^{{(0)}^{T}}\left|\nabla\mathbf{x}^{{(0)}^{*}}\right|^{2}}+\left[\mathbf{q}^{-}+2C_{0}\nabla\mathbf{x}^{(0)}\bar{\mathbf{G}}^{{(0)}^{T}}\widehat{\mathbf{G}}^{(0)}\mathbf{k}\right]\cdot\frac{\nabla\mathbf{x}^{{(0)}^{*}}}{{\left|\nabla\mathbf{x}^{{(0)}^{*}}\right|^{2}}}.$
(C.11)

##  Appendix D Appendix: General Description on Compound Matrix Method

Consider a two-point boundary value problem expressed in first order ordinary
differential equations

$\displaystyle\frac{d\mathbf{Y}}{dX}=\boldsymbol{\mathscr{A}}(\lambda,x)\mathbf{Y},\hskip
57.81621ptx\in(a,b)$ (D.1)

subjected to boundary condition

$\displaystyle\mathbf{BY}$ $\displaystyle=\mathbf{0},\hskip 72.26999ptx=a,$
$\displaystyle\mathbf{CY}$ $\displaystyle=\mathbf{0},\hskip 72.26999ptx=b,$
(D.2)

where $\lambda$ is the eigenvalue or critical buckling parameter, $\mathbf{Y}$
is $1\times 2n$ vector, $\boldsymbol{\mathscr{A}}$ is $2n\times 2n$ matrix and
$\mathbf{B}$ and $\mathbf{C}$ both are $n\times 2n$ full rank matrices i.e.,
$n$ boundary conditions are given at $x=a,b$. Assume

$\displaystyle\bigg{\\{}\mathbf{y}^{(1)}(\lambda,x),\mathbf{y}^{(2)}(\lambda,x),...,\mathbf{y}^{(n)}(\lambda,x)\bigg{\\}},$
(D.3)

is set of $n$ linearly independent solution (D.1) which satisfy the boundary
condition at $x=0$ and the general solution of (D.1) can be written as the
linear combination of its independent solution

$\displaystyle\mathbf{y}(\lambda,x)=\sum_{j=1}^{n}k_{j}\mathbf{y}^{j},$ (D.4)

where $k_{1},k_{2},...k_{n}$ are the constants. Solution matrix $\mathbf{M}$
to be $2n\times n$ is define whose $j$th column is $\mathbf{y}^{(j)}$ as
$\mathbf{M}=[\mathbf{y}^{(1)},\mathbf{y}^{(2)},...\mathbf{y}^{(n)}]$, then
(D.1) in terms of $\mathbf{M}$ is given

$\displaystyle\frac{d\mathbf{M}}{dx}=[\boldsymbol{\mathscr{A}}\mathbf{y}^{(1)},\boldsymbol{\mathscr{A}}\mathbf{y}^{(2)},...,\boldsymbol{\mathscr{A}}\mathbf{y}^{(n)}]=\boldsymbol{\mathscr{A}}\mathbf{M}.$
(D.5)

The compound variables are defined as minors of $\mathbf{M}$ and denoted as
$\Phi_{1},\Phi_{2},...$ and those are $(^{2n}C_{n})$ in numbers. For an
instance consider fourth order ODE ($n=2$), then the solution matrix is

$\displaystyle\mathbf{M}=\begin{bmatrix}y_{1}^{(1)}&y_{1}^{(2)}\\\
y_{2}^{(1)}&y_{2}^{(2)}\\\ y_{3}^{(1)}&y_{3}^{(2)}\\\
y_{4}^{(1)}&y_{4}^{(2)}\end{bmatrix},$ (D.6)

and 6 minors of $\mathbf{M}$

$\displaystyle\Phi_{1}$
$\displaystyle=(1,2)=\begin{vmatrix}y_{1}^{(1)}&y_{1}^{(2)}\\\
y_{2}^{(1)}&y_{2}^{(2)}\end{vmatrix},\hskip
36.135pt\Phi_{2}=(1,3)=\begin{vmatrix}y_{1}^{(1)}&y_{1}^{(2)}\\\
y_{3}^{(1)}&y_{3}^{(2)}\end{vmatrix},$ $\displaystyle\Phi_{3}$
$\displaystyle=(1,4),~{}~{}\Phi_{4}=(2,3),~{}~{}\Phi_{5}=(2,4),~{}~{}\Phi_{6}=(3,4).$
(D.7)

Using (D.1), the system of first order differential equation in terms of
compound variable is

$\displaystyle\Phi_{1}^{\prime}=\begin{vmatrix}y_{1}^{(1)}&y_{1}^{(2)}\\\
y_{2}^{(1)}&y_{2}^{(2)}\end{vmatrix}^{\prime}$
$\displaystyle=\begin{vmatrix}{y_{1}^{(1)}}^{\prime}&{y_{1}^{(2)}}^{\prime}\\\
{y_{2}^{(1)}}&{y_{2}^{(2)}}\end{vmatrix}+\begin{vmatrix}{y_{1}^{(1)}}&{y_{1}^{(2)}}\\\
{y_{2}^{(1)}}^{\prime}&{y_{2}^{(2)}}^{\prime}\end{vmatrix},$
$\displaystyle=\begin{vmatrix}\sum_{j=1}^{4}\mathscr{A}_{1j}y_{j}^{(1)}&\sum_{j=1}^{4}\mathscr{A}_{1j}y_{j}^{(2)}\\\
{y_{2}^{(1)}}&{y_{2}^{(2)}}\end{vmatrix}+\begin{vmatrix}y_{1}^{(1)}&y_{1}^{(2)}\\\
\sum_{j=1}^{4}\mathscr{A}_{2j}y_{j}^{(1)}&\sum_{j=1}^{4}\mathscr{A}_{2j}y_{j}^{(2)}\end{vmatrix},$
$\displaystyle=\mathscr{A}_{11}\Phi_{1}-\mathscr{A}_{13}\Phi_{4}-\mathscr{A}_{14}\Phi_{5}+\mathscr{A}_{22}\Phi_{1}+\mathscr{A}_{23}\Phi_{2}+\mathscr{A}_{24}\Phi_{3}.$
(D.8)

The system is now converted into $(^{2n}C_{n})$ ordinary differential
equations which is in the form of

$\displaystyle\mathbf{\Phi^{\prime}}=\boldsymbol{\mathscr{A}}^{*}(\lambda,x)\mathbf{\Phi},\hskip
72.26999ptx\in(a,b).$ (D.9)

Subjected to initial condition at $x=0$,

$\displaystyle\mathbf{\Phi}(a)=[\Phi_{1},~{}\Phi_{2},~{}\Phi_{3},~{}\Phi_{4},~{}\Phi_{5},~{}\Phi_{6}]$
(D.10)

Now, if we consider sixth order ODE system ($n=3$) then the solution matrix is
defined as

$\displaystyle\mathbf{M}=\begin{bmatrix}y_{1}^{(1)}&y_{1}^{(2)}&y_{1}^{(3)}\\\
y_{2}^{(1)}&y_{2}^{(2)}&y_{2}^{(3)}\\\ y_{3}^{(1)}&y_{3}^{(2)}&y_{3}^{(3)}\\\
y_{4}^{(1)}&y_{4}^{(2)}&y_{4}^{(3)}\\\ y_{5}^{(1)}&y_{5}^{(2)}&y_{5}^{(3)}\\\
y_{6}^{(1)}&y_{6}^{(2)}&y_{6}^{(3)}\end{bmatrix},$ (D.11)

where 20 minors of $\mathbf{M}$ are

$\displaystyle\Phi_{1}=(1,2,3),~{}~{}~{}\Phi_{2}=(1,2,4),~{}~{}~{}\Phi_{3}=(1,2,5),~{}~{}~{}\Phi_{4}=(1,2,6),~{}~{}~{}\Phi_{5}=(1,3,4),$
$\displaystyle\Phi_{6}=(1,3,5),~{}~{}~{}\Phi_{7}=(1,3,6),~{}~{}~{}\Phi_{8}=(1,4,5),~{}~{}~{}\Phi_{9}=(1,4,6),~{}~{}~{}\Phi_{10}=(1,5,6),$
$\displaystyle\Phi_{11}=(2,3,4),~{}~{}\Phi_{12}=(2,3,5),~{}~{}\Phi_{13}=(2,3,6),~{}~{}\Phi_{14}=(2,4,5),~{}~{}\Phi_{15}=(2,4,6),$
$\displaystyle\Phi_{16}=(2,5,6),~{}~{}\Phi_{17}=(3,4,5),~{}~{}\Phi_{18}=(3,4,6),~{}~{}\Phi_{19}=(3,5,6),~{}~{}\Phi_{20}=(4,5,6).$

The system of equation (D.9) is now numerically integrated using initial
condition (D.10) which produces the solution $y^{(j)}$ at $x=b$

$\displaystyle\mathbf{Cy}=\mathbf{C}\sum_{j=1}^{n}k_{j}\mathbf{y}^{(j)}(b)=\mathbf{CMk}=\mathbf{0}.$
(D.12)

For existence of non-trivial solution of differential equation (D.9)

$\displaystyle\det(\mathbf{CM})=0.$ (D.13)

##  Appendix E Appendix: Expression for unknown variable for radial and
circumferential growth

Using Eq. (C.11) we obtain the expression for $p^{(0)}$ as

$\displaystyle
p^{(0)}=\frac{2C_{0}\lambda^{4}}{\left|\nabla{\mathbf{x}^{(0)}}^{**}\right|^{2}},$
(E.1)

where $\nabla{\mathbf{x}^{(0)}}^{**}=\nabla{x}_{11}\mathbf{e}_{1}+\nabla
x_{22}\mathbf{e}_{2}+\nabla x_{33}\mathbf{e}_{3}$ and

$\displaystyle\nabla{x}_{11}=\frac{r^{(0)}}{R}\left[\frac{\partial\theta}{\partial
R}\frac{\partial
z^{(0)}}{\partial\Theta}-\frac{\partial\theta}{\partial\Theta}\frac{\partial
z^{(0)}}{\partial R}\right],\nabla{x}_{22}=\frac{1}{R}\left[\frac{\partial
r^{(0)}}{\partial\Theta}\frac{\partial z^{(0)}}{\partial R}-\frac{\partial
r^{(0)}}{\partial R}\frac{\partial z^{(0)}}{\partial\Theta}\right],$
$\displaystyle\nabla{x}_{33}=\frac{r^{(0)}}{R}\left[\frac{\partial
r^{(0)}}{\partial
R}\frac{\partial\theta}{\partial\Theta}-\frac{\partial\theta}{\partial
R}\frac{\partial r^{(0)}}{\partial\Theta}\right].$

On substituting the $p^{(0)}$ in (C.10) we obtain

$\displaystyle r^{(1)}=\displaystyle\frac{p^{(0)}\nabla
x_{11}}{2C_{0}\lambda^{2}},\quad\theta^{(1)}=\frac{p^{(0)}\nabla
x_{22}}{2C_{0}\lambda^{2}r^{(0)}},\quad z^{(1)}=\frac{p^{(0)}\nabla
x_{33}}{2C_{0}\lambda^{2}}.$ (E.2)


* Eilers and Marx (1996) authorEilers, P.H.C, and authorMarx, B.D. (year1996) titleFlexible smoothing with B-splines and penalties, journalStatist. Sci. volume11 pages89–102.
  * Eubank (1999) authorEubank, R.L. (year1999) titleNonparametric Regression and Spline Smoothing, 2nd ed., publisherCRC Press, New York.
  * Fan and Li (2001) authorFan, J., and authorLi, R. (year2001) titleVariable Selection via Nonconcave Penalized Likelihood and its Oracle Properties, journalJ. Amer. Statist. Assoc. volume96 pages1348–1360.
  * Ghement et al. (2008) authorGhement, I.R., authorRuiz, M., and authorZamar, R. (year2008) titleRobust estimation of error scale in nonparametric regression models, journalJ. Statist. Plann. and Inference volume138 pages3200–3216.
  * Ghosh and Basu (2013) authorGhosh, A., and authorBasu, A. (year2013) titleRobust estimation for independent non-homogeneous observation using density power divergence with applications to linear regression, journalElectron. J. Stat. volume7 pages2420–2456.
  * Ghosh and Basu (2015) authorGhosh, A., and authorBasu, A. (year2015) titleRobust estimation for non-homogeneous data and the selection of the optimal tuning parameter: the density power divergence approach, journalJ. Appl. Stat. volume42 pages2056-–2072.
  * Ghosh and Basu (2016) authorGhosh, A., and authorBasu, A. (year2016) titleRobust estimation in generalized linear models: the density power divergence approach, journalTEST volume25 pages269–290.
  * Green and Silverman (1994) authorGreen, P.J., and authorSilverman, B.W. (year1994) titleNonparametric Regression and Generalized Linear Models, publisherChapman & Hall, London.
  * Gu (2013) authorGu, C. (year2013) titleSmoothing Spline ANOVA Models, 2nd ed., publisherSpringer, New York.
  * Hampel et al. (2011) authorHampel, F.R., authorRonchetti, E.M., authorRousseeuw, P.J., and authorStahel, W.A. (year2011) titleRobust Statistics: The Approach Based on Influence Functions, publisherWiley, New York.
  * Hastie and Tibshirani (1990) authorHastie, T.J., and authorTibshirani, R.J. (year1990) titleGeneralized Additive Models, publisherChapman & Hall, Suffolk.
  * Hastie et al. (2009) authorHastie, T.J., authorTibshirani, R.J., and authorFriedman, J. (year2009) titleThe Elements of Statistical Learning: Data Mining, Inference, and Prediction, 2nd ed., publisherSpringer, New York.
  * Huber and Ronchetti (2009) authorHuber, P.J., and authorRonchetti, E.M. (year2009) titleRobust Statistics, 2nd ed., publisherWiley, Hoboken, NJ.
  * Kalogridis (2021) authorKalogridis, I. (year2021) titleAsymptotics for M-type smoothing splines with non-smooth objective functions, journalTEST volume31 pages373–389.
  * Kalogridis and Van Aelst (2021) authorKalogridis, I., and authorVan Aelst, S. (year2021) titleRobust penalized spline estimation with difference penalties, journalEcon. Statist., volumeappeared online.
  * Kauermann et al. (2009) authorKauermann, G., authorKrivobokova, T., and authorFahrmeir, L. (year2009) titleSome asymptotic results on generalized penalized spline smoothing, journalJ. R. Stat. Soc. Ser. B. Stat. Methodol. volume71 pages487–503.
  * Kim and Gu (2004) authorKim, Y.-J., and authorGu, C. (year2004) titleSmoothing spline Gaussian regression: more scalable computation via efficient approximation, journalJ. R. Stat. Soc. Ser. B. Stat. Methodol. volume66 pages337–356.
  * Mammen and van de Geer (1997) authorMammen, E., and authorvan de Geer, S. (year1997) titlePenalized Quasi-Likelihood Estimation in Partial Linear Models, journalAnn. Statist. volume25 pages1014–1035.
  * Maronna et al. (2019) authorMaronna, R.A., authorMartin, D., authorSalibián-Barrera, M. and authorYohai, V.J. (year2019) titleRobust Statistics: Theory and Methods, 2nd ed., publisherWiley, Chichester.
  * McCullagh and Nelder (1983) authorMcCullagh, P., and authorNelder, J.A. (year1983) titleGeneralized Linear Models, publisherChapman & Hall, London.
  * Nocedal and Wright (2006) authorNocedal, J., and authorWright, S.J. (year2006) titleNumerical Optimization, 2nd ed., publisherSpringer, New York.
  * O’Sullivan (1986) authorO’Sullivan, F. (year1986) titleA statistical perspective of ill-posed problems, journalStatist. Sci. volume1 pages502–518.
  * Pearce and Wand (2006) authorPearce, N.D., and authorWand, M.P. (year2006) titlePenalized Splines and Reproducing Kernel Methods, journalAmer. Statist. volume60 pages233–240.
  * Ruppert et al. (2003) authorRuppert, D., authorWand, M.P., and authorCarroll, R.J. (year2003) titleSemiparametric regression, publisherCambridge, NY.
  * Rynne and Youngston (2008) authorRynne, B., and authorYoungston, M.A. (year2008) titleLinear functional analysis, publisherSpringer, London.
  * Scott (2001) authorScott, D.W. (year2001) titleParametric Statistical Modeling by Minimum Integrated Square Error, journalTechnometrics volume43 pages274–285.
  * van de Geer (2000) authorvan de Geer, S. (year2000) titleEmpirical Processes in M-Estimation, publisherCambridge University Press, New York, NY.
  * van der Vaart (1998) authorvan der Vaart, A.W. (year1998) titleAsymptotic Statistics, publisherCambridge University Press, New York, NY.
  * Wahba (1990) authorWahba, G. (year1990) titleSpline models for observational data, publisherSiam, Philadelphia, Pen.
  * Wang et al. (2013) authorWang, X., authorJiang, Y., authorHuang, M., and authorZhang, H. (year2013) titleRobust Variable Selection With Exponential Squared Loss, journalJ. Amer. Statist. Assoc. volume108 pages632–643.
  * Wong et al. (2014) authorWong, R.K.W., authorYao, F., and authorLee, T.C.M. (year2014) titleRobust Estimation for Generalized Additive Models, journalJ. Comput. Graph. Statist. volume23 pages270–289.
  * Warwick and Jones (2005) authorWarwick, J., and authorJones, M.C. (year2005) titleChoosing a robustness tuning parameter, journalJ. Stat. Comput. Simul. volume75 pages581–588.
  * Wood (2017) authorWood, S.N. (year2017) titleGeneralized Additive Models, 2nd ed., publisherCRC Press, Boca Raton, FL.
  * Xiao (2019) authorXiao, L. (year2019) titleAsymptotic theory of penalized splines, journalElectron. J. Stat. volume13 pages747–794.
  * Zhou et al. (1998) authorZhou, S., authorShen, W., and authorWolfe, D.A. (year1998) titleLocal Asymptotics for Regression Splines and Confidence Regions, journalAnn. Statist. volume26 pages1760–1782.
  * Zou and Yuan (2008) authorZou, H., and authorYuan, M. (year2008) titleComposite quantile regression and the oracle model selection theory, journalAnn. Statist. volume36 pages1108–1126.