A visco-poroelastic theory for polymeric gels

Xiao Wang, Wei Hong


A polymeric gel can imbibe solvent and swell. Besides the dilatational mode of deformation, which involves long-range solvent migration, a gel may also undergo volume-conserved deformation. For a macroscopic gel with covalent cross-links, the volume-conserved deformation is usually much faster. However, these two modes are coupled for deformation at the microscopic level and for gels containing physical cross-links or large solvent molecules. In this paper, we seek to formulate a unified theoretical framework for the transient behaviour of polymeric gels to account for both solvent migration and viscoelastic deformation. Under this framework, we further develop a simple material model, and implement it into a finite-element code for numerical calculations. By simultaneously tracking the solvent migration and motion of polymer network, we evolve the inhomogeneous fields of stress and chemical potential. Several initial-boundary-value problems are solved as illustrative examples. For macroscopic gels with low viscosity, the time scales for viscoelasticity and poroelasticity are separated, and the long-term behaviour is just as that predicted by a poroelastic model. For structures or processes involving sizes comparable to the intrinsic length of a material, the viscoelasticity and poroelasticity must be considered simultaneously, especially when studying impact responses.

1. Introduction

A polymeric gel is often described conceptually as a cross-linked macromolecular network filled with species of small molecules, including the solvent and other mobile solutes. Because of this microstructure, a gel is capable of two basic modes of deformation (Hong et al. 2008). First, the stretching of polymer chains accompanied by local rearrangement of small molecules allows a gel to change its shape rapidly, but maintain a constant volume. The gel behaves just like an incompressible elastomer in this mode. The second mode, which involves long-range transportation of the small molecules, is usually slow and size-dependent. The deformation in the second mode allows the gel to swell or shrink in volume. A general inhomogeneous deformation of a gel usually contains both modes. The separate time scales between the two modes allow one to focus more on the dynamics of solvent transportation by modelling a polymeric gel as a poroelastic solid and assuming the deformation in the first mode to be instantaneous (Hirotsu & Onuki 1989; Hong et al. 2008; Hu et al. 2010).

Such an idealized microstructural picture, as well as the approximate treatment that neglects the kinetics of the first mode, has attracted much attention from theoretical aspects, owing to its simplicity (Li & Tanaka 1990; Birgersson 2008; Hong et al. 2008; Quesada-Perez et al. 2011; Cai & Suo 2012). Its practical value has also been realized from the experiments that use the model to extract material properties from simple tests (Lin & Hu 2006; Galli et al. 2009; Cai et al. 2010; Li et al. 2012). However, by neglecting the kinetics of the first mode of deformation, the simplified model also clearly has its limit. There are possible scenarios in which such a simplification is not applicable. The first is the characterization of a local mechanical property of a gel or that of a micro-gel particle, at length scale below micrometres, such as nano-indentation (Kalcioglu et al. 2012) and micro-rheology (Dutta & Belfort 2007; Chen et al. 2010; Patre & Toomey 2010). In this case, the short length for solvent migration reduces the swelling time to a value comparable to the characteristic time of viscoelasticity, and the deformation appears to be dominantly viscoelastic (Chen et al. 2010). Similarly, a gel may have a relatively long viscoelastic relaxation time comparable to that needed for solvent migration through a relatively long distance. This is the case when a gel’s polymer network contains a large amount of reversible cross-links that are capable of dynamic dissociation and re-association, including gels with only physical cross-links such as most biological gels (Roberts et al. 2011), and those with both physical and covalent cross-links such as the semi-interpenetrating network gel (Liu & Chan-Park 2009; Sun et al. submitted). In both cases, the second mode of deformation is no longer the only rate-limiting process, and it is necessary to consider the coupled visco-poroelastic deformation of a gel.

From the mechanics perspective, the study of coupled deformation and interstitial fluid migration originated from the general theory of stress-assisted diffusion in solids developed for flow through porous media (e.g. Gibbs 1878; Biot 1941, 1955, 1956; Rice & Cleary 1976). The theory of the viscoelastic deformation of a fluid-infiltrated polymer was developed to study the migration of a solvent at the transition of a dry polymer to a gel, upon absorption of moisture, known as case II diffusion (Weitsman 1987, 1990; Govindjee & Simo 1993; McBribe et al. 2011). Following this general framework, in this paper, we seek to incorporate viscoelasticity to the theory of poroelastic polymeric gels (Hong et al. 2008; Zhang et al. 2009; Kang & Huang 2010). Viscoelasticity has always been a phenomenological approach for a category of rate-dependent transient processes. While the actual physical processes are, in general, complicated (Kavanagh & Ross-Murphy 1998), and specialized models should be selected for specific materials, here we will simply use a Kelvin–Voigt-type rheology model (Reddy 2008) for demonstration. This simple model is suitable for capturing the fluidic nature of a permanently cross-linked gel in resisting instantaneous loads, together with its recoverable deformation and poroelastic solvent migration in micromechanical tests (Chen et al. 2010). In addition, this model could also be used to study gels with significant but recoverable viscoelastic strains, such as a semi-interpenetrating gel (i.e. a double-network gel with a covalent network and a physical network; Hyland et al. 2012; Sun et al. submitted). More complicated models requiring the description of sophisticated processes, such as the change of the network connectivity in a physical gel (An et al. 2010), are beyond the scope of this paper, but can in general be incorporated into the current theoretical framework.

This paper is organized as follows. Section 2 describes the theoretical framework established from non-equilibrium thermodynamics, accounting for the dissipation due to viscoelasticity and poroelasticity. In §3, the theory is specialized with the Flory–Huggins model for solvent–network interaction, and a Newtonian viscosity for volume-conserved viscous dissipation. A weak form is derived for finite-element implementation. In §4, we give numerical examples of our model, including the behaviour of a gel in confined and unconfined uniaxial deformation and compression by a cylindrical indenter. Specifically, we seek to elucidate the separation of viscoelastic and poroelastic processes through changes in length scales (Kalcioglu et al. 2012; Mohan et al. submitted).

2. Non-equilibrium thermodynamics

A gel is an aggregate of cross-linked polymer chains and mobile molecules. In this paper, we assume that at least some of the cross-links are permanent, so that part of the polymer network remains to be a continuum throughout the entire process. For simplicity, we will only look at simple gels with only one species of mobile molecules: the solvent. To trace the deformation, let us imagine a field of markers attached to the permanently cross-linked part of the polymer network. Each marker is distinguished by its coordinates X in the reference state at time 0. At time t, marker X moves together with the host polymer particle to a location with coordinates x(X,t). The state of deformation of the polymer network is characterized by the deformation gradient Embedded Image2.1

As sketched in figure 1, there are two ways of doing work to a gel: by applying a field of mechanical forces, and by exchanging solvent molecules with external buffers attached to the gel. In this paper, it is assumed that both means of doing work take place at the gel surface. If needed, the body forces or body sources of a solvent can be easily added to the theory. Following the usual approach of the finite-deformation field theory, we use the Lagrange frame and define all field variables with respect to the geometry in the reference state. Let dV be an element of volume, and NdA be an element of surface, where N is the unit outward-pointing normal vector and dA is the area of the element. We further denote the Helmholtz free energy of the gel in a volume element as UdV, the force acting on a surface element as tdA, and the solvent flux going through a surface element as IdA. The changing rate of the total free energy of the system, including the potential of the forces and buffers, reads Embedded Image2.2where Embedded Image is the chemical potential of the solvent molecules in the buffers, Embedded Image is the time derivative of free-energy density, and v=∂x(X,t)/∂t is the velocity. In general, the free energy of the gel is a function of both the deformation gradient F and the nominal concentration of solvent C. Following the usual approach and assuming molecular incompressibility, we assume that the solvent concentration C is related to the volumetric expansion det F, through a dimensionless function g: Embedded Image2.3where Ω is the volume occupied by each molecule in pure solvent.

Figure 1.

Sketch of a polymeric gel under combined chemo-mechanical loads. (Online version in colour.)

Under such an assumption, the free energy is only a function of the deformation gradient, U(F), and its changing rate could be written as Embedded Image2.4

The second law of thermodynamics dictates that the total free energy of a closed system keeps constant or decreases in any physical process. The energy decrease is due to dissipation in irreversible processes. In this paper, we assume that the free energy in a gel is dissipated mainly through two processes: solvent migration and viscous deformation. The rate of dissipation due to solvent migration should be a function of the diffusion flux J, and that, due to viscous deformation, a function of the velocity gradient Embedded Image. To ensure positive-definiteness of the dissipation rates, we write them into quadratic forms. Specifically, the dissipation rate due to solvent migration is Embedded Image2.5and that due to viscous deformation is Embedded Image2.6where Rs and Rv are generalized viscosity tensors, both being symmetric and positive-definite. The two generalized viscosity tensors are functions of state variables as well as their rates of change, especially when describing general nonlinear kinetic processes. Allowing state-dependent viscosity tensors, the formalisms in equations (2.5) and (2.6) are actually quite general.

Our assumption that the energy is dissipated only through solvent migration and viscous deformation leads to Embedded Image2.7for an arbitrary physical process. Substituting equations (2.5) and (2.6) into (2.7), and further combining them with equation (2.4), we arrive at Embedded Image2.8Equation (2.8) holds true for arbitrary admissible fields of velocity and flux, and governs the evolution of the system. Mathematically, equation (2.8) is equivalent to the variational problem of minimizing the following energy functional: Embedded Image2.9With the dissipative terms included, the energy functional Π serves conveniently in handling additional constraints. The variational problem can be regarded as an application of the principal of minimum energy dissipation (Onsager 1931). It should be noted that the existence of the energy functional is not guaranteed for general dissipative processes, and the specific functional form in equation (2.9) is due to the assumed quadratic dissipation-rate expressions in equations (2.5) and (2.6).

The conservation of solvent molecules requires that the time derivative of solvent concentration be related to the diffusion flux as Embedded Image2.10Together with the molecular incompressibility assumption equation (2.3), equation (2.10) prescribes a constraint between velocity and diffusion flux: Embedded Image2.11with g being the derivative of g and H the inverse of deformation gradient tensor F. The constraint can be enforced by adding to the energy functional a term including a Lagrange multiplier μ(X,t): Embedded Image2.12With the aid of integration by parts and the divergence theorem, the variation of the functional Embedded Image is reduced to Embedded Image2.13where Embedded Image2.14To minimize the functional Embedded Image, equation (2.13) needs to hold true for arbitrary test fields δv(X), δJ(X) and δμ(X). Equivalently, we have differential equations Embedded Image2.15and Embedded Image2.16in addition to constraint (2.11) in the bulk, and Embedded Image2.17and Embedded Image2.18on a surface. The mobility tensor Ms is the inverse of Rs. By comparing the field equations and boundary conditions in equations (2.15)–(2.18) to the formulation of field theories from different approaches (Hong et al. 2008), it could be recognized that the tensor s(X,t) given by equation (2.14) is the nominal stress, and the Lagrange multiplier μ(X,t) introduced to enforce solvent conservation is just the chemical potential. Equations (2.11), (2.15) and (2.16) constitutes an initial-boundary-value problem for the unknown fields x(X,t) and μ(X,t). Solutions could be determined upon prescription of proper given-displacement and given-flux boundary conditions, in place of the natural boundary conditions in equations (2.17) and (2.18), wherever applicable.

3. Specific material model

To solve the initial boundary value problem set by equations (2.11), (2.15) and (2.16), we need explicit forms for the free-energy function U(F), and the generalized viscosity tensors Rs and Rv. In this paper, we adopt a simple but widely used model for U(F). The free energy of the gel is divided into two packages, the elastic energy of stretching a polymer network Un, and the energy of mixing the solvent with polymer chains Um:U=Un+Um (Flory & Rehner 1943). Considering the initial stretch λ0 in the reference state, we further use the neo-Hookean form for the free energy of stretching (Flory 1953): Embedded Image3.1and the Flory–Huggins model for the energy of mixing (Huggins 1941; Flory 1942): Embedded Image3.2Here N is the number of polymer chains per unit volume in the reference state, kT is the temperature in the unit of energy and χ is the Flory–Huggins parameter for the enthalpy of mixing.

For simplicity, we assume the dissipative properties of the material to be isotropic in the current state for both solvent migration and viscous deformation. Let dissipation rate Ws to be dependent on the true solvent flux j, i.e. the number of solvent molecules going through unit current area of polymer network per unit time, Embedded Image3.3Similarly, we assume the viscous dissipation rate Wv to be dependent on the deviatoric part of the strain rate tensor Embedded Image as Embedded Image3.4with ∂vi/∂xj being the velocity gradient in the current configuration.

Here m and η are positive parameters, the physical meanings of which will be identified in the discussion follows. By writing Embedded Image and Embedded Image in these specific forms, we have assumed the dissipative properties to be independent of the deformation history. Using the geometric relations between the current state and the reference state, we have ji=FiKJK/det F and ∂vi/∂xj=HKjvi/∂XK, and the viscosity tensors Embedded Image and Embedded Image.

Furthermore, we neglect the mixing volume between a solution and polymer, and let Ω be the volume of each solvent molecule in a gel and in a liquid. Taking the dry state of the polymer network to be the reference, we can relate the nominal solvent concentration C to the volumetric expansion det F as (Hong et al. 2009) Embedded Image3.5

With the specific forms of the free-energy function and the dissipation rates, equation (2.14) is now specialized to Embedded Image3.6

Using the relation between the true and nominal fluxes, we rewrite equation (2.16) as Embedded Image3.7Clearly, equation (3.7) is the linear kinetic law for solvent migration, and parameter m is the mobility. Recalling the relation between true stress and nominal stress, σij=siKFjK/det F, we can also write equation (3.6) in terms of true quantities defined in the current geometry: Embedded Image3.8It is clear that the first term on the right-hand side of equation (3.8) is the stress contribution from elastically stretching the polymer network, the second term is that from mixing solvent and network (i.e. the swelling stresses), the third term is the viscous stress and the last term is the contribution from chemical potential, i.e. the osmotic stress. The physical significance of coefficient η can be readily interpreted as the viscosity of the gel in the current state. The assumptions we have made on energy dissipation, equation (3.4), is equivalent to the stress–strain-rate relation of a Newtonian fluid.

In the limiting case when η=0, the stress given by equation (3.6) or (3.8) reduces to that of a non-viscous gel (Hong et al. 2009). In general, the viscosity of a gel is dependent on the solvent concentration (Busfield et al. 2000; Tsunoda et al. 2000; Tang et al. 2007). For simplicity, we will only consider the case of constant viscosity in the current paper.

The mobility of solvent migration is related to the diffusion constant D through the well-known Einstein relation (Feynman et al. 1963) m=cD/kT, in which c is the true solvent concentration, c=C/det F. The expression for the diffusion flux of solvent, equation (2.16), is specialized to Embedded Image3.9

Like all viscoelastic materials, a viscoelastic gel carries an intrinsic time scale Embedded Image3.10the characteristic time for viscoelastic relaxation in the absence of solvent molecules. The combination of t0 and the diffusion constant D gives rise to an intrinsic length scale Embedded Image3.11

In contrast to size-independent viscoelastic deformation, diffusion processes are size-dependent. L0 is a length at which the two processes take comparable periods of time. At a characteristic length much smaller than L0, solvent diffusion equilibrates much faster than viscous deformation, and the evolution process of a gel is limited by the viscosity. At a characteristic length much larger than L0, the overall evolution process is limited by diffusion. Using the Einstein–Stokes equation, we estimate L0 approximately as Embedded Image for a highly swollen gel with permanent cross-links, of which the viscosity is essentially that of solvent. Here a is the size of a solvent molecule. Taking common values of for polymeric gels, L0 is one to two orders of magnitude larger than molecular size. For physical gels whose reversible cross-links also contribute to the viscosity, L0 could be much larger.

During a non-equilibrium process, the characteristic length of diffusion grows with time. Therefore, in the absence of any geometric dimension, a system may still undergo a transition from viscoelasticity-dominated deformation to diffusion-dominated deformation. The characteristic time for such a transition is given by t0. For example, shortly after a large piece of gel is brought to contact with a buffer of a liquid solvent of different chemical potential, or an impact load is applied to a gel, at time tt0 the response is dominantly viscoelastic, while at tt0 the behaviour is dominantly poroelastic.

4. Dimensionless equations and finite-element implementation

By normalizing all stresses by NkT, all lengths by L0, all fluxes by L0/Ωt0 and time by t0, we arrive at the dimensionless form of the governing equations (2.11), (2.14)–(2.16): Embedded Image4.1 Embedded Image4.2 Embedded Image4.3 and Embedded Image4.4Here, the hatted symbols represent the dimensionless forms of the corresponding variables, e.g. Embedded Image.

To enable finite-element formulation, we submit the kinetic law for solvent diffusion, equation (3.9) back into the functional weak form, equation (2.13), and insist local equilibrium of solvent, Embedded Image, on the surface next to a buffer. In terms of dimensionless variables, the reduced weak form reads Embedded Image4.5The weak form in equation (4.5), together with the dimensionless expression for nominal stress, equation (4.1), is implemented into a finite-element code with the commercial software COMSOL Multiphysics v. 4.2a. In all calculations, implicit time discretization and a direct linear solver are used. The spatial interpolation of the displacement field is one order higher than that of the chemical potential, in order to match the element order between stress and chemical potential. Numerically, the coupled visco-poroelastic model is found to be more stable than transient models with only poroelasticity, partly due to the added damping in both volumetric and deviatoric modes of deformation.

5. Numerical examples

In this section, we will demonstrate the application of the model and the finite-element implementation through several numerical examples involving time-dependent swelling and deformation, which are commonly seen in experiments of gels. We seek to illustrate the process of coupled viscoelastic deformation and solvent migration, and the numerical results are tested against other methods and experimental observations. For simplicity, we will limit the discussion to two-dimensional plane-strain problems, while the procedure could be easily extended to more general three-dimensional processes.

(a) Uniaxial deswelling and creep

Let us consider the one-dimensional problem of an infinitely large layer of gel under uniaxial deformation while immersed in a solvent bath, as sketched in figure 2. The dry polymer layer has thickness Embedded Image. It first swells freely and reaches equilibrium with external solvent of chemical potential Embedded Image, and is then bound to a rigid and impermeable substrate. The same free-swelling state is taken to be the initial state for all following problems. Equating stress to 0 in equation (4.1), we can obtain the initial swelling stretch λ0 from: Embedded Image5.1We use representative values of material parameters from typical experimental measurements (Flory 1953): =10−3 and χ=0.1, which lead to an initial stretch λ0=3.4. A uniform load s is applied on the top surface through a permeable plate. We set up the coordinate system such that the X1 and X2 axes are in-plane with the gel film on the gel–substrate interface, and the X3 axis is in the through-thickness direction pointing towards the gel surface. As a result of the lateral constraint, the lateral stretches are constant, λ1=λ2=λ0. The vertical stretch Embedded Image evolves with time as the solvent migrates through the top surface of the gel. Setting Embedded Image, equations (4.1)–(4.4) reduce to: Embedded Image5.2 Embedded Image5.3The boundary conditions are the vanishing chemical potential at the top surface Embedded Image and impermeability Embedded Image at the bottom. In the case when viscoelasticity is neglected, the last term on the right-hand side of equation (5.2) is omitted, and the reduced one-dimensional initial-boundary-value problem can be easily solved with the shooting method using a finite difference scheme (Hong et al. 2008). As a verification of the finite-element implementation, we solve the same problem with a two-dimensional plane-strain finite-element model by imposing symmetry boundary conditions on the two vertical sides.

Figure 2.

Permeable through-thickness compression of a large gel film bound to a rigid impermeable surface. The nominal stress applied is s/NkT=1. The through-thickness stretch λ3 decreases with time as solvent molecules are squeezed out of the gel. The results from one-dimensional finite difference calculation for a poroelastic model (solid curve, finite difference) are compared with the results from the current model (circles, finite element). (Online version in colour.)

Figure 2 plots the through-thickness stretch Embedded Image for a gel subject to a load of normalized magnitude Embedded Image. The finite-element results are compared against the finite difference calculation for the corresponding poroelastic model. At short time limit when load has just been applied, the solvent has no time to migrate, and stretch is unchanged λ3=λ0. The chemical potential increases to a higher value as a result of applied load, which drives the solvent to diffuse out of the gel through the top surface. Starting from the top surface, stretch ratio of the bulk gel gradually reduces to a much lower value, and equilibriums towards the long time limit. It is noted that as a consequence of the added viscosity, the stretching at the top surface does not reach equilibrium instantaneously as in the pure poroelastic model, but decreases its stretch gradually with a time comparable to t0, as shown by the open circles in figure 2.

Now let us turn to the unconstrained uniaxial deformation of a gel block. The block has a dimension Embedded Image when dry, but is initially swollen and has reached equilibrium with the surrounding solvent. During deformation, the surface of the block is permeable to solvent molecules. At Embedded Image, a constant nominal stress is applied on the top surface. The bottom surface is constrained from vertical displacement, and the lateral sides are taken to be traction-free. Snapshots of the deformation patterns of a gel of size Embedded Image are shown in figure 3. Two stages can be identified during the deformation process: the viscoelastic creep, and the poroelastic creep as a result of solvent migration. In the first stage, Embedded Image, the deformation is homogeneous and the uniform field of chemical potential decreases with the applied mechanical load. The volume of the gel is conserved in this stage, and the gel behaves like an incompressible material. The solvent migration takes place mostly in the second stage, starting from the thin surface layer of the gel. Owing to solvent migration, the deformation of the gel is no longer uniform, and the volume of the gel increases in the second stage. After a long enough time, when the global equilibrium has been reached, the deformation becomes uniform again. To show the creep deformation, we plot in figure 4a the overall stretch relative to the initial state as a function of the normalized time, for gels of various sizes. As shown by figure 4a, the distinction between the two stages is clearer for a gel with larger dimensionless size Embedded Image, i.e. a gel of lower viscosity or larger physical size. In other words, the viscoelastic and poroelastic timescales are separated (Mohan et al. submitted). In the case when only the long-term behaviour is of concern, one may neglect the viscous effect of a macroscopic gel.

Figure 3.

The inhomogeneous field of chemical potential and deformation patterns of a gel block under uniaxial tension. The gel is unconstrained laterally. Snapshots are taken at dimensionless time t/t0 of (a) 0, (b) 103, (c) 105, (d) 107. (Online version in colour.)

Figure 4.

Creep strain in a uniaxially stretched gel plotted as a function of the dimensionless time normalized by (a) the viscoelastic relaxation time and (b) the characteristic diffusion time (Online version in colour.)

While the viscoelastic creep is almost independent of the specimen size, as shown by figure 4a, the poroelastic deformation in the second stage is highly size-dependent. To show the size dependency, we plot the same creep curves in figure 4b by normalizing time Embedded Image with Embedded Image. In the renormalized plot, the creep strains in the poroelastic stage all collapse into one curve. As deformation is limited by the diffusion of solvent into the gel, the characteristic time of the second stage is given by L2/D. On the other hand, for a sample smaller than the intrinsic length scale, the poroelastic and viscoelastic stages overlap with each other, and it appears as if the material has only one relaxation mode.

(b) Indenting a visco-poroelastic gel

Indentation has frequently been used as an experimental method for measuring the elastic and kinetic properties of gels (Cai et al. 2010; Hu et al. 2010; Kalcioglu et al. 2012). Following Hu et al. (2010), we consider a two-dimensional plane-strain problem, in which a gel is indented by a cylindrical indenter, as shown in figure 5. In a relatively short time, the indenter is pressed into the gel to a depth h. With the indenter held in position, the solvent molecules gradually migrate and relax the indentation force. The actual boundary condition for solvent at the gel-indenter interface can be complicated. If the gel is immersed in solvent, the free surface can be treated as a boundary of constant solvent chemical potential before contact; after contact, the interface becomes impermeable to the normal diffusion flux, while the diffusion along the interface may still be present. Here for simplicity, we will only examine the two limiting cases: (i) completely non-interacting (impermeable) boundary and (ii) strongly interacting boundary with very fast reaction or transportation that equilibrates much faster than diffusion, so that a local equilibrium with constant chemical potential can be assumed.

Figure 5.

Sketch of a gel loaded by a cylindrical indenter. (Online version in colour.)

The normalized indentation force per unit length as a function of normalized time is plotted in figure 6. Relatively small indenter sizes are chosen to show both viscoelastic and poroelastic relaxations. To reduce the number of length scales involved, we use proportional dimensions and always indent to one-tenth of the indenter radius, so that both the contact geometry and the diffusion path scale with h. Clearly, the viscoelastic relaxation is independent of the indentation depth, while the poroelastic regime of the relaxation curve is size-dependent. From the short time limit after the viscoelastic deformation to the long time limit, the gel relaxes through redistribution of solvent. The solvent migration is caused by the non-uniform field of chemical potential induced by the stress concentration beneath the indenter, as shown by figure 7a. When the contact surface is taken to be permeable, a significant amount of solvent flux also goes through the surface, causing the gel to relax faster. Owing to the release of solvent from the permeable surface, the indentation force on it is also more relaxed at the long time limit, as shown by figure 6b.

Figure 6.

Relaxation of the force on a cylindrical indenter indenting a large gel. The line force values are normalized by indentation depth h and the modulus of the gel at the free swelling state, E=NkT/λ0. The force is plotted with respect to the time normalized by (a) the viscoelastic relaxation time and (b) the characteristic diffusion time. Results from three relative indenter sizes are compared. (Online version in colour.)

Figure 7.

(a) The inhomogeneous field of chemical potential μ/kT in a gel under a cylindrical indenter, at time tD/h2=10. The arrows indicate solvent migration flux in the gel. (b) Volumetric swelling ratio with respect to the initial state after poroelastic relaxation. (Online version in colour.)

With time rescaled by h2/D, the poroelastic parts of the relaxation curves all collapse to one curve, as in figure 6b. Unlike other modes of deformation (e.g. uniaxial stretch), the relaxation of a gel under indentation exhibits a poroelasticity-dominant stage even for very small indenters or very shallow indentations. This is due to the localized hydrostatic compression during indentation, which can never be relaxed solely by volume-conserved viscous creep. The presence of the poroelastic-dominant stage helps measurements to be carried out with micro-indenters to reduce the experiment times (Kalcioglu et al. 2012). Both the scaling relation and the shape of the relaxation curves during the poroelastic stage agree with experimental observations (Hu et al. 2010; Mohan et al. submitted). The viscoelastic property of the gel would determine its response shortly after the load. For example, the value of the higher contact force shortly after the load is dependent on the competition between the loading rate and the viscous relaxation. A higher loading rate will result in a higher contact force, and the resisting force to an instantaneous deformation is infinity. The viscous response has seldom been reported in experiments, partly due to the small characteristic time for the viscoelasticity of common polymeric gels.

6. Conclusions

The deformation in a gel may be divided into two fundamental modes: the volume-conserved deviatoric deformation, and the dilatational deformation involving solvent migration. While the second mode is often taken to be the limiting process, the energy dissipation and time needed for the first mode is often neglected. For gels with physical cross-links or large solvent molecules, the viscoelastic deformation of the first mode may consume significant time and shall not be neglected. In this paper, a theory for the visco-poroelastic deformation of polymeric gels is formulated. The coupled viscoelastic deformation and solvent transportation processes in polymeric gels are modelled under a continuum framework. For demonstration purposes, a specific material model is further developed using the Flory–Huggins model for swelling, linear kinetics for solvent diffusion and the Kelvin–Voigt-type model for viscoelasticity. Such a material model is suitable for gels with recoverable viscoelastic deformation, i.e. with at least a fraction of permanent cross-links. The specialized weak form has been implemented into a finite-element code via COMSOL Multiphysics v. 4.2a. We illustrate the application of the model through several numerical examples, including the two-stage relaxation of a visco-poroelastic gel under constrained and unconstrained uniaxial deformation, and the indentation of a visco-poroelastic gel using a cylindrical indenter. Problems of various length scales have been calculated to demonstrate the interplay between viscoelastic relaxation and poroelastic solvent migration. The theory is generally applicable to highly viscous polymeric gels, but it is also useful in dealing with the impact response of gels with low viscosity, or the behaviour of microgels.


W.H. acknowledges the support from the National Science Foundation through grant no. CMMI-0900342.

  • Received June 27, 2012.
  • Accepted August 2, 2012.


View Abstract