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 swelling-induced 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.
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 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 , 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 , 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 , and subsequently applied in . 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 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.
2. Swelling-induced deformation processes in polymer gels
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 , the gel body is identified with a region of the three-dimensional Euclidean space , denoted as the reference configuration, with boundary of unit normal m. If a time interval is chosen, the balance equations of forces and solvent are 2.1 and 2.2 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 across its boundary , 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 . 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 of the body; on the other hand, the reference configuration 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 is the steady state attained from through the homogeneous and isotropic map ; hence, Fo=∇fo=λo I, where λo is a uniform stretch and . 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 2.3 with 2.4 moreover, it is assumed that 2.5
Here, G=NκBT is the shear modulus, N is the number of polymeric chains per unit dry volume, T is the temperature, Ω is the solvent molar volume, 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 , it takes the form 2.6 with cdo the solvent concentration per unit dry volume, at the stress-free, swollen state (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  (see also [6,14]), have the following constitutive equations for the reference stress S, the chemical potential μ and the reference solvent flux h: 2.7 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 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 at F=I. We obtain the swollen volume and the stretching from as Jo=1+Ωcdo and , the solvent concentration per unit of reference volume as co=cdo/Jo, the pressure as po=G/λo and the chemical potential as . Under these conditions, , the balance of solvent yields and the balance of forces is identically 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 is assumed to satisfy the balance of forces and the local thermodynamic equilibrium 2.10 the difference between μ(X) and μext determines the driving force of the swelling (shrinking) process. For homogeneous materials, it holds that 2.11 with 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 .
In practice, –10−2 and χ∼0.1–0.5 for good solvents . With and χ varying in such ranges, the equilibrium stretch (that is, ) and equation (2.12) can be approximated by estimating the leading order term in the asymptotic expansion up to O. The following approximate shear modulus versus swelling ratio relation holds: 2.13 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 2.14 Equation (2.14) has already been introduced in  in terms of the volume fraction of the polymer , denoted as the equilibrium volume fraction; in , the approximate swelling ratio was denoted as the hydrogel deformation constant (see also ).
3. Bending deformations from free-swelling solutions
We considered a slightly pre-swollen gel bar , 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 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 . Equations (2.1)–(2.9) were recast in a weak form, using a mixed method to impose both the volumetric constraint on and the boundary conditions on ; thus, the problem was reformulated as follows: find u, p, c and a boundary concentration cb such that, for any , , and , with compatible with the boundary conditions, and on , it holds that 3.1 and 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 , which, in turn, is the steady state attained by 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 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. Moreover, we also found that 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 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 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 beam-like bodies , 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 two-dimensional 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 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 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 4.3 and 4.4 where, because of the material incompressibility, 3Gt and 3Gb identify the corresponding Young moduli.
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 4.5 and 4.6 Equations (4.5) and (4.6) deliver the following linear systems of two equations in Λ0 and : 4.7 and 4.8 with the pairs (Gt,λot) and (Gb,λob) satisfying the local thermodynamic equilibrium for the beam-like layers 4.9 and 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 4.11 4.12 being and 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 three-dimensional stress diffusion model and the structural simplified model looks excellent. Moreover, the longitudinal stretch Λ0 is almost unaffected by the aspect ratio, whereas κ scales almost linearly with it.
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.
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 can be obtained in the simpler form 4.14 with 4.15 4.16 The dashed lines in figure 3 show the patterns of (figure 3a) and (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 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, 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 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  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.
In figure 5, we plot (figure 5a) and (figure 5b) as functions of the moduli 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 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 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 : 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 , the curvature reduces even if the bottom beam is becoming softer and softer as Δo is decreasing. As , 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.
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 , the swelling of a soft thin layer bonded to a fixed rigid substrate 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 , 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 , 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.
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 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 5.2 and 5.3 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 5.4 By plugging equations (5.4) into the balance equations of force and torque, we obtained the following nonlinear system in the pair (Gt,Gb) 5.5 and 5.6 or, equivalently, 5.7 and 5.8 by combining linearly equations (5.5) and (5.6), and setting 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 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
with δ=β2/α2 and ϵ=α1/β1 the two smallness parameters depending on β, h/b, Λo and κ. With this, we employed the expansions in (5.10)1 and in (5.10)2, obtaining the asymptotic system 5.11 and 5.12 which can be solved to give 5.13 and 5.14 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 an iterative approach. We linearized equation (5.7) about the solution G(k)b at the kth iteration to obtain the new approximate solution G(k+1)b, and we employed the asymptotic relation for Gt to update its estimate using the value of G(k+1)b just found, 5.15 and 5.16 As initial values of the iterative process, we choose and 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.
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.
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’.
↵1 We thank the anonymous reviewer who suggested this situation to us.
↵2 For the sake of simplicity, we limit our analysis to bilayered gel beams; the extension to multi-layered gel beams is straightforward.
↵3 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).
↵4 Indeed, for the symmetry of the problem, if the pair (Λo,κ) is known for a beam with and shear moduli , , then the pair (Λ0,−κ) corresponds to a beam with ht=β h and shear moduli , and .
↵5 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.
- Received June 13, 2014.
- Accepted August 12, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.