## Abstract

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
2.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 d*V* be an element of volume, and **N**d*A* be an element of surface, where **N** is the unit outward-pointing normal vector and d*A* is the area of the element. We further denote the Helmholtz free energy of the gel in a volume element as *U*d*V*, the force acting on a surface element as **t**d*A*, and the solvent flux going through a surface element as *I*d*A*. The changing rate of the total free energy of the system, including the potential of the forces and buffers, reads
2.2where is the chemical potential of the solvent molecules in the buffers, 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*:
2.3where *Ω* is the volume occupied by each molecule in pure solvent.

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
2.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 . To ensure positive-definiteness of the dissipation rates, we write them into quadratic forms. Specifically, the dissipation rate due to solvent migration is
2.5and that due to viscous deformation is
2.6where **R**^{s} and **R**^{v} 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
2.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
2.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:
2.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
2.10Together with the molecular incompressibility assumption equation (2.3), equation (2.10) prescribes a constraint between velocity and diffusion flux:
2.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*):
2.12With the aid of integration by parts and the divergence theorem, the variation of the functional is reduced to
2.13where
2.14To minimize the functional , equation (2.13) needs to hold true for arbitrary test fields *δ***v**(**X**), *δ***J**(**X**) and *δμ*(**X**). Equivalently, we have differential equations
2.15and
2.16in addition to constraint (2.11) in the bulk, and
2.17and
2.18on a surface. The mobility tensor **M**^{s} is the inverse of **R**^{s}. 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 **R**^{s} and **R**^{v}. 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 *U*_{n}, and the energy of mixing the solvent with polymer chains U_{m}:U=U_{n}+U_{m} (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):
3.1and the Flory–Huggins model for the energy of mixing (Huggins 1941; Flory 1942):
3.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 *W*_{s} 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,
3.3Similarly, we assume the viscous dissipation rate *W*_{v} to be dependent on the deviatoric part of the strain rate tensor as
3.4with ∂*v*_{i}/∂*x*_{j} 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 and 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 *j*_{i}=*F*_{iK}*J*_{K}/det **F** and ∂*v*_{i}/∂*x*_{j}=*H*_{Kj}∂*v*_{i}/∂*X*_{K}, and the viscosity tensors and .

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)
3.5

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

Using the relation between the true and nominal fluxes, we rewrite equation (2.16) as
3.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}=*s*_{iK}*F*_{jK}/det **F**, we can also write equation (3.6) in terms of true quantities defined in the current geometry:
3.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
3.9

Like all viscoelastic materials, a viscoelastic gel carries an intrinsic time scale
3.10the characteristic time for viscoelastic relaxation in the absence of solvent molecules. The combination of *t*_{0} and the diffusion constant *D* gives rise to an intrinsic length scale
3.11

In contrast to size-independent viscoelastic deformation, diffusion processes are size-dependent. *L*_{0} is a length at which the two processes take comparable periods of time. At a characteristic length much smaller than *L*_{0}, 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 *L*_{0}, the overall evolution process is limited by diffusion. Using the Einstein–Stokes equation, we estimate *L*_{0} approximately as 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 *NΩ* for polymeric gels, *L*_{0} is one to two orders of magnitude larger than molecular size. For physical gels whose reversible cross-links also contribute to the viscosity, *L*_{0} 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 *t*_{0}. 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 *t*≪*t*_{0} the response is dominantly viscoelastic, while at *t*≫*t*_{0} the behaviour is dominantly poroelastic.

## 4. Dimensionless equations and finite-element implementation

By normalizing all stresses by *NkT*, all lengths by *L*_{0}, all fluxes by *L*_{0}/*Ωt*_{0} and time by *t*_{0}, we arrive at the dimensionless form of the governing equations (2.11), (2.14)–(2.16):
4.1
4.2
4.3
and
4.4Here, the hatted symbols represent the dimensionless forms of the corresponding variables, e.g. .

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, , on the surface next to a buffer. In terms of dimensionless variables, the reduced weak form reads 4.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 . It first swells freely and reaches equilibrium with external solvent of chemical potential , 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:
5.1We use representative values of material parameters from typical experimental measurements (Flory 1953): *NΩ*=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 *X*_{1} and *X*_{2} axes are in-plane with the gel film on the gel–substrate interface, and the *X*_{3} 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 evolves with time as the solvent migrates through the top surface of the gel. Setting , equations (4.1)–(4.4) reduce to:
5.2
5.3The boundary conditions are the vanishing chemical potential at the top surface and impermeability 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 plots the through-thickness stretch for a gel subject to a load of normalized magnitude . 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 *t*_{0}, 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 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 , 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 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, , 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 4*a* the overall stretch relative to the initial state as a function of the normalized time, for gels of various sizes. As shown by figure 4*a*, the distinction between the two stages is clearer for a gel with larger dimensionless size , 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.

While the viscoelastic creep is almost independent of the specimen size, as shown by figure 4*a*, the poroelastic deformation in the second stage is highly size-dependent. To show the size dependency, we plot the same creep curves in figure 4*b* by normalizing time with . 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 *L*^{2}/*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.

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 7*a*. 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 6*b*.

With time rescaled by *h*^{2}/*D*, the poroelastic parts of the relaxation curves all collapse to one curve, as in figure 6*b*. 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.

## Acknowledgements

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

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

- This journal is © 2012 The Royal Society