On a consistent finite-strain plate theory of growth

On a consistent finite-strain plate theory of growth

Accepted Manuscript On a consistent finite-strain plate theory of growth Jiong Wang, David Steigmann, Fan-Fan Wang, Hui-Hui Dai PII: DOI: Reference: ...

3MB Sizes 39 Downloads 72 Views

Accepted Manuscript

On a consistent finite-strain plate theory of growth Jiong Wang, David Steigmann, Fan-Fan Wang, Hui-Hui Dai PII: DOI: Reference:

S0022-5096(17)30459-3 10.1016/j.jmps.2017.10.017 MPS 3211

To appear in:

Journal of the Mechanics and Physics of Solids

Received date: Revised date: Accepted date:

1 June 2017 10 October 2017 26 October 2017

Please cite this article as: Jiong Wang, David Steigmann, Fan-Fan Wang, Hui-Hui Dai, On a consistent finite-strain plate theory of growth, Journal of the Mechanics and Physics of Solids (2017), doi: 10.1016/j.jmps.2017.10.017

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

On a consistent finite-strain plate theory of growth Jiong Wanga , David Steigmannb , Fan-Fan Wangc , Hui-Hui Daid,∗ a

AN US

CR IP T

School of Civil Engineering and Transportation, South China University of Technology, 510640 Guangzhou, Guangdong, China b Department of Mechanical Engineering, University of California, Berkeley, CA 94720, USA c Department of Mathematics,East China University of Science and Technology, Shanghai 200237, China d Department of Mathematics, City University of Hong Kong, 83 Tat Chee Avenue, Kowloon Tong, Hong Kong

Abstract

AC

CE

PT

ED

M

In this paper, a consistent finite-strain plate theory for growth-induced large deformations is developed. The three-dimensional (3D) governing system of the plate model is formulated through the variational approach, which is composed of the mechanical equilibrium equation and the constraint equation of incompressibility. Then, series expansions of the unknown functions in terms of the thickness variable are adopted. By using the 3D equilibrium equations and the surface boundary conditions, recursion relations for the expansion coefficients are successfully established. As a result, a 2D vector plate equation with three unknowns is obtained and the associated edge boundary conditions are proposed. It can be verified that the plate equation ensures the required asymptotic order for all the terms in the variations of the total energy functional. The weak formulation of the plate equation has also been derived for future numerical calculations. As applications of the plate theory, two examples regarding the growth-induced deformations and instabilities in thin hyperelastic plates are studied. Some analytical results are obtained in these examples, which can be used to describe the large deformations and reveal the bifurcation properties of the thin plates. Furthermore, the results obtained from the current plate theory are compared with those ∗

Corresponding author. Tel.: +852 34428660; fax: +852 34420250. Email addresses: [email protected] (Jiong Wang), [email protected] (David Steigmann), [email protected] (Fan-Fan Wang), [email protected] (Hui-Hui Dai) Preprint submitted to Journal of the Mechanics and Physics of Solids

October 27, 2017

ACCEPTED MANUSCRIPT

obtained from the classical F¨ oppl-von K´ arm´ an plate theory, from which the efficiencies and advantages of the current plate theory can be demonstrated.

CR IP T

Keywords: finite elasticity, growth, consistent plate theory, bifurcation analysis, analytical solutions 1. Introduction

AC

CE

PT

ED

M

AN US

Pattern formations are commonly observed from the morphogenesis of soft biological tissues during the growth (or atrophy) process, which have attracted extensive research interests in the fields of biology (Thompson, 1942; Nath et al., 2003; Coen et al., 2004), biochemistry (Givnish, 1987; Green, 1996), biophysics (Newell and Shipman, 2005; Marder and Papanicolaou, 2006) and biomechanics (Fung, 1990; Taber, 1995; Humphrey, 2002; Menzel and Kuhl, 2012). Besides the factors of genetic and chemistry, it is now understood that the mechanical effects play an important role in the pattern formations (Skalak et al., 1996; Lubarda and Hoger, 2002; Ben Amar and Goriely, 2005; Li et al., 2012; Rausch and Kuhl, 2013). In fact, growth (or atrophy) can induce the changes of body mass and the appearances of residual stress in biological tissues, which together with the external loads and environmental constraints result in large deformations of the tissues. Accompanying the large deformations, mechanical instabilities can be triggered, then various patterns are observed from the morphogenesis of the tissues. The patterns associated with the mechanical instabilities usually exhibits three phenotypes: wrinkling, folding and creasing (Jin et al., 2011; Li et al., 2012). At the continuum level, nonlinear elasticity theory (cf. Ogden, 1984) provides a nature framework for modeling the growth-induced instability in soft biological tissues. Following the approach of Rodriguez et al. (1994), the total deformation gradient tensor is usually decomposed into the multiplication of an elastic deformation tensor (representing the elastic response to the residual stress and the external loads) and a growth tensor (representing the changes of body mass). By adopting suitable material models (e.g., neo-Hookean, Mooney-Rivlin, Ogden) and through a conventional derivation procedure, the constitutive equation of the nominal stress tensor can be obtained, which satisfies the mechanical equilibrium equation and suitable boundary conditions. Furthermore, as soft biological materials are generally 2

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

assumed to be volumetric incompressible, another constraint equation for the elastic deformation tensor needs to be proposed. The mechanical equilibrium equation and the constraint equation formulate the governing system of the model. To study the mechanical instability phenomena in biological tissues, the bifurcation properties of the governing system need to be investigated. In the existing works, both the linear bifurcation and the post-bulking bifurcation analyses have been conducted to study the instability phenomena in soft materials (cf. Ben Amar and Goriely, 2005; Huang et al., 2005; Audoly and Boudaoud, 2008; Li et al., 2011; Cao and Hutchinson, 2012). In the current paper, we mainly focus on the samples with thin plate forms, which are commonly observed in biological organs or tissues, e.g., leaves, petals, wings of insects and human skin. To study the mechanical behaviors of plate samples made of biological (or general soft) materials, one usually needs to adopt a suitable plate model. Many of the previous works were conducted based on the well known F¨ oppl-von K´ arm´ an (FvK) plate theory. For example, Dervaux et al. (2009) derived a generalized FvK equation system based on the growth theory and through the variational approach, which was successfully applied to predict the bulking patterns in some typical examples. Cai et al. (2011) studied the periodic patterns of a nonlinear von K´ arm´ an plate bonded to a linear elastic foundation. For the different pattern modes, the energy states were analyzed by making use of an analytical upper-bound method and the finite element analysis. Several examples of continuous transitions between the different modes were also discussed. Budday et al. (2014) studied the development of the surface morphologies of brains during the growth process, where the cortex is viewed as a morphogenetically growing outer layer and the subcortex is viewed as a strain-driven growing inner core. An analytical model based on the FvK theory was established, which can be used to predict the critical conditions at the onset of folding. Huang et al. (2016) studied the surface wrinkling of a film-substrate system subject to uniaxial or biaxial compressions. The FvK plate equations were adopted to study the bulking of the thin film, which were solved by using both an analytical method and a semi-implicit Fourier spectral numerical method. Based on the obtained analytical and numerical solutions, the dependence of the wrinkling patterns on the mechanical property of the substrate and the external compressive strains were systematically investigated. Besides the works based on the FvK plate theory, some other approaches have also been proposed to study the mechanics of growing plates or shells. 3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

In Papastavrou et al. (2013) and Holland et al. (2013), biological surfaces were treated as growing membranes with zero thickness and was equipped with its own potential energy. By virtue of a surface projection operator, the governing system within a two-dimensional (2D) setting was established, which was further implemented into a finite element framework. The obtained numerical solutions can properly simulate the morphological changes of the growing surfaces during the onset of buckling and beyond. In Liu et al. (2013), an inhomogeneous field theory of hydrogel deformation was introduced, which was applied to describe the deformation patterns of leaves and surfaces of fruits during the growing or drying processes. Rausch and Kuhl (2014) proposed a computational model for finite membrane growth based on the Kirchhoff shell theory, which contains the constitutive equations of collagenous tissues with a pronounced microstructural direction and the equations for stretch-driven membrane growth. After the computational implementation of this model, some typical examples of growing biological membranes under different loading conditions were studied. Despite the existences of the above works and many other insightful results obtained, a number of issues in the mechanics of pattern formations in plate biological samples have not been resolved at a satisfactory level. In our opinion, the research in this area should be further developed in the following two aspects: 1) formulation of a comprehensive plate theory within the framework of finite stain elasticity (the FvK plate theory is for small strains) incorporating both bending and stretching (not either bending or stretching) by taking into account the growth effects, elastic incompressibility and certain consistency requirements (with no ad hoc hypotheses) ; 2) more analytical results from the post bifurcation analyses are desired to explain the key features of the pattern formations in thin growing samples. The development of plate theories has a long history and the literature in this field is extremely rich (cf. Timoshenko and Woinowsky-Krieger, 1959; Reddy, 2007). Some classical plate theories, e.g., the Kirchhoff-Love theory, the F¨ oppl-von K´ arm´ an (FvK) theory and the Mindlin-Reissner theory, were proposed based on a priori hypotheses, which may not be consistent with the three-dimensional (3D) formulation, and are not suitable for relatively thick plates or for finite-strain problems. Some other plate theories were established within the framework of finite strain elasticity. Usually, the kinematic quantities are first expanded along the thickness direction, and then the governing equations are derived from the 2D variational (or virtual work) principle by integrating over the thickness and conducting a truncation. De4

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

pending on the different levels of approximations, different kinds of plate theories can be obtained (Meroueh, 1986; Steigmann, 2007). In Steigmann (2013a,b), for the traction-free case, optimum truncated plate/shell energy functions incorporating both stretching and bending effects were established, from which the Koiter’s plate and shell models can be obtained with further approximations. In Dai and Song (2014), a new finite-strain plate theory was proposed for compressible hyperelastic materials, which can achieve the term-wise consistency with the 3D stationary potential energy to a required order. In a recent work of Wang et al. (2016), the plate theory proposed in Dai and Song (2014) was extended for incompressible hyperelastic materials. As the constraint equation of incompressibility and some additional variables need to be considered, the derivation procedure of the plate theory for incompressible materials becomes much more complicated. In the current work, the plate theory reported in Wang et al. (2016) will be further developed to study growth-induced large deformations in thin hyperelastic plates. To consider the growth effect, a growth tensor is incorporated in the total deformation gradient tensor. The 3D governing system for the mechanical response of the plate sample can be formulated through the variational approach, which is composed of the mechanical equilibrium equation and the constraint equation of elastic incompressibility. To derive a 2D vector plate equation from the 3D governing system, series expansions of the unknown functions in terms of the thickness variable are adopted. By using the bottom surface traction condition and the governing differential equations, we successfully derive the recursion relations for the expansion coefficients. By further considering top surface traction conditions, a 2D vector plate equation with 3 unknowns is obtained. To complete the plate equation system, suitable edge boundary conditions are also proposed. It can be verified that the plate equation ensures the required asymptotic order (the 4th order of the thickness) for all the terms in the variations of the total energy functional. The weak formulation of the plate equation are also derived for future numerical calculations. Recently, Steigmann (2015) developed a plate theory for materially uniform thin film for the traction-free case by deriving the optimum truncated plate strain energy. The current paper considers the presence of tractions at both top and bottom surfaces and the plate theory is derived by a different method, which relies on the establishment of the recursion relations. To show the efficiency of the current plate theory, two examples are studied. One is the large bending deformation of a thin plate induced by a gra5

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN US

CR IP T

dient growth field and the other is the growth-induced instability of a thin plate resting on a Winkler foundation. For simplicity, we only consider the uniaxial growth case and the plane-strain deformation of the plates. In the first example, the analytical solutions to the plate equations are derived and compared with the exact solutions of the original governing system, which show very good accuracies. In the second example, the bifurcation properties of the plate equation system are investigated through both the linear bifurcation and the post-bifurcation analyses. The critical bifurcation points corresponding to the nontrivial solutions with the different modes are determined. The amplitudes of the leading-order nontrivial solutions can also be derived, which is compared with the numerical solutions to show their accuracies. To demonstrate the advantages of the current plate theory, it was also compared with the classical FvK plate theory. In the uniaxial growth case, we show that the FvK-type plate equations can be derived from the current plate theory by adopted some extra scaling relations and assumptions. Some further investigations reveal that the FvK plate theory is not suitable for studying the first example. While, for the second example, the results obtained from the two plate theories are close to each other. This paper is organized as follows. In section 2, the 3D governing system of the plate model is established through the variational approach. In section 3, the 2D vector plate equation is derived and some proper edge boundary conditions are proposed. The consistency of the plate equation with the variational formulation is also verified. In section 4, the associated weak formulation of the 2D plate equations is set up. In section 5, two examples are studied and some analytical solutions are obtained, based on which the large deformations and the bifurcation properties of the thin plates are investigated. The relationships between the current plate theory and the Fvk plate theory are also analyzed. Finally, some conclusions are drawn. 2. The variational principle and the 3D governing system

AC

We consider a thin plate with constant thickness, which occupies the region Ω × [0, 2h] in its reference configuration κ. The total deformation of the plate comes from both the growth process and the elastic deformation (induced by external loads or boundary restrictions). Accompanying the deformation, the plate is transformed into its current configuration κt . The coordinates of a material point in κ and κt are denoted by X = (r, Z) and

6

ACCEPTED MANUSCRIPT

x, respectively. The total deformation gradient tensor is then given by F=

∂x ∂x ∂x ∂x = + ⊗ k = ∇x + ⊗ k, ∂X ∂r ∂Z ∂Z

(1)

F = AG,

CR IP T

where ‘∇, is the in-plane 2D gradient and k is the unit normal vector of the reference top surface Ω. Following the basic assumption of elastic growth theory (Rodriguez et al., 1994; Skalak et al., 1996; Ben Amar and Goriely, 2005), the total deformation gradient tensor F can be decomposed into (2)

AN US

where A is the elastic deformation tensor and G is the growth tensor. Generally, biological materials are elastically incompressible, thus we have the following constraint equation R(F, G) = R0 (FG−1 ) = R0 (A) = DetA − 1 = 0,

(3)

M

where Det(·) denotes the determinant of a tensor. Furthermore, we suppose that the material has the elastic strain-energy function φ(F, G) = JG φ0 (FG−1 ) = JG φ0 (A),

(4)

AC

CE

PT

ED

where JG = DetG. As the rate of deformation of the growth process is very small compared with that of the elastic deformation (Ben Amar and Goriely, 2005; Dervaux et al., 2009), in the following analyses the distribution of the growth tensor G in the plate is assumed to be given and does not change. The lateral surface of the plate, which is denoted as ∂Ω, is composed of the position boundary ∂Ω0 and the traction boundary ∂Ωq . For the case of dead-loading, we denote q± as the applied traction on the top and bottom ˜ as the applied traction on ∂Ωq . If neglecting the surfaces of the plate and q body force, the potential energy of the plate ψ is given by Z Z 2h ¯ ¯ ¯ ψ =φ − V , φ = JG φ0 (FG−1 )dZdr, Ω 0 Z − V¯ = q (r) · x(r, 0) + q+ (r) · x(r, 2h)dr (5) Ω Z Z 2h ˜ (s, Z) · x(s, Z)dZds, + q ∂Ωq

0

7

ACCEPTED MANUSCRIPT

By further adopting the Lagrangian multiplier p(X), the generalized potential energy functional of the plate can be written as Z Z 2h (6) p(X)JG R0 (FG−1 )dZdr. Ψ =ψ − 0

CR IP T



∂Ωq



+

Z

∂Ω0

+

0



Z

∂Ωq

DivS · δxdZdr −

Z





ST k + q− (r) · δx(r, 0)dr

(7)

ST k − q+ (r) · δx(r, 2h)dr Z 2h ST N · δx(s, Z)dZds Z

M

+

ZΩ

0 2h

0

2h



˜ (s, Z) · δx(s, Z)dZds, ST N − q

ED

=−

Z Z

AN US

Next, we shall calculate the variations of Ψ with respect to the independent variables x and p to derive the 3D governing system. First, the variation of Ψ with respect to x yields that   Z Z 2h  Z Z 2h  δΨ ∂φ0 ∂R0 tr JG tr pJG = · δF dZdr − · δF dZdr δx ∂F ∂F Ω 0 Ω 0 Z − q− (r) · δx(r, 0) + q+ (r) · δx(r, 2h)dr Ω Z Z 2h ˜ (s, Z) · δx(s, Z)dZds − q

0

CE

PT

where S = JG [∂φ0 /∂F − p(∂R0 /∂F)] is the nominal stress tensor. By virtue of (4), S can be further expressed as   ∂φ0 (FG−1 ) ∂R0 (FG−1 ) −p S(x, p) =JG ∂F ∂F   (8) ∂φ0 (A) ∂R0 (A) −1 =JG G −p . ∂A ∂A

AC

By considering the arbitrariness of δx in (7), we obtain the following PDE system DivS = 0, in Ω × [0, 2h], ST k|Z=0 = −q− (r), ST k|Z=2h = q+ (r), x = b(s, Z), on ∂Ω0 × [0, 2h],

ST N = q(s, Z),

on ∂Ωq × [0, 2h]. 8

on Ω,

(9)

ACCEPTED MANUSCRIPT

(10)

CR IP T

Besides that, the variation of Ψ with respect to p yields that Z Z 2h δΨ =− JG R0 (FG−1 )δpdZdr, δp Ω 0

from which one just obtains the constraint equation (3). Equations (3) and (9) formulate the governing system for the mechanical behavior of the thin hyperelastic plate in 3D case. 3. The vector plate equation

AN US

3.1. Derivation of the vector plate equation The governing equations (3) and (9) contain two independent unknowns x and p. To make approximations to eliminate the variable Z and obtain a 2D formulation, we consider the following series expansions of x and p in terms of Z about the bottom surface (assuming sufficient smoothness)

(11)

M

1 1 x(X) =x(0) (r) + Zx(1) (r) + Z 2 x(2) (r) + Z 3 x(3) (r) 2 6 1 4 (4) + Z x (r) + O(Z 5 ), 24

PT

ED

1 1 p(X) =p(0) (r) + Zp(1) (r) + Z 2 p(2) (r) + Z 3 p(3) (r) 2 6 (12) 1 4 (4) + Z p (r) + O(Z 5 ), 24 (n) n where (·) = ∂ (·)/∂Z n |Z=0 (n = 0, 1, · · · , 4). Accordingly, the deformation gradient tensors F, A, G−1 and JG G−1 have the expansions

(13)

AC

CE

1 1 F = F(0) (r) + ZF(1) (r) + Z 2 F(2) (r) + Z 3 F(3) (r) + O(Z 4 ), 2 6 1 1 A = A(0) (r) + ZA(1) (r) + Z 2 A(2) (r) + Z 3 A(3) (r) + O(Z 4 ), 2 6 1 −1 (0) (1) 2 (2) ¯ (r) + Z G ¯ (r) + Z G ¯ (r) + 1 Z 3 G ¯ (3) (r) + O(Z 4 ), G =G 2 6 1 −1 (0) (1) 2 (2) ˆ (r) + Z G ˆ (r) + Z G ˆ (r) + 1 Z 3 G ˆ (3) (r) + O(Z 4 ). JG G = G 2 6

By substituting (11) into (1) and comparing with (13)1 , one obtains F(n) = ∇x(n) + x(n+1) ⊗ k, 9

n = 0, 1, 2, 3.

(14)

ACCEPTED MANUSCRIPT

CR IP T

Further substituting (13)1−3 into (2) and considering the coefficients of Z n , we have ¯ (0) , A(0) =F(0) G ¯ (0) + F(0) G ¯ (1) , A(1) =F(1) G (15) ¯ (0) + 2F(1) G ¯ (1) + F(0) G ¯ (2) , A(2) =F(2) G ¯ (2) . ¯ (2) + F(0) G ¯ (1) + 3F(1) G ¯ (0) + 3F(2) G A(3) =F(3) G Similar as the deformation gradient tensors, the nominal stress tensor S can also be expanded as follows.

(16)

AN US

1 1 S(X) = S(0) (r) + ZS(1) (r) + Z 2 S(2) (r) + Z 3 S(3) (r) + O(Z 4 ). 2 6

×

M

On the other hand, from (8) and (13)2,4 , it can be derived that   ∂R0 (A) ∂φ0 (A) −1 −p S =JG G ∂A ∂A   1 1 (0) (1) 2 3 (2) (3) 4 ˆ + ZG ˆ + Z G ˆ + Z G ˆ + O(Z ) = G 2 6 ( A(0) (A(0) ) + A(1) (A(0) )[A − A(0) ]

(17)

AC

CE

PT

ED

1 + A(2) (A(0) )[A − A(0) , A − A(0) ] 2 1 + A(3) (A(0) )[A − A(0) , A − A(0) , A − A(0) ] 6   1 2 (2) 1 3 (3) (0) (1) 4 − p + Zp + Z p + Z p + O(Z ) 2 6 n (0) (0) (1) (0) × R (A ) + R (A )[A − A(0) ] 1 + R(2) (A(0) )[A − A(0) , A − A(0) ] 2 ) o 1 (3) (0) + R (A )[A − A(0) , A − A(0) , A − A(0) ] , 6

where

A(i) (A(0) ) =

∂ i+1 R0 ∂ i+1 φ0 (i) (0) | R (A ) = | (0) , (0) , A=A ∂Ai+1 ∂Ai+1 A=A 10

i = 0, 1, 2.

ACCEPTED MANUSCRIPT

AN US

CR IP T

The elastic moduli A(i) (i = 0, 1, 2) can be calculated with any given strain energy function φ0 . The explicit expressions of R(i) (i = 0, 1, 2) can be found in Appendix A. Comparison of (16) and (17) yields that   ˆ (0) A(0) − p(0) R(0) , S(0) (r) =G   ˆ (0) A(1) [A(1) ] − p(0) R(1) [A(1) ] − p(1) R(0) S(1) (r) =G   ˆ (1) A(0) − p(0) R(0) , +G  (2) (0) ˆ S (r) =G A(1) [A(2) ] + A(2) [A(1) , A(1) ] − p(0) R(1) [A(2) ]  − p(0) R(2) [A(1) , A(1) ] − 2p(1) R(1) [A(1) ] − p(2) R(0)   ˆ (1) A(1) [A(1) ] − p(0) R(1) [A(1) ] − p(1) R(0) +G   ˆ (2) A(0) − p(0) R(0) , +G  ˆ (0) A(1) [A(3) ] + 3A(2) [A(1) , A(2) ] + A(3) [A(1) , A(1) , A(1) ] S(3) (r) =G − p(0) (R(1) [A(3) ] + 3R(2) [A(1) , A(2) ] + R(3) [A(1) , A(1) , A(1) ])

AC

CE

PT

ED

M

 − 3p(1) (R(1) [A(2) ] + R(2) [A(1) , A(1) ]) − 3p(2) R(1) [A(1) ] − p(3) R(0)  ˆ (1) A(1) [A(2) ] + A(2) [A(1) , A(1) ] − p(0) R(1) [A(2) ] +G  (0) (2) (1) (1) (1) (1) (1) (2) (0) − p R [A , A ] − 2p R [A ] − p R   ˆ (2) A(1) [A(1) ] − p(0) R(1) [A(1) ] − p(1) R(0) +G   ˆ (3) A(0) − p(0) R(0) . +G (18) (i) From (14), (15) and (18), it can be observed that the dependence of S on x(i+1) (i = 1, 2, 3) is linearly algebraic. Due to the series expansions (11) and (12), the unknown functions in the governing system become x(n) (n = 0, · · · , 4) and p(n) (n = 0, · · · , 3). Next, we need to formulate a closed system for these 19 unknowns. First, from the bottom traction condition (9)2 , we obtain n oT  (0) T (0) (0) (0) (0) (0) ¯ (0) (1) (0) T ˆ ¯ S k = G (A − p R )(∇x G + x ⊗ (G ) k) k (19) − = −q (r). 11

ACCEPTED MANUSCRIPT

Equation (19) provides three algebraic equations for the unknowns x(1) and p(0) . Next, from the mechanical field equation (9)1 and by considering the series expansion of S, one obtains T

k = 0,

n = 0, 1, 2.

(20)

CR IP T

∇ · S(n) + S(n+1)

Equation (20) provides three algebraic equations for the unknowns x(n+2) and p(n+1) (n = 0, 1, 2). Further substituting (13)2 into the constraint equation (3) and considering the coefficients of Z n , we can derive the following equations

AN US

R0 (A(0) ) = 0, R(0) [A(1) ] = 0,

R(0) [A(2) ] + R(1) [A(1) , A(1) ] = 0,

(21)

R(0) [A(3) ] + 3R(1) [A(1) , A(2) ] + R(2) [A(1) , A(1) , A(1) ] = 0.

T

k+2h S(1)

T

k+2h2 S(2)

ED

S(0)

M

The equations in (21) provide the additional equations for the unknowns x(n+1) (n = 0, 1, 2, 3). Finally, from the top traction condition (9)3 , we obtain T

T 4 k+ h3 S(3) k+O(h4 ) = q+ (r). (22) 3

AC

CE

PT

This equation contains all the unknowns x(n) (n = 0, 1, · · · , 4) and p(n) (n = 0, 1, 2, 3). In summary, (19)-(22) formulate a closed system for the 19 unknowns (n) x (n = 0, · · · , 4) and p(n) (n = 0, · · · , 3). Notice that (20) and (21)(2−4) are linearly algebraic equations with respect to the highest orders of x(n+1) and p(n) (n = 1, 2, 3). By solving these equations, we can obtain the expressions of x(2) -x(4) and p(1) -p(3) in terms of x(0) , x(1) and p(0) . In fact, the explicit expressions of x(4) and p(3) are not needed in the current work. Thus, only the expressions of x(2) , p(1) and x(3) , p(2) are given below.  −1 ¯ (0) − B−1 f (2) ⊗ (G ¯ (0) )T k x(2) = − B−1 f (2) − D(0) R(0) ∇x(1) G  ¯ (1) B−1 (R(0) )T (G ˆ (0) )T k, +F(0) G   −1 ¯ (0) − B−1 f (2) ⊗ (G ¯ (0) )T k + F(0) G ¯ (1) . p(1) = − D(0) R(0) ∇x(1) G 12

(23)

ACCEPTED MANUSCRIPT

CR IP T

ˆ (0) )T k x(3) = − B−1 f (3) + p(2) B−1 (R(0) )T (G T (0) T  ˆ ) k + p(1) B−1 (R(0) )T (G ˆ (1) )T k, (G + 2p(1) B−1 R(1) A(1)    (2) (0) −1 ¯ (0) )T k p =−D R(1) [A(1) , A(1) ] + R(0) −B−1 f (3) ⊗ (G i h T (0) T ˆ ) k ⊗ (G ¯ (0) )T k + 2p(1) R(0) B−1 R(1) [A(1) ] (G h i ˆ (1) )T k ⊗ (G ¯ (0) )T k + p(1) R(0) B−1 (R(0) )T (G    (2) (0) (0) ¯ (2) (1) ¯ (1) (0) ¯ , ∇x G + 2F G + F G +R

AN US

where

(24)

¯ (0) )T k)α ((G ¯ (0) )T k)β , (B)ij =J (0) (A(1) − p(0) R(1) )αiβj ((G   ¯ (0) )T k ⊗ (G ¯ (0) )T k , D(0) =R(0) J (0) B−1 (R(0) )T (G  ¯ (0) + F(0) G ¯ (1) ] T (G ˆ (0) )T k f (2) =∇ · S(0) + (A(1) − p(0) R(1) )[∇x(1) G ˆ (1) )T k, + (A(0) − p(0) R(0) )T (G

ED

M

f (3) =∇ · S(1)   ¯ (0) + 2F(1) G ¯ (1) + F(0) G ¯ (2) ] T (G ˆ (0) )T k + A(1) − p(0) R(1) [∇x(2) G   T (0) T ˆ ) k + A(2) − p(0) R(2) [A(1) , A(1) ] (G   T (1) T  ˆ ) k + A(0) − p(0) R(0) T (G ˆ (2) )T k, + A(1) − p(0) R(1) [A(1) ] (G

AC

CE

PT

and J (0) = JG |Z=0 . Here, we have assumed the strong ellipticity of the constitutive relations, which guarantees the invertibility of B. The explicit expressions of x(1) and p(0) in terms of x(0) should be derived based on (19) and (21)1 . However, as these two equations are nonlinear algebraic equations, they can only be solved with the given constitutive form of φ0 (A). For incompressible neo-Hookean materials, the expressions of x(1) and p(0) in terms of x(0) have been given in Appendix B. By subtracting the top and bottom traction conditions (19) and (22), the following 2D vector plate equation can be obtained

where

¯= 1 S 2h

¯ = −¯ ∇·S q, Z

0

2h

2 SdZ = S(0) + hS(1) + h2 S(2) + O(h3 ), 3 13

(25)

¯= q

q+ + q− . 2h

ACCEPTED MANUSCRIPT

CR IP T

Equation (25) involves quantities up to x(3) and p(2) . After substituting the expressions of x(n) (n = 1, 2, 3) and p(n) (n = 0, 1, 2) into (25), it becomes a fourth-order differential equation for x(0) . Remark: Formally, the 2D vector plate equation (25) attains the order of h2 . However, in the case of linear elasticity, we find that the plate equation can only attain the order of h. Some discussions regarding this issue have been given in Appendix C. The readers can also refer to Song and Dai (2016) for a more throughout investigation on this issue.

PT

ED

M

AN US

3.2. Edge boundary conditions and consistency examination To ensure the completeness of the plate equation system, suitable boundary conditions also need to be proposed on the edge of the plate area ∂Ω. As the plate equation is a fourth-order differential equation, on both the position boundary ∂Ω0 and the traction boundary ∂Ωq , two conditions regarding x(0) or its derivatives are required. In the authors’ previous works (Dai and Song, 2014; Wang et al., 2016), both the position and the traction edge boundary conditions have been proposed by considering the original 3D lateral surface conditions, which will also be adopted in the current work. For the purpose of being self-contained, these boundary conditions are simply introduced below: Case 1. Position boundary conditions Suppose that on ∂Ω0 × [0, 2h] the position b is prescribed, then the following two boundary conditions can be proposed Z 2h Z 2h 1 1 (0) (0) ¯ on ∂Ω0 ¯= xdZ = bdZ = b, x = b (s), x 2h 0 2h 0 (26) 2 (2) 1 2 (3) 1 ¯ (0) (0) (1) 3 (0) ⇔ x = b (s), x + hx + h x + O(h ) = (b − b ), 3 3 h

AC

CE

where b(0) = b|Z=0 and a bar over a quantity represents the through-thickness average. By substituting the expressions of x(n) (n = 1, 2, 3) into (26)2 , a relation involving the third-order derivatives of x(0) can be obtained. Case 2. Traction boundary conditions Suppose that on ∂Ωq × [0, 2h] the traction q is specified, then the following

14

ACCEPTED MANUSCRIPT

two boundary conditions can be proposed T

(27)

CR IP T

S(0) N = q(0) , Z 2h Z 2h 1 1 T T ¯ ¯ on ∂Ωq , S NdZ = qdZ = q S N= 2h 0 2h 0 T  2 2 (2) (0) T (0) (0) (1) 3 ¯, ⇔S N=q , S + hS + h S + O(h ) N = q 3

AC

CE

PT

ED

M

AN US

¯ is the averaged traction along the thickness of the where q(0) = q|Z=0 and q plate. We remark that (27)1 can be replaced by specifying the moment about the middle line, see Dai and Song (2014). Equation (25) together with the edge boundary condition (26) or (27) formulate the 2D governing system of our current plate theory. In Dai and Song (2014), a consistency criterion has been proposed for plate theories, which requires that for all loadings that satisfy some smooth requirements, each term in the first variation of the energy functional should be of a required asymptotic order (say, O(h4 )). According to this criterion, the consistency of our current plate theory can be examined by analyzing the asymptotic order of each term in the variations (7) and (10). In fact, the asymptotic orders of the terms in (7) and (10) can be checked through the similar methods as that introduced in Wang et al. (2016). For the first term on the r.h.s. of (7), we can consider the series expansion of DivS in terms of Z. As the first three terms in the series expansion (cf. (20)) have been used together with (21)(2−4) to obtain the recursion relations of x(i) (i = 2, 3, 4) and p(i) (i = 1, 2, 3), we have DivS = O(Z 3 ). Thus, the first term is of O(h4 ). The asymptotic orders of the second and the third terms on the right hand side (r.h.s.) of (7) can be analyzed by considering the top and bottom traction conditions (19) and (22), which also attain O(h4 ). To analyze the fourth term on the r.h.s. of (7), we rewrite the integrand into the following form   Z 2h Z 2h 1 2 (2) (0) (1) T (0) T dZ S N · δxdZ = S N · δx + Zδx + Z δx 2 0 0 Z 2h  T (28) + ZS(1) N · δx(0) + Zδx(1) dZ 0

+

Z

0

2h

1 2 (2) T Z S N · δx(0) dZ + O(h4 ). 2 15

ACCEPTED MANUSCRIPT

AN US

CR IP T

From the position boundary conditions (26), it can be obtained that δx(0) = 0, δx(1) = O(h) and δx(1) + 2/3hδx(2) = O(h2 ). With these results, it is easy to check that all the three terms on the r.h.s. of (28) are of O(h4 ). Thus, the fourth term in (7) also satisfies the consistency condition. For the fifth term on the r.h.s. of (7), we can also expand the terms in the integrand in terms of Z and then use the traction boundary conditions (27) to analyze its asymptotic order (cf. (33) in Wang et al. (2016)). It can be verified that this term is also of O(h4 ). For the variation (10), the integrand R0 (FG−1 ) = R0 (A) has been expanded in terms of Z and the coefficients of Z i (i = 0, 1, 2, 3) are set to be zero in (21) to derive the recursion relations of x(i+1) (i = 0, 1, 2, 3). Thus the variation (10) also attains O(h4 ). In summary, all the terms in the variations (7) and (10) attain O(h4 ). Therefore, our current plate theory satisfies the consistency criterion introduced in Dai and Song (2014). 4. The associated weak formulation

ED

M

To be applicable for numerical calculations, the weak formulation of the plate equation (25) also needs to be derived. For that purpose, we multiply both sides of (25) with the term ξ = δx(0) and calculate the integrations over the region Ω, which yields that Z Z  ¯ ¯ · ξdr ⇒ ∇ · S · ξdr = − q ZΩ Z Ω Z (29)  T ¯ N · ξds − S ¯ : ∇ξdr = − q ¯ · ξ dr. S

PT

∂Ω





AC

CE

It should be noted that the weak formulation (29) involves the third-order derivative of x(0) , which originates from the terms F(2) and p(2) in S(2) . Next, we intend to eliminate these third-order derivative terms. Based on (14), (15)3 and (24), and through some calculations, we obtain that (2) (2) (2) (2) p(2) = p1 + p2 , A(2) = A1 + A2 . (30) (2)

(2)

The terms p1 and A1 do not contain the third-order derivative of x(0) , thus (2) their lengthy expressions are omitted here for brevity. The terms p2 and

16

ACCEPTED MANUSCRIPT

(2)

A2 have the following expressions       (2) (0) −1 ¯ (0) )T k ¯ (0) − R(0) B−1 ∇ · S(1) ⊗ (G p2 = − D R(0) ∇x(2) G

CR IP T

h i   (1) (2) ¯ (0) T ˆ (0) T (0) T (0) (1) −1 ¯ ) k , [∇x G ] (G ) k ⊗ (G A −p R B −R  ¯ (0) − B−1 ∇ · S(1) ⊗ (G ¯ (0) )T k =∇x(2) G   ¯ (0) ] T (G ˆ (0) )T k ⊗ (G ¯ (0) )T k − B−1 A(1) − p(0) R(1) [∇x(2) G (2) ˆ (0) )T k ⊗ (G ¯ (0) )T k. + p B−1 (R(0) )T (G (0)

(2)

A2

2

AN US

Similarly, based on (18) and (25), the following results can be derived ¯ : ∇ξ = W1 (∇∇x(0) , ∇ξ) + W2 (∇∇∇x(0) , ∇ξ). S where

(31)

PT

ED

M

 W1 = S(0) + hS(1) : ∇ξ o 2 ˆ (0) n (1) (2) (2) + h2 G (A − p(0) R(1) )[A1 ] − p1 R(0) : ∇ξ 3 2 2 ˆ (0)  (2) (A − p(0) R(2) )[A(1) , A(1) ] − 2p(1) R(1) [A(1) ] : ∇ξ + hG 3 (32) 2 ˆ (1)  (1) + h2 G (A − p(0) R(1) )[A(1) ] − p(1) R(0) : ∇ξ 3 2 ˆ (2) (0) (A − p(0) R(0) ) : ∇ξ, + h2 G 3 o 2 2 ˆ (0) n (1) (2) (2) W2 = h G (A − p(0) R(1) )[A2 ] − p2 R(0) : ∇ξ. 3 (2)

(2)

AC

CE

The expressions of A2 and p2 are substituted into (32)2 to eliminate the third-order derivative terms. Through some calculations, we obtain that

where

(1)

(2)

W2 = W2 + W2 ,  2  (1) W2 = h2 H0 : ∇x(2) + η · ∇ · S(1) , 3  2  (2) W2 = h2 I0 : ∇x(2) + ζ · ∇ · S(1) . 3 17

(33)

(34)

ACCEPTED MANUSCRIPT

and

CR IP T

¯ (0) (A(1) − p(0) R(1) )[∇ξ G ˆ (0) + η ⊗ (G ˆ (0) )T k], H0 =G n oT ˆ (0) ] (G ¯ (0) )T k, η = − B−1 (A(1) − p(0) R(1) )[∇ξ G  (0) (0) −1 ˆ (0) + η ⊗ (G ˆ (0) )T k] G ¯ R I0 =D(0) R(0) [∇ξ G o  (0) (1) (0) (1) −1 (0) T ˆ (0) T (0) T ¯ ¯ −G (A − p R )[B R (G ) k ⊗ (G ) k] ,

(35)

−1 ˆ (0) + η ⊗ (G ˆ (0) )T k]B−1 (R(0) )T (G ¯ (0) )T k. ζ = − D(0) R(0) [∇ξ G





where Z

Z

AN US

Based on the above analyses, the second term in (29) can be rewritten as Z Z Z Z (2) (1) ¯ (36) W1 dr + W2 dr + W2 dr, S : ∇ξdr = Ω



Z  (2) 2 2 2 2  (1) T  T = h H0 N · x ds + h S N · ηds 3 3 ∂Ω Ω ∂Ω Z 2 2 h (∇ · H0 ) · x(2) + S(1) : ∇η dr, − Ω 3 Z Z Z 2 2 T  (2) 2 2  (1) T  (2) W2 dr = h I0 N · x ds + h S N · ζds 3 ∂Ω Ω ∂Ω 3 Z 2 2 h (∇ · I0 ) · x(2) + S(1) : ∇ζ dr. − Ω 3

PT

ED

M

(1) W2 dr

AC

CE

By substituting (36) into (29), we obtain the following 2D weak formulation Z ¯ · ξ) dr (W1 + W3 − q Ω  Z    2 2 T T 2 T (2) 2 (1) ¯ N · ξ − h [H0 + I0 ] N · x − h S N · (η + ζ) ds, = S 3 3 ∂Ω (37) where 2  W3 = − h2 (∇ · (H0 + I0 )) · x(2) + S(1) : ∇(η + ζ) . 3 In fact, it can be proved that (cf. Wang et al., 2016) δx(1) = η + ζ,

δS(0) = H0 + I0 . 18

(38)

ACCEPTED MANUSCRIPT

AN US

CR IP T

By considering the edge boundary condition on ∂Ω, the 2D weak formulation given in (37) can be further simplified. In the following paragraphs, we rewrite (37) for two distinct situations, according to the different types of edge boundary conditions. Case 1. Edge position and traction in the 3D formulation are known In this case, on the right-hand side of (37), the integration on the position boundary ∂Ω0 is of O(h3 ) and is thus ignored. While on the traction T boundary ∂Ωq , it is easy to deduce that (H0 + I0 )T N = δ[S(0) N] = O(h2 ). Thus, the second term in the boundary integral can also be neglected. Then, the 2D weak formulation (37) can be reduced to Z ¯ · ξ) dr (W1 + W3 − q Ω  Z  2 2  (1) T  T ¯ = S N · ξ − h S N · (η + ζ) ds (39) 3 ∂Ωq Z = q0 · ξ − 2m0 · (η + ζ)ds, ∂Ωq

CE

PT

ED

M

The above weak formulation is suitable for finite element calculations with prescribed q0 and m0 . Case 2. Edge position and traction in the 3D formulation are unknown In practical situations, one usually does not know the exact traction distributions or displacement distributions on the edge boundary conditions. In these cases, the so-called natural boundary conditions can be proposed according to the weak formulation. To this end, we shall rewrite the boundary integral in (37) by using ξ,N , i.e., the normal derivative of ξ. First, we introduce a third-order tensor M in the following form  2  T (40) − h2 S(1) N · (η + ζ) + [H0 + I0 ]T N · x(2) = (M[N])T : ∇ξ, 3

AC

and denote

−1 ˆ (0) )T k, V =D(0) (R(0) )T (G −1  (0) (0) ˆ R W(0) =D(0) G ˆ (0) (A(1) − p(0) R(1) )[B−1 (R(0) )T (G ˆ (0) )T k ⊗ (G ¯ (0) )T k] , −G

A¯ =A(1) − p(0) R(1) .

19

ACCEPTED MANUSCRIPT

Next, we introduce the decomposition ∇ξ = ξ,S ⊗ T + ξ,N ⊗ N,

(41)

CR IP T

where T and ξ,S are respectively the unit tangential vector and tangent derivative. Substituting (40) and (41) into (37) and a simple integration by part (assuming that ∂Ω is smooth) leads to Z ¯ · ξ) dr (W1 + W3 − q ZΩ (42)  T ¯ = S N − (M[N, T]),S · ξ + M[N, N] · ξ,N ds, ∂Ω

AN US

where

2 n T M[N, N] = h2 BT1 B−1 S(1) N + BT1 B−1 B1 x(2) − B2 x(2) 3   T + (B−1 S(1) N) · V − W(0) x(2) · N

¯ (0) )T N − BT B−1 (R(0) )T (G ¯ (0) )T k × (R(0) )T (G 1

ED

M

2 2 n T −1 (1) T M[N, T] = h B3 B S N + BT3 B−1 B1 x(2) − B4 x(2) 3   T + (B−1 S(1) N) · V − W(0) x(2) · N

¯ (0) )T N − BT B−1 (R(0) )T (G ¯ (0) )T k × (R(0) )T (G 3

 ¯ (0) )T k (B1 )ij = J (0) A¯ijkl (G i ¯ (0) )T N (B2 )ij = J (0) A¯ijkl (G i ¯ (0) )T k (B3 )ij = J (0) A¯ijkl (G i (0) T (0) ¯ ¯ (B4 )ij = J Aijkl (G ) T i

PT

CE

and

 ¯ (0) )T N , (G k ¯ (0) )T N , (G k ¯ (0) )T T , (G k (0) T ¯ (G ) N .

o

,

o ,

k

AC

In the standard plate theory, the terms before ξ and ξ,N are considered as the generalized average traction and generalized bending moment, respectively. Formally, if we regard W := W1 + W3 as the variation of the plate stress work due to the virtual displacement ξ, the weak formulation (42) can be rewritten as Z Z Z ¯ · ξdr + ˆ (s) · ξ + m ˆ 0 (s) · ξ,N ds, Wdr = q q (43) Ω



∂Ω

20

ACCEPTED MANUSCRIPT

CR IP T

ˆ and m ˆ 0 are respectively the applied generalized traction and bending where q moment at the edge. Based on (42) or (43), one can propose the boundary conditions for various practical cases, e.g., the clamped edge, the pinned edge, the simply-supported edge (see the discussions given in Dai and Song (2014)). 5. Examples: growth-induced deformations and instabilities in hyperelastic plates

CE

PT

ED

M

AN US

As applications of the current plate theory, two examples will be studied in this section: (1) the large bending deformation of a thin plate induced by a gradient growth field; (2) the growth-induced instability of a thin plate resting on a Winkler foundation. The aim of this section is twofold. First, we hope to provide some analytical (exact or asymptotic) solutions for the growth-induced deformations in hyperelastic plates. To our knowledge, those kind of analytical results are rare in the literature. Second, the current plate theory will be compared with the classical FvK-type plate theory, from which the efficiency of the current plate theory can be revealed. The reference configuration of the plate and the associated Cartesian coordinate system are shown in Fig. 1. Without loss of generality, we take the half-length L = 1, then h just represents the thickness-length ratio of the plate. The external loads q− and q+ are applied on the lower (Z = 0) and the upper (Z = 2h) surfaces of the plate, respectively. Certain boundary conditions will also be proposed on the two ends X = ±1. For simplicity, we shall only consider the case of uniaxial growth (along the X-axis) and suppose that the plate only undergoes plane-strain deformation (along the Y -axis). So, the plate theory in fact becomes a beam theory in these examples. The growth tensor in the plate is denoted as G = diag(λ(X, Z), 1, 1). For this plane-strain problem, the series expansions (11) and (12) can be written as

AC

1 1 1 x = x0 (X) + Zx1 (X) + Z 2 x2 (X) + Z 3 x3 (X) + Z 4 x4 (X) + O(Z 5 ), 2 6 24 y = Y, 1 1 1 z = z0 (X) + Zz1 (X) + Z 2 z2 (X) + Z 3 z3 (X) + Z 4 z4 (X) + O(Z 5 ), 2 6 24 1 1 1 p = p0 (X) + Zp1 (X) + Z 2 p2 (X) + Z 3 p3 (X) + Z 4 p4 (X) + O(Z 5 ). 2 6 24

(44)

Accordingly, the deformation tensors F(n) and A(n) (n = 0, 1, 2, 3) in the 21

CR IP T

ACCEPTED MANUSCRIPT

(a)

(b)

Figure 1: (a) The reference configuration of the thin plate in 3D setting; (b) illustrations of the external loads and the uniaxial growth in the plane-strain problem.

AN US

series expansions (13)1 and (13)2 can be obtained. Furthermore, we consider the following expansion of the growth function 1 1 1 λ = λ0 (X) + Zλ1 (X) + Z 2 λ2 (X) + Z 3 λ3 (X) + Z 4 λ4 (X) + O(Z 5 ), (45) 2 6 24

M

To derive some concrete results, we suppose that the plate is made of incompressible neo-Hookean materials, which has the following elastic strainenergy function   φ0 (A) = C0 tr(AAT ) − 3 . (46)

ED

Then, the nominal stress tensor can be calculated through    ∂R0 (A) ∂φ0 (A) −1 −p = JG G−1 2C0 AT − pA−1 . S = JG G ∂A ∂A

(47)

AC

CE

PT

Based on (44)-(47), the series expansions of S and R0 (A) in terms of Z can be derived. Through the manipulations introduced in section 3.1, we can obtain the expressions of x1 -x3 , z1 -z3 and p0 -p2 in terms of x0 and z0 (cf. Appendix D). By substituting these expressions into the 2D vector plate equation (25), the following two equations for x0 and z0 can be derived  2C0 02 0 0 02 − 3λ40 x00 (x02 0 + z0 )λ0 − x0 (x0 + 2 02 02 3 λ0 (x0 + z0 )  02 00 0 0 00 02 3 00 5 02 λ0 (x02 0 + z0 ) x0 + λ0 ((3x0 − z0 )x0 + 4x0 z0 z0 ) (3) 0 ¯ 1 (x00 , · · · , x(3) hW 0 , z0 , · · · , z0 , λ0 , λ1 ) (4) 0 ¯ 2 (x0 , · · · , x(4) h2 W 0 , z0 , · · · , z0 , λ0 , λ1 , λ2 ) = 0, 0

q¯1 + + + +

22

z002 )3 λ00 (48)

ACCEPTED MANUSCRIPT

 0 2C0 02 4 02 02 2 0 z ((−x02 0 − z0 )(3λ0 + (x0 + z0 ) )λ0 02 3 2 02 λ0 (x0 + z0 ) 0  02 3 00 02 02 4λ50 x00 x000 ) + λ0 (−λ40 (x02 0 − 3z0 ) + (x0 + z0 ) )z0 (3) 0 ¯ 3 (x00 , · · · , x(3) hW 0 , z0 , · · · , z0 , λ0 , λ1 ) (4) 0 ¯ 4 (x00 , · · · , x(4) h2 W 0 , z0 , · · · , z0 , λ0 , λ1 , λ2 ) = 0,

q¯3 + +

CR IP T

+

(49)

x0 (X) = X + U0 (X),

AN US

+ ¯ given where q¯1 and q¯3 are the two components of the effective external load q ¯ 1 -W ¯ 4 are omitted here for brevity (the in (25). The lengthy expressions of W interested readers can contact the authors for the explicit expressions of these terms). In subsections 5.1 and 5.2, the plate equations (48)-(49) will be solved with certain external loads and end boundary conditions. Before that, we shall give some discussions on the relationship between (48)-(49) and the classical FvK-type plate equations (cf. Dervaux and Ben Amar, 2008; Dervaux et al., 2009; Dervaux and Ben Amar, 2010). First, we write z0 (X) = W0 (X),

λ0 = 1 + g1 (X),

(50)

PT

ED

M

where U0 and W0 are respectively the tangential and transverse displacements and g1 is referred as the growth rate. To derive the FvK-type plate equations, some scaling relations need to be adopted (cf. Dervaux et al., 2009): (1) h . W0  1, U0 ∼ W02 ; (2) g1 ∼ W02 . Since the tangential displacement U0 is much smaller that the transverse displacement W0 , the shear load q¯1 should be small and can be neglected. Besides that, the Z-dependence of the growth function is usually eliminated, which implies that λi = 0 (i = 1, 2, · · · ) in (45). With these scalings and simplifications, we substitute (50) into (48)-(49). By neglecting the terms higher than O(h3 , h2 W0 , hW02 , W03 ), the following two equations can be obtained (3)

AC

CE

8C0 (U000 + W00 W000 − g10 ) − 8C0 hW0 = 0, (51)   3 2 q¯3 + 8C0 W00 U000 + U00 W000 + W00 W000 − g10 W00 − g1 W000 2 (52)   8 (3) (4) 2 00 2 0 − 8C0 h W0 + W0 W0 − C0 h W0 = 0. 3 On the other hand, we consider the {11}-component of the second PiolaKirchhoff stress tensor T = SF−T , which is denoted as σ. Still by using (50) and neglecting the high-order terms, we obtain   1 02 0 σ = 8C0 U0 + W0 − g1 − 8C0 hW000 . (53) 2 23

ACCEPTED MANUSCRIPT

By using (53), (51) and (52) can be rewritten as dσ = 0, dX

8 d (4) C0 h2 W0 − (σW00 ) = q¯3 . 3 dX

(54)

AN US

CR IP T

Through some comparisons, it can be found that the two equations given in (54) are just the FvK-type plate equations in the uniaxial growth case (cf. Eqs. (3.1)-(3.2) in Dervaux and Ben Amar (2010)). Thus, here we have shown that with some scalings and simplifications, the plate equations (48) and (49) can be reduced into the classical FvK-type plate equations (in uniaxial growth case). In our future work, the relationship between our current plate theory and the FvK-type plate theory in general growth cases will be further investigated.

ED

M

5.1. Large bending deformation of a plate induced by a gradient growth field In the first example, we consider the large bending deformation of a thin plate induced by a gradient growth field. As shown in Fig. 2a, the two ends of the plate are free and there is no external load applied on the upper or lower surface of the plate. To remove the freedom of rigid body motion, we suppose the position of the middle point on the lower surface is fixed and the horizontal movement of the middle point on the upper surface is restricted. In the uniaxial growth case, the growth field λ(X, Z) in the plate is supposed to be λ(X, Z) = λ0 + λ1 Z, (55)

AC

CE

PT

where λ0 and λ1 are two parameters independent of X. Due to this gradient growth field, the plate will undergo bending deformations. In fact, for this plane strain problem, we can directly formulate the 2D governing system, including the mechanical field equation and the boundary conditions. The exact solution to the governing system can also be derived (cf. Wang et al. (2017) for the derivation procedure). If we use the Cartesian coordinate system (X, Z) for a reference point and the cylindrical coordinate system (r, θ) for a current point, the exact solution can be written as r(X, Z) =

λ0 + Z, λ1

θ(X, Z) = λ1 X,

p(X, Z) = 2C0 ,

(56)

where the geometrical representations of r and θ are shown in Fig. 2b. To test the correctness of our plate theory of growth, we shall also study this problem by using the plate equations (48) and (49) together with the 24

CR IP T

ACCEPTED MANUSCRIPT

(b)

(a)

AN US

Figure 2: (a) The reference configuration of the thin plate with gradient growth field; (b) illustration of the bending deformation of the plate.

0

ED

M

free end boundary conditions, and then make comparisons of the solutions. As the upper and lower surfaces of the plate are traction-free, we have q¯1 = 0 and q¯3 = 0. Then, (48) and (49) can be integrated once to yield that   λ40 2C0 0 − γ1 + x 1 − 02 λ0 0 (x0 + z002 )2  0 2C0 h 02 2 4 02 02 2 + 2 02 x0 (−λ1 (x02 0 + z0 ) (3λ0 + (x0 + z0 ) ) (57) 02 4 λ0 (x0 + z0 )  02 2 00 2 02 02 4 02 02 2 00 + 4λ20 z00 (λ40 + (x02 0 + z0 ) )x0 ) − 2λ0 (x0 − z0 )(λ0 + (x0 + z0 ) )z0 (3) 0 ¯ 5 (x0 , · · · , x(3) + h2 W 0 , z , · · · , z0 , λ0 , λ1 ) = 0, 0

CE

PT

  2C0 0 λ40 − γ2 + z 1 − 02 λ0 0 (x0 + z002 )2  2C0 h 02 2 6 02 02 00 0 0 00 + 2 02 − 3λ40 λ1 z00 (x02 0 + z0 ) − 2λ0 ((x0 − z0 )x0 + 2x0 z0 z0 ) (58) 02 4 λ0 (x0 + z0 )  02 4 2 02 02 2 02 02 00 0 0 00 − λ1 z00 (x02 0 + z0 ) − 2λ0 (x0 + z0 ) ((x0 − z0 )x0 + 2x0 z0 z0 ) (3) 0 ¯ 6 (x00 , · · · , x(3) + h2 W 0 , z0 , · · · , z0 , λ0 , λ1 ) = 0,

AC

¯ 5 and W ¯ 6 are also omitted for brevity. In where the lengthy expressions of W these two equations, γ1 and γ2 represent the average axial and shear stresses on the cross-section (perpendicular to the X-axis) of the plate. As no external loads are applied on the plate, γ1 and γ2 should be zero. On the two ends of the plate, we adopt the traction boundary conditions (27) with q(0) = 0 and ¯ = 0. Through some calculations, we found that (27)2 just have the same q 25

ACCEPTED MANUSCRIPT

CR IP T

expressions as the plate equations, which are satisfied automatically. From (27)1 , the following results can be obtained   2C0 0 λ40 |X=±1 = 0, x 1 − 02 λ0 0 (x0 + z002 )2   (59) 2C0 0 λ40 z 1 − 02 |X=±1 = 0. λ0 0 (x0 + z002 )2 Besides the boundary conditions (59), to fully determine the unknowns x0 and z0 , we also needs to consider the following restrictions to remove the rigid body motion z0 (0) = 0, 1 x0 (0) + (2h)x1 (0) + (2h)2 x2 (0) = 0. 2

AN US

x0 (0) = 0,

(60)

ED

M

To derive the analytical solution, we first consider the leading-order terms (i.e., the terms of O(h0 )) in the plate equations (57) and (58). It can be found that to ensure the vanishing of the leading order terms, x0 and z0 can be chosen as Z x0 (X) = λ0 cos (a(X) + hb(X)) dX + c1 , Z (61) z0 (X) = λ0 sin (a(X) + hb(X)) dX + c2 ,

CE

PT

where c1 and c2 are the integration constants and a(X) and b(X) are two functions to be determined. With the expressions of x0 and z0 given in (61), the boundary conditions (59) are automatically satisfied. By substituting (61) into (57)-(58) and through some calculations, it can be obtained that a(X) = −λ1 X + c3 ,

b(X) = c4 ,

(62)

AC

where c3 and c4 are another two constants. The values of c1 -c4 can be determined by using the restriction equations given in (60), which are given by λ0 (63) c1 = 0, c2 = − , c3 + hc4 = 0. λ1 Then, from (61)-(63), we obtain that x0 (X) =

λ0 sin(λ1 X), λ1

z0 (X) = 26

λ0 (cos(λ1 X) − 1). λ1

(64)

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 3: Configurations of the lower surfaces of the thin plates with the growth parameters λ0 = 1 and λ1 = π/5, π/3, π/2, 2π/3 and π.

AC

CE

PT

ED

M

By comparing (64) with (56), it can be found that the solution obtained from the plate equations are just the exact solution evaluated on the lower surface Z = 0. In fact, by substituting (64) into the expressions of x1 -x3 , z1 -z3 and p0 -p2 (cf. Appendix D), the exact solution given in (56) can be recovered. This recovery of the exact solution in the plane-strain setting provides a validation of the present plate theory of growth. Based on the solution given in (64), the deformations of the thin plate corresponding to the different values of λ0 and λ1 can be described. As an example, we choose λ0 = 1. For some selected values of λ1 , the configurations of the lower surface of the plate are shown in Fig. 3. It can be seen that due to the gradient growth field, the plate undergoes large bending deformations. Corresponding to the different values of λ1 , the configurations of the lower surface are some circular segments with different radii. When λ1 = π, a closed circle is formed. Some further discussions on the bending deformations of a thin plate induced by gradient growth fields can be found in Wang et al. (2017). In our opinion, the Fvk plate theory is not applicable to the current example. First, the scaling relations proposed in the Fvk plate theory does not hold. For example, from Fig. 3, we can find that U0 ∼ W0  h, g1 ∼ W0 ∼ 1 (cf. the definitions of U0 , W0 and g1 given in (50)). These relations are totally different from those proposed in the Fvk-type plate theory. Second, from (64) it can be seen that the strains of the plate can be relatively large (say > 50%). In this case, the generalized Hooke’s law adopted in the Fvktype plate theory does not hold. Third, in the work of Dervaux et al. (2009), 27

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 4: The reference configuration of a thin plate resting on a Winkler foundation and subject to the growth-induced compression.

the Z-dependence of the growth tensor components are eliminated in deriving the FvK-type plate equations. Thus, to study the deformation of a thin plate induced by a gradient growth field, their Fvk-type plate theory is not applicable.

AC

CE

PT

ED

M

5.2. Growth-induced instability of a plate resting on a Winkler foundation In the second example, we study the growth-induced instability of a thin plate resting on a Winkler foundation. The reference configuration of the plate and the associated Cartesian coordinate system are shown in Fig. 4, where the upper surface is denoted as Z = 0 and the lower surface is denoted as Z = 2h. The plate undergoes uniaxial growth along the X-axis and the growth function is assumed to be a constant λ(X, Z) = λ0 . At the two ends of the plate, the horizontal displacements are blocked by some smooth planes. Thus, due to the growth effect, the plate is subjected to axial compression. This example can be viewed as a prototype of many biological tissues (Li et al., 2012) and it has been investigated in many works (Moulton et al., 2013; Budday et al., 2014; Huang et al., 2016). Here, our derived plate (beam) theory of growth will be applied and an analytical post-bifurcation analysis will be provided. The upper surface of the plate is supposed to be traction-free. On the lower surface, the Winkler foundation provides the load q+ = (0, 0, −K0 λ0 W ), where K0 is the elastic constant of the foundation and W is the Z-component of the displacement. λ0 in q+ represents the fact that the load is applied on the grown configuration but not the reference configuration of the plate. 28

ACCEPTED MANUSCRIPT

¯ are given by Then, the components of the effective load q q¯1 = 0, 1 K0 λ0 [z(X, 2h) − 2h] 2h   1 1 1 2 3 = − K0 λ0 z0 + 2hz1 + (2h) z2 + (2h) z3 − 2h . 2h 2 6

(65)

CR IP T

q¯3 = −

By substituting (65) into (48) and (49), the two plate equations can be rewritten as

(3)

AN US

02 00 0 0 00 2C0 x000 λ20 (2C0 (3x02 0 − z0 )x0 + 8C0 x0 z0 z0 ) + 02 3 λ20 (x02 0 + z0 )   1 002 02 0 002 002 +h 8C0 (z003 (x002 0 − z0 ) + 3x0 z0 (−x0 + z0 ) 02 λ0 (x0 + z002 )3 (3)

00 00 0 + x00 z002 (−6x000 z000 + z00 x0 ) + x03 0 (2x0 z0 + z0 x0 ))   λ30 002 04 04 (3) 8C0 (Z003 (x002 + 4C0 (−x0 + z0 )z0 + 02 0 − 3z0 ) 02 5 (x0 + z0 ) (3)

ED

M

0 002 002 + x00 z002 (−12x000 z000 + z00 x0 ) + x02 0 z0 (−7x0 + 5z0 )   03 00 00 0 (3) 04 04 (3) + x0 (4x0 z0 + z0 x0 )) + 4C0 (−x0 + z0 )z0

(66)

AC

CE

PT

(4) 0 ¯ 7 (x0 , · · · , x(4) + h2 W 0 , z0 , · · · , z0 , λ0 ) = 0, 0   02 00 z0 λ0 x00 2C0 z000 2C0 λ20 (4x00 z00 x000 − (x02 0 − 3z0 )z0 ) K0 1 − − 02 + + 02 3 2h x + z 02 λ20 (x02 0 + z0 )   2 0 0 0 00 0 02 λ (3x z x − (x0 − 2z002 )z000 ) z00 (x00 x000 + z00 z000 )  + h K0 0 0 0 002 + 02 (x0 + z002 )3 λ20 (x02 0 + z0 ) 4C0 λ30  03 00 00 0 (3) 0 00 00 04 (3) − 02 − 24x02 0 z0 x0 z0 + x0 x0 + z0 (8x0 z0 − z0 x0 ) (x0 + z002 )5  (67) 002 0 (3) 03 002 002 0 (3) + 2x00 z002 (5x002 − 7z + z z ) + 2x (−3x + z + z z ) 0 0 0 0 0 0 0 0 0  4C0 0 (3) 03 00 00 0 00 00 04 (3) + − 12x02 0 z0 x0 z0 + x0 x0 + z0 (4x0 z0 − z0 x0 ) 02 3 λ0 (x02 + z ) 0 0  0 02 002 002 0 (3) 03 002 002 0 (3) + 2x0 z0 (3x0 − 3z0 + z0 z0 ) + 2x0 (−x0 + z0 + z0 z0 ) (4) 0 ¯ 8 (x00 , · · · , x(4) + h2 W 0 , z0 , · · · , z0 , λ0 ) = 0.

29

ACCEPTED MANUSCRIPT

¯ 7 and W ¯ 8 are also omitted for brevity. In Here, the lengthy expressions of W fact, (66) can be integrated once to yield that h 4C (2x0 z 0 x00 + (−x02 + z 02 )z 00 ) 2C0 x00 2C0 λ20 x00 0 0 0 0 0 0 0 + h − 02 2 02 2 2 02 ) + z λ0 (x0 + z0 ) λ0 (x02 0 0 02 00 i ) )z + z λ30 (8C0 x00 z00 x000 + 4C0 (−x02 0 0 0 + 02 4 ) + z (x02 0 0 (2) (3) (2) (3) 0 2 ¯ 0 + h W3 (x0 , x0 , x0 , z0 , z0 , z0 ) = 0,

CR IP T

− γ0 +

(68)

AN US

where the integration constant γ0 represents the average axial stress, and the ¯ 3 is omitted. lengthy expression of W To conduct further analyses, we need to propose suitable boundary conditions on the two ends of the plate X = −1, 1. As the displacements of the two ends of the plate are blocked along the X-axis, the following position boundary conditions can be proposed x|X=−1,Z=0 = x0 |X=−1 = −1,

x|X=1,Z=0 = x0 |X=1 = 1,

(69)1

ED

M

  Z 2h 1 2h2 h3 x(X, Z)dZ|X=−1 = x0 + hx1 + x2 + x3 |X=−1 = −1, 2h 0 3 3 (69)2   Z 2h 2 3 2h h 1 x(X, Z)dZ|X=1 = x0 + hx1 + x2 + x3 |X=1L = 1. 2h 0 3 3

Besides that, we propose the lubricated boundary conditions on the two ends of the plate, which are given by

PT

(0)

Σ13 |X=−1,Z=0 = Σ13 |X=−1 = 0,

(0)

Σ13 |X=1,Z=0 = Σ13 |X=1 = 0,

AC

CE

  Z 2h 2h2 (2) 1 (0) (1) Σ13 (X, Z)dZ|X=−1 = Σ13 + hΣ13 + Σ |x=−1 = 0, 2h 0 3 13   Z 2h 1 2h2 (2) (0) (1) Σ13 (X, Z)dZ|X=1 = Σ13 + hΣ13 + Σ |x=1 = 0. 2h 0 3 13

(69)3

(69)4

With the expressions of x(n) , z (n) and pn−1 (n ≥ 1) given in Appendix D, the boundary conditions can also be expressed in terms of x0 and z0 . 5.2.1. Linear bifurcation analysis First, we conduct a linear bifurcation analysis on the plate equations (67)(68) subject to the boundary conditions (69). For convenience, we take the 30

ACCEPTED MANUSCRIPT

transformation K0 = α1 C0 . The homogeneous deformation of the plate can be represented by x0 (X) = X, z0 (X) = E0 . (70)

γ0 = −

CR IP T

Substituting (70) into the expressions of x(n) , z (n) and pn−1 (n ≥ 1), it can be found that the only nonzero quantities are z1 = λ0 and p0 = 2C0 λ20 . It can also be verified that the homogeneous solutions (70) satisfy the boundary conditions given in (69). By substituting (70) into the plate equations (67) and (68), we obtain 2C0 (λ40 − 1) , λ20

E0 = −2h(λ0 − 1).

(71)

z0 (X) = −2h(λ0 − 1) + ∆W (X),

(72)

M

x0 (X) = X + ∆U (X),

AN US

In fact, E0 just ensures that the Z-component of the displacement of the lower surface equals zero, i.e., z0 + 2hz1 − 2h = 0. In addition to the homogeneous solutions, the plate equations may also have nontrivial solutions. Next, we conduct a linear bifurcation analysis to determine the critical parameter values at which the nontrivial solutions exist. For that purpose, we set

AC

CE

PT

ED

where  is a small quantity identifying the order of variations. Substituting (72) into (67)-(68) and dropping the nonlinear terms, we obtain   4h(1 + λ40 ) 1 4 2 2 + 6λ40 0 4 00 ∆U − ∆W + h −3 − 4 − 4λ0 ∆U (3) = 0, λ20 λ0 3 λ0 (73) 4 α1 ∆W (2 − (2 + hα1 )λ0 ) − + α1 λ0 ∆U 0 + ∆W 00 2h λ20 (74) 2h(6 + hα1 + 2(3 + hα1 )λ40 ) 4 2 (3) 4 (4) − ∆U + h (3 + λ0 )∆W = 0. 3λ0 3 From (73), it can be obtained that ∆U 0 =

2hλ0 (1 + λ40 ) ∆W 00 + O(h3 ). 1 + 3λ40

(75)

As the terms of O(h3 ) have been dropped in deriving the plate equations (67) and (68), there is no need to keep the term of O(h3 ) in (75). By substituting (75) into (74), the following single equation for ∆W can be obtained ψ0 ∆W + ψ2 ∆W 00 + ψ4 ∆W (4) = 0, 31

(76)

ACCEPTED MANUSCRIPT

where α1 , 2h (λ4 − 1) [2 + (6 + hα1 )λ40 ] , ψ2 = 0 λ20 + 3λ60  4h2  ψ4 = 3 + hα1 + (2 + 3hα1 )λ40 + (3 + 2hα1 )λ80 . 4 3 + 9λ0

CR IP T

ψ0 =

AN US

Equation (76) is a fourth-order ODE for ∆W , which requires four boundary conditions to be proposed at X = ±1. From (69) and by using (75), we can propose the following boundary conditions ∆W 0 (−1) = ∆W 0 (1) = 0,

∆W (3) (−1) = ∆W (3) (1) = 0.

(77)

The characteristic equation of (76) is

ψ0 + ψ2 η 2 + ψ4 η 4 = 0.

(78)

M

Through a conventional linear bifurcation analysis, it can be found that the nontrivial solutions of the ODE system (76)-(77) can only exist when (78) has the purely imaginary roots nπ i, 2

ED

η=

n = 1, 2, 3, · · · ,

(79)

CE

PT

where n is denoted as the mode number and i is the imaginary unit. For n = 2k and 2k − 1, the nontrivial solutions are given by ( ∆W (X) = Acos(kπX), if n = 2k, (80) ∆W (X) = Asin((k − 1/2)πX), if n = 2k − 1, k = 1, 2, 3, · · · .

AC

On the other hand, by substituting (79) into (78), we obtain the following equation involving the mode number n and the other parameters ψ0 −

ψ2 2 2 ψ4 4 4 n π + n π = 0, 4 16

(81)

By specifying the value of α1 in (81), we can obtain the relationship between λ0 and h for each mode number n. For example, when α1 = 0.2, the curves corresponding to these relationships (1 ≤ n ≤ 17) are shown in Fig. 5. 32

AN US

CR IP T

ACCEPTED MANUSCRIPT

(a)

(b)

Figure 5: The bifurcation curves corresponding to the modes n = 1, · · · , 17: (a) α1 = 0.2; (b) enlarged part of (a) near λ1 = 1.0 and h = 0.

AC

CE

PT

ED

M

It can be seen that for each mode n, the relationship between λ0 and h is non-monotonic. If the thickness-length ratio h is relatively large (say, h > 0.1), the mode n = 1 or 2 will be the first bifurcation mode and the corresponding bifurcation parameter λ0 has a large value (cf. Fig. 5a). If the thickness-length ratio h becomes smaller, the first bifurcation mode number increases and the corresponding value of λ0 decreases (cf. Fig. 5b). Thus, the nontrivial solutions corresponding to a large mode number is more possible to appear in thinner plates. Furthermore, we analyze the dependence of the critical mode nc (i.e., the mode with the minimum value of λ0 ) on the thickness-length ratio h. For that purpose, we assume that the definition domain of nc can be extended to the field of real numbers. For specified values of α1 and h, (81) just provides an implicit relationship between λ0 and n, from which we can define the function λ0 (n). By solving the equation ∂λ0 (n)/∂n = 0, we obtain the following relationship between the critical mode nc and the minimum bifurcation value (m) λ0 : v u (m) 4 (m) 4 (m) 8 3(2 + (6 + hα1 )λ0 )(−1 − 2λ0 + 3λ0 ) 1u t . nc = π 2h2 (λ(m) 2 + 3λ(m) 6 )(3 + hα + (2 + 3hα )λ(m) 4 + (3 + 2hα )λ(m) 8 ) 0

0

1

1

1

0

(m)

By substituting (82) into (81), we obtain the relationship between λ0 33

0

(82) and h.

CR IP T

ACCEPTED MANUSCRIPT

(b)

AN US

(a)

(m)

Figure 6: Dependencies of the minimum bifurcation parameter λ1 nc on the thickness-length ratio h (α1 = 0.2, α2 = 0).

and the critical mode

ED

M

Further from (82), the dependence of nc on h can also be determined. In Fig. (m) 6, we plot the h-λ0 and h-nc curves corresponding to the critical modes. From Fig. 6, it can be found that if h tends to zero, the minimum bifurcation (m) parameter λ0 tends to 1.0, while the critical mode number nc increases rapidly. This result further reveals that for thinner plates the nontrivial states corresponding to large mode numbers are more possible to appear.

AC

CE

PT

5.2.2. Post-bifurcation analysis To reveal the post-bifurcation behavior of the plate subject to growthinduced compression, we further conduct a nonlinear post-bifurcation analysis and construct the approximate analytical solutions to the plate equation system. Compared with the linear bifurcation analysis, the major advantage of the post-bifurcation analysis is that the amplitudes of the leading-order nontrivial solutions can be determined from the amplitude equations and then the bifurcation types can be classified. For the convenience of the readers, we put the lengthy derivation procedure of the amplitude equations in Appendix E and only introduce the key results here. First, to proceed with the perturbation method, we introduce a small parameter  and seek the solutions in the following perturbation expansion

34

ACCEPTED MANUSCRIPT

forms:

(2)

(3)

(83)

CR IP T

λ0 = λc + δλ, δλ = 2 ∆λ, x0 = xh + U1 + 2 U2 + 3 U3 + · · · , z0 = zh + W1 + 2 W2 + 3 W3 + · · · , γ0 = γh + 2 C0 γ0 + 3 C0 γ0 + · · · ,

ED

M

AN US

where xh , zh and γh are the homogeneous solutions to the plate equation system (which have been given in (70) and (71)) and λc is the bifurcation value of λ0 corresponding to any given mode. The above perturbation expansions imply thatpthe amplitude of x0 − xh is of the same order as z0 − zh and is of the order |λ0 − λc |. By substituting (83) into (67) and (68) and equating the coefficients of i  (i = 1, 2, 3) to be zero, we can obtain three ODE systems for the unknowns U1 -U3 and W1 -W3 (cf. (E1)-(E3) in Appendix E). The associated boundary conditions can also be derived from (69) (cf. (E4)). Next, we need to derive the post-bifurcation solutions from these ODE systems. At O(1 ), the ODE system is in fact the same as the linearized system (73)-(74). Thus, through some conventional manipulations, the following leading-order nontrivial solutions W1 and U10 can be obtained ( W1 (X) = 2R1 cos(ξ1 X + θ1 ) + 2R2 cos(ξ2 X + θ2 ), (84) U10 (X) = 2f (ξ1 )R1 cos(ξ1 X + θ1 ) + 2f (ξ2 )R2 cos(ξ2 X + θ2 ),

CE

PT

where R1 and R2 are the amplitude parameters to be determined, f is a function depending on h and λc (see (E6)), and ξi and θi (i = 1, 2) are some constants that can be determined from the boundary conditions (see (E7) and (E8)). At O(2 ), based on the obtained ODE system and through some calculations, we can derive the following asymptotic expressions of W2 and U20

AC

¯ 1 cos(ξ1 X + θ¯1 ) + 2R ¯ 2 cos(ξ2 X + θ¯2 ) + Hλ ∆λ + Hγ γ0(2) W2 =2R + 2H1 R12 cos(2(ξ1 X + θ1 )) + 2H2 R22 cos(2(ξ2 X + θ2 )) + 2H3 R1 R2 cos((ξ1 + ξ2 )X + θ1 + θ2 ) + 2H4 R1 R2 cos((ξ1 − ξ2 )X + θ1 − θ2 ) + H5 R12 + H6 R22 ,

35

(85)

ACCEPTED MANUSCRIPT

¯ 1 cos(ξ1 X + θ¯1 ) + 2f (ξ2 )R ¯ 2 cos(ξ2 X + θ¯2 ) + Lλ ∆λ U20 =2f (ξ1 )R (2)

(86)

CR IP T

+ Lγ γ0 + 2L1 R12 cos(2(ξ1 X + θ1 )) + 2L2 R22 cos(2(ξ2 X + θ2 )) + 2L3 R1 R2 cos((ξ1 + ξ2 )X + θ1 + θ2 ) + 2L4 R1 R2 cos((ξ1 − ξ2 )X + θ1 − θ2 ) + L5 R12 + L6 R22 ,

M

AN US

¯1, R ¯ 2 , θ¯1 , θ¯2 and γ0(2) are some unknown constants to be determined, where R and Hi and Li (i = 1, · · · , 6, λ, γ) are some constants depending on h, α1 , λc , ξ1 and ξ2 . Further from the boundary conditions, the values of θ¯1 , θ¯2 and (2) γ0 can also be determined (see (E11)-(E13)). At O(3 ), from the ODE system, we can first derive the asymptotic expression of U30 in terms of W3 , and then obtain one single equation for W3 (see (E14)). By further considering the solvability condition of this equation subject to the boundary conditions, we can obtain the following amplitude equations for the leading-order nontrivial solutions i h (2) R1 sλ ∆λ + sγ γ0 + s1 R12 + s2 R22 = 0, h i (87) (2) 2 2 R2 tλ ∆λ + tγ γ0 + t1 R2 + t2 R1 = 0.

PT

ED

To ensure the existence of nontrivial solutions, at least one of R1 and R2 should be nonzero. In the case R1 6= 0 and R2 = 0, from (87)1 and the (2) expression of γ0 , the following result for the amplitude parameter R1 can be obtained sλ L γ − sγ L λ . (88) R12 = r1 ∆λ, r1 = sγ L 5 − s1 L γ

CE

With the expressions of R1 , θ1 and ξ1 , the leading-order nontrivial solutions can be obtained as ( z0 =zh + 2 (r1 δλ)1/2 cos(k1 πX + k¯1 π), (89)1 x0 =x0 + 2f (k1 π) (r1 δλ)1/2 cos(k1 πX + k¯1 π), 0

h

z0 =zh + 2 (r1 δλ)1/2 cos((k1 − 1/2)πX + (k¯1 − 1/2)π),

AC

(

x00 =x0h + 2f ((k1 − 1/2)π) (r1 δλ)1/2 cos((k1 − 1/2)πX + (k¯1 − 1/2)π), (89)2 which are the analytical expressions for describing the post-bifurcation states of the plate. Similarly, when R1 = 0 and R2 6= 0, the approximate solutions of z0 and x00 can also be obtained. 36

ACCEPTED MANUSCRIPT

In the case that both R1 and R2 are not equal to zero, by solving the two equations in (87), we obtain R22 = r2 ∆λ,

(90)

CR IP T

R12 = r1 ∆λ, where

Lλ (sγ t1 − s2 tγ ) + Lγ (−sλ t1 + s2 tλ ) + L6 (sλ tγ − sγ tλ ) , Lγ (s1 t1 − s2 t2 ) + L6 (sγ t2 − s1 tγ ) + L5 (−sγ t1 + s2 tγ ) Lλ (−sγ t2 + s1 tγ ) + Lγ (sλ t2 − s1 tλ ) + L5 (−sλ tγ + sγ tλ ) . r2 = Lγ (s1 t1 − s2 t2 ) + L6 (sγ t2 − s1 tγ ) + L5 (−sγ t1 + s2 tγ )

r1 =

(91)

M

AN US

With the expressions of Ri , θi and ξi (i = 1, 2), we can obtain the expressions of the leading-order nontrivial solutions z0 and x00 . For example, the solutions corresponding to ξ1 = k1 π and ξ2 = k2 π (cf. (E7) and (E8) for the values of ξ1 and ξ2 ) are given by  z0 =zh + 2 (r1 δλ)1/2 cos(k1 πX + k¯1 π)      + 2 (r2 δλ)1/2 cos(k2 πX + k¯2 π), (92) 1/2   x00 =x0h + 2f (k1 π) (r1 δλ) cos(k1 πX + k¯1 π)    + 2f (k2 π) (r2 δλ)1/2 cos(k2 πX + k¯2 π).

AC

CE

PT

ED

Next, we give some further analyses on the properties of the post-bifurcation solutions. As an example, we choose h = 0.015 and α1 = 0.2. It has been found that with these parameters the post-bifurcation solutions only exist corresponding to the mode numbers 1 ≤ n ≤ 29. For the case R1 6= 0 and R2 = 0, we show the critical bifurcation parameters λc of the different modes in Fig. 7a. It can be seen that for the first mode n = 1, the value of λc is relatively large. The minimum value of λc is attained at the mode n = 7, which is consistent with the result shown in Fig. 5b. By using (88), the amplitudes of the post-bifurcation solutions can be defined by ρ = ± max |z − zh | = ±2 (r1 δλ)1/2 .

(93)

For the purpose of illustration, we plot the bifurcation diagrams of modes n = 6, 7 in Fig. 7b, where the λ0 -axis represents the trivial solution. It can be seen that both of the two modes exhibit the supercritical bifurcation property. In fact, we have found that the modes 1 ≤ n ≤ 23 have the supercritical bifurcation property and the modes 24 ≤ n ≤ 29 have the subcritical 37

ACCEPTED MANUSCRIPT

  

 

CR IP T

 0





 

    

   

  



 

(b)

(a)

AN US

Figure 7: (a) Relationship between the bifurcation modes and the critical bifurcation values of λ0 ; (b) The bifurcation diagrams for modes n = 6 and 7 (the parameters are h = 0.015 and α1 = 0.2).

ψ0 −

ED

M

bifurcation property. It appears that the existence of subcritical bifurcation modes for this type of problems has not been reported in literature. If both R1 and R2 are not equal to zero, we need to require that (81) is satisfied by two distinct modes n1 and n2 . In this case, the two modes have the same growth bifurcation parameter λc . This case can only happen for some special values of h and α1 such that the following two equations are satisfied simultaneously ψ2 2 2 ψ4 4 4 n π + n1 π = 0, 4 1 16

ψ0 −

ψ2 2 2 ψ4 4 4 n π + n2 π = 0. 4 2 16

(94)

CE

PT

Besides that, there exist another two post-bifurcation solutions corresponding to R1 6= 0 and R2 = 0, and R1 = 0 and R2 6= 0. For single-mode solutions, the definition of the amplitude has been given in (93), while for two-mode solutions, the amplitude can be defined by

AC

ρ = ± max |z − zh | = ±(2 (r1 δλ)1/2 + 2 (r2 δλ)1/2 ).

(95)

For the purpose of illustration, we consider the two modes n1 = 6 and n2 = 7. The parameter α1 is set to be α1 = 0.2. By solving the two equations in (94), we get that h = 0.015072 and λc = 1.016319. The evolutions of the amplitudes ρ accompanying the variation of λ1 for both the single-mode and the two-mode solutions are shown in Fig. 8. It can be seen that all of these solutions have the supercritical bifurcation property. 38

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 8: The bifurcation diagrams of the single-mode and the two-mode post-bifurcation solutions with the mode numbers n1 = 6 and n2 = 7, where α1 = 0.2, h = 0.015072 and λc = 1.016319.

AC

CE

PT

ED

M

5.2.3. Numerical verification To verify the analytical results obtained from the post-bifurcation analysis, we shall derive the numerical solutions to the plate equation (66) and (67) in this subsection. From the derivation procedure of the amplitude equations (cf. Appendix E), it can be seen that the terms higher than O(h) have been dropped in Ui0 (i = 1, 2, 3). This manipulation implies that the term of O(h2 ) in (66) (or (68)) are not taken into account in the analysis (it has also been pointed out in Appendix C that the term of O(h2 ) in (66) should be dropped to ensure the consistency of the plate equations). Thus, in the numerical calculations, the term of O(h2 ) in (66) should also be dropped. However, it was found that without this term the Jacobian matrix in the numerical iteration procedure becomes singular matrix. To overcome this problem, we multiply a small quantity d0 = 0.1 onto this term in the numerical calculations. The ODE package ‘bvp6c’ in Matlab is used to solve the system (66)-(67) subject to the boundary conditions (69). To use this package, (66) and (67) need to be written into a first-order ODE system

where

dY = f (X, Y), dX  T (3) (3) Y = (y1 , y2 · · · , y8 )T = x0 , x00 , x000 , x0 , z0 , z00 , z000 , z0 . 39

(96)

0.015

3

0.01

2

0.005

1

0

0

-0.005

-1

-0.01

-2

-0.015 -1

-0.5

0

0.5

×10 -3

-3 -1

1

X

-0.5

CR IP T

x(X)-X

z(X)

ACCEPTED MANUSCRIPT

0

0.5

1

X

(a)

(b)

AN US

Figure 9: Comparisons of the numerical solutions (curves) and the analytical solutions (dots) of the plate equation system with the mode number n = 7 (h = 0.015, α1 = 0.2, λc = 1.01626, δλ = 0.005): (a) z(X); (b) U (X) = x(X) − X.

M

The boundary conditions (69) can be rewritten as ( y1 (−1) = −1, y1 (1) = 1, y3 (−1) = 0, y3 (1) = 0, y6 (−1) = 0, y6 (1) = 0, y8 (−1) = 0, y8 (1) = 0.

(97)

AC

CE

PT

ED

To derive some concrete numerical results, we still choose the parameters h = 0.015 and α1 = 0.2. For simplicity, only the single-mode post-bifurcation solutions will be considered. The analytical solutions (89) are adopted as the initial guesses for the numerical calculations. As an example, we choose the mode number n = 7, i.e., k1 = 3 in (89)2 . The critical bifurcation value of this mode is λc = 1.01626, which is the minimum bifurcation value among the different modes. By setting the incremental growth parameter δλ = 0.005 and solving the ODE system (96)-(97) with the ‘bvp6c’ package, we obtain the numerical values of x0 and z0 during the interval −1 ≤ X ≤ 1, which are shown in Fig. 9. The analytical results of x0 and z0 are also plotted in Fig. 9 for comparison. It can be seen that there are very good agreements. The bifurcation diagrams obtained from the numerical and the analytical solutions are plotted in Fig. 10, which also show very good agreements. The results show in Fig. 9 and Fig. 10 demonstrate the validity of the analytic results obtained from the post-bifurcation analysis. As mentioned before, to derive the numerical solutions, the term of O(h2 ) in (66) is not dropped but is multiplied with a small quantity d0 = 0.1. To justify this manipulation, we need to verify that addition of this term has no 40

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 10: The bifurcation diagrams obtained from the numerical solutions (dots) and the analytical solutions (curves) for mode n = 7 (h = 0.015, α1 = 0.2, α2 = 0, λc = 1.01626).

ED

M

much influence on the solutions. For that purpose, we show that the value of the term of O(h2 ) (with the multiplication of d0 ) is much smaller than the leading-order term in (66). The comparison for the solution in Fig. 9 is shown in Fig. 11. It can be seen that the term of O(h2 ) is indeed much smaller than the leading-order term. Thus, one can expect that, keeping the last term in (49) with a multiplication of a factor 0.1 for the ease of numerical computations does not make much difference to the case of dropping it.

AC

CE

PT

5.2.4. Comparisons with the FvK plate theory From the results shown in Figs. 7 and 9, it can be seen that the scaling relations adopted in the FvK plate theory, e.g., h . W0  1, U0 ∼ W02 and g1 ∼ W02 , are approximately satisfied. Thus, the current example can also be solved by using the FvK-type plate equations (54). Here, we aim to compare the results obtained from the different plate theories and reveal their consistency. ¯ into (54), the FvKBy substituting the specific forms of λ(X, Z) and q type plate equations can be rewritten as (3)

8C0 U000 + 8C0 W00 W000 − 8C0 hW0

= 0,

(98)1

α1 W 0 + α1 (U00 + W002 − g1 ) + 8W00 U000 + 4(−2g1 + 2U00 + 3W002 )W000 2h (98)2 h i 16 (3) (3) (4) + h −W000 (α1 + 16W000 ) − 8(U0 + 2W00 W0 ) + h2 W0 = 0, 3 −

41

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 11: Comparison of the leading-order term (the solid curve) and the term of O(h2 ) (the dashed curve) in (66).

M

where we have made the transformations λ0 = 1 + g1 and K0 = α1 C0 . In fact, (98)1 can be integrated once to yield (53). On the two ends of the plate, we can still propose the boundary conditions (69), which can be simplified into (3)

(3)

ED

U0 (−1) = U0 (1) = 0, W00 (−1) = W00 (1) = 0, W0 (−1) = W0 (1) = 0, (99) The plate equation system (98)-(99) have the following homogeneous solutions U0 (X) = 0, W0 (X) = −2hg1 , σ = −8C0 g1 . (100)

PT

It can be checked that these solutions are same as those given in (70) and (71). For linear bifurcation analysis, we take

CE

U0 (X) = ∆U (X),

W0 (X) = −2hg1 + ∆W (X).

(101)

AC

Substituting (101) into (98) and dropping the nonlinear terms (where (98)1 needs to be integrated once), we obtain



8∆U 0 − 8h∆W 00 = 0,

(102)1

16 α1 ∆W + α1 ∆U 0 − (8g1 + hα1 )∆W 00 − 8h∆U (3) + h2 ∆W (4) = 0. 2h 3 (102)2 42

ACCEPTED MANUSCRIPT

Table 1: Comparisons of the critical bifurcation values of g1 obtained from the two plate theories (with parameters h = 0.015 and α1 = 0.2). n=7

n=6

n=8

n=5

n=9

FvK plate theory

0.01596

0.01604

0.01712

0.01814

0.01916

Current plate theory

0.01626

0.01630

0.01751

0.01839

0.01969

CR IP T

Mode number

From (102)1 , we get that ∆U 0 = h∆W 00 . Substituting this expression into (102)2 , we have

AN US

8 α1 ∆W + 8g1 ∆W 00 + h2 ∆W (4) = 0. 2h 3

(103)

AC

CE

PT

ED

M

The bifurcation property of (103) can be discussed through the conventional approach. It can be found that the nontrivial solutions just have the same forms as those given in (80), where the mode number n satisfies the following equation 8 α1 − 8g1 n2 π 2 + h2 n4 π 4 = 0. (104) 2h 3 Based on (104) and with given values of h and α1 , the critical bifurcation value of g1 corresponding to the different mode number n can be determined. As an example, we choose h = 0.015 and α1 = 0.2. In Table 1, the bifurcation values of g1 for the first five modes (i.e., the modes with the minimum bifurcation values of g1 ) are listed. For the purpose of comparison, we also list the critical values of g1 = λ0 − 1 obtained from our consistent plate theory. It can be seen that the results obtained from the two plate theories are very close to each other, which reveals the consistency of the two plate theories in studying the current problem. Furthermore, the plate equation system (98)-(99) are also solved numerically by using the package bvp6c in Matlab. The procedure to derive the numerical solutions are similar as that introduced in subsection 5.2.3. Especially, the analytical solutions (89) are still adopted as the initial guesses to generate the nontrivial solutions. As an example, we choose h = 0.015 and α1 = 0.2. Then, for each mode number n, the numerical solutions to system (98)-(99) can be obtained, which are similar as those shown in Fig. 9. Besides that, accompanying the variation of g1 , the bifurcation diagrams of the nontrivial solutions can also be plotted. In Fig. 12, we choose the mode n = 7 and plot the bifurcation diagram obtained from the numerical 43

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 12: The bifurcation diagrams obtained from the numerical solutions of the FvK plate theory (dots) and the analytical solutions obtained in subsection 5.2.2 (curves) for mode n = 7 (h = 0.015, α1 = 0.2).

M

solutions of (98)-(99), which is compared with the diagram of the analytical solutions obtained in subsection 5.2.2. It can also be found that the results obtained from the two plate theories are close to each other. 6. Conclusions

AC

CE

PT

ED

In this paper, we developed a novel plate theory for large deformations of thin hyperelastic plates induced by growth. This plate theory was established within the framework of nonlinear elasticity, taking into account the growth effects and the elastic incompressibility. Furthermore, it was verified that the plate equations of this theory ensure the asymptotic order of O(h4 ) for all the terms in the variations of the total energy functional. For future numerical calculations, the weak formulation of the plate equations was also derived. As typical applications of this plate theory, we studied two examples in uniaxial growth case. One is the large bending deformation of a thin plate induced by a gradient growth field, and the other is the growth-induced instability of a thin plate resting on a Winkler foundation. In the first example, the analytical solution to the plate equations was derived, which actually recovers the exact solution of the original governing system in the plane strain setting. In the second example, both the linear bifurcation and the post bifurcation analyses were conducted, based on which the critical bifurcation points of the nontrivial solutions and the amplitudes of the leading-order 44

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

CR IP T

nontrivial solutions can be derived. The influences of the different parameters on the bifurcation properties of the plate equations were also studied. The analytical solutions obtained from the post bifurcation analysis were further compared with the numerical solutions, from which good agreements were found. The current plate theory was also compared with the classical FvK plate theory to reveal its advantages. In the uniaxial growth case, we showed that the FvK-type plate equations can be derived from the current plate theory by adopting some extra scalings and assumptions. For the first example studied in section 5, the Fvk-type plate theory is not applicable, because the scaling relations cannot be satisfied and the strain values exceed the small strain range. For the second example, based on the FvK-type plate equations, the bifurcation points of the nontrivial solutions with different modes were determined and the post-bifurcation solutions were derived through numerical calculations. It was found that the results obtained from the two plate theories are very close to each other. These investigations demonstrate the advantages of our current plate theory. For example, the finite-strain, bending and stretching effects are all taken into consideration and no ad hoc hypotheses are involved. It would be expected that the current plate theory of growth could have broad applications. In our future work, the current plate theory will be applied to study the growth-induced pattern formations in hyperelastic plates with other loading and boundary conditions. The comparisons of the current plate theory and other classical plate theories will also be made systematically. Furthermore, a consistent finite-strain shell theory of growth will be developed to investigate the mechanical behaviors of soft biological tissues with other geometrical shapes.

CE

Acknowledgements

AC

This paper was supported by a GRF grant (Project No.: CityU 11303015) from the Research Grants Council of Hong Kong SAR, China. Jiong Wang is supported by a grant from the National Nature Science Foundation of China (Project No.: 11602086) and the Guangdong Nature Science Funds for Distinguished Young Scholar (Project No.: 2015A030306009). David Steigmann gratefully acknowledges the support of the US NSF through grant CMMI 1538228.

45

ACCEPTED MANUSCRIPT

Appendix A. The explicit expressions of R(i)

CR IP T

In the current work, the constraint function R(A) = Det(A) − 1 plays the similar role as the strain energy function φ(A). Thus, we need to consider the derivatives of R(A) with respect to A. By using the formula   d dA −1 (Det(A)) = Det(A)tr A , (A1) dτ dτ the following results can be obtained

∂R(A) = Det(A)A−1 = Det(A)A−1 ab ea ⊗ eb , ∂A  ∂ 2 R(A) −1 −1 −1 = Det(A) A−1 R(1) (A) = ea ⊗ eb ⊗ ec ⊗ ed , ab Acd − Aad Acb 2 ∂A ∂ 3 R(A) (A2) R(2) (A) = ∂A3  −1 −1 −1 −1 −1 −1 −1 −1 =Det(A) A−1 ab Acd Aef − Aaf Acd Aeb − Aab Acf Aed −1 −1 −1 −1 −1 −1 −1 −1 − A−1 ad Acb Aef + Aaf Acb Aed − Aad Acf Aeb

AN US

R(0) (A) =

M

× ea ⊗ eb ⊗ ec ⊗ ed ⊗ ee ⊗ ef ,

ED

For incompressible materials, we have Det(A) = 1, thus the results given in (A2) can be further simplified.

PT

Appendix B. Expressions of x(1) and p(0) for incompressible neoHookean materials

AC

CE

In this appendix, based on (19) and (21)1 , we shall derive the explicit expressions of x(1) and p(0) in terms of x(0) for incompressible neo-Hookean materials. First, we choose two unit normal vectors e1 and e2 from the plate plane, such that {e1 , e2 , k} form an orthonormal basis. One material point in the reference and the current configurations of the plate are denoted as X = (r, Z) = (X1 , X2 , Z) and x = (x1 , x2 , x3 ), respectively. With respect to the basis {e1 , e2 , k}, the in-plane deformation gradient (0) ∇x can be written as ∇x(0) =

∂xi ei ⊗ ej , ∂Xj

i = 1, 2, 3, j = 1, 2,

46

(B1)

ACCEPTED MANUSCRIPT

AN US

CR IP T

where e3 = k. For convenience, we denote (∇x(0) )∗ as     ∂x3 ∂x2 ∂x3 ∂x1 ∂x1 ∂x3 ∂x2 ∂x3 (0) ∗ − e1 − − e2 (∇x ) = ∂X1 ∂X2 ∂X1 ∂X2 ∂X1 ∂X2 ∂X1 ∂X2   ∂x1 ∂x2 ∂x2 ∂x1 + − k. ∂X1 ∂X2 ∂X1 ∂X2 (B2) The constitutive relations of the incompressible neo-Hookean materials have been given in (46) and (47). By substituting these constitutive relations into (19) and (21)1 , and considering the expression of (∇x(0) )∗ , the following two equations can be obtained (0) T 2 (1) ¯ ) k x = −q− − 2C0 ∇x(0) G ¯ (0) (G ˆ (0) )T k + p(0) (∇x(0) )∗ , 2C0 J (0) (G (B3) ¯ (0) )−1 , x(1) · (∇x(0) )∗ = Det(G (B4) where J (0) = JG |Z=0 . From (B3) and (B4), it is easy to obtain that

(0) T 2 ¯ ) k 2C0 J (0) (G = ¯ (0) ) |(∇x(0) )∗ |2 Det(G   (∇x(0) )∗ − (0) ¯ (0) ˆ (0) T + q + 2C0 ∇x G (G ) k 2. |(∇x(0) )∗ |

PT

ED

p(0)

¯ (0) (G ˆ (0) )T k + p(0) (∇x(0) )∗ −q− − 2C0 ∇x(0) G , ¯ (0) )T k 2 2C0 J (0) (G

(B5)

M

x(1) =

(B6)

Appendix C. Discussion on the asymptotic order of the 2D vector plate equation

AC

CE

Formally, the 2D vector plate equation (25) attains the accuracy of O(h2 ). However, some further analyses need to be conducted to check the accuracy and to ensure that a consistent equation system has been derived. First, equation (25) can be rewritten in the following component form ! (2) (1) X ∂Σ(0) ∂Σ 2 ∂Σ αi + h αi + h2 αi = −¯ qi , i = 1, 2, ∂X ∂Xα 3 ∂Xα α α=1,2 ! (C1) (1) (2) X ∂Σ(0) ∂Σ 2 ∂Σ α3 + h α3 + h2 α3 = −¯ q3 , ∂X ∂Xα 3 ∂Xα α α=1,2 47

ACCEPTED MANUSCRIPT

Besides that, the component form of (19) is given by (0)

Σ3i = −qi− ,

i = 1, 2, 3. (0)

(C2) (0)

CR IP T

In the framework of linear elasticity, we have Σαi = Σiα , thus (C1)2 is in fact equivalent to ! ! X ∂q − X ∂Σ(1) 2 ∂Σ(1) 1 α3 α + h α3 = − q¯3 + , (C3) ∂X 3 ∂X h ∂X α α α α=1,2 α=1,2

M

where I3 = e3 ⊗ e3 .

AN US

It can be seen that (C3) only provides a plate equation with O(h) in the case of linear elasticity. Thus, the term of O(h2 ) in (C1)1 should also be dropped to ensure the consistency of the plate equations. In summary, we have the following plate system of order O(h) in the case of linear elasticity   2 2 (2) (0) (1) ∇ · S + hS + h S I3 = −¯ q, (C4) 3

Appendix D. Some iterative relations

AC

CE

PT

ED

In section 5, the constitutive relations of incompressible neo-Hookean materials are adopted. Besides that, we consider the case of uniaxial growth and suppose that the plate only undergoes plane-strain deformation. In this case, the iterative expressions of x1 -x3 , z1 -z3 and p0 -p2 in terms of x0 and z0

48

ACCEPTED MANUSCRIPT

are given by λ0 x00 2C0 λ20 λ0 z00 , z = , p = , 1 0 02 02 02 x02 x02 x02 0 + z0 0 + z0 0 + z0  2 0 02 1 z1 λ0 (p0 x0 + 2C0 (−2λ1 x00 + λ0 x01 )z00 ) x2 = 02 ) + z 2C0 λ30 (x02 0 0 02 0 0 0 0 − x1 λ20 (2C0 λ1 (x02 0 − z0 ) + z0 (p0 x0 + 2C0 λ0 z1 ))  0 00 0 00 02 0 0 + 2C0 x00 (x02 0 λ0 + z0 λ0 − λ0 x0 x0 − λ0 z0 z0 ) ,  2 0 1 02 z2 = z1 λ0 (x0 (−2C0 λ0 x01 + p00 z00 ) + 2C0 λ1 (x02 0 − z0 )) 02 3 02 2C0 λ0 (x0 + z0 ) + x1 λ20 (−z00 (4C0 λ1 x00 + p00 z00 ) + 2C0 λ0 x00 z10 )  0 02 0 0 00 0 00 + 2C0 z00 (x02 0 λ0 + z0 λ0 − λ0 x0 x0 − λ0 z0 z0 ) ,  1 p1 = z1 λ0 (4C0 λ1 x00 − 2C0 λ0 x01 + p00 z00 ) 02 λ0 (x02 + z ) 0 0  + x1 λ0 (p00 x00 − 4C0 λ1 z00 + 2C0 λ0 z10 ) − 2C0 z00 x000 + 2C0 x00 z000 , 

M

1 3 02 3 0 02 − 4C0 x2 λ30 λ1 x02 0 − 2C0 x1 λ0 λ2 x0 + z2 λ0 p0 x0 4 02 2C0 λ0 (x0 + z002 ) 3 0 0 2 2 0 0 3 0 0 + z1 λ30 p01 x02 0 − 8C0 z2 λ0 λ1 x0 z0 + 4C0 z1 λ0 λ1 x0 z0 − 4C0 z1 λ0 λ2 x0 z0 − x2 λ30 p00 x00 z00 − x1 λ30 p01 x00 z00 + 4C0 z2 λ40 x01 z00 − 4C0 z1 λ30 λ1 x01 + p1 λ30 x00 x01 z00 + 2C0 z1 λ40 x02 z00 + 4C0 x2 λ30 λ1 z002 − 4C0 x1 λ20 λ21 z002 0 4 0 0 3 0 0 + 2C0 x1 λ30 λ2 z002 − p1 λ30 x02 0 z1 − 4C0 x2 λ0 z0 z1 + 4C0 x1 λ0 λ1 z0 z1 0 02 0 0 0 02 0 − 2C0 x1 λ40 z00 z20 − 4C0 λ1 x03 0 λ0 + 2C0 λ0 x0 x1 λ0 − 4C0 λ1 x0 z0 λ0 0 0 02 0 02 00 + 2C0 λ0 x00 z00 z10 λ00 + 2C0 λ0 x03 0 λ1 + 2C0 λ0 x0 z0 λ1 + 2C0 λ0 λ1 x0 x0  00 0 0 00 2 0 0 00 − 2C0 λ20 x02 0 x1 + 2C0 λ0 λ1 x0 z0 z0 − 2C0 λ0 x0 z0 z1 , 1

AC

2C0 λ40 (x02 0 − + − + − + +



2 2 02 3 02 4C0 z2 λ30 λ1 x02 0 − 4C0 z1 λ0 λ1 x0 + 2C0 z1 λ0 λ2 x0 + z002 ) 4C0 z2 λ40 x00 x01 + 4C0 z1 λ30 λ1 x00 x01 − 2C0 z1 λ40 x00 x02 − 8C0 x2 λ30 λ1 x00 z00 4C0 x1 λ20 λ21 x00 z00 − 4C0 x1 λ30 λ2 x00 z00 + z2 λ30 p00 x00 z00 + z1 λ30 p01 x00 z00 4C0 z2 λ30 λ1 z002 − 2C0 z1 λ30 λ2 z002 − x2 λ30 p00 z002 − x1 λ30 p01 z002 + p1 λ30 x01 z002 4C0 x2 λ40 x00 z10 − 4C0 x1 λ30 λ1 x00 z10 − p1 λ30 x00 z00 z10 + 2C0 x1 λ40 x00 z20 0 0 02 0 0 0 0 0 0 03 0 4C0 λ1 x02 0 z0 λ0 + 2C0 λ0 x0 x1 z0 λ0 − 4C0 λ1 z0 λ0 + 2C0 λ0 z0 z1 λ0 0 0 03 0 0 0 00 2 0 0 00 2C0 λ0 x02 0 z0 λ1 + 2C0 λ0 z0 λ1 + 2C0 λ0 λ1 x0 z0 x0 − 2C0 λ0 x0 z0 x1  2C0 λ0 λ1 z002 z000 − 2C0 λ20 z002 z100 ,

CE

z3 =

PT

ED

x3 =

AN US

CR IP T

x1 = −

49

ACCEPTED MANUSCRIPT

 1 8C0 z2 λ20 λ1 x00 − 4C0 z1 λ0 λ21 x00 + 4C0 z1 λ20 λ2 x00 + x2 λ20 p00 x00 02 + z0 ) + x1 λ20 p01 x00 − 4C0 z2 λ30 x01 + 4C0 z1 λ20 λ1 x01 − p1 λ20 x00 x01 − 2C0 z1 λ30 x02 − 8C0 x2 λ20 λ1 z00 + 4C0 x1 λ0 λ21 z00 − 4C0 x1 λ20 λ2 z00 + z2 λ20 p00 z00 + z1 λ20 p01 z00 + 4C0 x2 λ30 z10 − 4C0 x1 λ20 λ1 z10 − p1 λ20 z00 z10 + 2C0 x1 λ30 z20 + 2C0 x01 z00 λ00  − 2C0 x00 z10 λ00 + 2C0 λ1 z00 x000 − 2C0 λ0 z00 x001 − 2C0 λ1 x00 z000 + 2C0 λ0 x00 z100 .

λ20 (x02 0

CR IP T

p2 =

Appendix E. Derivation of the amplitude equations in the postbifurcation analyses

AN US

In this appendix, we introduce the derivation procedure of the amplitude equations in the post-bifurcation analyses. The perturbation expansions of the unknowns x0 , z0 and the parameters λ1 , γ0 have been given in (83). By substituting (83) into (67) and (68) and equating the coefficients of i (i = 0, 1, 2, · · · ) to be zero, we can obtain a series of ODE systems. The equations corresponding to 0 are automatically satisfied. From the coefficients of i (i = 1, 2, 3), we have (3)

PT

ED

M

2(1 + 3λ4c )U10 4h(1 + λ4c )W100 4h2 (1 + 3λ4c + 4λ8c )U1 (E1)1 − = 0, − λ2c λc 3λ4c   (3) 4 00 3 )U + 4(1 + λ α W h λ 4 00 1 c c 1 1 α1 W 1 2(λc − 1)W1 − + λc α1 U10 − − 2h λ2c λc (E1)2  2  2h (3) (4) = 0, − (1 + 2λ4c )α1 U1 − 2λc (3 + λ4c )W1 3λc −4(1 + λ4c )∆λ + 2λc (U20 + λ4c (−6U102 + 3U20 + 2W102 )) λ3c   4h 4 0 00 4 00 4 0 00 + (2(1 + λc )W1 U1 + 2(1 + 3λc )U1 W1 − (1 + λc )W2 ) λc + h2 F1 (W1 , U1 , U2 ) = 0, (2)

AC

CE

− γ0 +

50

(E2)1

ACCEPTED MANUSCRIPT

 α1 W2 2(1 − λ4c )W200 + + λ − α1 U102 + α1 (U20 + W102 ) c 2h λ2c  h + 8λc (W10 U100 + U10 W100 ) + 2 λc (8(1 + 3λ4c )U1002 − 8(1 + λ4c )W1002 λc

− α1 ∆λ −

(3)

(3)

(E2)2

CR IP T

− λ3c α1 W200 + 4U10 (λ3c α1 W100 + 2(1 + 3λ4c )U1 ) − 4(1 + λ4c )U2 )  (3) + W10 ((1 + 3λ4c )α1 U100 − 8λc (1 + λ4c )W1 )

AN US

+ h2 F2 (W1 , U1 , W2 , U2 ) = 0. 2 (3) − γ0 + 3 10λ5c U103 − 2U10 ((1 − 3λ4c )∆λ + λ5c (6U20 + 5W102 )) λc  4 0 5 0 0 + λc (1 + 3λc )U3 + 4λc W1 W2 + hF3 (W1 , U1 , W2 , U2 , W3 )

(E3)1

M

+ h2 F4 (W1 , U1 , W2 , U2 , U3 ) = 0,   α1 W3 4(1 + λ4c )W100 0 − + ∆λ α1 U1 − + λc α1 U103 + 2W10 W20 2h λ3c  − U10 (2U20 + 3W102 ) + U30 + 4λ2c − 10U10 W10 U100 + 2W20 U100 + 2W10 U200  2(1 − λ4c )W300 − 5U102 W100 + (2U20 + 3W102 )W100 + 8λ2c U10 W200 + λ2c 2 + hF5 (W1 , U1 , W2 , U2 , W3 , U3 ) + h F6 (W1 , U1 , W2 , U2 , W3 , U3 ) = 0,

CE

PT

ED

(E3)2 where the lengthy expressions of F1 -F6 are omitted for brevity. By further substituting (83) into the boundary conditions (69) and considering the coefficients of i (i = 1, 2, 3), we have ( (3) (3) Wi0 (−1) = 0, Wi0 (1) = 0, Wi (−1) = 0, Wi (1) = 0, (E4) Ui (−1) = 0, Ui (1) = 0, Ui00 (−1) = 0, Ui00 (1) = 0, i = 1, 2, 3.

AC

Next, we shall derive the post-bifurcation solutions by solving the equations (E1)-(E3) subject to the boundary conditions (E4). At O(1 ), the system (E1) is in fact same as the linearized system (73)(74). With the boundary conditions (E4)i=1 and through the conventional derivations, we obtain the following expressions of W1 and U10 ( W1 (X) = 2R1 cos(ξ1 X + θ1 ) + 2R2 cos(ξ2 X + θ2 ), (E5) U10 (X) = 2f (ξ1 )R1 cos(ξ1 X + θ1 ) + 2f (ξ2 )R2 cos(ξ2 X + θ2 ), 51

ACCEPTED MANUSCRIPT

where R1 and R2 are the amplitude parameters, ξi and θi (i = 1, 2) are some constants to be determined and the function f is defined by 2hλc (1 + λ4c )ξ 2 . 1 + 3λ4c

(E6)

CR IP T

f (ξ) = −

AC

CE

PT

ED

M

AN US

To obtain nontrivial solutions, at least one of R1 and R2 should be nonzero. Therefore, we consider two cases: (1) R1 6= 0 and R2 = 0 (the similar results can also be obtained for R1 = 0 and R2 6= 0); (2) both R1 and R2 are nonzero. By substituting (E5) into the boundary conditions (E4)i=1 and through some simple analyses, it can be obtained that   ξ1 = k1 π, θ1 = k¯1 π, k1 = 1, 2, 3, · · · , k¯1 = 0, ±1, ±2, · · · , ¯  ξ1 = 2k1 − 1 π, θ1 = 2k1 − 1 π, k1 = 1, 2, 3, · · · , k¯1 = 0, ±1, ±2, · · · 2 2 (E7) in the first case and   ξ1 = k1 π, θ1 = k¯1 π, k1 = 1, 2, 3, · · · , k¯1 = 0, ±1, ±2, · · · , ¯  ξ1 = 2k1 − 1 π, θ1 = 2k1 − 1 π, k1 = 1, 2, 3, · · · , k¯1 = 0, ±1, ±2, · · · , 2 2 (E8)1  ¯ ¯ ξ = k π, θ = k π, k = 1, 2, 3, · · · , k = 0, ±1, ±2, · · · ,  2 2 2 2 2 2 ¯  ξ2 = 2k2 − 1 π, θ2 = 2k2 − 1 π, k2 = 1, 2, 3, · · · , k¯2 = 0, ±1, ±2, · · · , 2 2 (E8)2 in the second case. We note that ξi (i = 1, 2) are the imaginary parts of the roots of (78). The values of the amplitude parameters R1 and R2 will be determined later. At O(2 ), the asymptotic expression of U20 can be derived from (E2)1 , in which the terms of order higher than O(h) have been dropped (analog to the case of O()). By substituting U20 into (E2)2 and using the expressions of U1 and W1 obtained in the previous step, we obtain the following equation for W2 (4)

ψ4 W2 + ψ2 W200 + ψ0 W2 (2)

=cλ ∆λ + cγ γ0 + 2c1 R12 cos(2(ξ1 X + θ1 )) + 2c2 R22 cos(2(ξ2 X + θ2 )) + 2c3 R1 R2 cos((ξ1 + ξ2 )X + θ1 + θ2 ) + 2c4 R1 R2 cos((ξ1 − ξ2 )X + θ1 − θ2 ) + c5 R12 + c6 R22 , 52

(E9)

ACCEPTED MANUSCRIPT

where ci (i = 1, · · · , 6, λ, γ) are some constants related to h, α1 , λc , ξ1 and ξ2 . The general solution of (E9) can be written as

CR IP T

¯ 1 cos(ξ1 X + θ¯1 ) + 2R ¯ 2 cos(ξ2 X + θ¯2 ) + Hλ ∆λ + Hγ γ0(2) W2 =2R + 2H1 R12 cos(2(ξ1 X + θ1 )) + 2H2 R22 cos(2(ξ2 X + θ2 )) + 2H3 R1 R2 cos((ξ1 + ξ2 )X + θ1 + θ2 ) + 2H4 R1 R2 cos((ξ1 − ξ2 )X + θ1 − θ2 ) + H5 R12 + H6 R22 ,

(E10)1

AN US

¯1, R ¯ 2 , θ¯1 , θ¯2 and γ0(2) are some unknown constants to be determined where R and Hi (i = 1, · · · , 6, λ, γ) are constants related to h, α1 , λc , ξ1 and ξ2 . Also, we can get the following expression of U20 ¯ 1 cos(ξ1 X + θ¯1 ) + 2f (ξ2 )R ¯ 2 cos(ξ2 X + θ¯2 ) + Lλ ∆λ U20 =2f (ξ1 )R (2)

+ Lγ γ0 + 2L1 R12 cos(2(ξ1 X + θ1 )) + 2L2 R22 cos(2(ξ2 X + θ2 )) (E10)2 + 2L3 R1 R2 cos((ξ1 + ξ2 )X + θ1 + θ2 ) + 2L4 R1 R2 cos((ξ1 − ξ2 )X + θ1 − θ2 ) + L5 R12 + L6 R22 ,

ED

M

where Li (i = 1, · · · , 6, λ, γ) are constants related to h, α1 , λc , ξ1 and ξ2 . The boundary conditions (E4)i=2 should also be satisfied by W2 and U2 . By substituting (E10) into (E4)i=2 and through some analyses, it can be obtained that   sin ξ1 − θ¯1 = 0, sin ξ1 + θ¯1 = 0,   (E11) sin ξ2 − θ¯2 = 0, sin ξ2 + θ¯2 = 0. (2)

PT

Furthermore, the following expression of γ0 can be derived (2)

CE

γ0 = −(Lλ ∆λ + L5 R12 )/Lγ

(E12)

AC

in the case R1 6= 0 and R2 = 0, and (2)

γ0 = −(Lλ ∆λ + L5 R12 + L6 R22 )/Lγ

(E13)

in the case R1 6= 0 and R2 6= 0. At O(3 ), the asymptotic expression of U30 can be derived from (E3)1 , where the terms of order higher than O(h) have also been dropped. By substituting this expression into (E3)2 , we obtain the following single equation

53

ACCEPTED MANUSCRIPT

for W3 (4)

ψ4 W3 + ψ2 W300 + ψ0 W3 (3)

AN US

CR IP T

=dγ γ0 + 2d1 R13 cos(3(ξ1 X + θ1 )) + 2d2 R23 cos(3(ξ2 X + θ2 )) + 2d3 R12 R2 cos((2ξ1 + ξ2 )X + 2θ1 + θ2 ) + 2d4 R1 R22 cos((ξ1 + 2ξ2 )X + θ1 + 2θ2 ) + 2d5 R12 R2 cos((2ξ1 − ξ2 )X + 2θ1 − θ2 ) + 2d6 R1 R22 cos((ξ1 − 2ξ2 )X + θ1 − 2θ2 ) i h  (2) (E14) + 2 sλ ∆λ + sγ γ0 R1 + s1 R1 3 + s2 R1 R22 cos(ξ1 X + θ1 )  i h (2) + 2 tλ ∆λ + tγ γ0 R2 + t1 R2 3 + t2 R2 R12 cos(ξ2 X + θ2 )  ¯ 1 d7 R1 cos(2ξ1 X + θ¯1 + θ1 ) + d9 R2 cos((ξ1 + ξ2 )X + θ¯1 + θ2 ) + 2R  +d9 R2 cos((ξ1 − ξ2 )X + θ¯1 + θ2 ) + d10 R1 cos(θ¯1 − θ1 )  ¯ 2 d11 R1 cos((ξ1 + ξ2 )X + θ¯2 + θ1 ) + d12 R1 cos((−ξ1 + ξ2 )X + 2R  +θ¯2 − θ1 ) + d13 R2 cos(2ξ2 X + θ¯2 + θ2 ) + d14 R2 cos(θ¯2 − θ2 ) ,

PT

ED

M

where di (i = 1, · · · , 14, γ), si and ti (i = λ, γ, 1, 2) are some constants related to h, α1 , λc , ξ1 and ξ2 . Next, we shall derive the amplitude equations for R1 and R2 by considering the solvability condition of (E14) subject to the boundary conditions (E4)i=3 . d4 d2 For convenience, we denote the operator L = ψ4 dX 4 + ψ2 dX 2 + ψ0 . In the function space W = {w|w ∈ C 4 [−1, 1], w0 (−1) = w0 (1) = 0, w(3) (−1) = ∗ w(3) (1) = 0}, we gi = hf, L∗ gi, R 1define the adjugate operator L of L as hLf, where hf, gi = −1 f gdX. In fact, it is easy to check that L∗ = L. Thus, the solution of equation L∗ p = 0 is just given by

CE

˜ 1 cos(ξ1 X + θ˜1 ) + 2R ˜ 2 cos(ξ2 X + θ˜2 ). p(X) = 2R

(E15)

AC

R1 Consequently, the solvability condition for (E14) requires that −1 pFr dX = 0, where Fr denotes the right-hand side of (E14). Through some calculations, it can be found that to ensure the solvability condition, the coefficients of cos(ξ1 X + θ1 ) and cos(ξ2 X + θ2 ) in (E14) should vanish, which yield that h i (2) R1 sλ ∆λ + sγ γ0 + s1 R12 + s2 R22 = 0, h i (E16) (2) R2 tλ ∆λ + tγ γ0 + t1 R22 + t2 R12 = 0. 54

ACCEPTED MANUSCRIPT

CR IP T

In the case R1 6= 0 and R2 = 0, by substituting (E12) into (E16)1 , the following result for the amplitude parameter R1 can be obtained sλ L γ − sγ L λ . (E17) R12 = r1 ∆λ, r1 = sγ L 5 − s1 L γ

AN US

With the expressions of θ1 , ξ1 and R1 , we obtain the following expressions of the leading-order nontrivial solutions ( z0 =zh + 2 (r1 δλ)1/2 cos(k1 πX + k¯1 π), (E18)1 x00 =x0h + 2f (k1 π) (r1 δλ)1/2 cos(k1 πX + k¯1 π), ( z0 =zh + 2 (r1 δλ)1/2 cos((k1 − 1/2)πX + (k¯1 − 1/2)π), x00 =x0h + 2f ((k1 − 1/2)π) (r1 δλ)1/2 cos((k1 − 1/2)πX + (k¯1 − 1/2)π), (E18)2 which are the analytical expressions for describing the post-bifurcation states of the plate. Similarly, when R1 = 0 and R2 6= 0, the approximate solutions of z0 and x00 can also be obtained. In the case that both R1 and R2 are not equal to zero, by substituting (E13) into (E16) and solving the two equations, we obtain

M

R12 = r1 ∆λ, where

R22 = r2 ∆λ,

(E20)

PT

ED

Lλ (sγ t1 − s2 tγ ) + Lγ (−sλ t1 + s2 tλ ) + L6 (sλ tγ − sγ tλ ) , Lγ (s1 t1 − s2 t2 ) + L6 (sγ t2 − s1 tγ ) + L5 (−sγ t1 + s2 tγ ) Lλ (−sγ t2 + s1 tγ ) + Lγ (sλ t2 − s1 tλ ) + L5 (−sλ tγ + sγ tλ ) . r2 = Lγ (s1 t1 − s2 t2 ) + L6 (sγ t2 − s1 tγ ) + L5 (−sγ t1 + s2 tγ ) r1 =

(E19)

AC

CE

With the expressions of θi , ξi and Ri (i = 1, 2), we can obtain the expressions of the leading-order nontrivial solutions z0 and x00 . For example, the solutions corresponding to ξ1 = k1 π and ξ2 = k2 π are given by  z0 =zh + 2 (r1 δλ)1/2 cos(k1 πX + k¯1 π)      + 2 (r2 δλ)1/2 cos(k2 πX + k¯2 π), (E21) 1/2 0 0 ¯1 π)  x =x + 2f (k π) (r δλ) cos(k πX + k  1 1 1 0 h    + 2f (k2 π) (r2 δλ)1/2 cos(k2 πX + k¯2 π).

Remark: The interested readers can contact the authors for the explicit expressions of the constants ci , Hi , Li (i = 1, · · · , 6, λ, γ), di (i = 1, · · · , 14, γ), si and ti (i = λ, γ, 1, 2). 55

ACCEPTED MANUSCRIPT

References References

CR IP T

Audoly, B., Boudaoud, A., 2008. Buckling of a stiff film bound to a compliant substrate: part I-III. J. Mech. Phys. Solids 56, 2401-2458. Ben Amar, M., Goriely, A., Growth and instability in elastic tissues. J. Mech. Phys. Solids 53, 2284-2319.

Budday, S., Steinmann, P., Kuhl, E., 2014. The role of mechanics during brain development. J. Mech. Phys. Solids 72, 75-92.

AN US

Cai, S., Breid, D., Crosby, A. J., Suo, Z., Hutchinson, J. W., 2011. Periodic patterns and energy states of buckled films on compliant substrates. J. Mech. Phys. Solids 59, 1094-1114.

Cao, Y.P., Hutchinson J. W., From wrinkles to creases in elastomers: the instability and imperfection-sensitivity of wrinkling, Proc. R. Soc. A 468, 94-115.

ED

M

Coen, E., Rolland-Lagan, A. G., Matthews, M., Bangham, J. A., Prusinkiewicz, P., 2004. The genetics of geometry. Proc. Natl. Acad. Sci. 101, 4728-4735. Dai, H.-H., Song, Z.L., 2014. On a consistent finite-strain plate theory based on a 3-D energy principle. Proc. R. Soc. A. 470, 20140494.

PT

Dervaux, J., Ben Amar, M., 2008. Morphogenesis of Growing Soft Tissues. Phys. Rev. Lett. 101, 068101.

CE

Dervaux, J., Ciarletta, P., Ben Amar, M., 2009. Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the F¨ oppl-von K´ arm´ an limit. J. Mech. Phys. Solids 57, 458-471.

AC

Dervaux, J., Ben Amar, M., 2010. Localized growth of layered tissues. IMA J. Appl. Math. 75, 571-580. Fung, Y. C., 1990. Biomechanics: Motion, Flow, Stress, and Growth. Springer, New York.

56

ACCEPTED MANUSCRIPT

Givnish, T. J., 1987. Comparative studies of leaf form: assessing the relative roles of selective pressures and phylogenetic constraints. New Phytologist 106, 131-160.

CR IP T

Green, P. B., 1996. Transductions to generate plant form and pattern: an essay on cause and effect. Ann. Botany 78, 269-281. Holland, M. A., Kosmata, T., Goriely, A., Kuhl, E., 2013. On the mechanics of thin films and growing surfaces. Math. Mech. Solids 18, 561-575. Huang, Z. Y., Hong, W., Suo, Z., 2005. Nonlinear analyses of wrinkles in a film bonded to a compliant substrate. J. Mech. Phys. Solids 53, 2101-2118.

AN US

Huang, X., Li, B., Hong, W., Cao, Y. P., Feng, X. Q., 2016. Effects of tensionCcompression asymmetry on the surface wrinkling of filmCsubstrate systems. J. Mech. Phys. Solids 94, 88-104.

Humphrey, J. D., 2002. Continuum biomechanics of soft biological tissues. Proc. R. Soc. Lond. A 459, 3-46.

M

Jin, L. H., Cai, S.Q., Suo, Z.G., 2011. Creases in soft tissues generated by growth. Europhysics Letters 95, 64002.

ED

Li, B., Cao, Y. P., Feng, X. Q., Gao, H. J., 2011. Surface wrinkling of mucosa induced by volumetric growth: theory, simulation and experiment. J. Mech. Phys. Solids 59, 758-774.

PT

Li, B., Cao, Y.P., Feng, X.Q., Gao, H.J., 2012. Mechanics of morphological instabilities and surface wrinkling in soft materials: a review. Soft Matter 8, 5728-5745.

CE

Liu, Z.S., Swaddiwudhipongc, S., Hong, W., 2013. Pattern formation in plants via instability theory of hydrogels. Soft Matter 9, 577-587.

AC

Lubarda, V. A., Hoger, A., 2002. On the mechanics of solids with a growing mass. International Journal of Solids and Structures 39, 4627-4664. Marder, M., Papanicolaou, N., 2006. Geometry and elasticity of strips and flowers. Journal of Statistical Physics 125, 1065-1092. Menzel, A., Kuhl, E., 2012. Frontiers in growth and remodeling. Mech. Res. Commun. 42, 1-14. 57

ACCEPTED MANUSCRIPT

Meroueh, K., 1986. On a formulation of a nonlinear theory of plates and shells with applications. Comput. Struct. 24, 691-705.

CR IP T

Moulton, D. E., Lessinnes, T., Goriely, A., 2013. Morphoelastic rods. Part I: A single growing elastic rod. J. Mech. Phys. Solids 61, 398-427.

Nath, U., Crawford, B. C. W., Carpenter R., Coen E., 2003. Genetic control of surface curvature. Science 299, 1404-1407.

Newell, A. C., Shipman, P. D., 2005. Plants and fibonacci. J. Stat. Phys. 121, 937-968.

AN US

Ogden, R. W., 1984. Non-linear Elastic Deformation. Dover, Newyork.

Papastavrou, A., Steinmann, P., Kuhl, E., 2013. On the mechanics of continua with boundary energies and growing surfaces. J. Mech. Phys. Solids 61, 1446-1463.

M

Rausch, M. K., Kuhl, E., 2013. On the effect of prestrain and residual stress in thin biological membranes. J. Mech. Phys. Solids 61, 1955-1969. Rausch, M. K., Kuhl, E., 2014. On the mechanics of growing thin biological membranes. J. Mech. Phys. Solids 63, 128-140.

ED

Reddy, J. N., 2007. Theory and analysis of elastic plates and shells. CRC Press, Taylor and Francis.

PT

Rodriguez, A. K., Hoger, A., McCulloch, A., 1994. Stress-dependent finite growth in soft elastic tissue. J. Biomech. 27, 455-467.

CE

Skalak,R., Zargaryan, S., Jain, R. K., Netti, P. A., Hoger, A., 1996. Compatibility and the genesis of residual stress by volumetric growth. J. Math. Biol. 34, 889-914.

AC

Song, Z.L., Dai, H.-H., 2016. On a consistent dynamic finite-Strain plate theory and its linearization. J. Elast. 125, 149-183. Steigmann, D. J., 2007. Thin-plate theory for large elastic deformations. Int. J. Non-Linear Mech. 42, 233-240. Steigmann, D. J., 2013a. A well-posed finite-strain model for thin elastic sheets with bending stiffness. Math. Mech. Solids 18, 103-112. 58

ACCEPTED MANUSCRIPT

Steigmann, D. J., 2013b. Koiters shell theory from the perspective of threedimensional nonlinear elasticity. J. Elast. 111, 91-107.

CR IP T

Steigmann, D. J., 2015. Mechanics of materially uniform thin films. Math. Mech. Solids 20, 309-326. Taber, L., 1995. Biomechanics of growth, remodeling and morphogenesis. Appl. Mech. Rev. 48, 487-545.

Thompson D. W., 1942. On growth and form. Cambridge and New York, University Press and Macmillan, 1942.

AN US

Timoshenko, S., Woinowsky-Krieger, S., 1959. Theory of plates and shells. McGrawHill, New York.

Wang, J., Dai, H.-H., Song, Z.L., 2016. On a consistent finite-strain plate theory for incompressible hyperelastic materials. Inter. J. Solids Struct. 78-79, 101-109.

AC

CE

PT

ED

M

Wang, J., Wang, Q.Y., Dai, H.-H., Exact solution for growth-induced large bending deformation of a hyperelastic plate. arXiv:1710.03120v1 [physics.class-ph] (http://arxiv.org/abs/1710.03120).

59