Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

Swelling-induced and controlled curving in layered gel beams A. Lucantonio, P. Nardinocchi and M. Pezzulla Proc. R. Soc. A 2014 470, 20140467, published 10 September 2014

References

This article cites 18 articles, 2 of which can be accessed free

Subject collections

Articles on similar topics can be found in the following collections

http://rspa.royalsocietypublishing.org/content/470/2171/2014046 7.full.html#ref-list-1

computational mechanics (7 articles) materials science (161 articles) mechanics (156 articles)

Email alerting service

Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here

To subscribe to Proc. R. Soc. A go to: http://rspa.royalsocietypublishing.org/subscriptions

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

Swelling-induced and controlled curving in layered gel beams rspa.royalsocietypublishing.org

A. Lucantonio1 , P. Nardinocchi2 and M. Pezzulla2 1 SISSA, International School for Advanced Studies, via Bonomea

Research Cite this article: Lucantonio A, Nardinocchi P, Pezzulla M. 2014 Swelling-induced and controlled curving in layered gel beams. Proc. R. Soc. A 470: 20140467. http://dx.doi.org/10.1098/rspa.2014.0467 Received: 13 June 2014 Accepted: 12 August 2014

Subject Areas: materials science, mechanics, computational mechanics Keywords: morphing beams, swelling-induced bending, shape control Author for correspondence: P. Nardinocchi e-mail: [email protected]

265, Trieste 34136, Italy 2 Dipartimento di Ingegneria Strutturale e Geotecnica, Sapienza Università di Roma, via Eudossiana 18, Roma 00184, Italy We describe swelling-driven curving in originally straight and non-homogeneous beams. We present and verify a structural model of swollen beams, based on a new point of view adopted to describe swelling-induced deformation processes in bilayered gel beams, that is based on the split of the swellinginduced deformation of the beam at equilibrium into two components, both depending on the elastic properties of the gel. The method allows us to: (i) determine beam stretching and curving, once assigned the characteristics of the solvent bath and of the non-homogeneous beam, and (ii) estimate the characteristics of non-homogeneous flat gel beams in such a way as to obtain, under free-swelling conditions, three-dimensional shapes. The study was pursued by means of analytical, semi-analytical and numerical tools; excellent agreement of the outcomes of the different techniques was found, thus confirming the strength of the method.

1. Introduction Soft active materials have been largely employed to realize actuators in which deformations and displacements are triggered through a wide range of external stimuli (electric field, pH, temperature and solvent absorption). The effectiveness of these actuators critically depends on the capability of achieving prescribed changes in the shape and size of the system. In particular, in gel-based actuators, the shape and the size of the structures are related to the spatial distribution of solvent inside the gel and to the magnitude of the solvent uptake. Currently, several approaches to the shape control of swellable materials are being pursued, which involve

2014 The Author(s) Published by the Royal Society. All rights reserved.

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

We thank the anonymous reviewer who suggested this situation to us.

...................................................

1

2

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

materials in the form of thin non-homogeneous sheets. As was recently demonstrated [1,2], when chemically distinct regions are present within a gel structure, free-swelling can also give rise to deformation patterns that are non-homogeneous and resemble well-known macroscopic structural deformations. In [1], an approach termed halftone gel lithography was used to photopattern a thick flat disc made of photo-cross-linkable copolymers. The patterned disc was homogeneous along the (small) thickness but non-homogeneous on the two-dimensional base disc; therefore, its free swelling gives rise to different three-dimensional shapes, depending on the two-dimensional inhomogeneity. In [2], a non-homogeneous gel system was presented based on the coupling of two different kinds of gels, with different elastic and swelling properties; the global system corresponds to a thin fibred plate-like structure which is homogeneous across its thickness and presents a non-homogeneous material structure on the two-dimensional domain. A key issue described in [1,3], which is still to be resolved, concerns the estimation of the characteristics of non-homogeneous flat gel systems in such a way as to obtain, under free-swelling conditions, three-dimensional shapes; and then, the capability to realize, through appropriate techniques, non-homogeneous flat gel systems with the estimated characteristics. We concentrated our investigations on layered gel beams, which can be seen as the simplest prototype of non-homogeneous gel beams. Indeed, layered gel beams are only non-homogeneous along their thickness and therefore represent special cases of non-homogeneous gel systems with beam-like homogeneous components. When embedded in a solvent bath, each homogeneous part would swell homogeneously, if it was free from the rest of the system; if the temperature and the chemical potential of the solvent bath were fixed, each part would have an equilibrium uniform swelling ratio that depended only on its own shear modulus. When the assembly of two homogeneous beams is considered, because of the possibly different shear moduli and, consequently, different equilibrium swelling ratios, a non-uniform equilibrium swelling is expected which results in global bending of the system. Within the context of continuum mechanics, swelling-induced deformation processes can be studied through the so-called stress-diffusion models [4,5]. We have at our disposal a fully nonlinear three-dimensional stress-diffusion model of swelling-induced large deformations in polymer gels; it was presented and discussed in [6], and subsequently applied in [7]. Our model can be implemented to find free-swelling non-homogeneous solutions of the stress-diffusion problem corresponding to non-homogeneous gels. In particular, in the case of the bilayered beam, we show how, depending on the mismatch between the shear moduli of the two beam-like parts (hence, on the corresponding mismatch between the swelling equilibrium ratios) and on the characteristics of the solvent bath, different swollen curved shapes can be obtained from a straight and dry configuration. We also show how, even if the two beam-like parts have the same shear modulus, a bending can be obtained starting from a straight beam if the enthalpy of the mixing parameter is different for the two parts.1 However, apart from the previous results, the key issue we stress here concerns the proposal of a new idea that makes the modelling of swelling-driven large deformations in non-homogeneous thin structures effective and aims to establish predictive structure–function relationships for bending bilayer gel-based composites. We assumed that, within each homogeneous beam-like part of the system, the bending deformation can be multiplicatively split into the uniform free-swelling stretch that would take place if the part was free from the rest of the beam (an ideal stretch) and a further elastic component. The uniform free-swelling stretches of each homogeneous beam are determined from the mechano-chemical equilibrium equations, which may be written within the Flory–Rehner stress-diffusion model, and depend on the shear modulus of the part, the chemical potential of the bath and the Flory parameter. The elastic deformations result from the multiplicative decomposition once the compatibility of the global bending deformation has been ensured; they deliver internal stresses within the gel beam corresponding, in the absence of external forces, to null forces and torques on each cross section of the beam. Importantly, we note that both the components of the deformation have an elastic

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

Swelling-induced deformation processes can be fully described within a stress-diffusion continuum theory based on the balance equations for forces and solvent, on the thermodynamics inequalities restricting the class of admissible constitutive prescriptions, and on the choice of a free energy which accounts for both the elastic and mixing contributions. Following the study in [6], the gel body is identified with a region B of the three-dimensional Euclidean space E, denoted as the reference configuration, with boundary ∂B of unit normal m. If a time interval I ⊂ R is chosen, the balance equations of forces and solvent are div S = 0

and c˙ = −div h , on B × I

(2.1)

− h · m = q , on ∂B × I,

(2.2)

and Sm=t

and

with S and h as the stress in the reference state (a.k.a. the first Piola–Kirchhoff stress) and the solvent flux, t and q as the boundary traction and the rate of solvent transported into B across its boundary ∂B, and c as the solvent concentration field per unit reference volume. State variables of the problem are c and the displacement field u from the reference configuration B. The constitutive equations for the stress S, the chemical potential μ and the solvent flux h come from thermodynamical issues, once the free energy is fixed, which we assume in the Flory–Rehner form [11,12]. As is well known, the Flory–Rehner free energy is defined with respect to a dry configuration Bd of the body; on the other hand, the reference configuration B is assumed to be not completely dry, it being mandatory to introduce a non-dry reference state when numerical simulations are involved, owing to the singularity of the Flory–Rehner energy in the dry state. In particular, we assume that the given reference state B is the steady state attained from Bd through the homogeneous and isotropic map fo : Bd → B; hence, Fo = ∇fo = λo I, where λo is a uniform stretch and Jo = det Fo = λ3o . With F = I + ∇u, and the solvent concentration cd = Jo c per unit dry volume, the Flory–Rehner free energy density φd per unit dry volume takes the form φd = φˆ d (Fd , cd ) =

1 RT G (Fd · Fd − 3) + h(cd ), 2 Ω

with

 h(cd ) = Ω cd log

Ω cd 1 + Ω cd

 + χ (T)

Fd = F Fo ,

Ω cd ; 1 + Ω cd

(2.3)

(2.4)

moreover, it is assumed that 1 + Ω cd = Jd = det Fd .

(2.5)

...................................................

2. Swelling-induced deformation processes in polymer gels

3

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

nature, which is different from what happens when a similar multiplicative decomposition is used to describe growth or thermal processes [8–10]; hence, both the components of the multiplicative decomposition depend on the elastic properties of the gel. We have shown that our simple bending model predicts material responsiveness in excellent agreement with the outcomes of our pilot study. In particular, the dependence of the beam curvature on the ratios between the thicknesses and between the shear moduli of the two homogeneous beam-like parts was explicitly derived. Interesting limit cases are identified and discussed, corresponding to a nearly homogeneous gel beam and a very thin gel beam placed on a stiffer beam-like substrate. The same model allowed us to deal with the inverse problem of sizing the components of the gel beam in such a way as to obtain a target shape of the bilayer system. Owing to the intrinsic nonlinearity of the mechano-chemical equilibrium equations delivering the free-swelling components of the bending deformation, the inverse problem did not admit an explicit solution. However, the key system allowed for a numerical analysis, and an asymptotic approximation which can be found through semi-analytical procedures. Distinguished and interesting characteristics of the problem are found and discussed.

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

(2.6)

with cdo the solvent concentration per unit dry volume, at the stress-free, swollen state B (the same concentration may be measured per unit reference volume as co = cdo /Jo ). We assume that the theory is consistent with a thermodynamic-like inequality and, following the approach proposed in [13] (see also [6,14]), have the following constitutive equations for the reference stress S, the chemical potential μ and the reference solvent flux h: G ˆ S = S(F) − p F = FFo FTo − p F , Jo μ = μ(c) ˆ + pΩ = and

F = (det F)F−T ,

RT  h (c) + p Ω, Ω

h (c) =

ˆ p, c) = − D c C−1 ∇(μ(c) ˆ + p Ω), h = h(F, RT



dh dcd

J = det F,

(2.7)

 cd =Jo c

C = FT F.

(2.8) (2.9)

The field p represents the osmotic pressure within the polymer, when the external pressure is assumed to be zero, and constitutively couples the stress S and the chemical potential μ, due to the incompressibility constraint (2.3)2 that can be alternatively written as J Jo = 1 + Ω Jo c. The free-swollen reference state B may be completely characterized in terms of the amount of solvent cdo absorbed, starting from the general equations (2.1) to (2.9) and assuming that the following conditions are realized: homogeneous solvent concentration cdo , uniform stretching, ˆ null loads, zero stress S(F) − pF = 0 at F = I. We obtain the swollen volume and the stretching 1/3 from Bd as Jo = 1 + Ωcdo and λo = Jo , the solvent concentration per unit of reference volume as ˆ o ) + po Ω. Under these co = cdo /Jo , the pressure as po = G/λo and the chemical potential as μo = μ(c ˆ p, c) = 0, the balance of solvent yields c˙ = 0 and the balance of forces is identically conditions, h(F, satisfied. Typically, we select a swelling stretch λo  1, only to overcome the singularity problem at λo = 1. Once the reference conditions have been defined, problems (2.1)–(2.9) admit a distinguished solution, still in the form of a free-swelling solution: the polymer is assumed to be immersed in a solvent bath of chemical potential μext = μo under load-free conditions and it can freely swell or shrink. The equilibrium solution (F∞ , c∞ ) is assumed to satisfy the balance of forces and the local thermodynamic equilibrium S(X) = 0 and

μ(X) = μext , ∀X ∈ B;

(2.10)

the difference between μ(X) and μext determines the driving force of the swelling (shrinking) process. For homogeneous materials, it holds that F∞ = λ∞ I and c∞ = with GΩ + RT λ∞

 log



λ3∞ − 1 λ3∞

 +

λ3∞ − 1 , Ω

1 λ3∞

+

χ (T) λ6∞

(2.11)  = μext ,

(2.12)

and λ∞ is the free-swelling equilibrium ratio. Once the temperature T and the external chemical potential μext have been fixed, the swelling ratio λ∞ cannot be explicitly evaluated in terms of the shear modulus G; on the contrary, the shear modulus versus swelling ratio relationship is directly expressed by equation (2.12) and is represented in figure 1a for μext = 0 and χ = 0.2, 0.3, 0.4 by the blue, orange and green lines, respectively. It is evident, and also well known, that a lower shear modulus G corresponds to a higher equilibrium swelling ratio λ∞ .

...................................................

Jo = 1 + Ω cdo ,

4

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

Here, G = NκB T is the shear modulus, N is the number of polymeric chains per unit dry volume, T is the temperature, Ω is the solvent molar volume, R is the universal gas constant and χ is the dimensionless measure of the enthalpy of mixing. Equation (2.5) expresses the incompressibility constraint of both polymer and solvent, which makes the change in volume of the gel equal to the corresponding volume of the absorbed solvent; at B, it takes the form

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

(a)

(b)

500

5

G (kPa)

x3

200

ht

c

0

2

x2 hb

x1

100 3

4 l•

5

L /2

b

Figure 1. (a) Solid and dashed lines show the relationship between the shear modulus G and the swelling ratio λ∞ at the mechano-chemical equilibrium as prescribed by equation (2.12) and by the asymptotic approximation (2.13), respectively, for different values of χ (0.2, 0.3, 0.4). (b) The two beam-like homogeneous parts of the bilayer system are shown in white and grey. The beam thickness h is the sum of the thickness of the top beam ht and the thickness of the bottom one hb ; due to symmetry, only a half beam is represented. (Online version in colour.)

In practice, GΩ/R T ∼ 10−5 –10−2 and χ ∼ 0.1–0.5 for good solvents [15]. With GΩ/RT and χ varying in such ranges, the equilibrium stretch λ∞ 1 (that is, 1/λ∞ 1) and equation (2.12) can be approximated by estimating the leading order term in the asymptotic expansion up to O(1/λ8∞ ). The following approximate shear modulus versus swelling ratio relation holds:   GΩ 1 (2.13) ∼ − χ λ−5 ∞. RT 2 The dashed lines in figure 1a show the pattern of G versus λ∞ when the equilibrium equation (2.12) is replaced by the asymptotic equilibrium approximation (2.13). It is worth noting that equation (2.13) allows for an explicit representation of the equilibrium swelling ratio λ∞ in terms of the shear modulus as   RT 1/2 − χ 1/5 . (2.14) λ∞ = Ω G Equation (2.14) has already been introduced in [5] in terms of the volume fraction of the polymer φeq = 1/λ3∞ , denoted as the equilibrium volume fraction; in [16], the approximate swelling ratio λ∞ was denoted as the hydrogel deformation constant (see also [17]).

3. Bending deformations from free-swelling solutions We considered a slightly pre-swollen gel bar B, made by the assembly of two different gel beams.2 We assumed that the top (t) beam has a shear modulus Gt smaller than the shear modulus Gb of the bottom (b) beam and set α = Gt /Gb ≤ 1. When immersed in a solvent bath with assigned chemical potential μext , the top layer of the gel beam would swell isotropically, if free, with a homogeneous equilibrium free-swelling ratio λot given by equation (2.12)) with G = Gt ; as α ≤ 1, λot is larger than the equilibrium swelling ratio λob (which can be derived from the same equation (2.12) with G = Gb ) of the bottom beam (as figure 1a shows). Actually, as the two homogeneous beam-like parts realize an assembly (the non-homogeneous gel beam), a macroscopic bending is realized; an extension of the gel beam axis is expected, too. We assumed the beam had width b and thickness h = ht + hb , with ht = β h and hb = (1 − β)h being the thicknesses of the homogeneous top and bottom beam-like parts, respectively, and the parameter β ranging between 0 and 1. Let us note that, for β = 0, the gel beam is identified by the bottom part 2 For the sake of simplicity, we limit our analysis to bilayered gel beams; the extension to multi-layered gel beams is straightforward.

...................................................

300

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

400

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

 0=

B

ˆ p, c) · ∇ c˜ dV h(F,

 and

0=

∂B

(μ(c) ˆ + p Ω − μext )˜cb dA.

(3.2)

Equations (3.1)1,2 express the weak form of the balance of forces and of the incompressibility constraint, respectively; equation (3.2)1 represents the weak form of the balance of solvent under steady-state conditions and equation (3.2)2 expresses the condition that all the boundary is exposed to a pure solvent bath, and the chemical equilibrium is maintained there at all times. The values for the initial conditions are in accordance with the assumption that the initial state is represented by the reference configuration B, which, in turn, is the steady state attained by Bd along a free-swelling process. For fixed Gt and Gb , a parametric analysis was implemented with 0 < β < 1 and h/b going from 1 to 10, that is, for each value of the parameters β and h/b, equations (3.1)–(3.2) were solved together with equations (2.7)–(2.9), with the temperature, the solvent molar volume and the dimensionless measure of the enthalpy of mixing fixed as T = 293 K, Ω = 6.023 × 10−5 m3 mol−1 and χ = 0.2. For any values of the parameters, a uniform bending solution was obtained, corresponding to a plane bending occurring in the flexure plane x2 = 0 (figure 1b). Any solution delivers uniform beam axis curvature κ and longitudinal axis stretch Λo , evaluated as κ=

w (1 + u ) − w u ((1 + u )2 + w2 )3/2

and Λo = (C11 (x1 , 0, 0))1/2 ,

with C11 = Ce1 · e1 ,

(3.3)

respectively, being u(x1 ) = u(x1 , 0, 0) · e1 , w(x1 ) = u(x1 , 0, 0) · e3 , and v(x1 ) = u(x1 , 0, 0) · e2 ≡ 0. 1/2 Moreover, we also found that ∂ C11 /∂ x3 was almost a constant field (apart from boundary layers), that is, cross sections stay flat as standard in a bending process. Figure 2 shows the outcomes of our numerical experiment corresponding to Gt  15 kPa and Gb  45 kPa (α  0.3) and aspect ratios h/b = 1, 2, 4. The free-swelling ratios λot and λob are expected to be 3.8 and 3.1, respectively, with a free-swelling mismatch o = λot − λob = 0.7. We show the pattern of the longitudinal axis stretch Λo and the curvature κ versus β. As expected, for β  0, that is, for a homogeneous gel beam with shear modulus Gb , the longitudinal stretch of the beam axis is equal to the free-swelling ratio λob and the beam curvature is zero. Likewise, for β  1, that

...................................................

and

6

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

and is homogeneous with shear modulus G = Gb , whereas for β = 1 the gel beam is identified by the top part and is homogeneous with shear modulus Gt . The length of the gel beam is denoted as L and fixed to L = 10 h; we also fix the width b of the gel beam as b = 1 mm and investigated beams with different aspect ratios h/b. As a reference configuration for the numerical analysis, we used an almost dry straight configuration of the gel beam corresponding to swelling ratios of the two components around 1.001, with a solvent bath chemical potential μo < 0, and compatible with the shear moduli of the two homogeneous beam-like parts according to the free-swelling equilibrium equation (2.12). In particular, for Gt  15 kPa and Gb  45 kPa (α = 0.3), we obtained μo  −11 × 103 J mol−1 . The beam was constrained so as to allow the free-swelling change in shape; hence, because of the symmetry of the problem, we hampered the longitudinal displacement over the entire middle cross section, the out-of-plane displacement along the line spanned by the coordinates (0, 0, x3 ) for −h/2 < x3 < h/2, and the transverse displacement at the point (0, 0, −h/2). The remaining boundary was assumed to be traction-free (t = 0). We allowed the sample to freely swell in a solvent bath whose chemical potential μext was controlled and set to 0. The model was implemented and subsequently solved numerically in the finite element software COMSOL Multiphysics [6]. Equations (2.1)–(2.9) were recast in a weak form, using a mixed method to impose both the volumetric constraint on B and the boundary conditions on ∂B; thus, the problem was reformulated as follows: find u, p, c and a boundary concentration cb ˜ p˜ , c˜ and c˜ b , with u˜ compatible with the boundary conditions, and c˜ = c˜ b on such that, for any u, ∂B, it holds that     1 + ΩJo c ˆ J− p˜ dV (3.1) −(S(F) − p F∗ ) · ∇ u˜ dV and 0 = 0= Jo B B

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

(a)

(b)

3.8

k (m–1)

L0

3.2

60 40 20

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

b

0.6

0.8

1.0

b

Figure 2. Longitudinal axis stretch Λo (a) and beam axis curvature κ (b) versus β for α = 0.3 and h/b = 1 (circles), 2 (triangles), 4 (squares). (Online version in colour.)

(a)

(b)

3.8

100 80 k (m–1)

L0

3.6 3.4 3.2 3.0

60 40 20

0

0.2

0.4

0.6

0.8

1.0

0

b

0.2

0.4

0.6

0.8

1.0

b

Figure 3. Solid lines show the patterns of Λ0 (a) and κ (b) versus β, for different aspect ratios h/b (= 1, 2, 4) of the cross section of the beam, with μext , Gt and Gb fixed. At β = 0 and β = 1, we have a homogeneous gel beam (as h = hb and h = ht , respectively); for such special β values, Λo is determined by λob and λot , respectively, whereas κ is null. Dashed lines show the corresponding patterns when the asymptotic approximation of the mechano-chemical equilibrium is implemented. (Online version in colour.)

is, for a homogeneous gel beam with shear modulus Gt , the longitudinal stretch is equal to the free-swelling ratio λot and the beam curvature is still zero. For β in between 0 and 1, the gel beam is non-homogeneous, and free-swelling determines a longitudinal stretch in between λot and λob and a beam curvature which attains its maximum value for β  0.6.

4. Structural deformations in bilayer gel beams It would have been interesting to study the change in beam curvature due to changes in the moduli ratio α; the goal cannot be straightforwardly realized starting from the computational model. Indeed, changing α means changing at least one of the shear moduli Gt and Gb ; hence, this would mean looking for a new almost dry and straight reference configuration, which, all in all, would be a time-consuming procedure. On the other hand, the outcomes of the numerical

...................................................

3.4

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

80

3.6

3.0

7

100

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

(4.1)

where Λ0 and κ are the possibly large longitudinal stretch and curvature of the beam axis and κ > 0.4 It is assumed that the top and bottom homogeneous beam-like parts suffered elastic longitudinal deformations λet and λeb such that λet = λλ−1 ot

and λeb = λλ−1 ob .

(4.2)

Then, assuming zero out-of-plane stresses, the corresponding longitudinal stresses σt and σb on the cross sections of the top and bottom layers, respectively, are evaluated as σt (x1 , x3 ) = 3Gt (λ(x1 , x3 )λ−1 ot − 1),



 h h − β h < x3 < 2 2

(4.3)

and σb (x1 , x3 ) = 3Gb (λ(x1 , x3 )λ−1 ob − 1),

h − < x3 < 2



 h −βh , 2

(4.4)

where, because of the material incompressibility, 3Gt and 3Gb identify the corresponding Young moduli.

It should be λ(x1 , x3 ) = Λ0 (x1 )(1 + y3 κ(x1 )), where y3 is the coordinate along the beam thickness corresponding to x3 in the deformation process; we assumed that y3 = Λo x3 (isotropic stretch).

3

4

¯ t, Indeed, for the symmetry of the problem, if the pair (Λo , κ) is known for a beam with ht = β¯ h and shear moduli Gt = G

¯ b , then the pair (Λ0 , −κ) corresponds to a beam with ht = β h and shear moduli Gt = G ¯ b , Gb = G ¯ t and β = 1 − β. ¯ Gb = G

...................................................

λ(x1 , x3 ) = Λ0 (x1 )(1 + x3 Λ0 (x1 )κ(x1 )),

8

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

experiments appeared so simple (plane bending, flat cross sections, uniform curvatures) that we looked for an alternative approach to the problem which allowed us to easily establish predictive structure–function relationships. We propose a new point of view in the analysis of free-swelling solutions of special gel beamlike bodies B, with homogeneous beam-like parts, which may be easily extended to the analysis of free-swelling-induced deformation patterns in multi-layered beams such as in plate-like layered systems. We assumed that, within each homogeneous beam-like part of the system, the bending deformation can be multiplicatively split into the uniform free-swelling ratio that would take place if the part were free from the rest of the beam and a further elastic component. The uniform free-swelling components of each part are determined from the appropriate mechano-chemical equilibrium equations which can be written within the Flory–Rehner stress-diffusion model in terms of the shear modulus, the chemical potential of the bath and the Flory parameter. The elastic deformations result from the multiplicative decomposition once the global compatibility of the bending deformation has been ensured; they deliver internal stresses within the gel beam corresponding, in the absence of external forces, to null forces and torques on each cross section of the beam. Both the deformation components have an elastic nature, and hence depend on the shear moduli of the materials. Owing to the numerical evidence of the plane bending processes, we moved to a twodimensional scenario: given the coordinate system (x1 , x2 , x3 ) in the reference configuration, we assumed that the flexure plane x2 = 0 is a symmetry plane (see figure 1b). Bending deformations were then modelled assuming that (i) there exists a swelling ratio Λ0 (x1 ) > 0 in between λot and λob that determines an isotropic stretch at x1 and (ii) the beam’s cross sections stay flat. With this, the measure λ(x1 , x3 ) of the longitudinal deformation that takes place within the gel beam can be represented as3

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

and

 h/2−βh b

−h/2

h/2−βh

x3 σb (x1 , x3 ) dx3 + b

 h/2 h/2−βh

x3 σt (x1 , x3 ) dx3 = 0.

(4.6)

Equations (4.5) and (4.6) deliver the following linear systems of two equations in Λ0 and Λ1 = κΛ20 :     Gt Gb Gt Gb Λ0 + S t Λ1 = At Gt + Ab Gb At + Ab + Sb (4.7) λot λob λot λob and

 St

Gt Gb + Sb λot λob



  Gt Gb Λ0 + I t Λ1 = St Gt + Sb Gb , + Ib λot λob

(4.8)

with the pairs (Gt , λot ) and (Gb , λob ) satisfying the local thermodynamic equilibrium for the beamlike layers     λ3ot − 1 Gt Ω χ (T) 1 + R T log (4.9) μˇ t (Gt , λot ) = + 3 + 6 = μext λot λ3ot λot λot and μˇ b (Gb , λob ) =

Gb Ω + RT λob



 log

λ3ob − 1 λ3ob

 +

1 λ3ob

+

χ (T) λ6ob

 = μext .

(4.10)

Moreover, the pairs (At , Ab ), (St , Sb ) and (It , Ib ) measure the area and the static and inertia moments with respect to the x2 -axis of the homogeneous top and bottom beams, respectively. For fixed b, external chemical potential μext and Flory parameter χ , the pair (Λ0 , Λ1 ) may be explicitly expressed in terms of the free-swelling ratios λot and λob and β through equations (4.7)–(4.10) as 2 6 6 Λ0 (λot , λob , β) = (β 4 λot λ12 ob f (λot ) − β(β − 1) (a(β)λob + d(β)λot )λot λob f (λot )f (λob ) 2 −1 + (β − 1)4 λ12 ot λob f (λob ) ) F(λot , λob , β) ,

(4.11)

Λ1 (λot , λob , β) = 6β(1 − β)(λot − λob )λ6ot λ6ob f (λot )f (λob )h−1 F(λot , λob , β)−1 ,

(4.12)

being 2 6 6 4 12 2 F(λot , λob ) = β 4 λ12 ob f (λot ) − c(β)λob λot f (λot )f (λob ) + (β − 1) λot f (λob )

and

  1 f (λo ) = χ + λ3o + λ6o log 1 − 3 , λo

(4.13)

with a(β) = 3 + β(β − 3), c(β) = 2β(β − 1)(2 + β(β − 1)) and d(β) = β 2 + β + 1. Following the previous preliminary assumptions, we suppose that α = Gt /Gb ≤ 1, hence o ≥ 0, and investigate equations (4.11) and (4.12). The solid lines in figure 3 show the patterns of Λ0 (figure 3a) and κ (figure 3b) versus β, for λot = 3.8 and λob = 3.1, that is, from equations (4.9) and (4.10) with μext = 0, for Gt  15 kPa and Gb  45 kPa (hence, α = 0.3), and for an aspect ratio h/b = 1 (blue circles), 2 (green triangles), 4 (orange squares). The corresponding coloured symbols identify the longitudinal stretch and the curvature of the beam axis, as derived from the numerical simulations; the agreement between the outcomes of the fully nonlinear threedimensional stress diffusion model and the structural simplified model looks excellent. Moreover,

...................................................

−h/2

9

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

We looked for free-swelling solutions to the gel beam problem with uniform Λ0 and κ; under no external load, the total force and the total moment on the gel beam have to be identically null, that is, it holds that  h/2  h/2−β h σb (x1 , x3 ) dx3 + b σt (x1 , x3 ) dx3 = 0 (4.5) b

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

(a)

(b)

3.2

k (m–1)

L0

2.9

30 20 10

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

b

0.6

0.8

1.0

b

Figure 4. Longitudinal axis stretch Λo (a) and beam axis curvature κ (b) versus β for α = 0.3 and h/b = 1 (circles), 2 (triangles), 4 (squares); the shear modulus is the same in the two homogeneous beams Gt = Gb = 44 KPa, whereas χt = 0.2 and χb = 0.3. Solid and dashed lines show the corresponding outcomes of the structural analytical model when equation (2.12) and equation (2.14) are used, respectively. (Online version in colour.)

the longitudinal stretch Λ0 is almost unaffected by the aspect ratio, whereas κ scales almost linearly with it. Remark. We also show how, even if the two beam-like parts have the same shear modulus, a bending can be obtained starting from a straight beam if the enthalpy of the mixing parameter is different for the two parts. Figure 4 shows the longitudinal axis stretch Λo (a) and beam axis curvature κ (b) versus β for α = 0.3 and h/b = 1 (circles), 2 (triangles), 4 (squares), for a bilayered beam when the enthalpy of mixing parameters are χt = 0.2 and χb = 0.3 for the top and bottom layer, respectively. Solid and dashed lines show the corresponding outcomes of the structural analytical model when equation (2.12) and equation (2.14) are used, respectively.

(a) Asymptotics When equations (4.9) and (4.10) are replaced by their asymptotic versions (as shown for the exact and the corresponding asymptotic versions (2.12) and (2.13)), an approximated solution (Λ˜ 0 , Λ˜ 1 ) can be obtained in the simpler form Λ˜ 0 (α, β, λot ) = fo (α, β) λot

and Λ˜ 1 (α, β, Gt ) = f1 (α, β)

o (α, Gt ) , h

(4.14)

with fo (α, β) = and

f1 (α, β) =

β 4 − β(β − 1) (a(β) + d(β)α −1/5 )α −1 + (β − 1)4 α −11/5 β 4 − c(β)α −6/5 + (β − 1)4 α −12/5 β4

6β(1 − β)α −6/5 . − c(β)α −6/5 + (β − 1)4 α −12/5

(4.15) (4.16)

The dashed lines in figure 3 show the patterns of Λ˜ 0 (figure 3a) and κ˜ = Λ˜ 1 /Λ˜ 20 (figure 3b) versus β, determined by equations (4.15) and (4.16). As can be appreciated, the exact and the approximate solutions derived from the asymptotic equilibrium are in extremely good agreement. We also note that the approximate longitudinal stretch Λ˜ o can be expressed as a fraction of the larger swelling ratio λot through a function fo , only depending on the moduli and thickness ratios α and β; moreover, it holds that fo (α, β) < 1 for 0 < α < 1 and 0 < β < 1, and fo (α, 1) ≡ 1. Likewise, Λ˜ 1 can be represented as a fraction of the ratio o /h through a function f1 , only depending on α and β, and with f1 (α, β) > 0 for 0 < α < 1 and 0 < β < 1, and f1 (0, β) = f1 (1, β) ≡ 0; in this sense, with

...................................................

3.0

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

40

3.1

2.8

10

50

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

(a)

(b) 4

k~ (m–1)

~

L0

2 1

0

60 40

b

20

0.5 a

1.0

0

0.5 a

1.0

Figure 5. Solid, dashed and dashed–dotted lines show the patterns of Λ˜ 0 (a) and κ˜ (b) versus α, for an aspect ratio h/b = 2 and for different thickness ratios β = 0.2, 0.5, 0.8, respectively. (a) For β → 1, the longitudinal axis stretch Λ˜ 0 always moves towards λot regardless of the value of α, as the beam is almost determined by the top part. (b) At α = 0 and α = 1, the curvature of the beam κ˜ is always equal to zero, due to either the infinite stiffness of the bottom layer that hampers the bending (α = 0) or the homogeneity of the beam that freely swells without bending (α = 1). (Online version in colour.)

reference to a thermal analogy, o /h may be thought of as the curvature induced by the mismatch between the free-swelling ratios of the top and bottom beams. An analogous thermal analogy motivated the equation derived in [18] for the swelling-induced curvature of the pine cone model; therein, the theory of bimetallic thermostats was used and swelling deformations were viewed as hygrometric expansions, which depended linearly on the change in humidity of the environment through hygrometric coefficients that are unaffected by the elastic moduli of the materials. In our case, the free-swelling ratios λot and λob strictly depend on the shear moduli of the gel (see equation (2.14)), so determining a more complex dependence of the curvature κ on the elastic characteristics of the two layers which are also behind the swelling ratio mismatch. ˜ β, Gt ) (figure 5b) as functions of the moduli In figure 5, we plot Λ˜ o (α, β, λot ) (figure 5a) and κ(α, ratio α for different values of β between 0 and 1, for a fixed aspect ratio h/b = 2 and shear modulus Gt = 15 kPa of the top layer.5 We note that, for β → 1, the longitudinal axis stretch Λ˜ o always moves towards λot regardless of the value of α, as the beam is almost determined by the top part. For β → 0, the beam is almost determined by the bottom part and Λ˜ o is identified with λot only at α = 1, that is, when Gb = Gt /α = Gt and the beam is homogeneous. In figure 5b, we note that the curvature of the beam is equal to zero at α = 0 and α = 1, for any β. This is due to either the infinite stiffness of the bottom layer that hampers the bending (α = 0) or the homogeneity of the beam that freely swells without bending (α = 1). For any α between 0 and 1, the curvature is different from zero and strongly depends on β. In particular, let us focus on the dashed–dotted line corresponding to β → 1: for α  0, we obtain a large curvature as the higher stiffness of the bottom component increases o and the higher thickness of the top beam determines the bending; moreover, as α → 1, the curvature reduces even if the bottom beam is becoming softer and softer as o is decreasing. As β → 0, the pattern of curvature is smoothed out and attains its maximum at values of α that are inversely proportional to β (see, for example, the solid line in figure 5b); it appears that o has a weaker influence on the beam’s curvature. Remark. This last situation—a thin soft layer bonded to a thick rigid one—might resemble other cases studied in the literature, concerning the deformation modes of a soft thin layer bonded to a more rigid substrate. In [19], the swelling of a soft thin layer bonded to a fixed rigid substrate In our case, o is an elastic deformation and depends on the elastic properties of the beams; hence, we always need to fix two parameters such as Gt and α to determine it.

5

...................................................

80 b

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

3

11

100

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

(b)

0

–50 –50

0 x1 + u (mm)

50

12

50

0

–50 –50

0 x1 + u (mm)

50

Figure 6. Solid lines represent the two-dimensional view of the top, middle and bottom surfaces for α = 0.03 (a) and α = 0.3 (b), β = 0.05, and h/b = 2. A wrinkling pattern is clearly evidenced, involving only the top beam, at α = 0.03. (Online version in colour.)

was studied; it was shown that wrinkling of the thin layer is an admissible deformation mode of the system, which did not bend due to the constraints of the substrate. A similar problem can also be found in [20], where bending of a homogeneous beam swollen by a solvent droplet put on the top surface is studied; it was shown that, at equilibrium, the beam recovers the straight initial configuration, whereas, during the transient, bending can only occur if the thickness of the beam is below the so-called elastoswelling length scale. Otherwise, wrinkling of the swollen top surface appears. In our opinion, we might have an analogous situation: for α  0 and β  0, the larger swelling of the top beam cannot drive a global bending and wrinkling might occur, only involving the top soft thin layer. With the aim of investigating the problem further, we ran a numerical free-swelling experiment considering a beam with a very thin top component (β = 0.05) much softer than the bottom one (Gt = 1.5 kPa and α = 0.03) and all the other parameters fixed (as in the previous conditions). We looked for the bending deformation process; figure 5b suggests that, for β → 0, the curvature is almost null and varies slightly with α. Interestingly, a wrinkling pattern was evidenced involving only the top surface of the beam. Figure 6a shows that the beam is only slightly bent: the bottom insert corresponds to the three-dimensional deformed configuration with the colour code denoting the vertical displacement w; green, black and red curves correspond to a view of the deformed shapes of the bottom, middle and top surfaces of the beam in the flexure plane; top inserts show the wrinkling pattern of the top surface. Therefore, the beam presents a richer deformation pattern, with the global bending almost null and the top layer undergoing a bending with small wavelengths (wrinkling). Figure 6b shows the same results obtained for α = 0.3. As expected, the curvature is very close to the one realized for α = 0.03 but no wrinkling pattern is now evidenced.

5. Shaping problem in gel beams Shape control of swelling materials is a recent central topic in the literature dealing with soft active materials; it typically involves materials in the form of thin non-homogeneous sheets [1–3]. With reference to the bilayer gel beam, the shaping problem consists in the determination of the geometrical and material characteristics of the homogeneous beam-like parts, allowing for a fixed uniform longitudinal axis stretch and beam axis curvature. We propose a solution to the problem, set within the context of the structural model discussed in §4.

...................................................

x3 + w (mm)

50

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

x3 + w (mm)

(a)

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

and

λob = λˇ ob (Gt , Gb , Λ0 , Λ1 )

(5.1)

was plugged into the chemical equilibrium equations (4.9) and (4.10) for the two homogeneous beam-like parts of the system. The new equations μ¯ t (Gt , Gb , Λ0 , Λ1 ) = μˇ t (Gt , λot (Gt , Gb , Λ0 , Λ1 )) = μext

(5.2)

μ¯ b (Gt , Gb , Λ0 , Λ1 ) = μˇ b (Gb , λob (Gt , Gb , Λ0 , Λ1 )) = μext

(5.3)

and

were highly nonlinear in the pair (Gt , Gb ) and an exact explicit solution cannot be derived. However, the system allowed for a numerical analysis and its asymptotic counterpart can be solved through semi-analytical procedures. Distinguishing and interesting characteristics of the problem were found and are discussed below.

(a) The asymptotic shaping problem With reference to the asymptotic solution (2.14) of the equilibrium equation for μext = 0, we set   1/5 RT 1 1 1 1 = 3 −χ and λob = with . (5.4) λot = Γ Ω 2 Γ (3Gt )1/5 Γ (3Gb )1/5 By plugging equations (5.4) into the balance equations of force and torque, we obtained the following nonlinear system in the pair (Gt , Gb ) 3At Gt + 3Ab Gb − Γ (St Λ1 + At Λ0 )(3Gt )6/5 − Γ (Sb Λ1 + Ab Λ0 )(3Gb )6/5 = 0

(5.5)

3St Gt + 3Sb Gb − Γ (It Λ1 + St Λ0 )(3Gt )6/5 − Γ (Ib Λ1 + Sb Λ0 )(3Gb )6/5 = 0,

(5.6)

and

or, equivalently, 3α1 Gt + 3β1 Gb + γ1 (3Gb )6/5 = 0

(5.7)

3α2 Gt + 3β2 Gb + γ2 (3Gt )6/5 = 0,

(5.8)

and

by combining linearly equations (5.5) and (5.6), and setting

and

Γ α1 =

St At − , It Λ1 + St Λ0 St Λ1 + At Λ0

Γ α2 =

St At − , Ib Λ1 + Sb Λ0 Sb Λ1 + Ab Λ0

Γ β1 =

Sb Ab − , It Λ1 + St Λ0 St Λ1 + At Λ0

Γ β2 =

Sb Ab − Ib Λ1 + Sb Λ0 Sb Λ1 + Ab Λ0

γ2 =

St Λ1 + At Λ0 It Λ1 + St Λ0 − . Sb Λ1 + Ab Λ0 Ib Λ1 + Sb Λ0

γ1 =

Sb Λ1 + Ab Λo Ib Λ1 + Sb Λo − , St Λ1 + At Λ0 It Λ1 + St Λ0

(5.9)

In the following, it is assumed, without loss in generality, that Λ1 > 0. A numerical analysis of the key system (5.2) and (5.3) shows that there exists a critical value κ c of the curvature, depending on β, h/b and Λ0 , above which the system admits two real distinct solutions. These are represented in figure 7a for β = 0.5, h/b = 2 and Λ0 = 3.386, in a semi-log graph: solid and dashed lines identify Gb and Gt , respectively (since α < 1 for hypothesis, Gt < Gb ); each colour identifies a couple (Gt ,Gb ). The critical curvature κ c is around 50 m−1 . Consequently, a target curvature κ > κ c can be realized by choosing between two different couples (Gt , Gb ), depicted in red and green (the values of Gb for the two different couples are very close to each other and cannot be distinguished

...................................................

λot = λˇ to (Gt , Gb , Λ0 , Λ1 )

13

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

We solved equations (4.7)–(4.10) with respect to the pair (Gt , Gb ) of shear moduli for different κ, with a fixed aspect ratio h/b and thickness ratio β, and a selected longitudinal stretch Λo . In a first step, we solved the balance equations (4.7) and (4.8) of force and torque with respect to the pair (λto , λob ); then, the solution

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

(a)

(b)

14

25 R

Gt (kPa)

shear moduli (kPa)

20

1

15 10

1 k 0

5 0

100

50 k (m–1)

0

50 Gb (kPa)

50

100

Figure 7. (a) Solid and dashed lines identify Gb and Gt , respectively, corresponding to β = 0.5, h/b = 2 and Λo = 3.386; the critical curvature κ c is around 50 m−1 . (b) Behaviour of the solutions of the shaping problem in the plane (Gt , Gb ). Grey lines are κ-isolines for the key system (5.2) and (5.3), whereas symbols represent the semi-analytical solutions. (Online version in colour.)

in the plot). A target curvature κ < κ c can be realized by choosing among four different couples (Gt , Gb ), depicted in blue, orange, red and green. Again, the values of Gb are very close to one another and cannot be distinguished graphically. Figure 7b shows the behaviour of the solutions of the system in the plane (Gt , Gb ) as κ varies. Grey lines are κ-isolines for the key system (5.2) and (5.3) corresponding to κ = 40 m−1 (light grey) and κ = 60 m−1 (dark grey). The solid coloured lines represent the traces of the four solutions (Gt , Gb ) as κ varies. The grey dashed–dotted line corresponds to Gt = Gb : all the solutions stay below it since α < 1.

(b) Semi-analytical solutions We manipulated the system (5.7) and (5.8) to obtain the new equations6 3 Gt + 3Gb +

γ1 (3Gb )6/5 = 0 β1

and 3Gt + 3δ Gb +

γ2 (3Gt )6/5 = 0, α2

(5.10)

with δ = β2 /α2 and  = α1 /β1 the two smallness parameters depending on β, h/b, Λo and κ. With this, we employed the expansions Gb ∼ G0b + G1b in (5.10)1 and Gt ∼ G0t + δG1t in (5.10)2 , obtaining the asymptotic system   1 α2 5 + 5δGb (5.11) Gt ∼ − 3 γ2 and Gb ∼ −

1 3



β1 γ1

5

+ 5Gt ,

(5.12)

which can be solved to give Gt ∼

1 5δ(β1 /γ1 )5 + (α2 /γ2 )5 3 25δ − 1

(5.13)

Gb ∼

1 5(α2 /γ2 )5 + (β1 /γ1 )5 . 3 25δ − 1

(5.14)

and

We divided equation (5.7) by β1 and equation (5.8) by α2 which are both different from 0 within the ranges considered for Λo and κ (with β < 0.9).

6

...................................................

k

R

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

R

102

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

, and we employed the asymptotic relation for Gt to

(k+1) Gb

update its estimate using the value of just found,   γ1 6 (k+1) (k) (k) (k) (k+1) (k) (Gb )6/5 + (Gb 31/5 − Gb )(Gb )1/5 + Gb + Gt = 0 β1 5

(5.15)

and (k+1)

Gt

=−

1 3



α2 γ2

5

(k+1)

+ 5δGb (0)

.

(0)

(5.16) (0)

As initial values of the iterative process, we choose Gb = −Gt and Gt computed from (5.13). We represent the outcomes of just one iteration by orange triangles (Gb ) and circles (Gt ) in figure 7a; the larger value of Gt among the smallest was found (the orange circles overwrite the dashed orange line) together with the solution for Gb . Symbols in figure 7b show the two semi-analytical solutions. It is worth noting that the above analysis, which has been illustrated and discussed for specific values of Λo , β and h/b, is strong enough to hold when different values are considered. Indeed, we verified that, even if the number of solutions varies, nevertheless the semi-analytical solutions still hold.

6. Conclusion We described swelling-induced curving in originally straight non-homogeneous gel beams immersed in a solvent bath; in particular, we focused the present work on bilayered gel beams, which can be considered as a simple non-trivial prototype of a non-homogeneous along the thickness gel beam, composed of only two homogeneous beam-like parts. We presented and verified a structural model of swollen beams that was based on a new point of view adopted to describe swelling-induced deformation processes in bilayered gel beams, that is, the split of the swelling-induced deformation of the beam at equilibrium into two components. We proposed two simple formulae for beam stretching and curvature, which depend on the solvent bath and on key material parameters such as the ratios between shear moduli and layer thicknesses. The formulae work for extremely large curving and stretching of the gel beam. This study was pursued through analytical, semi-analytical and numerical tools; an excellent agreement of the outcomes of the different techniques was found, thus revealing the strength of the method. In particular, the analysis showed interesting behaviours of the system, such as the absence of bending and rising of wrinkling that take place for special values of the moduli and thickness ratios, which will be investigated in the future. Funding statement. This work was supported by MIUR (the Italian Ministry of University and Research) through Project PRIN2009 ‘Mathematics and Mechanics of Biological Assemblies and Soft Tissues’.

References 1. Kim J, Hanna JA, Byun M, Santangelo CD, Hayward RC. 2012 Designing responsive buckled surfaces by halftone gel lithography. Science 335, 1201–1205. (doi:10.1126/science.1215309)

...................................................

(k+1)

obtain the new approximate solution Gb

15

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

For fixed β, h/b and Λo , the fields Gt and Gb defined by equations (5.13) and (5.14) are represented as functions of κ in figure 7a by blue triangles (Gb ) and circles (Gt ), respectively. It is evident that the asymptotic procedure cannot find the solutions with smaller Gt (orange, red and green dashed lines); indeed, for small Gt , the first two terms in equation (5.10)1 may be comparable, and the successive asymptotic solution does not hold any more. To find the other non-trivial roots of the original system, when Gt Gb , and κ < κ c we used (k) an iterative approach. We linearized equation (5.7) about the solution Gb at the kth iteration to

Downloaded from rspa.royalsocietypublishing.org on October 31, 2014

16 ...................................................

rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20140467

2. Wu ZL, Moshe M, Greener J, Therien-Aubin H, Nie Z, Sharon E, Kumacheva E. 2013 Threedimensional shape transformations of hydrogel sheets induced by small-scale modulation of internal stresses. Nat. Commun. 4, 1586. (doi:10.1038/ncomms2549) 3. Dias MA, Hanna JA, Santangelo CD. 2011 Programmed buckling by controlled lateral swelling in a thin elastic sheet. Phys. Rev. E 84, 036603. (doi:10.1103/PhysRevE.84.036603) 4. Hong W, Zhao X, Zhou J, Suo Z. 2008 A theory of coupled diffusion and large deformation in polymeric gels. J. Mech. Phys. Solids 56, 1779–1793. (doi:10.1016/j.jmps.2007.11.010) 5. Doi M. 2009 Gel dynamics. J. Phys. Soc. Jpn 78, 052001. (doi:10.1143/JPSJ.78.052001) 6. Lucantonio A, Nardinocchi P, Teresi L. 2013 Transient analysis of swelling-induced large deformations in polymer gels. J. Mech. Phys. Solids 61, 205–218. (doi:10.1016/j.jmps. 2012.07.010) 7. Lucantonio A, Roché M, Nardinocchi P, Stone HA. 2014 Buckling dynamics of a solventstimulated stretched elastomeric sheet. Soft Matter 10, 2800–2804. (doi:10.1039/c3sm52941j) 8. Rodriguez EK, Hoger A, McCulloch AD. 1994 Stress-dependent finite growth in soft elastic tissues. J. Biomech. 27, 455–457. (doi:10.1016/0021-9290(94)90021-3) 9. Di Carlo A, Quiligotti S. 2002 Growth and balance. Mech. Res. Commun. 29, 449–456. (doi:10.1016/S0093-6413(02)00297-5) 10. Nardinocchi P, Teresi L. 2007 On the active response of soft living tissues. J. Elast. 88, 27–39. (doi:10.1007/s10659-007-9111-7) 11. Flory PJ, Rehner J. 1943 Statistical mechanics of cross-linked polymer networks. I. Rubberlike elasticity. J. Chem. Phys. 11, 512–520. (doi:10.1063/1.1723791) 12. Flory PJ, Rehner J. 1943 Statistical mechanics of cross-linked polymer networks. I. Swelling. J. Chem. Phys. 11, 521–526. (doi:10.1063/1.1723792) 13. Coleman BD, Noll W. 1963 The thermodynamics of elastic materials with heat conduction and viscosity. Arch. Ration. Mech. Anal. 13, 167–178. (doi:10.1007/BF01262690) 14. Gurtin ME, Fried E, Anand L. 2010 The mechanics and thermodynamics of continua. Cambridge, UK: Cambridge University Press. 15. Mark JE. 1999 Polymer data handbook. New York, NY: Oxford University Press. 16. Dai H, Song Z. 2011 Some analytical formulas for the equilibrium states of a swollen hydrogel shell. Soft Matter 7, 8473–8483. (doi:10.1039/C1SM05425B) 17. Lucantonio A, Nardinocchi P, Stone HA. 2014 Swelling dynamics of a thin elastomeric sheet under uniaxial pre-stretch. J. Appl. Phys. 115, 083505. (doi:10.1063/1.4866576) 18. Reyssat E, Mahadevan L. 2009 Hygromorphs: from pine cones to biomimetic bilayers. J. R. Soc. Interface 6, 951–957. (doi:10.1098/rsif.2009.0184) 19. Kang MK, Huang R. 2010 Swell-induced surface instability of confined hydrogel layers on substrates. J. Mech. Phys. Solids 58, 1582–1598. (doi:10.1016/j.jmps.2010.07.008) 20. Pandey A, Holmes DP. 2013 Swelling-induced deformations: a materials-defined transition from macroscale to microscale deformations. Soft Matter. 9, 5524–5528. (doi:10.1039/ c3sm00135k)