Swelling kinetics of constrained hydrogel spheres

Abstract

A cross-linked polymer network immersed in a solvent will absorb molecules from its surroundings, leading to transient swelling. Under the constraint of a semi-permeable membrane, the system will swell less and generate a larger internal pressure in return, a system rarely analyzed to date. We use a nonlinear poroelastic theory to model the kinetics of swelling under mechanical constraint. We find the simulation results agree well with our experimental data using hydrogel beads made of a mixture of 3-sulfopropyl acrylate potassium salt and acrylamide, bathed in water. Understanding and predicting the response speed and the actuation stress developed during the swelling of constrained hydrogels can guide the design of polymer-based soft actuators with unusually high strength.

---

1. Introduction

A three-dimensional polymer network can be created by cross-linking long polymers. Immersed in a solvent, this network may absorb molecules from its surroundings, thus creating a mixture called a polymer gel. If the solvent is water, it is referred to as a hydrogel. Bathed in a solvent and under no external constraint, the gel will swell by absorbing new solvent molecules until it reaches equilibrium.

Gels are widely used in micromechanical systems as autonomous flow controllers,¹ pH sensors² or actuators;³ in biology for injury repair⁴ or as carriers for drug delivery⁵ and in soft robotics as stimuli-responsive muscle-like systems.⁶ They are biocompatible, similar to living tissues and adjustable in shapes and chemical properties.⁷ A major flaw exists, however, namely their slow response speed which was shown to be enhanced using electro-osmosis.³

In the field of hydrogel-based soft robotics, there is a growing need for actuators that can generate substantial forces while maintaining precise control, similar to what is achievable with rigid robots. In our prior study,³ we introduced a method to harness large swelling forces using a constrained system consisting of a hydrogel wrapped with a stiff semi-permeable membrane. However, we still lack mathematical understanding of how stress and strain develop in these constrained hydrogels.

The intrinsic swelling properties of gels have been extensively studied. Historically, one of the first attempts was conducted by Tanaka et al.,^8,9 where they considered only the polymer network which swells under the effect of its own osmotic pressure. The effect of the fluid is modeled through the addition of a friction factor in the mechanical balance. Another approach based on Biot's work on soil consolidation,¹⁰ and further developed by Scherer,¹¹ treats the gel as a continuum with the pore pressure as a state variable. Biot also later treated the same problem with a mixture theory approach and showed that these approaches are equivalent.¹² Hui and Muralidharan¹³ showed that extending Tanaka's formalism using a momentum balance equation for the fluid leads to Biot's mixture theory. Nonetheless, all these theories are built within the restrictive framework of linear elasticity, whose assumption of small deformations hardly holds for a swelling ratio of several tens or even hundreds reached by some hydrogels including polyacrylamide-based ones.^14,15

From a thermodynamic point of view, two processes are dominant in the energetic balance of the swelling. The absorption of water molecules makes the hydrogel swell and the polymer network is stretched in return. When an energetic equilibrium is reached between the energetically favourable absorption and the energetically unfavourable stretching, the swelling stops. An important aspect of hydrogels is that, as they swell, their pores expand. Hence, as the size of the pores is used to characterise the permeability, it should increase throughout the process, drastically affecting the dynamics of the swelling. Hong et al.¹⁶ built a theory for large deformations using a thermodynamic approach, based on the works of Flory¹⁷ for the mixing part, Huggins¹⁸ for the stretching part, and Grattoni¹⁹ for the deformation-dependent permeability, among others.

This study aims at analyzing the kinetics of the constrained swelling of a hydrogel coated with a stiff semi-permeable membrane (Fig. 1). In our system, swelling is solely driven by osmosis without external stimuli. When the polymer network absorbs new water molecules through the membrane, the resulting swelling stress of the inner hydrogel is applied on the membrane, which in turn acts as a constraining stress on the gel. Therefore, the constrained hydrogel will only swell slightly, and an important internal pressure will develop. This idea can be used to build hydrogel-based soft actuators, taking advantage of the aforementioned benefits of hydrogels and combining them with unusually high actuating strengths.³ Understanding its physics can allow manufacturers to tailor the strength of their actuators based on their needs and predict its response time, and designers to increase even more their strength or enhance their swelling speed.

Fig. 1 Schematics of free and constrained swelling of the hydrogel. A dry polymer network swells as it absorbs new solvent molecules, forming a polymer gel. Solvent molecules fill and expand the pores of the polymer network. With a stiff semi-permeable membrane coating on the hydrogel, the constrained hydrogel swells less than in the free case.

First, we will tackle the problem of the kinetics of free swelling, namely the swelling of the hydrogel without the membrane. Then, we will consider the same question with a constrained hydrogel, referring to the swelling sphere composed of an inner core of bare hydrogel and a thin outer layer of stiff semi-permeable membrane (Fig. 1). The addition of a membrane restricts the swelling of the hydrogel, leading to smaller pore sizes within the hydrogel, which results in a reduced final swelling ratio and a slower swelling rate. On the other hand, an important internal pressure driven by the osmotic mixing potential of the polymer can grow dramatically, which can be of interest in many applications of soft actuators.

In the following, we begin with the theoretical analysis. Then, further insights into the internal behaviour during the swelling process are presented through simulations. The swelling tests were conducted on both free and constrained cases (Fig. 2) and the simulated results will be compared with experimental data. Then, the experimentally validated model will be used to perform simulations with other values for the different parameters to understand in more depth their importance. Finally, a comparison with the linear theory of Tanaka et al. is provided,^8,9 assuming a linear relationship between stress and strain, will be discussed in view of small deformations in the constrained case.

Fig. 2 Sequential images of swelling hydrogel beads for the (a) free and (b) constrained case. The constrained hydrogel refers to the addition of a layer of membrane that coats the hydrogel. Pictures are taken at logarithmically spaced time intervals and from a top view. Surface instabilities, which we neglected for both cases, can be seen in the beginning of the experiments. Each grid measures 5 mm.

2. Theoretical analysis

This section presents a combination of the theoretical frameworks proposed by Engelsberg et al.¹⁴ and Bertrand et al.,¹⁵ which are to be applied to our experimental systems of hydrogel sphere swelling. We note that the models themselves are based on the fundamental work of Hong et al.¹⁶

2.1 Geometry and microscopic incompressibility

We consider the problem of the swelling of a hydrogel sphere, first freely swelling in water and then constrained by a membrane. Let's consider a spherical sphere of initial radius a, the dry state being our reference state. In this study, we ignore the effects of surface instabilities that may arise during the swelling²⁰ and assume a spherical symmetry during the whole swelling process.

If e⃑_r is the unit vector in the radial direction, R⃑ = Re⃑_r is the position vector in the reference state and r⃑(R⃑, t) = r(R, t)e⃑_r is the position vector in the current state at time t. We define F = ∂r⃑/∂R⃑ as the gradient deformation tensor, thus assumed diagonal with the following eigenvalues, which are the stretches in the three principal directions,

λ_r = ∂r/∂R, λ_θ = λ_ϕ = r/R, (1)

and the following Jacobian determinant, which is the swelling ratio,

J = λ_rλ_θλ_ϕ = λ_rλ_θ². (2)

We state that the increase in volume of the system (polymer and fluid) can only be due to the addition of water molecules from outside the system. This microscopic incompressibility condition leads to

1 + ΩC = det(F) = J = 1/φ, (3)

where Ω = 2.99 × 10⁻²⁹ m³ is the volume of one water molecule, C is the number of water molecules per unit volume of the dry polymer and φ is the volume fraction of the polymer in the mixture, which consists of a cross-linked polymer and the absorbed solvent molecules. We note that J(φ = 1) = 1 in the dry state (reference state). The volume fraction of fluid is thus φ_f = 1 − φ, which ensures that φ + φ_f = 1.

2.2 Free energy density

The fundamental principle behind the swelling of a gel in water is based on two mechanisms. The first one is the spontaneous absorption of water molecules from the solvent to the gel driven by an entropic and/or enthalpic affinity between the polymer and water molecules, a process referred to as the mixing. Due to the additional water molecules, the gel will build an increasing pore pressure which will tend to stretch the cross-linked polymer network and raise the free energy of the system, a process referred to as the stretching. An equilibrium state will be reached when the energetically favourable mixing is exactly compensated by the energetically unfavourable stretching.

The basic assumption of the ideal elastomeric gel theory^16,21 is that the mixing and the stretching part can be totally separated in the computation of the free energy density W. The stretching part will only depend on the three eigenvalues of the deformation tensor F, and the mixing part will only depend on the amount of additional water molecules C:

W(F, C) = W_stretch(F) + W_mix(C). (4)

W _stretch(F) models the elastic contribution to the free energy of an arbitrarily deformed network of cross-linked polymers, ignoring the influence of the fluid on the network. The term is usually derived using a Gaussian-chain elastic model,¹⁸ which yields

where N is the number of chains per unit volume of the dry polymer, k_B is the Boltzmann constant and T is the temperature. This model is usually chosen because it yields satisfactory results, for the swelling ratios involved, while staying rather simple. Experimentally, N can be computed from the value of the shear modulus G of the swollen gel through

W _mix(C) models the contribution to the free energy of the mixing of polymers and water molecules, ignoring the elastic response of the polymer network. The term is usually derived from the Flory–Huggins lattice theory,¹⁷ which yields

where α = Ω_p/Ω is the ratio of the volume of a polymer molecule to that of a solvent molecule and χ is called the Flory–Huggins interaction parameter, or simply the mixing parameter. The first two terms are entropic in nature and the term containing χ is enthalpic.²² If χ < 0, the enthalpic contribution of the mixing is energetically favourable, but it is not if χ > 0. In a first approximation we can consider it as a constant because it is sufficient to grasp the overall dynamics of the swelling.²¹ But experimental studies showed that the mixing parameter should depend on the temperature and the volume fraction of the polymer.^23,24 Small improvements have been made^22,25,26 and a systematic correction has been proposed by Bawendi et al.²⁷ Their simplest extension,²⁸ even within a mean-field approximation, yields

χ = χ₀ + χ₁φ(1 − φ), (8)

with χ₀ = zε/(2k_BT) and χ₁ = zε²/[4(k_BT)²]. Here, z = χ₀²/χ₁ is the coordination number used in the lattice theory and ε is the free energy change when creating a polymer–solvent contact. The influence of the χ₁ term is largest in experiments with a small swelling ratio. While in the free swelling experiments, the polymer volume fraction usually drops below 1%,^14,15 our constrained experiment forces the gel to swell less. Therefore the effects of the χ₁ term might be of importance in our study. This is the reason why this form of the mixing parameter will be used in the present study. Experimentally, χ₀ and χ₁ will be used as fitting parameters, but their value can be controlled by their link to the coordination number z.

2.3 Stress tensor and chemical potential

The condition of microscopic incompressibility leads us to introduce a Lagrange multiplier p in the expression of the free energy density

W(F, C) = W_stretch(F) + W_mix(C) + p[1 + ΩC − det(F)]. (9)

This allows us to compute the nominal stress s_i (i = r, θ or ϕ) and the chemical potential μ as s_i = (∂W/∂λ_i)_T,C,λj (with j ≠ i) and μ = (∂W/∂C)_T,λi (for all i):The true stress σ_i is then computed from the nominal stress as σ_i = s_iλ_i/J. We note that an ionic contribution to the chemical potential should be added. In the case of swelling in pure water, it reduces²¹ to a −k_BTεφ term in the expression of μ. ε and 1/α thus play the same role, numerically speaking. Bertrand et al.¹⁵ neglected ε and used α = 250. Engelsberg et al.¹⁴ neglected 1/α and used ε = 0.008, which is numerically equivalent to a fictitious α = 125. For the sake of simplicity in the formulation, we use α = 300 and neglect ε.

The elastic contribution to the true stress, , and the Lagrange multiplier, p, thus relate to the total stress following the usual form of Biot's theory of poroelasticity: . Therefore, we interpret p as the pore pressure, i.e. the pressure of the liquid inside of the pores of the polymer network. Similarly, the chemical potential is interpreted as the sum of the osmotic pressure, Π, being the mixing contribution, and the pore pressure, p, being the mechanical contribution: μ/Ω = Π + p. This yields

2.4 Kinetic law

We consider a model where the kinetic behaviour of water molecules in the gel is driven by the gradient of chemical potential. Indeed, we study a process driven by a combination of a chemical part (diffusion, represented by the osmotic pressure and following Fick's law) and a mechanical part (flow through porous media, represented by the pore pressure and following Darcy's law) combined in a single term: the chemical potential. Those two laws lead to an equivalent description of the water molecules' displacement in the gel.²⁹ The flux of water molecules thus assumes the form of Darcy's law with a deformation-dependent permeability which accounts for the fact that the mean cross section of the pores of the polymer network increases with the swelling.

The conservation of water molecules yields

with j⃑_R the nominal flux of water molecules and the divergence being taken with respect to the coordinate in the reference state.¹⁶ The true deformation-dependent Darcy's law (with j⃑_r the true flux of solvent) iswhere c is the true concentration of water molecules (c = C/J). We note that we retrieve the formal Darcy's formulation by setting the effective permeability D(φ) = k(φ)k_BT/(Ωcη), with the hydraulic permeability k(φ) related to the mean cross section of the pores and η the viscosity of the solvent, and by identifying μ/Ω to p. The deformation-dependent Darcy's law takes the following form with respect to the coordinates of the reference state:The deformation-dependent permeability takes the usual form¹⁹where the permeability parameter D₀ will be used as a fitting parameter. Engelsberg et al.,¹⁴ Bertrand et al.¹⁵ and Grattoni et al.¹⁹ used values of 1.75, 1.5 and 1.85, respectively for the deformation-dependence parameter β. For our case, the value of β will be fixed at 1.5.

Inputting eqn (3), (14) and (15) into eqn (12) and taking only the radial component, we obtain the time evolution equation for J:

2.5 Mechanical equilibrium

We consider a situation with no external volumetric forces. Thus, the mechanical equilibrium requires that ∇_r·σ = 0 inside the gel. Switching to the reference state and using eqn (11), we obtain

2.6 Governing equation

Inputting eqn (17) into eqn (16), we then proceed to non-dimensionalize the governing equation, introducing the non-dimensionalized variables ρ = R/a and τ = Dt/a². Note that the value of the effective permeability D(J) will be different for the two cases: D should be smaller in the constrained case than in the free case. That being said, the spatio-temporal evolution of the swelling ratio is ruled by the following equationWe relate λ_θ to J by integrating eqn (2) using eqn (1) and obtain the following closing condition:

2.7 Boundary conditions

As will be explained in Section 4.2, the membrane of a constrained gel will be manufactured as a strengthened hydrogel. But theoretically speaking, we only consider its induced stress onto the surface of the core hydrogel as a boundary condition, i.e. we do not aim at solving for the swelling of the membrane and we consider only the mechanical aspect of the swollen membrane. Thus, for both cases we want the boundary of the gel to be in equilibrium with the surrounding fluid, namely water, so that the chemical potential on the boundary μ_b is expressed as μ_b(ρ = 1, τ) = 0.

• In the case of free swelling, where the gel sphere is free to swell in water, we impose that the radial component of the true stress, σ_r, is equal to 0. This leads, using eqn (11) and only on the boundary, to

• In the case of constrained swelling, where the gel is coated by a semi-permeable membrane and immersed in water, we compute the radial stress–radial stretch relationship for the case of an incompressible Neo–Hookean thin shell,³⁰ deriving the stress–stretch relationship from the incompressible simplification of eqn (5). This leads, using eqn (11) and only on the boundary, to:

where G_m is the shear modulus of the membrane and h is its thickness. The condition of no flux of solvent molecules at the center of the sphere is satisfied by imposing (∂J/∂ρ)(ρ = 0, τ) = (∂λ_θ/∂ρ)(ρ = 0, τ) = 0.

2.8 Numerical scheme

To integrate the governing equation, we use the Crank–Nicolson discretization method for the temporal derivative with a first-order finite difference discretization scheme for the spatial derivative. This method being semi-implicit and the governing equation being non-linear, we solve the problem by, at each time step, computing the values of J and λ_θ at every spatial step at once, using the Newton–Raphson method to find the root of the residue function. This method is not subject to any stability condition, so we choose the number of time steps and space steps as respectively N_τ = 1000 and N_ρ = 500, which were found to yield sufficient accuracy in a reasonable computing time. For numerical reasons, the hydrogel is considered as not perfectly dried at t = 0, namely λ = 1.007 and J = 1.02. Furthermore, to let the swelling start, eqn (15) is replaced by D(φ) = D₀(J − 0.9)J^β−1, which is of negligible importance as J starts increasing.

3. Simulation

The governing equation is integrated to understand in more depth the swelling process in both cases. We proceed to plot several variables against the non-dimensional radius ρ (ρ = 0 refers to the center of the sphere and ρ = 1 to its boundary, at all time). We plot them for both the free and the most constrained case of our experiments (G_m = 2.24 MPa), at different time steps, logarithmically spaced between t = 0 and t ⋍ 20 h. Parameters used in this simulation are either fixed, computed from experiments or fitted to experimental data, as will be explained in Section 6.

Fig. 3 shows that the hydrogel swells from its boundary to its center, solvent diffusion being driven by the chemical potential gradient. The gel's radius increases by a factor 3.18 in the free case and by a factor 1.56 in the constrained case (Fig. 3a). The final state in the free case is almost fully swollen with a final value of the fluid volume fraction φ_f ⋍ 0.97 while φ_f ⋍ 0.74 in the constrained case, showing the non-negligible presence of polymers in the final mixture (Fig. 3b). The chemical potential rises until it reaches a value of equilibrium almost equal to that of the solvent in the case of high swelling ratio, and lower than the said value in the constrained case to account for the non-negligible polymer presence in the mixture (Fig. 3c). Hence, there is a jump in chemical potential across the boundary of the gel. It rises rapidly near the boundary, reaching almost immediately its final value, and more slowly inside due to the slow diffusion-like process.

Fig. 3 Evolution of state descriptive values against radius at logarithmically-spaced time steps. (a) Azimuthal stretch, or radius over the initial radius, λ_θ, (b) fluid volume fraction, φ_f and (c) normalized chemical potential, [small mu, Greek, tilde] = μ/(ΩNk_BT). (d) Aforementioned logarithmically-spaced time steps for both cases. Non-dimensionalized time steps are computed through τ = Dt/a², see Section 6 for values of a and D.

Fig. 4 underlines the mechanical properties of the hydrogel and especially the balance between elastic stresses and the pore pressure. The elastic stresses are always and everywhere tensile, because the polymer network has to be stretched in all directions throughout the whole swelling process to leave space for the solvent molecule to flow in the mixture (Fig. 4a and b). Solvent permeating through the boundary of the gel either stays near the boundary to inflate this region or diffuses toward the center. In the free swelling case, the fluid volume fraction increases sharply near the boundary leading to the fast growth of the outer pores and to the comparatively slow diffusion of the solvent to the inner regions. This is because pore sizes which directly affect the flux of the solvent vary across different regions (eqn (14)). The rapid increase in fluid volume fraction in the outer region indicates that the solvent tends to fill the larger pores near the boundary before gradually flowing inward. In the constrained case, the same phenomenon is witnessed, with values of the elastic radial and azimuthal stresses slightly higher but of the same order, but the final state is reached at a lower fluid volume fraction, i.e. at a lower swelling ratio due to the presence of the membrane. As the outer part of the gel rapidly reaches its final state in the constrained case, water molecules have time to flow inward and the process is closer to a uniform volumetric expansion.

Fig. 4 Evolution of true stresses and chemical potential components against radius at logarithmically-spaced time steps. (a) Elastic radial stress , (b) elastic azimuthal stress , (c) pore pressure p and (d) normalized osmotic pressure [capital Pi, Greek, tilde] = Π/(Nk_BT). Legends for the line colors follow those of Fig. 3.

Mechanical balance is ensured by the pore pressure, which grows largely in the constrained case and remains negligible in the free case where the polymer network can be largely stretched and its pores expanded, making room for the newly absorbed solvent. Hence, as the volume of the pores increases, their pressure decreases leading to a comparatively low pore pressure. Whereas in the constrained case, the stiffness of the membrane forces the gel to swell less (the osmotic pressure thus remains non-negligible) and the pores to remain small (the pore pressure thus becomes large), as expected³ (Fig. 4c). The pore pressure drives the transport of solvent molecules from the boundary toward the center of the gel. So does the osmotic pressure, starting at a large negative value in the dry state for both cases because polymer and solvent molecules are strictly separated in space, i.e. are not mixed at all. Then mixing spontaneously occurs, lowering (in absolute values) the value of the osmotic pressure. In the free case, the final state is almost perfectly mixed and the osmotic pressure becomes very small. In the constrained case, a larger value of the osmotic pressure remains even in the final state because of the membrane preventing the gel from swelling completely (Fig. 4d). A qualitative comparison between the chemical potential in Fig. 3c and the osmotic pressure in Fig. 4d hints towards the osmotic pressure being the driving phenomenon of the swelling process (compared to the pore pressure).

Fig. 5 puts emphasis on the mechanical aspect of the actuation, studying the developed true stresses. The combination of the two contributions in Biot's expression of the total stresses in eqn (11) shows that the total radial and total azimuthal stresses get highly compressive in the constrained case due to the large pore pressure (Fig. 5a and b). Fig. 5c shows that on the boundary, the total stress remains equal to 0 in the free swelling case as expected, with the pore pressure compensating the elastic radial stress. The pore pressure decreases in the free case, as the swelling of the pores occurs. Whereas, the total stresses of the constrained hydrogel remain substantially higher due to the presence of the stiff membrane, and consequently the pore pressure in the constrained hydrogel becomes high (∼−1.25 MPa) compared to the negligible elastic stress of the polymer network (Fig. 5a and b). This high compressive pore pressure is the primary contributor to the total stress in the constrained hydrogel. In our previous work,³ the blocking stresses of bare and constrained hydrogel were experimentally measured to be respectively 27 kPa and 1.44 MPa (at a swelling ratio of 1.53 that is similar to our experimental condition 1.56). These experimental results match well with our simulation results. Both cases show a sharp increase of the elastic radial stress at the beginning of the swelling, when the outer layer is swelling while being stuck to the inner core not yet impregnated, followed by a slow decrease as the inner part becomes swollen. In both cases, the gel absorbs solvent molecules due to chemical affinity. Indeed, the study of the decomposition of the chemical potential in terms of osmotic pressure and pore pressure in eqn (11) shows that the driving phenomenon in both cases, even though the increase of the pore pressure in the constrained case helps the diffusion process by adding a mechanical forcing, is the entropic mixing part, i.e. the osmotic pressure (Fig. 5d).

Fig. 5 (a) and (b) Evolution of total stresses against radius at logarithmically-spaced time steps. (c) and (d) Evolution of the total radial stress and the normalized chemical potential on the boundary against time. Line colors follow those defined in Fig. 3.

4. Materials and methods

4.1 Materials

3-Sulfopropyl acrylate potassium salt (SPA; Sigma-Aldrich 251631) and acrylamide (AAm; Sigma-Aldrich A8887) were used as monomers for the core copolymer hydrogel network. N,N′-Methylenebisacrylamide (MBAAm; Sigma-Aldrich M7279) and 1-hydroxy-cyclohexyl-phenyl-ketone (Irgacure 184; Sigma-Aldrich 405612) were respectively used as a cross-linking agent and a photoinitiator. Poly(vinyl alcohol) (PVA, M_w 146 000–186 000, 99+% hydrolyzed; Sigma-Aldrich 363065) and Alginic acid (Alg, medium viscosity; Sigma-Aldrich A2033) were used as monomers for the shell hydrogel that is coated onto the core hydrogels. Calcium chloride (CaCl₂; Sigma-Aldrich C4901) and copper(II) chloride (CuCl₂; Sigma-Aldrich 222011) were used to physically crosslink the Alg chain. Ecoflex (00-10; Smooth-On Inc.) was used as a mold for the core hydrogel. Deionized water filtered with a Direct-Q®3 machine (Merck Millipore) was used to make the electrolyte solution. All chemical reagents were used without further purification.

4.2 Fabrication of hydrogels

Unless otherwise specified, the precursor solution of the core hydrogel P(AAm-co-SPA) was prepared with a molar ratio of 272 : 2727 : 10 : 1 = SPA : AAm : MBAAm : Irgacure 184, where the molar concentration of AAm and SPA monomer were respectively 0.272 M and 2.727 M. The copolymer was used to obtain appropriate swelling stress without rupturing the surrounding membrane, and therefore reliable swelling data. To create spherical hydrogel beads, the mold made of elastomer, Ecoflex™ (SMOOTH-ON 00-10), was used due to its softness and ease of use. The process began by mixing and degassing the A/B components of Ecoflex in a 1 : 1 ratio using a centrifugal mixer (THINKY, ARE-310). Subsequently, the spherical metal beads with a diameter of 4 cm were carefully embedded into the mixed solution, which was then cured into the oven at 100 °C for 1 h. After curing, the metal beads were carefully removed from the elastomer. To produce spherical hydrogels, the precursor solution was injected into the molds using a syringe and cured by exposing them to 365 nm UV irradiation for 10 min. The diameter of the core hydrogel sphere was fixed at 4 mm. To fabricate the core–shell hydrogel, first, the coating solution was prepared by dissolving Alg and PVA with a weight ratio of 2 and 3 wt% respectively in deionized water under stirring and heating (90 °C). After degassing by sonication for 1 h, a clear solution was obtained. The core hydrogel that contains multivalent ions can physically crosslink the Alg chain. To compare the effect of the modulus of the shell (G_m), two different cations were used to coat the core hydrogel. The core hydrogel was briefly immersed in 3.0 M CaCl₂ or 2.0 M CuCl₂ solution for 1 min and any excess solution on the surface was completely dried using an air blower. Subsequently, the spherical hydrogel was immersed into the coating solution (PVA/Alg) allowing the pre-gel to grow from the surface of the core hydrogel within a corresponding coating time. After 10 min of coating, the pre-coated hydrogel was thoroughly cleaned with deionized water to remove any excess of uncrosslinked polymers. To strengthen the shell, the pre-coated hydrogel was annealed in a 120 °C oven for 2 h. The densification of the polymer chain facilitates the additional hydrogen bonding, thereby enhancing the mechanical properties of the shell and preventing swelling in water. For further enhancement via salting-out, the annealed core–shell hydrogel was immersed into 2.0 M ZnSO₄ for 15 min, followed by thorough cleaning with deionized water.

4.3 Tensile tests

The square-shaped membranes were prepared following the above procedure. These membranes were grown on the hydrogels containing multivalent ions. Subsequently, the sheets were cut into dog-bone-shaped specimens with a gauge width of 2 mm. The thickness of each individual specimen was measured using a Vernier caliper. To gather force-displacement data, a universal testing machine (Instron 3343) with 1 kN load cell was employed, maintained at a constant strain rate of 10 mm min⁻¹. Stress–strain curves were generated by dividing the measured force by the initial gauge cross-sectional area and the measured displacement by the initial clamp distance.

4.4 Compressive tests

The disk-shaped core hydrogels were prepared by using an acrylic mold with a diameter of 10 mm and a thickness of 2 mm, followed by UV irradiation for 10 min. The swelling ratios of the specimens were controlled by making the gel swell in a specific volume of deionized water. The specimens were then kept inside a humid chamber for 24 h to reach the equilibrium state without evaporation. To determine the densities of the hydrogels, the weight of each measured sample was divided by the volume calculated from the diameter and thickness of the samples. Compressive force–displacement data were collected using a universal testing machine (Instron 3343) equipped with a 1 kN load cell. The tests were conducted at a constant strain rate of 1 mm min⁻¹. During the tests, the specimens were positioned between the stage and load cell, with mineral oil to fill the gap between them. This setup helped to prevent evaporation and minimize the friction between the specimens and the load cell. The shear modulus was calculated from the true stress–strain curve, where the measured force was divided by the instantaneous cross-sectional area, and the measured displacement was divided by the initial thickness.

4.5 Swelling rate measurement

Before measuring the swelling rate of the spherical core–shell hydrogels, all samples were completely dried at room temperature for a duration of 1 day to ensure consistent starting conditions. The swelling ratio-time curves were then generated by tracking changes in the diameters of the core–shell hydrogels over time. These measurements were facilitated by using a handheld digital microscope (Dinolite AM4815ZT) to capture the time-lapse images while the hydrogels swelled in a bath containing a substantial quantity of deionized water. The swelling ratio was determined by dividing the measured diameter at that specific swelling time by the initial diameter of the core–shell hydrogel.

5. Experimental results and validation

We conducted swelling tests on both unconstrained and constrained hydrogel spheres with an initial radius of a = 2 mm (Fig. 2). The results clearly demonstrate that constrained swelling leads to a slower swelling rate and a smaller final swelling ratio. In the free swelling case, the bare hydrogel reaches its equilibrium state in approximately 5 h, whereas in the constrained case, it continues to swell until approximately 15 h. Furthermore, the final radius of the bare hydrogel is approximately twice as large as that of the constrained hydrogel when measured at identical time intervals.

We first aim at fitting the final swelling ratio for the free swelling case. From our experimental measurements of shear modulus and using eqn (6) at T = 298 K and J = 32, the number of chains per unit volume of dry polymer is fixed at N = 1.849 × 10²⁵ m⁻³ for a shear modulus of G = 23.89 kPa. The value of G is in agreement with usual values from the literature.¹⁴ That of N is roughly one order of magnitude larger than usual values.^14–16 Indeed, the goal of the present study was not to reach a very high swelling ratio and we needed the membrane to withstand the stress induced by the swelling, thus leading to the choice of a denser hydrogel. Hence, the final swelling ratio for the free swelling case now only depends on χ₀ and χ₁, which are fitted to match the experimental results. It gives χ₀ = 0.49 and χ₁ = 0.04. We note that the coordination number z ⋍ 6 is reasonable with respect to that of a simple 3D square lattice. We also underline the fact that, even though the expression of the mixing parameter in this study is different with respect to the constant approximation classically used, the value of χ still lies within the generally accepted range 0.49–0.51 for polyacrylamide hydrogels, with the value of χ being generally higher in the more constrained case, as recommended.²⁴

We next consider the final swelling ratio for the constrained swelling case, with membranes of measured initial thickness h = 250 μm. From tensile tests performed on the membranes, we fitted the true stress–stretch curves to the theoretical expression σ_r(λ_θ) = G_m(λ_θ² − λ_θ⁻⁴). That way, we computed their shear moduli G_m, which are of the order of the MPa, depending on the membrane.

Lastly the permeability parameter D₀ is fitted to match the time scales of the experiments, which leads to D₀ = 5.76 × 10⁻¹¹ m² s⁻¹. The non-dimensionalization step done to obtain eqn (18) is conducted using the value of the effective permeability D = 1.02 × 10⁻⁸ m² s⁻¹ for the free case, which is of the same order of magnitude of values from other studies,^14,15 and D = 4.53 × 10⁻¹⁰ m² s⁻¹ or D = 3.19 × 10⁻¹⁰ m² s⁻¹ for the two constrained cases. These values correspond to their respective final values computed from eqn (15), when the mixture has reached its equilibrium state. We note that the value of D increases throughout the swelling process by roughly two orders of magnitude in the free case and only one in the constrained cases.

The simulation is then compared with the experimental data, as shown in Fig. 6a. We obtain a satisfying fit in terms of prediction of both the final swelling ratio and the kinetics of the two processes. The theoretical model then allows us to compute the actuating pressure inside of the gel (Fig. 6b). As expected, the more constrained the gel is, the larger internal pressure it builds. The negative values of σ_r indicate that the stress is compressive. In all cases, at early times the power law behaviour for the radius (λ_θ − 1) against time resembles that of a diffusive process with an exponent close to 0.5 (Fig. 6c). The seemingly increasing exponent (with respect to the constraint) may be purely artificial, as the fitting to a power law seems to get less and less relevant as the gel gets more constrained.

Fig. 6 (a) Comparison between experiments averaged over three to four samples (circles) and simulations (full lines). (b) Simulated total radial compressive stress produced on the boundary in the same three cases. (c) Strain against time in logarithmic scale. Dashed lines underline the power law behaviour at an early time and are fitted to the experimental values.

6. Parametric analysis

We now want to use our experimentally verified theoretical model to extend our understanding of gels with different parameters through simulations. We will investigate the importance of the shear modulus of the membrane G_m, the density of crosslinks N (which characterises the stiffness of the gel hence its ability to deform) and finally the length of polymers α and the mixing parameter χ₀ (which are related to the chemical affinity between the gel and the solvent). The permeability parameter β (which rules the effect of the swelling of the pores on the gels permeability) does not affect the final swelling ratio. From eqn (11), N characterises the elastic stress and the couple (α, χ₀) the osmotic pressure (at least at the two leading orders when J ≫ 1). We neglect here the effect of χ₁, which is of importance in thermo-sensitive problems²¹ for example, but not in our problem, as will be discussed in Section 7.1.

To evaluate the importance of each parameter, we plot the final azimuthal strain, which characterises the final radius, against the values of the parameters (Fig. 7a–d). As expected, the stiffer the membrane, the lower the final radius (Fig. 7a). N and α play important roles in the free case where the swelling depends on the parameters of the gel alone, but almost none in the constrained case where the membrane prevails (Fig. 7b and c). Finally the value of χ₀ seems to be non-negligible in the two cases, but less relevant than the two previous parameters in the free case (Fig. 7d). We then want to relate the slopes of the log–log curves to power laws derived from the constitutive equations (eqn (20) and (21)). At equilibrium, λ_r = λ_θ = λ_ϕ = λ = J^1/3. Fig. 7e shows that for high values of G_m, λ ∼ G_m^−0.102. Inputting it in σ_r ∼ G_m(λ² − λ⁻⁴) and neglecting the second term, we obtain a dependence, σ_r ∼ G_m^0.796, close to Fig. 7i. It shows that if the membrane gets stiffer, the final radius will decrease slightly but the stress produced will be much larger. Fig. 7f and g show that N and α play almost no role in the constrained case.

Fig. 7 (a)–(d) Effects of G_m, N, α and χ₀ (respectively) on the kinetics (stretch on the boundary against the dimensionless time) and (e)–(h) on the range of swelling (final stretch on the boundary against the value of the studied parameter, log scale). Here, colors are used as visual indicators of the strength of the membrane for each plot. In the first column (a), (e) and (i), the color goes from blue to red as the strength of the membrane increases. In the other three (b)–(d) and (f)–(h), the blue curves represent free swelling cases and red curves constrained swelling cases (with G_m = 2.24 MPa, our experimentally most constrained case). Crosses indicate the values used so far for the different parameters. The range for the values of G_m goes from low values where the swelling behaves almost freely (G_m ⋍ 1.0 kPa) to one order of magnitude more than our constrained cases (G_m ⋍ 10 MPa). Ranges for the values of N and α are centered on the previously fixed value. That for the values of χ₀ is the commonly accepted one (0.49 < χ₀ < 0.51). (i) Final radial stress on the boundary against the value of G_m. For other parameters, the value of the stress remains almost constant.

Now we evaluate which terms are dominant in eqn (21). Indeed the final stretch can be computed by equating the value of the chemical potential on the boundary to that of the solvent and solving for λ, namely: . Using J ≫ 1, we further develop in power series the expression of the first term of the osmotic pressure, ln((J −1)/J) = ln(1 − 1/J) ⋍ − 1/J − 0.5/J² − 0.33/J³ + …, leading to (with G_m = 0 for the free case):

To derive simple power laws, we wish to simplify eqn (22). Let us start with the constrained case for once. From Fig. 5c we can assume that the elastic stress can be neglected compared to the pore pressure. This is supported by Fig. 7f which underlines the marginal importance of the value of N in the constrained case. An analysis of the order of magnitude of the terms of Π shows that the most important term with current values of the parameters and at a typical value of λ ⋍ 1.55 is the first enthalpic mixing term (see also Fig. 7g and h). This leads to a power law of the form:For G_m (Fig. 7e) and T (not shown) respectively, the exponents derived (±0.125) slightly differs from the exponent fitted through simulations (⋍−0.102 and ⋍0.093), but remains fairly close. We attribute the larger discrepancy for the mixing part to the bold assumption to reduce the osmotic pressure to a single term. For the free case, G_m = 0 hence the elastic stress balances the pore pressure. The analysis of the order of magnitude of the contributions of Π with λ ⋍ 3.18 shows that the most important term is the entropic mixing contribution (see also Fig. 7g and h). This leads to a power law of the form:Here, the free case is more complex due to high deformations and the lower stiffness of the total system, and therefore the power law fails to capture the correct exponents. It only predicts a similar dependence of the final strain to N and α, despite their totally different nature.

7. Mixing parameter and linear theory

In this last section, we briefly discuss the importance of the deformation dependence of the mixing parameter and the possibility to make the model simpler without losing most of the simulation's accuracy.

7.1 Flory's constant approximation

Even though Orofino and Flory proposed a virial expansion²³ of Flory's mixing parameter to better fit their experimental results, Flory's first order approximation is still widely used. Fig. 8 shows the results obtained with constant χ. Some parameter values have been changed in this simulation to match the data, namely χ₁ = 0 and α = 220. χ₀ = χ has been kept to its previous value to remain in the 0.49–0.51 range. We note that it closely resembles our previous simulation, and thus conclude that, for the purpose of this study, a constant mixing parameter in the hyperlinear model is sufficient to simulate the process.

Fig. 8 Comparison between experiments (circles), simulation of the hyperlinear model with constant χ (dotted lines) and simulation of the linear model (dashed line). The blue set of data corresponds to the free swelling, and the orange (G_m = 1.14 MPa) and red (G_m = 2.24 MPa) to the constrained cases.

7.2 Tanaka and Fillmore formalism

Because in our study the hydrogels reach rather small values of swelling ratio, we conduct a comparison with the linear theory of Tanaka and Fillmore,⁹ referred to in the next lines as TF's model, using the same notations. This theory assumes a linear relationship between stress and strain. We solve the following problem with a simple first order explicit finite difference scheme:where u is the radial displacement of the gel from its final value, a₀ is the initial radius, a is the final radius, Δa₀ = a − a₀, D is a permeability coefficient, σ_rr is the radial stress and E_m is Young's modulus.

Some critical aspects of this model must be emphasized. First, the theoretical foundations of the model have been pointed out as somewhat artificial.³¹ Second, our experiments largely break the small deformation assumption of linear elasticity. Third, TF's model is not predictive: it uses the initial and final values of the radius to compute mechano-chemical properties of the gel such as its osmotic pressure. With that in mind, we allow ourselves to use TF's model in an exclusively practical way, with the only purpose of simulating the kinetics of the swelling, knowing the final state of our experiments and ignoring all physical considerations. This leads to dashed lines in Fig. 8, where again the final radius in the three cases are known and all other parameters have been fitted to match the experimental data. The simulations match strikingly well the experimental data. Also, the governing equation is much simpler, even allowing for an analytical solution for the free swelling case.⁹ This match to TF's model is only possible because the stretches involved in our experiments are small compared to the usual values for softer hydrogels and because the mechano-chemical aspects can be decently modelled with simple expressions. Larger swelling and finer aspects of hydrogels such as thermo-responsive behaviours or pH-driven phase transitions are absolutely impossible to model with TF's theory: one needs the hyperlinear framework.

8. Conclusion

We built a model to study the kinetics of swelling of constrained hydrogel spheres using a thermodynamic approach based on the theory of Hong et al.¹⁶ We added a constraint to the system to model the membrane around the gel, switching the key aspect of our study from swelling to actuation. That way, we managed to generate a radial stress of the order of the megapascal. Then, we used that model to simulate the process to obtain insights into finer aspects of the system's behaviour. We performed experiments to validate the model, and found that the simulation matched the experimental results satisfactorily well. We also complexified the expression of the mixing parameter χ and it turned out to be of small significance. Due to rather small deformations and simple behaviours, we finally conducted a simulation based on Tanaka et al.'s linear theory which, despite its lack of physical sense and its many fitting parameters, can predict reasonably well the process. The constrained gel, capable of generating unusually high actutation force and supporting large external load, could be employed in various fields, encompassing soft robotics, structural components, impact absorption systems, and implantable structures. Offering a predictive mechanical model for constrained hydrogels could provide guidance for the design of novel soft actuators in those application areas.

Author contributions

H.-Y. K. and J.-Y. S. conceived and supervised the research. All authors contributed to preparing the manuscripts. T. C. built the theoretical model, and performed the simulation and data analysis. H. N. designed the materials, conducted the experiments and collected data.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Research Foundation of Korea (grant no. 2018052541 and 2021017476) via SNU SOFT Foundry. Administrative support from SNU Institute of Engineering Research is acknowledged.

Notes and references

1. D. J. Beebe, J. S. Moore, J. M. Bauer, Q. Yu, R. H. Liu, C. Devadoss and B.-H. Jo, Nature, 2000, 404, 588–590.
2. R. Bashir, J. Z. Hilt, O. Elibol, A. Gupta and N. A. Peppas, Appl. Phys. Lett., 2002, 81, 3091–3093.
3. H. Na, Y.-W. Kang, C. S. Park, S. Jung, H.-Y. Kim and J.-Y. Sun, Science, 2022, 376, 301–307.
4. A. Jain, Y.-T. Kim, R. J. McKeon and R. V. Bellamkonda, Biomaterials, 2006, 27, 497–504.
5. Y. Xiao, W. Xu, Q. Zhu, B. Yan, D. Yang, J. Yang, X. He, S. Liang and X. Hu, Carbohydr. Polym., 2009, 77, 612–620.
6. M. Bassil, J. Davenas and M. E. L. Tahchi, Sens. Actuators, B, 2008, 134, 496–501.
7. K. Deligkaris, T. S. Tadele, W. Olthuis and A. van den Berg, Sens. Actuators, B, 2010, 147, 765–774.
8. T. Tanaka, L. O. Hocker and G. B. Benedek, J. Chem. Phys., 1973, 59, 5151–5159.
9. T. Tanaka and D. J. Fillmore, J. Chem. Phys., 1979, 70, 1214–1218.
10. M. A. Biot, J. Appl. Phys., 1941, 12, 155–164.
11. G. W. Scherer, J. Non-Cryst., 1989, 113, 107–118.
12. M. A. Biot, J. Appl. Phys., 1955, 26, 182–185.
13. C.-Y. Hui and V. Muralidharan, J. Chem. Phys., 2005, 123, 154905.
14. M. Engelsberg and W. Barros, Phys. Rev. E, 2013, 88, 062602.
15. T. Bertrand, J. Peixinho, S. Mukhopadhyay and C. W. MacMinn, Phys. Rev. Appl., 2016, 6, 064010.
16. W. Hong, X. Zhao, J. Zhou and Z. Suo, J. Mech. Phys. Solids, 2008, 56, 1779–1793.
17. P. J. Flory, J. Chem. Phys., 1942, 10, 51–61.
18. M. L. Huggins, J. Chem. Phys., 1941, 9, 440.
19. C. A. Grattoni, H. H. Al-Sharji, C. Yang, A. H. Muggeridge and R. W. Zimmerman, J. Colloid Interface Sci., 2001, 240, 601–607.
20. J. Dervaux, Y. Couder, M.-A. Guedeau-Boudeville and M. Ben Amar, Phys. Rev. Lett., 2011, 107, 018103.
21. M. Quesada-Pérez, J. A. Maroto-Centeno, J. Forcada and R. Hidalgo-Alvarez, Soft Matter, 2011, 7, 10536–10547.
22. B. Qiao and D. Zhao, J. Chem. Phys., 2004, 121, 4968–4973.
23. T. A. Orofino and P. J. Flory, J. Chem. Phys., 1957, 26, 1067–1076.
24. J. Li, Y. Hu, J. J. Vlassak and Z. Suo, Soft Matter, 2012, 8, 8121–8128.
25. R. Koningsveld and L. A. Kleintjens, Macromolecules, 1971, 4, 637–641.
26. P. L. Taylor, Y. Yu and X. Y. Wang, J. Chem. Phys., 1996, 105, 1237–1241.
27. M. G. Bawendi and K. F. Freed, J. Chem. Phys., 1988, 88, 2741–2756.
28. M. G. Bawendi, K. F. Freed and U. Mohanty, J. Chem. Phys., 1987, 87, 5534–5540.
29. S. S. L. Peppin, J. A. W. Elliott and M. G. Worster, Phys. Fluids, 2005, 17, 053301.
30. J. Kiendl, M.-C. Hsu, M. C. Wu and A. Reali, Comput. Methods Appl. Mech. Eng., 2015, 291, 280–303.
31. T. Komori and R. Sakamoto, Colloid Polym. Sci., 1989, 267, 179–183.

---

Footnote

† These authors contributed equally to this work.

---