An incremental plate theory for polymer gels in equilibrium

Mechanics Research Communications 96 (2019) 49–55

Contents lists available at ScienceDirect

Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom

An incremental plate theory for polymer gels in equilibrium Xiaoyi Chen a, Hui-Hui Dai b,∗ a b

Beijing Normal University-Hong Kong Baptist University United International College, 2000 Jintong Road, Tangjiawan, Zhuhai, Guangdong Province, China Department of Mathematics, City University of Hong Kong, 83 Tat Chee Avenue, Kowloon, Hong Kong SAR, China

a r t i c l e

i n f o

Article history: Received 31 January 2019 Revised 1 March 2019 Accepted 2 March 2019 Available online 4 March 2019 Keywords: Incremental plate theory Instability Polymer gel

a b s t r a c t Based on the finite-strain plate theory derived by one of the authors, an incremental plate theory for deformations superimposed on a general base state is formulated in this paper. The unknown functions of the incremental deformation are first expanded into Taylor series in terms of the thickness variable and then expanded around the base state by retaining third-order nonlinearity. From the field equations and the boundary conditions at the top and bottom surfaces, the recursive relations of the expansion coefficients as well as the incremental balance equations are derived. With the constitutive relation for swollen polymer gels in equilibrium, an incremental plate theory is obtained. As an application, this theory is used to study the incremental deformation of a polymer gel layer with a homogeneous base state in a plane strain setting. A linear bifurcation analysis is carried out, which gives the critical values of the external chemical potential and the mode number. The results agree with those obtained directly from a full two-dimensional analysis. Post-bifurcation is also conducted through a perturbation procedure, which reveals that the bifurcation is of supercritical type. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction In solid mechanics, instability problems are usually studied by the incremental theory for an incremental deformation superimposed on a given finite deformation or base state. The main mathematical framework of this method is described in books by Biot [1] and Ogden [2]. For the instability of thin structures, it is possible to use an incremental plate theory which combines the plate theory and the incremental theory. One type of plate theories [3–6] are based on ad hoc kinematic and other assumptions. One disadvantage is that these assumptions may not be valid in a general loading condition. An incremental plate theory [7,8] based on it naturally inherits this defect. Over the years, many efforts have been made to derive plate theories in a more consistent manner and great advances have been achieved in [9–11], which are based on the variational (weak) formulation. Recently, Dai and Song [12] derived a consistent finitestrain plate theory directly from the 3D differential formulation, which combines bending and stretching effects together. Wang et al. [13] adopted the procedure in [12] to derive a plate theory of growth and used it to study the instability of a biological-tissue layer induced by growth, and their results demonstrated the efficiency of the model. In this paper, by adopting the plate model in

∗

Corresponding author. E-mail address: [email protected] (H.-H. Dai).

https://doi.org/10.1016/j.mechrescom.2019.03.001 0093-6413/© 2019 Elsevier Ltd. All rights reserved.

[12], an incremental plate theory is systemically constructed. With the use of the constitutive relation for polymer gels in equilibrium, an incremental plate theory for polymer gels is presented. We point out that recently a seven-parameter large-strain poroelastic plate theory for polymer gels was developed in Lucantonio et al. [14] for a swelling process. In contrast, here we only consider equilibrium states. Nevertheless, the present plate theory has a few different features, including (i) no ad hoc assumption on the form of the unknowns (except the smoothness assumption); (ii) no higher-order stress resultants involved; (iii) only three differential equations with three unknowns; (iv) satisfying the 3D field equations and top/bottom traction conditions in a pointwise manner. The paper is arranged as follows. In Section 2, we briefly recall the plate theory in [12]. The constitutive relation for polymer gels in equilibrium is given in Section 3 according to the thermodynamic derivation of Lucantonio et al. [15]. A weakly nonlinear incremental plate theory is derived in Section 4, which is then specialized to a plane strain problem in Section 5. In Section 6, the framework in Section 5 is applied to study the instability of a polymer gel layer attached to a Winkler foundation, and both linear and post-bifurcation analyses are carried out analytically. In particular, for the linear bifurcation, good agreements are found with the results based on a full 2D setting, which gives a validation of this incremental plate theory. The post-bifurcation analysis reveals that the bifurcation is of supercritical type.

50

X. Chen and H.-H. Dai / Mechanics Research Communications 96 (2019) 49–55

2. Plate equations in Dai and Song [12]

[16–19], the free energy of the swollen polymer gel (taking the dry state as the reference state) adopts the form of

We consider the deformation of a thin plate subjected to surface tractions and edge traction/position constraints. The plate is occupying the region  =  × [0, H] (H is much smaller than the length scale of the bottom surface ) in the reference configuration and the position of the material point in this state is denoted as X. In the current configuration, the material point is transported to the position x(X) and the deformation gradient is F = ∂ x/∂ X and the nominal stress is denoted as S. Neglecting the body force, the governing equation and the boundary conditions at the bottom and top surfaces are



Div S = 0, in , ST (−e3 ) = q− , Z=0



ST e3 

Z=H

(1)

= q+ , on ,

where q+ and q− are respectively tractions on the top and bottom surfaces. Note the boundary conditions at lateral edges will be proposed under specific applications. The governing Eq. (1)1 are PDEs with variables X, Y and Z. By taking advantage of the small geometric size in the thickness, plate equations can be derived from (1), which are PDEs only with variables X and Y by the method introduced in [12]. Here, to be concise, the main results are briefly listed. Assuming sufficient smoothness, then one has the following Taylor expansions about Z = H:

x (X ) =

4  ( Z − H )k ( k ) x (X, Y ) + O((Z − H )5 ), k! k=0

3  ( Z − H )k ( k ) F (X ) = F (X, Y ) + O((Z − H )4 ), k! k=0

3  ( Z − H )k ( k ) S (X ) = S (X, Y ) + O((Z − H )4 ). k!

(2)

k=0

From the relations between x, F and S, we have

F(n ) = ∇ x(n ) + x(n+1)  e3 ,

n = 0, 1, 2, 3,

S ( 0 ) = S ( F ( 0 ) ) , S ( 1 ) = A 1 ( F ( 0 ) )[ F ( 1 ) ] , S ( 2 ) = A 1 ( F ( 0 ) )[ F ( 2 ) ] + A 2 ( F ( 0 ) )[ F ( 1 ) , F ( 1 ) ] , where ∇ x(n ) =

(3)

(n )

(∂ xi /∂ X j )ei  e j (n = 0, 1, 2, 3 ) (summation convention is applied; X1 , X2 denote X, Y respectively), and A 1 = ∂ S/∂ F, A 2 = ∂ 2 S/∂ F∂ F. Note S(3) is intermediate and its expression is not shown here. A closed governing system for x(0 ) − x(3 ) can be derived by employing the field equations and the boundary conditions at the top and bottom surfaces (cf. [12]), which are

 ∇ · S (0) + S(1)T e3 = 0, ∇ · S(1) + S(2)T e3 = 0, S(0)T e3 = q+ , 2 ∇ · S(0) − H2 S(1) + H6 S(2) + O(H 3 ) = −q¯ ,

where ∇ · S(k ) = ∂ Si(jk ) e j /∂ Xi , (k = 0, 1, 2 ) and q¯ = (q+ + q− )/H. x(3) ,

Eq. (4)1,2 are linear algebraic equations for and which can be used to eliminate these two unknown vectors. Eq. (4)3 is a nonlinear algebraic equation for x(1) . By substituting x(1 ) − x(3 ) to (4)4 , one vector equation only containing x(0) is obtained. The associated boundary conditions have been extensively discussed in [12] and here we omit the details. Equations in (4) are the plate equations which will be used in Section 4. 3. The constitutive relation for a polymer gel Consider a dry polymer network in the shape of  (see Section 2). It swells through the interaction with the external solvent and the swollen network is an elastomeric gel. By referring to

1 NkT [I1 − 3] − NkT ln J 2  νc χ νc  + kT c ln + , 1 + νc 1 + νc

(5)

where N is the number of polymer chains per unit dry volume, k is the Boltzmann’s constant, T is the absolute temperature, F is the deformation gradient, J = det F, I1 = tr(FT F ), χ is a dimensionless parameter describing the enthalpy of mixing, ν is the volume of one solvent molecule and c is the number of solvent molecules inside the polymer network per unit dry volume. Note in (5), the term −N kT ln J is included to describe the finite deformation of the polymer gel based on the derivations of Flory and Rehner [18] and Simo et al. [19]. Actually, as pointed out by Treloar [20], whether this term should be added or not is still inconclusive. In [16,18,19], this term is present. Based on (5) without the term −NkT ln J, Lucantonio et al. [15] derived a theoretical and computational model to describe the deformation of polymer gels by introducing an intermediate free swelling state as the reference configuration. The dry state, free swelling state and the current state are denoted respectively as B, Bm and Bt . The deformation gradient is given by F = Fm F0 , where F0 = λ0 I and Fm are respectively the deformation gradients from B to Bm and from Bm to Bt . In the sequel, the subscript “m” will be used to denote a quantity associated with Bm . The molecule chain is incompressible and this constraint is put as 1 = det Fm = Jm or equivalently J = 1 + ν c = λ30 . According to a thermodynamic derivation of Lucantonio et al. [15], the constitutive relations for the nominal stress and chemical potential are given by

⎧ ∂φm ⎪ ⎨Sm = ∂φm − pF−1 ( Sm )i j = − p(F−1 m , m )i j , ∂ Fm ∂ (Fm ) ji ⎪ ⎩μ = ∂φm + pν, ∂ cm

(6)

where φm = J0−1 φ , J0 = det F0 = λ30 , cm = J0−1 c and μ is the chemical potential of one solvent molecule inside the polymer network and p is the Lagrangian multiplier. The chemical potential of the external solvent molecule is denoted as μ0 . When the swelling reaches equilibrium, we have the relation μ = μ0 . This relation together with (6)3 (with the use of (5)) determines p as

p=

NkT

λ30

+



λ3 − 1 1 μ0 kT χ − ln 0 3 + 3 + 6 . ν ν λ0 λ0 λ0

(7)

The parameter λ0 is fixed in the free swelling state. By setting Fm = 1 and Sm = 0, we have

μ0 (4)

x(2)

φ ( F, c ) =

kT

= Nν

1

λ0

−

1



λ30

 λ30 − 1 1 χ + ln + 3 + 6 . λ30 λ0 λ0

(8)

The swelling development at steady state can be represented by λ0 or through (8) by μ0 /(kT). Later on, for convenience, μ0 /(kT) is taken as the driving parameter for the swelling states in equilibrium. By sketching the graph, it is easily found that λ0 is an increasing function of μ0 /(kT). Taking (7) into (6)1 , upon using the relation S = λ20 Sm , we obtain

S = S¯ +

μ0 ˜ kT



S,

S˜ = −



S¯ = N kT FT − F−1 +

kT

ν

kT

ν

JF−1 ,



J ln



J−1 1 χ + + 2 F−1 . J J J

(9)

In the next section, the incremental plate theory is constructed by taking the nominal stress in the form of S = S¯ + γ S˜ without

X. Chen and H.-H. Dai / Mechanics Research Communications 96 (2019) 49–55

specifying the expressions for S¯ , S˜ and γ . This may make the derivation applicable to other cases, such as taking γ = χ or even γ = 0 (which corresponds to a usual compressible hyperelastic material). Eq. (9)1 is a special case of this form for γ = μ0 /(kT ).

+ A 2 (Fb(0) , γb )[F˙ (0) , F˙ (2) ] +



1 kT J−1 NkT [I1 − 3 − 2 ln J ] + (J − 1 ) ln 2 ν J kT χ J − 1 μ + − ( J − 1 ), ν J ν

1 3 (0 ) A (Fb , γb )[F˙ (0) , F˙ (0) , F˙ (2) ] 2 1 + B 6 (Fb(0) , Fb(2) , γb )[F˙ (0) , F˙ (0) , F˙ (0) ] 6





+ γ˙ G3 (Fb(0) , Fb(2) ) + G 2 (Fb(0) , Fb(2) , γb )[F˙ (0) ] + C 1 (Fb(0) )[F˙ (2) ]



1 4 (0 ) (2 ) (0 ) (0 ) G (Fb , Fb )[F˙ , F˙ ] 2 + 2B 1 (F(0) , F(1) , γb )[F˙ (1) ] + B 7 (F(0) , F(1) , γb )[F˙ (0) ] + γ˙ C 2 (Fb(0) )[F˙ (0) , F˙ (2) ] +

(10)

and it can be verified that ∂ψ ∂ F = S according to (9). Thus, one may say the polymer gel in equilibrium can be modeled as a compressible hyperelastic material without an incompressibility constraint.

+ 2B ( 2

b b (0 ) (1 ) Fb , Fb ,

b

− Sb ,

(12)

0 ≤ k ≤ 2,

where S(k) and Sb(k ) can be obtained from (3) and (9). For a small increment, the quantities S˙ (0 ) − S˙ (2 ) are further expanded at the base state, and here we keep only the third-order weak nonlinearity. The results are

1 S˙ (0) = A 1 (Fb(0) , γb ) [F˙ (0) ] + A 2 (Fb(0) , γb ) [F˙ (0) , F˙ (0) ] 2 1 3 (0 ) ( 0 ) + A (Fb , γb ) [F˙ , F˙ (0) , F˙ (0) ] 6   1 + γ˙ G1 (Fb(0) ) + C 1 (Fb(0) )[F˙ (0) ] + C 2 (Fb(0) )[F˙ (0) , F˙ (0) ] , 2 (13) S˙ (1) = A 1 (Fb(0) , γb )[F˙ (1) ] + B 1 (Fb(0) , Fb(1) , γb )[F˙ (0) ] + A 2 (Fb(0) , γb )[F˙ (0) , F˙ (1) ] +

1 2 (0 ) (1 ) B (Fb , Fb , γb )[F˙ (0) , F˙ (0) ] 2

1 3 (0 ) A (Fb , γb )[F˙ (0) , F˙ (0) , F˙ (1) ] 2 1 + B 3 (Fb(0) , Fb(1) , γb )[F˙ (0) , F˙ (0) , F˙ (0) ] + γ˙ G2 (Fb(0) , Fb(1) ) 6   + γ˙ G 1 (F(0) , F(1) )[F˙ (0) ] + C 1 (F(0) )[F˙ (1) ]



b

+ γ˙ C 2 (Fb(0) )[F˙ (0) , F˙ (1) ] +

b



+ γ˙





1 6 (0 ) (1 ) (0 ) (0 ) G (Fb , Fb )[F˙ , F˙ ], 2

(15)

⎧ ∂ S˜ (0) ∂ 2 S˜ ⎪ ⎪ G1 (Fb(0) ) = S˜ (Fb(0) )), C 1 = (Fb ), C 2 (Fb(0) ) = 2 (Fb(0) ) ⎪ ⎪ ∂F ∂F ⎪ ⎪ (0 ) (1 ) (0 ) (1 ) (0 ) (1 ) 1 2 2 ⎪ B ( F , F , γ ) = A ( F , γ ) [ F ] , B ( F , F , γb ) ⎪ b b ⎪ b b b b b b ⎪ ⎪ (0 ) (1 ) 3 ⎪ = A ( F , γ ) [ F ] , ⎪ b b b ⎪ ⎪ (0 ) (1 ) (0 ) (1 ) (0 ) (1 ) ⎪ ⎪ B 3 (Fb , Fb , γb ) = A 4 (Fb , γb )[Fb ], G2 (Fb , Fb ) ⎪ ⎪ ⎪ ( 0 ) ( 1 ) ⎪ = C 1 ( F b )[ F b ] ⎪ ⎪ ⎪ ⎪ (0 ) (1 ) (0 ) (1 ) (0 ) (1 ) (0 ) (1 ) 1 ⎪ G (Fb , Fb ) = C 2 (Fb )[Fb ], G 3 (Fb , Fb ) = C 11 (Fb )[Fb ], ⎪ ⎪ ⎪ ⎪ 3˜ ⎪ ∂ S (0 ) (0 ) (0 ) (2 ) (0 ) (2 ) 11 4 2 ⎪ ⎪ ⎪C (Fb ) = ∂ F3 (Fb ), B (Fb , Fb , γb ) = A (Fb , γb )[Fb ], ⎪ ⎪ ⎪ ⎨ B 5 ( F ( 0 ) , F ( 2 ) , γ ) = A 3 ( F ( 0 ) , γ )[ F ( 2 ) ] , B 6 ( F ( 0 ) , F ( 2 ) , γ ) b

b

b

b

(0 )

b

b

(2 )

b

b

b

= A (Fb , γb )[Fb ], ⎪ ⎪ ⎪ (0 ) (2 ) 3 1 ⎪ G ( F , F ) = C (Fb(0) )[Fb(2) ], G 2 (Fb(0) , Fb(2) ) ⎪ b b ⎪ ⎪ ⎪ ⎪ = C 2 (Fb(0) )[Fb(2) ], ⎪ ⎪ ⎪ (0 ) (2 ) (0 ) (2 ) (0 ) (1 ) ⎪ 4 ⎪ G (Fb , Fb ) = C 11 (Fb )[Fb ], B 7 (Fb , Fb , γb ) ⎪ ⎪ ⎪ ⎪ = A 3 (Fb(0) , γb )[Fb(1) , Fb(1) ], ⎪ ⎪ ⎪ ⎪ ( 0 ) ( 1 ) (0 ) (1 ) (1 ) (0 ) (1 ) 8 ⎪ B (Fb , Fb , γb ) = A 4 (Fb , γb )[Fb , Fb ], B 9 (Fb , Fb , γb ) ⎪ ⎪ ⎪ ⎪ ⎪ = A 5 (Fb(0) , γb )[Fb(1) , Fb(1) ] ⎪ ⎪ ⎪ ( 0 ) ( 1 ⎪G4 = C 2 (F )[F ) , F(1) ], G 5 (F(0) , F(1) ) = C 11 (F(0) )[F(1) , F(1) ], ⎪ b b b b b b b b ⎪ ⎪ ⎪ 4 ⎪ ⎩G 6 (F(0) , F(1) ) = C 12 (F(0) )[F(1) , F(1) ], C 12 (F(0) ) = ∂ S˜ (F(0) ). b b b b b b ∂ F4 b 4

Based on the plate equations of (4), the incremental plate equations are obtained as





∇ · S˙ (0) + S˙ (1)T e3 = 0, ∇ · S˙ (1) + S˙ (2)T e3 = 0, S˙ (0)T e3 = q˙ + , 2 ∇ · S˙ (0) − H2 S˙ (1) + H6 S˙ (2) + O(H 3 ) = −q¯˙ . (17)

(14) S˙ (2) = A 1 (Fb(0) , γb )[F˙ (2) ] + B 4 (Fb(0) , Fb(2) , γb )[F˙ (0) ]

b

(16)

1 3 (0 ) (1 ) (0 ) (0 ) G (Fb , Fb )[F˙ , F˙ ] , 2

and

b

+ γ˙ 2G 3 (Fb(0) , Fb(1) )[F˙ (0) , F˙ (1) ] + G 5 (Fb(0) , Fb(1) )[F˙ (0) ]

+

b

b



+ γ˙ 2G 1 (Fb(0) , Fb(1) , γb )[F˙ (1) ] + C 2 (Fb(0) , γb )[F˙ (1) , F˙ (1) ]

where the tensor coefficients are defined by

⎧ ⎨x˙ (k) = x(k) − xb(k) , 0 ≤ k ≤ 4, F˙ (k ) = F(k ) − Fb(k ) = ∇ x˙ (k ) + x˙ (k+1)  e3 , 0 ≤ k ≤ 3, ⎩ ˙ (k ) (k ) (k ) =S

b

1 + B 9 (Fb(0) , Fb(1) , γb )[F˙ (0) , F˙ (0) , F˙ (0) ] 6 + A 3 (F(0) , γb )[F˙ (0) , F˙ (1) , F˙ (1) ] + γ˙ G4 (F(0) , F(1) )

(11)

where the subscript “b” denotes quantities in the base state and the dot denotes the incremental quantities from the base state to the current state. All quantities are expanded similarly as in Section 2, and we have the relations for the coefficients

S

b

γb )[F˙ (0) , F˙ (1) ]

1 8 (0 ) (1 ) (0 ) B (Fb , Fb , γb )[F˙ (0) , F˙ (0) ] + A 2 (Fb , γb )[F˙ (1) , F˙ (1) ] 2 + B 3 (F(0) , F(1) , γb )[F˙ (0) , F˙ (0) , F˙ (1) ]

At certain loading condition, the position of the gel is obtained (denoted as xb ) and this state is denoted as the base state. A small variation of the loading will cause an incremental displacement x˙ , and superimposing x˙ on xb gives the current position x = xb + x˙ . In the incremental theory, behaviors of the incremental quantities are studied. In this section, based on plate Eq. (4), the incremental plate theory is constructed. The incremental quantities are + q˙ − q˙ + = q+ − q+ , q˙ − = q− − q− , q¯˙ = q˙ + , γ˙ = γ − γb , H b b x˙ = x − xb , F˙ = F − Fb , S˙ = S − Sb ,

b



+

4. Incremental plate theory



1 5 (0 ) (2 ) B (Fb , Fb , γb )[F˙ (0) , F˙ (0) ] 2

+

Remark 1. By choosing a strain energy function as (cf. [21])

ψ=

51

Eq. (17) is a closed system for x˙ (0 ) − x˙ (3 ) . More specifically, Eq. (17) are linear algebraic equations for x˙ (2 ) and x˙ (3 ) respectively 1,2

and (17)3 is a nonlinear algebraic equation for x˙ (1 ) (which can be solved by perturbation method). By solving all these algebraic



52

X. Chen and H.-H. Dai / Mechanics Research Communications 96 (2019) 49–55

equations, (17)4 is reduced to a vector equation only containing x˙ (0 ) . In the following, the incremental system is linearized and more explicit results can be obtained. The linearization of S˙ (0 ) − S˙ (2 ) from 13–(15) leads to



= A 1 [F˙ (0) ] + γ˙ G1 , S˙ (1) = A 1 [F˙ (1) ] + B 1 [F˙ (0) ] + γ˙ G2 , = A 1 [F˙ (2) ] + 2B 1 [F˙ (1) ] + (B 4 + B 7 )[F˙ (0) ] + γ˙ (G3 + G4 ),

S˙ (0) S˙ (2)

(18)

The incremental (or linearized incremental) plate theory obtained here is applicable to a general base state whose deformation can be inhomogeneous. For simplicity, in the next section, we consider a homogeneous base state for a plane strain problem. 5. The plane strain problem with a homogeneous base state Consider a plane strain problem in the X − Z plane and also set the base state to be homogeneous, i,e.,

where for simplicity, the variables Fb(0 ) and γ b are omitted hence-

x = x(X, Z ), y = Y, z = z(X, Z ), Fb = Fb(0) , Fb(k ) = 0, k = 1, 2.

  ⎧ (1 ) x˙ = −(E1 )−1 D 1 [∇ x˙ (0) ] + γ˙ (G1 )T e3 − q˙ + , ⎪ ⎪  1  ⎪ ⎪ (2 ) 1 −1 1 T (1 ) 2 1 (1 ) ⎪ ⎪x˙ = −(E ) (D + (H ) )[∇ x˙ ] + (E + ∇ · H )[x˙ ] ⎪   ⎪ ⎪ ⎪ −(E1 )−1 (A 1 )T [∇ ∇ x˙ (0) ] + (D 2 + ∇ · A 1 )[∇ x˙ (0) ] ⎪ ⎪     ⎪ ⎪ ⎪ −(E1 )−1 ∇ · γ˙ G1 + γ˙ (G2 )T e3 , ⎪ ⎪  ⎪ ⎨x˙ (3) = −(E1 )−1 (D 1 + (H 1 )T )[∇ x˙ (2) ] + (2E2 + ∇ · H 1 )[x˙ (2) ]  ⎪ −(E1 )−1 (A 1 )T [∇ ∇ x˙ (1) ] + (2D 2 + ∇ · A 1 ⎪ ⎪  ⎪ ⎪ ⎪ + (H 2 )T )[∇ x˙ (1) ] ⎪ ⎪   ⎪ ⎪ ⎪ −(E1 )−1 (E3 + ∇ · H 2 )[x˙ (1) ] + (B 1 )T [∇ ∇ x˙ (0) ] ⎪ ⎪  ⎪ ⎪ ⎪ −(E1 )−1 (D 3 + ∇ · B 1 )[∇ x˙ (0) ] + γ˙ (G3 + G4 )T e3 ⎪ ⎪   ⎩ + ∇ · γ˙ G2 .

The incremental quantities in Section 4 become

forth. From Eq. (17)1−3 , x˙ (1 ) − x˙ (3 ) are solved explicitly as

F˙ (k ) =

−q¯˙ =

2

2

H H H (H 1 )T [∇ x˙ (3) ] + (∇ · H 1 )[x˙ (3) ] + (A 1 )T [∇ ∇ x˙ (2) ] 6 6 6   H H2 + − ( H 1 )T + (∇ · A 1 + 2(H 2 )T ) [∇ x˙ (2) ] 2 6



+



−

 +



H H2 (∇ · H 1 ) + (∇ · H 2 ) [x˙ (2) ] 2 3

1 2 (0 ) (0 ) 1 A [F˙ , F˙ ] + A 3 [F˙ (0) , F˙ (0) , F˙ (0) ] 2 6   1 2 (0 ) (0 ) 1 1 ˙ (0 ) ˙ + γ˙ G + C [F ] + C [F , F˙ ] , 2 1 S˙ (1) = A 1 [F˙ (1) ] + A 2 [F˙ (0) , F˙ (1) ] + A 3 [F˙ (0) , F˙ (0) , F˙ (1) ] 2   + γ˙ C 1 [F˙ (1) ] + C 2 [F˙ (0) , F˙ (1) ] , 1 3 (0 ) (0 ) (2 ) A [F˙ , F˙ , F˙ ] 2 + A 2 [F˙ (1) , F˙ (1) ] + A 3 [F˙ (0) , F˙ (1) , F˙ (1) ]

S˙ (2) = A 1 [F˙ (2) ] + A 2 [F˙ (0) , F˙ (2) ] +





H H2 − ( A 1 )T + (B 1 )T [∇ ∇ x˙ (1) ] 2 3





H2 (2∇ · B 1 + (H 3 )T ) [∇ x˙ (1) ] 6

  H H2 1 2 3 + ∇ ·H − ∇ ·H + ∇ · H [x˙ (1) ] 2 6   H H2 1 T 1 T 4 7 T + (A ) − (B ) + (B + B ) [∇ ∇ x˙ (0) ] 2 6   H H2 1 1 4 7 + ∇ ·A − ∇ ·B + ∇ · (B + B ) [∇ x˙ (0) ] 2 6   H 2 H2 3 1 4 + ∇ · γ˙ G − G + ( G + G ) + O ( H 3 ), 2

6

⎪ F˙ (k ) =  F˙ 1(k ) +  2 F˙ 2(k ) +  3 F˙ 3(k ) + · · · , 0 ≤ k ≤ 2, ⎪ ⎪ ⎪ ⎩ ˙ (k ) S =  S˙ 1(k ) +  2 S˙ 2(k ) +  3 S˙ 3(k ) + · · · , 0 ≤ k ≤ 2.

(24)

Note that the superscripts denote the order of (Z − H ) in the Taylor expansions while the subscripts denote the order of  in the asymptotic expansions. The incremental deformation is caused by γ˙ , thus we should have q˙ − = q˙ + = 0. However, for the later application in Section 6, due to the existence of a Winkler foundation, some adjustments are needed. We rewrite q˙ + and q˙ − as follows. 2 ± 3 ± ˙ ˙ q˙ ± =  q˙ ± 1 +  q2 +  q3 + · · · .

(20)

where Di1jk = A31i jk , Di2jk = B31i jk , Di3jk = (B 4 + B 7 )3i jk ,

E1i j = A31i3 j , E2i j = B31i3 j , E3i j = (B 4 + B 7 )3i3 j ,

(H 1 )Ti jk = A ji13k , Hi 2jk = Bi1j3k , (H 2 )Ti jk = B 1ji3k ,

Hi 3jk = (B 4 + B 7 )i j3k ,

(22)

∂ S˙ (0)T e1 ˙ (1)T ∂ S˙ (1)T e1 ˙ (2)T +S e3 = 0, +S e3 = 0, S˙ (0)T e3 = q˙ + , ∂X ∂X

T ∂ ˙ (0 ) H ˙ (1 ) H 2 ˙ (2 ) S − S + S e1 + O(H 3 ) = −q¯˙ . (23) ∂X 2 6

⎧ γ˙ =  2 γ , ⎪ ⎪ ⎪ ⎪ ⎨x˙ (k) =  x˙ (k) +  2 x˙ (k) +  3 x˙ (k) + · · · , 0 ≤ k ≤ 3, 1 2 3



1 T

Hi 1jk = Ai1j3k ,



+ γ˙ C 1 [F˙ (2) ] + C 2 [F˙ (0) , F˙ (2) ] + C 2 [F˙ (1) , F˙ (1) ] ,

From Eq. (23)1−2 , x˙ (2 ) − x˙ (3 ) can be readily worked out. However, since (23)3 is a nonlinear algebraic equation for x˙ (1 ) , the perturbation method is applied to solve it. Assume the following expansion forms

H + (H ) − (∇ · A 1 + (H 2 )T ) [∇ x˙ (1) ] 2 +

k = 0, 1, 2,

and the incremental equations are reduced to

becomes

4 2

 e1 + x˙ (k+1)  e3 ,

S˙ (0) = A 1 [F˙ (0) ] +

(19) Eq. (17)

∂ x˙ (k) ∂X

(21)

1 (H 3 )Ti jk = (B 4 + B 7 ) ji3k , (A 1 )Ti jkl = A jikl .

Note that through (19 )1−3 , x˙ (1 ) − x˙ (3 ) can all be represented by x˙ (0 ) . Then, Eq. (19)4 can be transformed to an equation only containing x˙ (0 ) .

(25)

Remark 2. In (24), we assume that the prescribed loading parameter γ˙ is of O( 2 ) as a priori. Then, we set the induced quantities x˙ (k ) , F˙ (k ) and S˙ (k ) to be O( ) so that γ˙ will not enter the leading order equations (see (29 ) j=1 ), which should be consistent with the governing equations at the bifurcation point in order to capture the bifurcated branch. Substituting the expansions 24–(25) into (22) and (23), the following relations are obtained

F˙ (jk ) =

∂ x˙ (jk) ∂X

(k+1 )

 e1 + x˙ j

(k ) (k )  e3 , S˙ 1 = A 1 [F˙ 1 ],

X. Chen and H.-H. Dai / Mechanics Research Communications 96 (2019) 49–55

0 ≤ k ≤ 2, 1 ≤ j ≤ 3,

6. Application

⎧ 1 2 (0 ) (0 ) (0 ) 1 (0 ) 1 ⎪ ⎪ ⎨S˙ 2 = A [F˙ 2 ] + A [F˙ 1 , F˙ 1 ] + γ G , ˙ (1 )

2

˙ (1 )

˙ (0 ) ˙ (1 )

S2 = A [F2 ] + A [F1 , F1 ], ⎪ ⎪ ⎩ (2 ) S2 = A 1 [F˙ 2(2) ] + A 2 ([F˙ 1(1) , F˙ 1(1) ] + [F˙ 1(0) , F˙ 1(2) ] ), ⎧ 1 ⎪ S˙ 3(0) = A 1 [F˙ 3(0) ] + A 2 [F˙ 1(0) , F˙ 2(0) ] + A 3 [F˙ 1(0) , F˙ 1(0) , F˙ 1(0) ] ⎪ ⎪ 6 ⎪ ⎪ ⎪ ⎪ + γ C 1 [F˙ 1(0) ], ⎪ ⎪ ⎪   ⎪ ⎪ S˙ (1) = A 1 [F˙ 3(1) ] + A 2 [F˙ 1(0) , F˙ 2(1) ] + [F˙ 2(0) , F˙ 1(1) ] ⎪ ⎪ ⎨ 3 1 + A 3 [F˙ 1(0) , F˙ 1(0) , F˙ 1(1) ] + γ C 1 [F˙ 1(1) ], ⎪ 2 ⎪ ⎪ ⎪ ˙ (2) = A 1 [F˙ (2) ] + A 2 (2[F˙ (1) , F˙ (1) ] + [F˙ (0) , F˙ (2) ] + [F˙ (0) , F˙ (2) ] ) ⎪ S ⎪ 3 3 1 2 2 1 1 2 ⎪ ⎪ ⎪ 1 ⎪ ( 0 ) ( 1 ) ( 1 ) ( 0 ) ( 0 ) ( 2 ) 3 ⎪ + A (2[F˙ 1 , F˙ 1 , F˙ 1 ] + [F˙ 1 , F˙ 1 , F˙ 1 ] ) ⎪ ⎪ 2 ⎪ ⎩ + γ C 1 [F˙ 1(2) ]. 1

2

(26) From the incremental governing Eq. (23), by taking the coefficient of  j ( j = 1, 2, 3) to be zero, we have

⎧ ( 0 )T ∂ S˙ e1 ∂ S˙ (j1)T e1 ⎪ ⎪ ⎪ j + S˙ (j1)T e3 = 0, + S˙ (j2)T e3 = 0, S˙ (j0)T e3 = q˙ +j , ⎨ ∂X ∂X

T ⎪ ⎪ ∂ ˙ (0 ) H ˙ (1 ) H 2 ˙ (2 ) ⎪ Sj − Sj + Sj e1 = −q¯˙ j + O(H 3 ), ⎩ ∂X 2 6

(27) where q¯˙ j = (q˙ +j + q˙ −j )/H. One can observe that now Eq. (27)3 are linear algebraic equations for x˙ (j1 ) . So, all the unknown vectors x˙ (jk ) (1 ≤ j, k ≤ 3 ) can be solved explicitly and the results are

⎧ (k ) 1 −1 k k ˙ (0 ) k ˙ ⎪ ⎨x1 = −(E ) H ∂ x1 /∂ X , x˙ 2(k ) = −(E1 )−1 Hk ∂ k x˙ 2(0) /∂ X k − (E1 )−1 pk (X ), ⎪ ⎩ (k ) x˙ 3 = −(E1 )−1 Hk ∂ k x˙ (0) /∂ X k − (E1 )−1 hk (X ),

(28)

and pk (X) is a function of {x˙ n(m ) |1 ≤ m ≤ k, 1 ≤ n ≤ 2}\{x˙ 2(k ) } and (k )

is a function of {x˙ n |1 ≤ m ≤ k, 1 ≤ n ≤ 3}\{x˙ 3 }. Taking (28) into (27)4 , the resulting governing equations for the leading terms x˙ (j0 ) ( j = 1, 2, 3 ) are obtained as

(29)

E24 = E16 − E15 (E1 )−1 E14 , E25 = E16 (E1 )−1 E14 + E15 (E1 )−1 E21 , E

26

= E (E ) 16

1 −1

E

21

+ E (E ) 15

1 −1

X=±L

= 0.

(31)

The base state is chosen to be the homogeneous state B∗ and the deformation gradient is

Fb = Fb(0) = e1  e1 + e2  e2 + λ3 e3  e3 .

(32)

The traction-free boundary condition on the top surface leads to

0 = λ3 −

1

λ3

+

1 Nν



  1 χ  −1 μ0 1 λ3 ln 1− +1+ λ − , λ3 λ3 3 kT Nν

(33)

which provides a relation between μ0 /(kT) and λ3 . The incremental deformation is superimposed on the above base state, and for this application, we have q˙ + = 0, q˙ − = −K¯ w˙ (X, 0 )e3 = −K¯ eT x˙ (X, 0 )e3 , where K¯ = K/NkT and q˙ + and q˙ − 3

are also scaled by NkT but retain the original notations.

The critical condition for the instability is obtained by conducting the linear bifurcation analysis. Setting γ˙ = μ˙ = 0 and applying (29 ) j=1 , we have x˙ 1(0) = u0 (X )e1 + w0 (X )e3

E ,

(34)

2  3 ¯ 3 + e9 H u 0 + e10 H w0 + O (H , K H ) = 0,

where P is the integration constant. The coefficients ei (i = 1, 2, · · · ) are related to material constants and the base state variable, whose expressions are omitted. It is easy to verify that (0 ) P = S˙ 11 −

H (1 ) H 2 (2 ) 1 S˙ + S˙ + O(H 3 ) = 2 11 6 11 H

 0

H

S˙ 11 (X, Z )d Z,

(35)

and thus P is the averaged resultant. From (31), we have the boundary conditions

u0 (±L ) = 0, u1 (±L ) = 0,

where





x · e1 |X=±L = X · e1 |X=±L , e3 · ST e1 

 e4 K¯ w0 /H + e5 K¯ u0 + e6 K¯ Hw0 + e7 K¯ H 2 u 0 + e8 w0

E17 = C 1 , E18 = C 1 , E19 = C 1 , E20 = C 1 ,

∂ 2 x˙ (j0) H 25 ∂ 3 x˙ (j0) H 2 26 ∂ 4 x˙ (j0) ∂ f j (X ) + E − E + 2 6 ∂X ∂X2 ∂X3 ∂X4 3 ˙ ¯ = −q j + O(H ),

(30)

where K is the stiffness of the Winkler foundation and w(X, Z) is the vertical displacement. For the edge conditions, we propose blocked and lubricated ends, which are put as

⎪ ⎩

ij 3i3 j 3i1 j 1i3 j 1i1 j ij ij ⎪ ⎪ 21 16 15 14 1 −1 14 22 15 ⎪ E = E − ( E + E )( E ) E , E = ( E + E14 )(E1 )−1 , ⎪ ⎪ ⎪ ⎩ 23 E = −E16 (E1 )−1 E14 − E22 E21 ,

E24

q+ = 0, q− = −Kw(X, 0 )e3 ,

⎧   2  3 ⎪ ⎨e1 u0 + e2 Hw0 + e3 H u0 + O(H ) = P,

⎧ 1 H = E14 , H2 = E21 , H3 = E23 , E14 = A31i1 j , E15 = A11i3 j , ⎪ ij ij ⎪ ⎪ ⎪ 16 1 ⎪ ⎪ ⎨Ei j = A1i1 j ,

(m )

Based on the results in Section 5, the instability of a swelling polymer gel layer is studied in a plane strain setting. The dry polymer network is of shape  = [−L, L] × [0, H] and is placed on an elastomer which is regarded as a Winkler foundation here. This initial state of the system is denoted as B0 . Imposing restrictions on the two ends (X = ±L ) and putting the system into a solvent, the polymer layer starts to swell. Due to the end constraint, stress accumulates during swelling, which may cause the onset of the bifurcation at certain critical moment. The states before and after the bifurcation are respectively denoted as B∗ and Bt . The tractions on the top and bottom surfaces are set to be

6.1. Linear bifurcation analysis

where ∂ k /∂ Xk denote the kth derivative with respect to X and

hk (X)

53



H 0

S˙ 13 (±L, Z )dZ = 0,

(36)

where x˙ · e1 = u0 (X ) + u1 (X )(Z − H ) + · · · can be further simplified as

u0 (±L ) = 0, w0 (±L ) = 0, w 0 (±L ) = 0.

(37)

23

and the vector function f1 (X ) = 0 and f2 (X) is a function of x˙ 1(0 )

and f3 (X) is a function of x˙ 1(0 ) and x˙ 2(0 ) . Eq. (29) will be used in the next section to study the instability problem of a swelling polymer gel layer.

By solving the eigenvalue problem (34) and (37), the critical condition is derived as (cf. [13])

ϕ1

 n π 4 2L

− ϕ2

 nπ 2 2L

+ ϕ3 = 0,

(38)

54

X. Chen and H.-H. Dai / Mechanics Research Communications 96 (2019) 49–55

where

ϕ1 = H (e10 − e2 e9 /e1 ), ϕ2 = H K¯ (e6 − e2 e5 /e1 ) + e8 , ϕ3 = K¯ e4 /H. 2

(39)

Given the mode number n, the critical value for λ3 can be worked out from (38) and this value is denoted as λc3 . Due to the relation of μ0 /(kT) and λ3 (see (33)), we also obtain μc /(kT). The branch with the smallest value of λc3 is the most unstable one and is most likely to take place, and the mode number of this branch is denoted as nm with the eigenvalue denoted by λm . If regarding n as 3 a continuous variable, then the value of nm can be obtained by ∂λc

further letting ∂ n3 = 0. Together with Eq. (38), nm and λm are ob3 tained as



 2L ϕ − 4ϕ1 ϕ3 λm = 0, nm = π 3 2 2



 ϕ2  . 2ϕ1  m λ

(40)

that, it retains high accuracy with explicit results for the critical loading. Also, one can conduct a post-bifurcation analysis analytically, which is given in the next subsection. 6.2. Post-bifurcation analysis near critical state The critical condition of λ3 (or μ0 /(kT)) is obtained in Section 6.1. In this subsection, near the critical state, the amplitude of the leading term is determined. The incremental loading is γ˙ = μ˙ =  2 μ and from (29), three sets of governing equations are obtained x˙ (j0 ) = u0 j (X )e1 + w0 j (X )e3 , ( j = 1, 2, 3 )

⎧   2  3 ⎪ ⎨e1 u01 + e2 Hw01 + e3 H u01 + O(H ) = P1 ,

 e4 K¯ w01 /H + e5 K¯ u01 + e6 K¯ Hw01 + e7 K¯ H 2 u 01 + e8 w01

⎪ ⎩

2  3 ¯ 3 + e9 H u 01 + e10 H w01 + O (H , K H ) = 0,

3

In Fig. 1, the critical values of λ3 (from (38)) are plotted at each mode number for different values of H and K¯ , see blue solid curves. The lowest points of the curves represent the most unstable mode and related quantities λm and nm are calculated through 3 (40). For comparison, the results from the classical linear bifurcation analysis (cf. [22]) based on a full two-dimensional plane strain theory are also shown, see the red dashed curves. Good agreements are found between the blue and red curves, which indicate the validness of the current incremental plate theory. Also, smaller H and K¯ yield better results (compare (a,b,d) to (c)), and this is consistent with the derivation of the incremental plate theory, which requires H and K¯ to be small, see remainders in (34). From the comparisons between (a), (b), (c) and (d), it can also be observed that smaller H makes the gel more unstable with a larger mode number, while a smaller K¯ (softer foundation) also leads to a more unstable gel but with a smaller mode number. The advantage of the present theory compared to the classical approach is

(41)

and

⎧   2  3 ⎪ ⎨e1 u02 + e2 Hw02 + e3 H u02 + F1 (X ) + O(H ) = P2 ,

 e4 K¯ w02 /H + e5 K¯ u02 + e6 K¯ Hw02 + e7 K¯ H 2 u 02 + e8 w02

⎪ ⎩

(42)

2  3 ¯ 3 + e9 H u 02 + e10 H w02 + F2 (X ) + O (H , K H ) = 0,

and

⎧   2  3 ⎪ ⎨e1 u03 + e2 Hw03 + e3 H u03 + F3 (X ) + O(H ) = P3 ,

 e4 K¯ w03 /H + e5 K¯ u03 + e6 K¯ Hw03 + e7 K¯ H 2 u 03 + e8 w03

⎪ ⎩

(43)

2  3 ¯ 3 + e9 H u 03 + e10 H w03 + F4 (X ) + O (H , K H ) = 0,

where Pi (i = 1, 2, 3 ) are integration constants, F1 (X), F2 (X) are functions of u01 , w01 and F3 (X), F4 (X) are functions of u01 , w01 , u02 and w02 . The boundary conditions are

u0 j (±L ) = 0, w0 j (±L ) = 0, w 0 j (±L ) = 0, j = 1, 2, 3.

(44)

The solutions for the system (41) and (44 ) j=1 are

w01 = R1 cos(θ1 + b1 X ) + R2 cos(θ2 + b2 X ), u01 = ϕ0 w01 ,

(45)

where ϕ0 = −He2 /e1 and R1 , R2 , θ 1 , θ 2 are arbitrary constants and the pure imaginary numbers b1 i, b2 i are roots of the characteristic equation ϕ1 x4 + ϕ2 x2 + ϕ3 = 0. By further solving (42) with (44 ) j=2 and applying the solvability condition for (43) with (44 ) j=3 , the amplitude equations for R1 and R2 can be obtained (cf. [13])

 



R 1 L μ μ ¯ + L p P2 + L1 R21 + L2 R22 = 0,





R 2 M μ μ ¯ + M p P2 + M1 R22 + M2 R21 = 0,

(46)

where Lμ , Lp , L1 , L2 , Mμ , Mp , M1 , M2 are constants related to coefficients ei (i = 1, 2, · · · ). The results of (46) can be divided into three cases, (i) R1 = 0, R2 = 0, (ii) R1 = 0, R2 = 0, and (iii) R1 = 0, R2 = 0. Take the first case for example, the amplitude equation computed from (46) is

R21 = r1 μ ¯,

r1 =

Hμ L p − H p Lμ , H p L1 − h13 L p

(47)

Then the leading order solutions for x˙ (0 ) are

 w01 =  u

01



μ˙ r1 cos (b1 X + θ1 ),  = −ϕ0 b21 μ˙ r1 cos (b1 X + θ1 ),

where μ˙ = μ/(kT ) Fig. 1. The critical value of λ3 at each mode number for different values of H and K¯ . Blue lines are from the incremental plate theory and red dashed lines are from a full 2D analysis. The lowest points of the curves represent the most unstable modes. The parameter values are L = 1.0, Nν = 10−4 , χ = 0.1 (cf. [21]), which are also used in Figs. 2 and 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

− μc /(kT )

⎧ nπ ⎨b1 = , θ1 = k 1 π , 2L ⎩b = nπ , 1 2L

θ1 = k 2 π +

(48)

and b1 and θ 1 are

n = 2, 4, · · · , k1 = 0, ±1, ±2, · · ·

π 2

.

, n = 1, 3, · · · , k2 = 0, ±1, ±2, · · · (49)

X. Chen and H.-H. Dai / Mechanics Research Communications 96 (2019) 49–55

55

Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mechrescom.2019.03. 001.

References Fig. 2. Bifurcation diagram (loading-amplitude) for the branch with nm in Fig. 1 (c and d). Parameters other than H and K¯ are given in Fig. 1. For the incremental loading, μ˙ ∈ [0, 0.1].

Fig. 3. The profiles for the post-bifurcation deformation from Z = 0 to Z = H with step 0.2H. For parameters, H = 10−2 and K¯ = 10−1 and other parameters are given in Fig. 1. The incremental loading is μ˙ = 0.1.



The amplitude is defined as ρ = ±max|w0 | = ± μ˙ r1 . The amplitudes for cases (ii) and (iii) can be similarly obtained. The bifurcation diagram is plotted in Fig. 2 for the branch with nm in Fig. 1. By varying the values of H and K¯ , it is found that the bifurcation type is always supercritical, which implies the bifurcation of the gel system is not imperfection sensitive. Also, we observe that the stiffness K¯ value does not affect much the amplitude for the same incremental of the chemical potential. The profile of the gel in post-bifurcation state is shown in Fig. 3 from Z = 0 to Z = H with step 0.2H. Acknowledgments This work is supported by the National Natural Science Foundation of China (Grant No. 11702027) and a GRF grant from Research Grants Council of the HKSAR, China (Project No.: CityU 1130241) and SRG grant from City University of Hong Kong (Project No.: 7004670).

[1] M.A. Biot, Mechanics of Incremental Deformation, Inc., John Wiley & Sons, 1965. [2] R.W. Ogden, Non-linear Elastic Deformations, Ellis Harwood Ltd., 1984. [3] G. Kirchhoff, Über das gleichgewicht und die bewegung einer elastischen scheibe, Journal für die reine und angewandte Mathematik 40 (1850) 51–88. [4] A.E.H. Love, The small free vibrations and deformation of a thin elastic shell, Phil. Tans. R. Soc. Lond. A 179 (1888) 491–546. [5] T. von Kármán, Festigkeitsprobleme im maschinenbau, in: In encyclopädie der mathematischen wissenschaften 4, Teubner Verlag, Lepzig, Germany, 1910, pp. 311–385. [6] R.D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic elastic plates, J. Appl. Mech 18 (1951) 31–38. [7] H.L. Cox, The Buckling of Plates and Shells, Pergamon-McMillan, New York, 1963. [8] A.M.A. van der Heijden, W. T. Koiter’s Elastic Stability of Solids and Structures, Cambridge University Press, Cambridge, 2009. [9] P.G. Ciarlet, A justification of the von kármán equations, Arch. Ration. Mech. Anal. 73 (1980) 349–389. [10] G. Friesecke, R.D. James, S. Müller, A hierarchy of plate models derived from nonlinear elasticity by gamma-convergence, Arch. Ration. Mech. Anal. 180 (2006) 183–236. [11] D.J. Steigmann, Thin-plate theory for large elastic deformations, Int. J. Non-Linear Mech. 42 (2007) 233–240. [12] H.H. Dai, Z. Song, On a consistent finite-strain plate theory based on three-dimensional energy principle, Proc. R. Soc. A 470 (2014) 2014049. [13] J. Wang, D. Steigmann, F.F. Wang, H.H. Dai, On a consistent finite-strain plate theory of growth, J. Mech. Phys. Solids 111 (2018) 184–214. [14] A. Lucantonio, G. Tomassetti, A. DeSimone, Large-strain poroelastic plate theory for polymer gels with applications to swelling-induced morphing of composite plates, Composites B 115 (2017) 330–340. [15] A. Lucantonio, P. Nardinocchi, L. Teresi, Transient analysis of swelling-induced large deformations in polymer gels, J. Mech. Phys. Solids 61 (2013) 205–218. [16] P.J. Flory, Thermodynamics of high polymer solutions, J. Chem. Phys. 10 (1942) 51. [17] M. Doi, Gel dynamics, J. Phys. Soc JPN 78 (2009) 052001. [18] P.J. Flory, J. Rehner, Statistical mechanics of cross-linked polymer networks I. Rubberlike elasticity, J. Chem. Phys. 11 (1943) 512. [19] J.C. Simo, K.S. Pister, Remarks on rate constitutive equations for finite deformation problems: computational implications, Comput. Methods Appl. Mech. Eng. 46 (1984) 201–215. North-Holland. [20] L.R.G. Treloar, The Physics of Rubber Elasticity, Clarendon press, Oxford, 2005. [21] W. Hong, Z. Liu, Z. Suo, Inhomogeneous swelling of a gel in equilibrium with a solvent and mechanical load, J. Mech. Phys. Solids 46 (2009) 3282–3289. [22] D.G. Roxburgh, R.W. Ogden, Stability and vibration of pre-stressed compressible elastic plates, Int. J. Eng. Sci. 32 (1994) 427–454.