1 Introduction
Hydrogels are polymer networks swollen in an aqueous solution. Natural entities, e.g. extracellular matrix or the vitreous body of the eye, can be understood within this concept. Similar, various synthetic hydrogels are exploited in commercial products due to their particular properties (e.g. disposable diapers and contact lenses), and have a significant potential for applications like smart valves Beebe et al. (2000), tissue engineering Chan and Mooney (2008); Annabi et al. (2014); Seliktar (2012), drug delivery systems Bysell and Malmsten (2006); Li and Mooney (2016); Culver et al. (2017), and biological sensors Tierney et al. (2008); Buenger et al. (2012). For many of these applications, the capability of hydrogels to swell or shrink as a response to stimuli (e.g. changes in temperature, mechanical forces, pH, salinity level, electric field, specific molecules recognized by included capture moieties) is exploited. During this swelling process, instabilities can occur at the surface of the gels causing wrinkling (i.e. global harmonic waves) and/or creasing (i.e. localized sharp folds) Tanaka et al. (1987); Chan and Crosby (2006); Trujillo et al. (2008); Chan et al. (2008); Guvendiren et al. (2009); Breid and Crosby (2009); Chung et al. (2009); Yang et al. (2010); Li et al. (2012); Chen et al. (2014); Zhou et al. (2017), as illustrated in Figure 1. While creasing mainly occurs as the first mode of instability when homogeneous gels are exposed to large swelling Trujillo et al. (2008); Li et al. (2012), wrinkling patterns are first and foremost found in gels having a gradient stiffness through the thickness, caused by a variation in the crosslinking density within the gel Guvendiren et al. (2009, 2010a, 2010b), or by the deposition of a thin and stiff film at the outer surface of the gel, possibly to alter properties like permeability, stability or biocompatibility Sultan and Boudaoud (2008); Prot et al. (2013); Sherstova et al. (2016). In addition, the impact of the mechanically layered or anisotropic character of natural entities on morphogenesis attains growing interest Liu et al. (2012); Li et al. (2012) with structure formation in aging skin as an example Limbert and Kuhl (2018).
A phenomenon closely related to swelling induced instability in layered hydrogels is the evolution of buckling patterns in layered polymer plates under mechanical loading Bowden et al. (1998); Groenewold (2001); Yang et al. (2010); Yin et al. (2018). Assuming linear elastic material behavior, analytical expressions can be found for the critical strain at the onset of buckling and the wavelength of the resulting wrinkling pattern based on the stiffness ratio between the film and the substrate and the thickness of the film Stafford et al. (2004). Along the same lines, linear perturbation analysis (LPA) frameworks have been developed for studying the onset of instability in layered hydrogel systems exposed to equilibrium swelling Wu et al. (2013, 2017), i.e. assuming a homogeneous field for the chemical potential in the gel at all times. While these analytical or semianalytical approaches are computationally efficient, they neglect the transient nature of the osmotic swelling process in hydrogels, giving a timedependent and inhomogeneous field of the chemical potential through the gel, which could possibly alter buckling initiation. Accounting for the effects of transient swelling in layered gels calls upon numerical procedures like the finite element method (FEM).
For finite element simulations of the equilibrium swelling of gels, numerous constitutive models are available in the literature, both general formulations defining the gel by its basic properties Hong et al. (2009); Kang and Huang (2010a) and formulations that relate the parameters in the simulations more explicitly to the network properties of the specific gels Marcombe et al. (2010); Prot et al. (2013). To also describe the transient nature of hydrogel swelling there are generally two methods that are used in the literature. The most available approach is to utilize the analogy between gel diffusion and heat transfer and implement the swelling process using the thermal modeling capabilities already available in commercial finite element codes Toh et al. (2013); Duan et al. (2013). The other and more general approach is to formulate diffusiondeformation specific elements Zhang et al. (2009); Chester et al. (2015); Bouklas et al. (2015), however, the implementation of such elements can be considered challenging and timeconsuming.
The studies in the literature using finite element simulations for predicting the onset of instability mainly focus on buckling caused by equilibrium swelling (i.e. omitting transient effects) using 2D Kang and Huang (2010b); Wu et al. (2013); Weiss et al. (2013); Budday et al. (2015) or 3D models Xu et al. (2014, 2015); Xu and PotierFerry (2016). Transient effects related to swelling induced buckling are less described. Bouklas et al. Bouklas et al. (2015) briefly discussed simulations of transient swelling, however, they did not demonstrate the difference between equilibrium and transient swelling for the initiation of buckling. Toh et al. Toh et al. (2015) only considered plates with homogeneous material properties, in addition, they used a multipoint constraint routine to force a sinusoidal buckling pattern and hence neglected creasing as a possible mode of instability. In their recent study, Yu et al. Yu et al. (2018) investigated buckling initiation in a 2D model of a thin constrained plate exposed to transient swelling, however, they did not consider gels with a stiffness gradient. Finally, Dortdivanlioglu and Linder Dortdivanlioglu and Linder (2019) studied the initiation of buckling in layered confined hydrogels during transient swelling, though, they did not discuss the dependence of the diffusion coefficient for the initiation of instability nor the effective swelling ratio at the onset of buckling. In addition, none of the mentioned studies discuss the onset of instability during transient swelling of softonhard layered gels.
Although buckling caused by mechanical compression or equilibrium swelling of soft layered materials is well studied using experimental, analytical, or numerical procedures, the effect of the transient nature of hydrogel swelling for the onset of buckling remains insufficiently understood. Hence, we address the effects of transient swelling on the critical swelling ratio, the time to buckling, and the wavelength of the resulting buckling pattern for confined hydrogels. Furthermore, we study both hardonsoft systems, with stiffness ratios between the film and the substrate in the moderate (20, 100) and low range (2, 5), and a softonhard system with a stiffness ratio of 0.5. Due to its importance for the presented simulation results, we give a detailed discussion on the introduction of surface imperfections, the quantification of onset of instability, and the use of a plane strain assumption. For the case of equilibrium swelling, we benchmark our model against an LPA framework available in the literature Wu et al. (2013), obtaining excellent correspondence. The results of this study are believed to be of importance for the development of models that can predict hydrogel swelling behavior in a general manner, as well as providing valuable insight into the nature of swelling induced buckling processes in layered hydrogels, which can affect guiding principles for the design of hydrogel systems.
The article is organized as follows: In the next section, the geometrical properties of the problem studied herein are described. In Section 3, the constitutive formulation used for transient hydrogel swelling is presented. Section 4 briefly describes the LPA approach used to predict the onset of buckling during equilibrium swelling as a benchmark for the finite element approach. Thereafter, in Section 5, the finite element modeling is described in detail, including the implementation of the constitutive model for the commercial finite element software Abaqus and a validation case. Section 6 presents and discusses the results obtained from the finite element procedure and compares them to the linear perturbation results for the case of homogeneous swelling. Finally, conclusions from the study are presented in Section 7.
2 Problem definition
To study the onset of swelling induced instability in constrained layered hydrogel plates we define a square plate as illustrated in Figure 2, having edges of length and a total thickness of . The upper part of the plate consists of a film with a contrast stiffness and a thickness . In the present study, the lengths of the plate are set to mm, the total plate thickness is set to mm, while the thickness of the film is set to . We introduce the film to plate thickness ratio , and note that a value of is used herein. The influence of the chosen dimensions on the obtained results is discussed in Section 6.6.
All lateral edges of the plate are constrained from inplane deformation, while there are no restrictions to swelling in the outofplane direction. At its upper surface, the plate is exposed to a solvent. Hence, a change in the chemical potential in the surrounding of the gel will cause an inhomogeneous and timedependent profile for the chemical potential through the thickness of the gel. This chemical potential can cause swelling and hence give the plate a new and timedependent total height of . The chemical potential profile through the thickness of the gel will depend on the diffusion coefficient of the solvent molecules. As the goal of this study is to investigate the effect of transient swelling on the onset of buckling, we study the described problem with various values for the diffusion coefficient of the solvent molecules. Calculations are performed for diffusion coefficients between m/s and m/s. While m/s represents the lower range of physical values for a gel swelling in water Caccavo et al. (2018), m/s is an unrealistic value, however, it is used to obtain buckling results when the chemical potential is approximately homogeneous through the gel.
3 Constitutive formulation
3.1 Equilibrium swelling
The constitutive modeling of hydrogel behavior applied herein is based on the work by Hong et al. Hong et al. (2008, 2009), Kang and Huang Kang and Huang (2010a), and Toh et al. Toh et al. (2013). The freeenergy function for the hydrogel is assumed to originate from the additive contributions of stretching of the polymer network and mixing of the polymer and the solvent molecules (Flory, 1953; Huggins, 1941; Flory, 1942)
(1) 
where is the number of polymeric chains per reference volume, is the temperature in the unit of energy,
is the deformation gradient tensor,
is the first invariant of the left CauchyGreen tensor , is the volume ratio of the material, is the volume per solvent molecule, is the nominal concentration of solvent molecules, and is the FloryHuggins parameter.By assuming that both the polymer network and the solvent molecules retain their volumes through the swelling process, we find that the volume increase of the gel only can come from an increase in the number of solvent molecules inside the gel, hence we can write
(2) 
Further, to enable implementation for the finite element method, we can introduce a new freeenergy function by the use of a Legendre transformation as
(3) 
and hence ensure the deformation gradient and the chemical potential to be the two independent variables of the model. The Cauchy stress tensor can then be obtained from the freeenergy potential function as
(4) 
where is the secondorder identity tensor.
Normalized quantities are applied in the following and are defined according to , , and .
3.2 Kinetics
In order to capture transient swelling in the constitutive model, the migration of solvent molecules into the hydrogel network must be accounted for. Here, we adopt the modeling procedure by Hong et al. Hong et al. (2008) and assume that the small solvent molecules diffuse in the hydrogel network. A kinetic law for the gel can be written as
(5) 
where
is the nominal flux vector for the solvent molecules,
denotes the position of the material points in the reference configuration, while is the mobility tensor.Further, the diffusion coefficient of the solvent molecules is taken as isotropic and independent of the deformation gradient F and the nominal concentration of solvent molecules Hong et al. (2008). The true flux vector is assumed to follow Fick’s first law Feynman et al. (1963)
(6) 
where is the position of the material points in the current configuration. The nominal and true flux vectors can be related through
(7) 
Combining Equations (2), (5), (6), and (7) we find that the mobility tensor can be written as
(8) 
It can be noted that the outlined material model neglects possible viscoelastic effects in the polymer network. The benefit of this approach is that the transient swelling is the only timedependent feature of the model, and it is thereby easy to isolate the transient effects. The impact of this simplification on the obtained results is discussed in Section 6.4.
4 Linear perturbation analysis
Based on the work by Wu et al. Wu et al. (2013) for a bilayer hydrogel we implemented an LPA framework for studying the onset of mechanical instability during equilibrium swelling (i.e. assuming a homogeneous chemical potential). The procedure is based on a 2D perturbation analysis where the perturbed deformation gradient is assumed to take the form
(9) 
letting be the outofplane deformation, while the assumed perturbations from the equilibrium state in the  and direction are given by and . The further steps of the method are thoroughly described in Wu et al. Wu et al. (2013) and its implementation culminates in solving the determinant of an 88 matrix equal to zero, where all nonzero elements are given between Equations (43) and (44) in Wu et al. Wu et al. (2013).
A stability plot, showing the estimated critical swelling ratio at the onset of instability for different wavelengths of the initial harmonic perturbation, is given in Figure
3. Here, denotes the stiffness ratio between the film and the substrate, , with and being the stiffness of the film and the substrate, respectively. For readability, only results with the stiffness of the substrate set to is shown, however, similar trends would be found for . The analysis is based on the plate height and film thickness as given in Section 2, while the plate width is assumed infinite in the LPA framework.For the softonhard system () we can see that the minimum point for the critical swelling ratio occurs at the zerowavelength limit. This indicates creasing as the dominating mode of instability and a similar response was found for softonhard systems with stiffness ratio values of 0.1 and 0.8 (results not shown here due to readability). For the hardonsoft systems () on the other hand, the minimum point of the critical swelling ratio is found at a defined wavelength, indicating wrinkling to be the dominating mode of instability. This difference between softonhard and hardonsoft systems is in accordance with previously reported results Wu et al. (2013).
From Figure 3 it can also be noted that for hardonsoft systems the critical swelling ratio increases as the stiffness ratio is reduced. Further, the instability curves for the hardonsoft gels display a wider plateau in buckling wavelength around the critical point as the stiffness ratio is increased. Hence, it can be expected that the critical wavelength obtained in experiments would be more sensitive to experimental uncertainties, like the initial surface geometry, as the stiffness ratio is increased. Along the same lines, we can expect that the critical wavelength obtained using numerical methods as FEM could be more sensitive to numerical approximations as the stiffness ratio is increased.
To study the effect of transient swelling for the initiation of buckling, we will make use of the finite element method as described in the following section. The swelling ratio and the wavelength obtained by the LPA will be used as a benchmark for the finite element simulations with a large diffusion coefficient (i.e. m/s), causing a nearly homogeneous chemical potential through the thickness of the gel, to mimic the assumptions of the LPA calculations.
5 FE modeling
5.1 Implementation and validation of the constitutive model
5.1.1 Implementation
To accommodate for finite element simulations of the defined swelling problem, the constitutive model described in Section 3.1 was implemented as a Fortran routine to be used with the commercial finite element software Abaqus/Standard. The equilibrium behavior of the model is implemented as a userdefined material, UMAT Abaqus (2014), defining the equations for the normalized stress and the Jacobian tangent stiffness. As the chemical potential is singular in the dry state of the hydrogel , an intermediate configuration is introduced. Consequentially, the full deformation tensor can be split in a multiplicative manner as . maps the reference configuration to the intermediate configuration, and maps the intermediate configuration to the current configuration, as illustrated in Figure 4. In the intermediate configuration, the gel is assumed stressfree and in a state of isotropic strain such that . In the present study, all finite element simulations start in the intermediate configuration with a homogeneous distribution of the normalized chemical potential in the gel and dimensions as given in Section 2. The value of is found by numerically solving at the start of the finite element simulation for the initial chemical potential and the material parameters of the hydrogel.
To account for the transient behavior of the gel swelling process, a coupled solvent diffusion and large deformation procedure is utilized. This procedure follows the main ideas presented by Toh et al. Toh et al. (2013) and a detailed description is included in A.
The implemented Fortran code and an example input file for Abaqus is made available as a Mendeley dataset linked to this work Ilseng and Prot .
5.1.2 Validation
To validate the implemented model, we study the case of 1D swelling of a homogeneous perfect plate. We start from an intermediate state where the gel is in a state of isotropic swelling so that it equilibrates at a homogeneous chemical potential of . The plate is then confined from inplane expansion but can still freely swell outofplane. The chemical potential at the surface of the plate is changed to and the gel gradually swells in the outofplane direction giving a deformation gradient with the nonzero elements and where . The swelling ratio of the gel can then be expressed as . As the swelling in the direction is free, the normalized Cauchy stress component must be zero throughout the gel. From Equation (4) we then get
(10) 
Equation 10 yields an expression for the normalized chemical potential that varies through the thickness of the gel as
(11) 
From the condition of molecular incompressibility, we can express the nominal concentration of solvent molecules as
(12) 
By combining Equations (5) and (8) the nominal flux in the direction can be written as
(13) 
where we can express
through use of the chain rule
(14) 
The conservation of solvent molecules requires that
(15) 
where can be obtained by the use of Equation (12)
(16) 
Finally, by inserting Equation (14) into Equation (13), we can with the use of Equation (16) write Equation (15) as
(17) 
5.2 Boundary conditions
To induce swelling in the hydrogel plate described in Section 2, the chemical potential is changed at the upper surface from to . To avoid an abrupt change in boundary condition in the finite element simulation, the chemical potential is changed during a smooth step with a duration of 0.2 seconds (assumed to be short compared to the time to buckling initiation for transient swelling). The potential at the exposed surface is then kept constant at zero and the hydrogel swelling is driven by the diffusion process. Through this process, the plate is constrained from inplane deformation at all lateral faces, while the lower surface is constrained from outofplane deformation. To comply with the LPA results and also accommodate for large swelling ratios, we want to mimic a nearly dry state at the start of the simulation, and hence set the initial homogeneous chemical potential to , giving an initial free swelling of .
5.3 Material parameters
The material parameters used in this study are given in Table 1. defines the normalized stiffness of the substrate, while the parameter defines the ratio between the stiffness in the film and the stiffness of the substrate according to . For the diffusion coefficient , a value in the range between 10 and 10 is reported as representative for the diffusion of water molecules in a gel (Hong et al., 2008; Caccavo et al., 2018). However, we perform simulations for every decade of in the range from to to study the change in the buckling behavior of hydrogels as we gradually go from a transient to an equilibrium swelling process. The set of parameters as given in Table 1 leads to a total of 120 simulations to cover the parameter space.
In addition, simulations were performed for and . These simulations confirm the general trends for softonhard systems as presented in the following. Hence, the results for these two stiffness ratios are not included in the current presentation to enhance readability.
0.01, 0.001  0.5, 2, 5, 20, 100  0.5   1 
5.4 Buckling initiation
To obtain a measure for the global swelling ratio at the point of buckling and beyond, a plane is optimized with respect to the position of the uppermost layer of nodes in the film as illustrated in Figure 6. The height of the plane, , is found by minimizing the residual given as the sum of the squared distances between the nodes in their current position and the plane, i.e. . The minimization procedure is performed using the SciPy package of Python. The global swelling ratio is defined as the plane height divided by the initial plate thickness mm.
There are multiple approaches that could be used to quantify the point where buckling initiates in a finite element simulation, e.g. abrupt changes in stress, strain energy or geometry. For practical applications, we consider the topology of the surface to be the most relevant measure. Hence, we define the point of onset of buckling by the use of a threshold value for the difference between the highest and lowest point on the surface, denoted , i.e. buckling occurs if . Obviously, the obtained results for the onset of buckling will be influenced by the choice of the buckling limit, and the value to be used for should depend on the application of the simulated gel. For the purpose of this study, we choose to use =. Although the absolute values for the swelling ratio and time to initiation of buckling would be altered with a different threshold value, all trends presented in Section 6 would be preserved if the buckling limit would be set to a lower value like .
5.5 Initial imperfection
To accommodate buckling, an initial imperfection is introduced in the upper surface of the stiff film. A Python script is used to change the position of all the surface nodes in the Abaqus input file. To ensure that the method of introducing this imperfection is not dominating the obtained results, we use two different modes of the initial imperfection. First, we introduce a single sinusoidal wave over the width of the gel
(18) 
where is the new position of the node , is the node coordinate in the perfect model, is the maximum perturbation given as a fraction of the initial height, and is the position of node . Second, we use a random perturbation given by
(19) 
where is a pseudorandom number. While the expression in Equation (18) gives a smooth surface, it assumes a specific shape of the surface of the gel. The random imperfection in Equation (19), on the other hand, avoids a priori assumptions of the shape of the initial surface but can produce highly uneven surfaces.
The effect of the initial imperfection size, , on the maximum height difference in the plate surface, , for slow and fast diffusion processes, and m/s, respectively, is shown in Figure 7. Clearly, a larger imperfection size would lead to faster growth of the height difference in the plane, however, the same asymptote is reached for all simulations where buckling is initiated. It can be seen that for the random initial imperfections, there is a significant increase in at the beginning of the simulations with a slow diffusion process. This effect arises as the random imperfection introduces peak and valley nodes in the mesh, where the peak nodes would be exposed to more solvent and hence swell faster than the valley points (similar to free swelling of a cube where the corners would swell faster than the center of the faces Zhang et al. (2009)). This initial effect is not seen for m/s, as this would resemble an equilibrium process with equal swelling ratios in both peak and valley points, nor for the sinusoidal imperfection as this surface will be initially smooth. To demonstrate the effect of the randomness introduced in the mesh by Equation (19), results from three subsequent random imperfections applied to an initially perfect mesh using are shown (i.e. overlapping blue curves in Figure 7). Finally, the plot also shows that the imperfection must be larger than a critical size to trigger buckling, as seen for the curve for a random imperfection with where no buckling is initiated.
In all further calculations, we set to be 0.01%, meaning that the plate height at each node will be in the range between 0.49995 mm and 0.50005 mm and that we initially have . It is important to note that the initial imperfection must be smaller than the threshold value used to quantify the onset of buckling. In the following, we perform simulations using either the harmonic imperfection (Equation (18)) or the random imperfection (Equation (19)).
5.6 Plane strain assumption
To capture both the gradient of the chemical potential through the thickness of the plate in the transient swelling analyses and the evolving inplane buckling pattern requires a relatively fine mesh, putting significant demands on computational resources.
However, it can be noted that for the problem at hand, a plane strain assumption would be correct up to the point of buckling. After buckling has initiated though, a plane strain assumption would restrain the model to a 1D buckling pattern, in contradiction with experimental and theoretical studies Cai et al. (2011) having shown that other buckling patterns, like a hexagonal or herringbone mode, would be energetically favorable.
To investigate if a 2D plane strain assumption would predict buckling at the same swelling ratio as a full 3D model, despite the difference in bucking mode, a simple case of equilibrium swelling with and was used. Due to the homogeneous chemical potential and the large wavelength of the resulting buckling pattern, a relatively coarse mesh could be used. The 2D and 3D problems were both discretized in the same manner through the thickness using fully integrated higher order elements. The obtained buckling patterns and the increase in during swelling are shown in Figures 8 and 9 respectively. The modes shown in Figure 8 are obtained from the first increment after initiation of buckling according to the threshold level . The initial buckling pattern in the 3D simulation can be seen to have a checkerboard pattern in the central region of the plate, this gradually expands to cover the plate before it would transition to a herringbone mode at larger swelling ratios Cai et al. (2011). For the plane strain simulations, the 2D nature of the model forces the plate to buckle in the less advantageous 1D buckling pattern more commonly observed when one of the inplane stresses dominates Breid and Crosby (2011). From the fringe plots, it can be seen that the value of is larger in the 3D simulation. However, from Figure 9 it is seen that the swelling ratio (i.e. the height of the optimized plane) at the onset of buckling is similar between the 2D and 3D simulations. Hence, to reduce the computational demands, a plane strain model is used in the further to study the point of onset of instability. This limits the study from considering detailed postbuckling analysis as the correct buckling pattern cannot be obtained.
5.7 Element type and size
For the further 2D simulations, the plate is discretized by fully integrated linear plane strain temperaturedisplacement elements (denoted CPE4T in Abaqus) where the temperature field is used to mimic the distribution of the normalized chemical potential through the plate (see Appendix A for the analogy between gel diffusion and heat transfer). Linear elements were preferred over secondorder elements as the first use a lumped heat capacity matrix and hence avoid spurious oscillations in the temperature field during changes in the boundary conditions Bergheau and Fortunier (2013). Thereby, with our use of the heat transfer and gel swelling analogy, oscillations in the chemical potential during transient swelling is avoided. See Bouklas et al. (2015); Dortdivanlioglu et al. (2018) and the references therein for a further discussion on oscillations in the chemical potential during finite element simulations of transient gel swelling.
The effect of the element size on the obtained results was studied by running simulations with , and m/s using between 1 and 30 square elements over the film thickness. To avoid a change in the surface topology as the mesh is refined, the harmonic imperfection (Equation (18)) was used. As the gradient in the chemical potential is largest close to the surface of the gel, the total size of the model was reduced by gradually increasing the mesh size going from the film and into the substrate as illustrated in the enlarged view in Figure 10. The resulting swelling ratio at the onset of instability as a function of the number of elements over the thickness of the stiff film is shown in Figure 11. The red squares refer to the swelling ratio, while the blue circles refer to the relative change in the swelling ratio compared with the previous (larger) element size. The same trend is observed with respect to the time to onset of instability. As the swelling ratio at the onset of buckling seems to be converged for a mesh with 30 square elements over the film thickness, this mesh size, leading to a total of 194 400 2D elements, is used in all further calculations. It is worth noting that for this element size, the random imperfection with set to 0.01% produces a maximum perturbation of the upper nodes of 3% of the initial element height. A sketch of the final mesh with a random imperfection can be seen in Figure 10.
6 Results and discussion
6.1 LPA vs FEM for equilibrium swelling
A comparison between the FEM and the LPA results for hardonsoft systems during equilibrium swelling are shown in Figure 12 where the left ordinate indicates the critical swelling ratio at buckling initiation and the right ordinate indicates the wavelength of the buckling pattern at initiation. The FEM results are shown for and m/s. For the critical swelling ratio, a nearly perfect correspondence can be observed between the two analysis methods. The predicted wavelength, on the other hand, shows a perfect correspondence between the two methods for and , for the two other stiffness ratios a reasonably good correspondence is seen. A reason for the discrepancy in the wavelength results can stem from the finite width of the FEM model, while the LPA method assumes an infinitely wide plate. However, the generally good correspondence between the FEM and LPA results indicates that the length of the plate in the FEM model is sufficient to yield results representative for infinitely wide plates. It can be noted that a comparison between FEM and LPA results using shows a similar correspondence, and were omitted from Figure 12 merely to enhance readability.
For the softonhard system () both the LPA and the FEM approaches predict creasing as the initial instability mode. For the swelling ratio at the onset of buckling, on the other hand, the FEM and LPA results deviate with about 8%. However, it is observed that for and the FEM results have not fully reached the equilibrium plateau level at the diffusion coefficient of 1 m/s (seen in Section 6.3, Figure 14), explaining a part of the discrepancy.
6.2 Stress and chemical potential profile at the initiation of buckling
The distribution of the normalized chemical potential, , and the normalized inplane stress, , through the thickness of the plate at the onset of buckling are shown in Figure 13 for all values of the diffusion coefficient using stiffness parameters of and (note that the presented trends are representative also for other parameter combinations). Here, is used to denote the position in the plate in the intermediate state (i.e. the starting point for the simulations) such that the curves are easily comparable. The data are extracted from the right edge of the plane strain plate, however, the profiles would be representative for the whole plate at the initiation of buckling.
For the normalized chemical potential shown in Figure 12(a), three distinct regions of the diffusion coefficient can be found. First, for all simulations with m/s, we see that the chemical potential at the onset of buckling is nearly homogeneous through the thickness of the plate and we have overlapping curves, i.e. equilibrium swelling is represented. We denote this region I. Second, we have a transition region in the range m/s m/s where the chemical potential through the plate is increasingly inhomogeneous as the diffusivity is reduced. We denote this region II. Finally, the results for a diffusion coefficient m/s produces overlapping curves as the profile of the chemical potential through the thickness at the onset of buckling is nearly independent of the diffusion coefficient. We denote this region III.
For the profiles of the normalized inplane stresses at the onset of buckling seen in Figure 12(b), we see the same three regions as for the normalized chemical profile, with overlapping curves for fast (I) and slow (III) diffusion and a transition zone inbetween (II). For the equilibrium swelling simulations (i.e. m/s), the chemical potential is homogeneous through the thickness of the plate, and consequentially the compressive stress changes through the plate only due to the stiffness difference between the film and the substrate.
For the slow diffusion processes (i.e. m/s), on the other hand, the compressive stress is gradually reduced (i.e. becomes less negative) going downwards through the film thickness. Then there is a discontinuous reduction as one goes from the stiff film to the soft substrate. In the substrate, there is a further gradual decrease towards zero as one approaches the bottom of the plate. It can also be seen that the maximum compressive stress, located at the top of the film, increase as the diffusion coefficient is reduced.
6.3 Swelling ratio at buckling initiation
Figure 14 shows the swelling ratio of the plate (i.e. the global amount of swelling ) at the onset of buckling for values of the diffusion coefficient in the range from to , between 100 and 0.5, and a normalized stiffness of the substrate of 0.001 or 0.01. The plotted data were obtained using a random initial imperfection in the upper surface of the plate (i.e. Equation 19), however, a similar plot would be produced if a harmonic initial imperfection was used (i.e. Equation 18). In accordance with the results for the normalized chemical potential and compressive stress profiles at the onset of buckling (discussed in Section 6.2) the swelling ratio at the onset of buckling can be divided into the three distinct regions I, II, and III. In region I, there is a relatively large and stable swelling ratio at buckling initiation that equals the equilibrium solution. Then, in region II there is a transition zone where the swelling ratio at buckling initiation is gradually reduced as the diffusivity is lowered. Finally, in region III a relatively low and stable swelling ratio at buckling initiation is obtained. In the plot, it can also be seen that the range of the diffusion coefficient in which the three regions occur depends on the value of , with the transitions between the regions moving towards higher values of as is reduced.
Further, the swelling ratio at the onset of instability can also be seen to depend on the absolute stiffness of the substrate and the film, and not only the ratio between the two. For the softonhard system (), a larger critical swelling ratio is obtained at high values of for the softer gels (i.e. ) compared to the stiffer gels (i.e. ). For hardonsoft gels () on the other hand, a larger swelling ratio at the onset of buckling is obtained for the stiffer gels (i.e. ) compared to the softer gels (i.e. ). This is especially seen in the results for combined with low values of , where a swelling ratio close to unity is present at the onset of buckling if , considerably below a swelling ratio of 1.2 as found for . A further discussion on the large difference between these parameter combinations is given in Section 6.5 discussing the obtained buckling profiles.
6.4 Time to buckling initiation
Figure 15 shows the time to the onset of instability on a loglog plot for the simulations where wrinkling was obtained as the mode of instability (see Section 6.5 for a discussion on buckling profiles). Simulations predicting creasing were excluded from the plot as the time to onset of creasing is known to show a strong mesh dependence Bouklas et al. (2015).
In Figure 15, it is seen that there is a negligible effect on the time to the onset of buckling in region I and II (i.e. m/s), however, reducing the diffusion coefficient even further (i.e. into the physically reasonable region for water molecules diffusing in a gel) gives a clear rise in the time needed for buckling to occur, with an increasing effect as is reduced for hardonsoft gels. For all parameter combinations, a linear loglog relation between the time to instability and the diffusion coefficient can be seen in region III, corresponding to the stable region for swelling at the onset of instability for low values of as shown in Figure 14. It is also seen that there is a clear dependence on the absolute value of the stiffness in the plate and the substrate and not only the ratio between the two, with longer times to buckling initiation for stiffer gels.
From the time to instability results it can be noted that the time used on the smooth increase of the chemical potential at the surface of the gel in the FEM model, i.e. 0.2 seconds, is negligible compared to the time to the onset of buckling for low values of the diffusivity. In addition, with the timescale needed to obtain buckling in region III, it is assumed that the contribution from viscoelastic effects in the polymer network would be negligible. For region I and II, on the other hand, further studies are needed to investigate how viscous effects might alter the initiation of buckling.
6.5 Buckling profiles
To study the dependence on the diffusion coefficient for the obtained buckling pattern we extract the wavelength found in the FEM simulations for all values of . Figure 16 shows the obtained wavelengths for , and . For and , no clear trend in the critical wavelength could be seen.
For the softonhard gel (), a critical wavelength of zero was obtained for all values of the diffusion coefficient as creasing was found to be the first mode of instability. This result was also found in the simulations with and . For the hardonsoft gels on the other hand ( and ), wrinkling was found as the first mode of instability in most of the simulations, with a tendency of a reduction in the critical wavelength as the diffusion coefficient was reduced. For and , the critical wavelength drops to zero when the diffusivity goes below m/s, as creasing (rather than wrinkling) was found as the critical mode of instability in these simulations. This is illustrated in Figure 17 showing the instability mode for the fast and the slow diffusion processes for and . It can be seen that for fast swelling, a structured wrinkling pattern is obtained where the size of the initial imperfection is negligible compared to the length scale of the wrinkles. For the slow swelling simulation, on the other hand, it can be noted that the outofplane deformation is limited to the upper part of the film and that a localized folding process has started at the surface of the gel. The location of the developing crease depends on the initial surface imperfection in the model. The results shown in Figures 16 and 17 were obtained with a random imperfection in the upper surface of the plate (i.e. Equation (19)), however, similar results were obtained using the harmonic imperfection (i.e. Equation (18)). In addition, the shift from wrinkling to creasing as the diffusion coefficient decrease from m/s to m/s was found to be insensitive of the mesh size (tested in the range between 20 and 50 elements over the film thickness).
The wrinkling to creasing transition for a hardonsoft gel clearly demonstrates how the timedependent nature of swelling can cause a discrepancy between the buckling mode obtained in experiments compared to that predicted by an equilibrium analysis. Further, the fact that and reach instability by creasing before wrinkling at low values of explains the large difference in the critical swelling ratio between and during slow swelling for seen in Figure 14.
The authors hypothesize that creasing can occur during transient swelling of hardonsoft gels with low film to substrate stiffness ratios due to a small difference between the strain energy in the creased and wrinkled configurations. For a sufficiently slow swelling process (i.e. a high ratio between the film thickness and the diffusion coefficient) in a sufficiently soft gel, instability will be triggered while the compressive stresses are focused in the upper part of the film, causing creasing to occur rather than wrinkling. However, further experimental and theoretical research is needed to fully understand the wrinkling to creasing transition during transient swelling for a hardonsoft system and to quantify the critical conditions for this transition to occur.
6.6 Dependence on initial dimensions
6.6.1 Preliminaries
For the length of the plate , the generally good agreement between the LPA results, assuming an infinitely wide plate, and the FEM results, giving the plate a finite length, indicates that the ratio H/L is small enough for the results herein to be representative for infinitely wide plates. To discuss the influence on the obtained results by the chosen values for the film vs plate ratio and total plate thickness , we will consider equilibrium and transient swelling separately.
6.6.2 Equilibrium swelling
For a discussion on the geometrical dependencies under equilibrium swelling, the LPA framework is an effective method. First, looking at the film vs plate thickness ratio , the dependence on the critical swelling ratio and the normalized critical wavelength are shown in Figures 17(a) and 17(b), respectively. The data in the plots were obtained using , although similar trends could be obtained with . Starting with the softonhard case (), both the critical swelling ratio and the normalized critical wavelength can be seen to be nearly independent of . For the hardonsoft cases () on the other hand, we see that near the film to plate height ratio of used herein the critical swelling ratio (Figure 17(a)) is relatively stable with respect to the film fraction. However, there is a clear tendency of an increasing critical swelling ratio as gets larger than 0.25, with a stronger dependence observed for lower stiffness ratios. The predicted normalized critical wavelength (Figure 17(b)) for the hardonsoft cases is changing significantly near , with increasing wavelengths as the film to plate thickness ratio is increased.
Second, for a change in the initial height of the plate (keeping the ratio constant), both the critical swelling ratio and the normalized critical wavelength (as presented in Figure 18) were found to be independent of under the assumption of equilibrium swelling of an infinitely wide plate.
6.6.3 Transient swelling
For a discussion on how the results for transient swelling depend on the initial thickness of the plate (keeping constant), it is worth noting that the defined problem introduces a relation between the length scale of the swelling plate and the diffusion coefficient of the gel material. This relation means that scaling with a factor would be equivalent to scaling the dimensions of the plate with a factor 1/. Hence, the results presented in the previous sections for varying diffusivity and constant dimensions could have been obtained with a constant diffusivity and varying the plate dimensions. This scaling is demonstrated in Figure 19, plotting the critical swelling ratio as a function of the initial plate thickness for a constant diffusivity of m/s.
It is important to note the difference in the dependence on the initial plate height between equilibrium and transient swelling. This difference highlights the importance of the length scale of a swelling hydrogel system for whether simple equilibrium analyses would yield precise predictions for the initiation of buckling, with improved precision of equilibrium estimates for thinner gels. In addition, the length scale dependence in transient swelling means that the regions I, II, and III and their respective diffusivity values as given in Section 6.2 are specific for the plate dimensions used.
A further discussion on how a change in the ratios and would influence the initiation of buckling during transient swelling is considered outofscope for the present work.
7 Conclusion
In this paper, we show that the onset of instability in hydrogels with gradient stiffness is highly influenced by the kinetic nature of the swelling process. Both the critical swelling ratio and the time to buckling initiation were found to increase as the diffusion coefficient of the material was reduced. In addition, these measures were found to depend on the absolute stiffness of the film and the substrate and not only the ratio between the two.
For the geometrical configuration studied herein, changing the diffusivity within the physical region for hydrogels in water would have a negligible effect on the swelling ratio at the onset of buckling, while the time to the onset of buckling could change significantly, with an increasing effect as the stiffness ratio between the stiff film and the soft substrate is reduced.
For softonhard gels, creasing was found as the first mode of instability independent of the diffusion coefficient. For hardonsoft gels on the other hand, the diffusivity of the gel material could alter the buckling pattern by reducing the wavelength as the diffusivity was reduced. This is most evident for the combination of material parameters and , where a stable wrinkling pattern was predicted by an equilibrium analysis, while creasing was triggered as the first mode of instability when using a physical value for the diffusion coefficient of the gel.
The results presented herein for a change in the diffusivity with constant plate dimensions could also be read as a change in plate dimensions (retaining the shape) with a constant value for the diffusivity. Whether simplified analyses based on an assumption of equilibrium swelling would be predictive for the onset of buckling in a layered hydrogel system would hence depend on the length scale of the problem.
For further work, we suggest that the transition from wrinkling to creasing in hardonsoft gels is studied both theoretically and experimentally. Further, viscoelastic effects of the polymer network should be included in the constitutive modeling to study the influence this will have on the initiation of buckling. In addition, analyses of the postbuckling behavior of confined swelling gels where the threedimensional nature of the buckling pattern is accounted for could increase the understanding of the fundamental mechanisms governing swelling induced buckling and pattern evolution.
Conflicts of interest
There are no conflicts to declare.
Acknowledgments
This work was supported by the Norwegian Research Council (Project no 240299/F20).
Appendix A Gel diffusion and heat transfer analogy
This appendix outlines the analogy between gel swelling and heat transfer, paving the way for implementing coupled solvent diffusion and large deformations in Abaqus using default temperaturedisplacement elements.
The energy balance in heat transfer as it is implemented in Abaqus reads
(20) 
where is the material density, is the specific heat capacity, is temperature, is the heat flux flowing into the body, the outward unit normal, and is the heat supplied internally into the body. and define the volume and the surface of the body in , respectively. The heat flux is given by
(21) 
where is the thermal conductivity of the material. By assuming the conservation of solvent molecules in a gel and its solution reads Chester and Anand (2010); Toh et al. (2013)
(22) 
while the flux of the solvent molecules in a gel was given in Equation (6). By comparing Equation (20) to Equation (22) and Equation (21) to Equation (6), we find that the equivalents of temperature, specific heat capacity, mass density, and conductivity can be identified as , , , and respectively. These values are calculated as internal variables in the Fortran code using the USDFLD subroutine in Abaqus. The internal heat source in Abaqus, , is set to zero.
To find an analytical expression for to implement as the specific heat capacity is nontrivial. A strategy suggested by Toh et al. Toh et al. (2013) is to derive an expression for from Equation (4), using Equation (2) to relate and
(23) 
where denotes the trace of the normalized Cauchy stress tensor and . Assuming , the righthand side of Equation (23) is differentiated with respect to by use of symbolic differentiation. Using , an expression for is obtained as
(24) 
Note that Equation (24) is slightly different from the expression given in Toh et al. Toh et al. (2013). The good benchmark results shown in Figure 5 indicates that Equation (24) yields good accuracy for the problem studied herein.
References
References
 6.144. Dassault Systèmes. Cited by: §5.1.1.
 25th Anniversary Article: Rational Design and Applications of Hydrogels in Regenerative Medicine. Advanced materials (Deerfield Beach, Fla.) 26, pp. 85–124. External Links: Document Cited by: §1.
 Functional hydrogel structures for autonomous flow control inside microfluidic channels. Nature 404, pp. 588–590. External Links: Document Cited by: §1.
 Finite Element Simulation of Heat Transfer. Wiley. Cited by: §5.7.
 A nonlinear, transient finite element method for coupled solvent diffusion and large deformation of hydrogels. Journal of the Mechanics and Physics of Solids 79, pp. 21–43. External Links: Document Cited by: §1, §1, §5.7, §6.4.
 Spontaneous formation of ordered structures in thin films of metals supported on an elastomeric polymer. Nature 393, pp. 146–149. External Links: Document Cited by: §1.
 Surface wrinkling behavior of finite circular plates. Soft Matter 5, pp. 425–431. External Links: Document Cited by: §1.
 Effect of stress state on wrinkle morphology. Soft Matter 7, pp. 4490–4496. External Links: Document Cited by: §5.6.
 Perioddoubling and periodtripling in growing bilayered systems. Philosophical Magazine 95, pp. 3208–3224. External Links: Document Cited by: §1.
 Hydrogels in sensing applications. Progress in Polymer Science 37, pp. 1678–1719. External Links: Document Cited by: §1.
 Visualizing the interaction between polyLlysine and poly(acrylic acid) microgels using microscopy techniques: Effect of electrostatics and peptide size. Langmuir 22, pp. 5476–5484. External Links: Document Cited by: §1.
 Hydrogels: experimental characterization and mathematical modelling of their mechanical and diffusive behaviour. Chem. Soc. Rev. 47, pp. 2357–2373. External Links: Document Cited by: §2, §5.3.
 Periodic patterns and energy states of buckled films on compliant substrates. Journal of the Mechanics and Physics of Solids 59, pp. 1094–1114. External Links: Document Cited by: §5.6, §5.6.
 Spontaneous formation of stable aligned wrinkling patterns. Soft Matter 2, pp. 324–328. External Links: Document Cited by: §1.
 Surface wrinkles for smart adhesion. Advanced Materials 20, pp. 711–716. External Links: Document Cited by: §1.
 New materials for tissue engineering: towards greater control over the biological response. Trends in Biotechnology 26, pp. 382–392. External Links: Document Cited by: §1.
 Stimuliresponsive buckling mechanics of polymer films. Journal of Polymer Science, Part B: Polymer Physics 52, pp. 1441–1461. External Links: Document Cited by: §1.
 A coupled theory of fluid permeation and large deformations for elastomeric materials. Journal of the Mechanics and Physics of Solids 58, pp. 1879–1906. External Links: Document Cited by: Appendix A.
 A finite element implementation of a coupled diffusiondeformation theory for elastomeric gels. International Journal of Solids and Structures 52, pp. 1–18. External Links: Document Cited by: §1.
 Diffusioncontrolled, selforganized growth of symmetric wrinkling patterns. Advanced Materials 21, pp. 1358–1362. External Links: Document Cited by: §1.
 AnalyteResponsive Hydrogels: Intelligent Materials for Biosensing and Drug Delivery. Accounts of Chemical Research 50, pp. 170–178. External Links: Document Cited by: §1.
 Mixed isogeometric analysis of strongly coupled diffusion in porous materials. International Journal for Numerical Methods in Engineering 114, pp. 28–46. External Links: Document Cited by: §5.7.
 Diffusiondriven swellinginduced instabilities of hydrogels. Journal of the Mechanics and Physics of Solids 125, pp. 38–52. External Links: Document Cited by: §1.
 Simulation of the Transient Behavior of Gels Based on an Analogy Between Diffusion and Heat Transfer. Journal of Applied Mechanics 80, pp. 41017. External Links: Document Cited by: §1.
 The Feynman lectures on physics. AddisonWesley Pub. Co.. Cited by: §3.2.
 Thermodynamics of High Polymer Solutions. The Journal of Chemical Physics 10, pp. 51–61. External Links: Document Cited by: §3.1.
 Principles of polymer chemistry. Cornell University Press. Cited by: §3.1.
 Wrinkling of plates coupled with soft elastic media. Physica A: Statistical Mechanics and its Applications 298, pp. 32–45. External Links: Document Cited by: §1.
 Kinetic study of swellinginduced surface pattern formation and ordering in hydrogel films with depthwise crosslinking gradient. Soft Matter 6, pp. 2044–2049. External Links: Document Cited by: §1.
 Solvent induced transition from wrinkles to creases in thin film gels with depthwise crosslinking gradients. Soft Matter 6, pp. 5795–5801. External Links: Document Cited by: §1.
 SwellingInduced surface patterns in hydrogels with gradient crosslinking density. Advanced Functional Materials 19, pp. 3038–3045. External Links: Document Cited by: §1.
 Inhomogeneous swelling of a gel in equilibrium with a solvent and mechanical load. International Journal of Solids and Structures 46, pp. 3282–3289. External Links: Document Cited by: §1, §3.1.
 A theory of coupled diffusion and large deformation in polymeric gels. Journal of the Mechanics and Physics of Solids 56, pp. 1779–1793. External Links: Document Cited by: §3.1, §3.2, §3.2, §5.3.
 Solutions of Long Chain Compounds. The Journal of Chemical Physics 9, pp. 440–440. External Links: Document Cited by: §3.1.
 [35] User subroutine for modeling transient swelling of hydrogels with Abaqus. Mendeley Data. External Links: Document Cited by: §5.1.1.
 A Variational Approach and Finite Element Implementation for Swelling of Polymeric Hydrogels Under Geometric Constraints. Journal of Applied Mechanics 77, pp. 61004. External Links: Document Cited by: §1, §3.1.
 Swellinduced surface instability of confined hydrogel layers on substrates. Journal of the Mechanics and Physics of Solids 58, pp. 1582–1598. External Links: Document Cited by: §1.
 Mechanics of morphological instabilities and surface wrinkling in soft materials: a review. Soft Matter 8, pp. 5728–5745. External Links: Document Cited by: §1.
 Designing hydrogels for controlled drug delivery. Nature Reviews Materials 1, pp. 16071. Cited by: §1.
 On skin microrelief and the emergence of expression microwrinkles. Soft Matter 14, pp. 1292–1300. External Links: Document Cited by: §1.
 Pattern formation in plants via instability theory of hydrogels. Soft Matter 9, pp. 577–587. External Links: Document Cited by: §1.
 A theory of constrained swelling of a pHsensitive hydrogel. Soft Matter 6, pp. 784–793. External Links: Document Cited by: §1.
 Swelling of a hemiellipsoidal ionic hydrogel for determination of material properties of deposited thin polymer films: an inverse finite element approach. Soft Matter 9, pp. 5815–5827. External Links: Document Cited by: §1, §1.
 Designing cellcompatible hydrogels for biomedical applications. Science 336, pp. 1124–1128. External Links: Document Cited by: §1.
 Nanoindentation and finite element modelling of chitosan–alginate multilayer coated hydrogels. Soft Matter 12, pp. 7338–7349. External Links: Document Cited by: §1.
 A bucklingbased metrology for measuring the elastic moduli of polymeric thin films. Nature Materials 3, pp. 545–550. External Links: Document Cited by: §1.
 The Buckling of a Swollen Thin Gel Layer Bound to a Compliant Substrate. Journal of Applied Mechanics 75, pp. 51002–51005. Cited by: §1.

Mechanical instability of gels at the phase transition
. Nature 325, pp. 796–798. External Links: Document Cited by: §1.  Determination of swelling of responsive gels with nanometer resolution. Fiberoptic based platform for hydrogels as signal transducers. Analytical Chemistry 80, pp. 5086–5093. External Links: Document Cited by: §1.
 Wrinkling of a Polymeric Gel During Transient Swelling. Journal of Applied Mechanics 82, pp. 061004. External Links: Document Cited by: §1.
 Inhomogeneous Large Deformation Kinetics of Polymeric Gels. International Journal of Applied Mechanics 05, pp. 1350001. External Links: Document Cited by: Appendix A, Appendix A, §1, §3.1, §5.1.1.
 Creasing instability of surfaceattached hydrogels. Soft Matter 4, pp. 564. External Links: Document Cited by: §1.
 Creases and wrinkles on the surface of a swollen gel. Journal of Applied Physics 114, pp. 073507. External Links: Document Cited by: §1.
 Swellinduced surface instability of hydrogel layers with material properties varying in thickness direction. International Journal of Solids and Structures 50, pp. 578–587. External Links: Document Cited by: §1, §1, §1, §4, §4.
 Onset of swellinduced surface instability of hydrogel layers with depthwise graded material properties. Mechanics of Materials 105, pp. 138–147. External Links: Document Cited by: §1.
 Instabilities in thin films on hyperelastic substrates by 3D finite elements. International Journal of Solids and Structures 6970, pp. 71–85. External Links: Document Cited by: §1.
 3D finite element modeling for instabilities in thin films on soft substrates. International Journal of Solids and Structures 51, pp. 3619–3632. External Links: Document Cited by: §1.
 A multiscale modeling framework for instabilities of film/substrate systems. Journal of the Mechanics and Physics of Solids 86, pp. 150–172. External Links: Document Cited by: §1.
 Harnessing surface wrinkle patterns in soft matter. Advanced Functional Materials 20, pp. 2550–2564. External Links: Document Cited by: §1, §1.
 Surface wrinkling of anisotropic films bonded on a compliant substrate. International Journal of Solids and Structures 141142, pp. 219–231. External Links: Document Cited by: §1.
 A threedimensional transient mixed hybrid finite element model for superabsorbent polymers with straindependent permeability. Soft Matter 14, pp. 3834–3848. External Links: Document Cited by: §1.
 A finite element method for transient analysis of concurrent large deformation and mass transport in gels. Journal of Applied Physics 105, pp. 093522. External Links: Document Cited by: §1, §5.5.
 Transition of surface–interface creasing in bilayer hydrogels. Soft Matter 13, pp. 6011–6020. External Links: Document Cited by: §1.
Comments
There are no comments yet.