1 Introduction

Elastomeric materials are three-dimensional networks chemically cross-linked by flexible, long polymeric molecules which are laterally attached to one another at occasional points along their length. The majority of elastomeric materials can imbibe large quantities of suitable solvent molecules without disrupting the essential skeletal network structure of the elastomer, forming aggregates known as polymeric gels capable of large and recoverable deformation. The dual characters of a solid and a liquid make the gel ideal for diverse applications: medical devices^[1], drug-delivery systems^[2], tissue engineering^[3], and actuators responsive to physiological cues ^[4-5].

The past century has witnessed numerous theories describing the behavior of diffusion and deformation in gel system. Gibbs and Bumstead^[6] originally formulated field theories of mass transport in elastic solids, assuming that the solid and the fluid have equilibrated. Biot^[7] combined the thermodynamic theory with Darcy's law to model the motion of a fluid in a porous elastic solid, which has been used to analyze phenomena ranging from compaction of soils to deformation of tissues. However, the works of Gibbs and Biot were not tailored for the polymeric gel. Flory and Rehner^[8] proposed a free energy function for the gel, giving an explicit form of the free energy for polymer gels. In their model, the effects of the entropy of stretching the network, the entropy of mixing the network polymers and the solvent molecules, and the enthalpy of mixing were taken into consideration.

In recent years, more complete coupled diffusion-deformation theories have emerged to describe the response of polymeric gels including swelling and drying, squeezing of fluid by applied the mechanical deformation, and forced permeation. Baek and Pence^[9] studied the inhomogeneous deformation of hydrogels in equilibrium under saturated and unsaturated conditions by taking advantage of the theory of Flory and Rehner. Hong et al.^[10] formulated a theory of coupled mass transport and large deformation following Flory and Rehner. Based on the microscopic mixed entropy and Flory-Rehner model, Zhang^[11] investigated the strain-stress relation for the macromolecular microsphere composite hydrogel. The Flory-Rehner model has also been widely used to describe various approximate behaviors of hydrogels by Cai and Suo^[12], Hong et al. ^[13-14], Liu et al.^[15], and so on.

Though classical, the Flory-Rehner model does not give any consideration to the chain extensibility, because it assumes that the end-to-end distance of a chain is smaller than the length of the fully stretched chain. What's more, the cross-linkages between the networks are assumed stable, which is not true in practice. So only when the deformation is small can the Flory-Rehner model successfully predict the mechanics behavior of polymeric gel.

Many attempts have been made to improve the free energy function of polymer gel. Based on the Mooney-Rivlin model of elasticity, Wineman and Rajagopal^[16] developed a specific free energy function for polymer gel. However, the Mooney-Rivlin model is only valid when the strains are less than 100%, because it is based on statistical and empirical arguments. Chester and Anand^[17] adopted the non-Gaussian model, also referred as "8-chain model", to account for the limited extensibility of the polymer chains. But this model does not take into account the effects of entanglements. Mergell and Everaers^[18] developed a tube model to mimic entanglement effects by introducing constraining potentials. There are some other models which describe the entanglement contribution to polymer network chains: constrained-junction model (Flory^[19] and Ronca and Allegra^[20]), primitive path model (Edwards^[21]), diffused-constraint model (Kloczkowski et al.^[22]).

Edwards and Vilgis^[23] presented the famous slip-link model based on the concept of entanglements, which successfully predicted the hardening of the rubber at a high deformation. Higgs and Gaylord^[24], Urayama^[25], and Meissner and Matejka^[26] compared the above-mentioned entanglement model with experimental data and found that the Edwards-Vilgis slip-link model can give the best agreement with experiment. Yan an Jin ^[27-29] and Yan et al.^[30] showed the significant influence of entanglements on the mechanical behavior of neutral polymer gels, and amphoteric pH-sensitive macro-and micro-hydrogels, and the numerical results are close to the experimental data.

To well understand the mechanical properties of gels, we must not only quantify how their structure and function change in response to stimuli, but also quantify their structure and function at a given time. Following Yan and Jin^[27] and Chester et al.^[31], we develop a hybrid free energy function of polymeric gels with the Edwards-Vilgis model and the Flory-Huggins solution theory, in which the entanglements will be taken into consideration. Based on the theory we will develop a finite element method, by which transport and deformation are solved concurrently. The method will be implemented in ABAQUS via a user element subroutine (UEL).

2 Free energy function based on slip-link model

As illustrated in Fig. 1, a fluid-free macroscopically homogeneous polymer gel body occupies a region ℜ at time t₀ in the three-dimensional Euclidean space E, which is often considered as a reference configuration. X is used to denote an arbitrary material point in ℜ. At time t, the dry gel imbibes solvent molecules and deforms to current configuration occupying region r(t)⊂ E. The mapping x=x(X, t) from reference configuration ℜ to current configuration r(t) is considered as a smooth bijection. The deformation gradient is given by

(1)

where I is the unit tensor, and u is the displacement from X to x.

We assume that the deformation gradient has a multiplicative decomposition,

(2)

where F^e represents the stretching and rotation of the swollen network structure, and F^s represents the local distortion due to swelling. The swelling is taken to be isotropic with F^s=λ^sI, where λ^s denotes the swelling stretch. We assume that the swelling stretch is given by

(3)

where Ω is the volume of a mole of fluid molecules, and c_R denotes the fluid concentration measured in moles of fluid per unit reference volume of the dry gels.

Following Flory and Rehner, stretching the network of the polymers and mixing the solvent molecules and the polymers are responsible for the free energy of polymeric gels. So the free energy of the gel takes the form

(4)

where W_mechanical(C, c_R) is the contribution to the free energy due to the deformation of polymer network, W_mixing(c_R) represents the free energy caused by mixing of the solvent with the polymer network, and C=F^TF is the right Cauchy-Green deformation tensor.

2.1 Free energy due to deformation

The free energy due to mechanical stretching is described in terms of the Flory-Rehner framework by Chester et al.^[31]. However, the Flory-Rehner free energy function is based on the simplest affine network model, which does not take entanglements into account. In fact, the real polymer networks have many chain entanglements. Following Yan and Jin^[27], the free energy of a gel due to mechanical stretching, which is derived from the Edwards-Vilgis slip-link model^[23], can be written as

(5)

Here, I_i(i=1, 2, 3) are principal invariants of the right Cauchy-Green deformation tensor C=F^TF; J is the determinant of deformation gradient F; N_s is the concentration of slip links with representative values N_s from 10²⁴ n/m³ to 10²⁷ n/m³; N_c is the concentration of crosslinks with representative values N_c from 10²⁴ n/m³ to 10²⁷ n/m³; α is the inextensibility parameter with representative values α from 0.05 to 0.20; T is the absolute temperature; k is the Boltzmann constant; and η is the slippage parameter.

Next, we suppose that W_mechanical also have an energetic component

(6)

which is a contribution meant to reflect the internal energy associated with volumetric mechanical deformation of the swollen elastomer. Here, J^e and J^s are the determinants of elastic and swelling distortions, namely, J^e=det F^e and J^s=det F^s. K is a bulk modulus of the gel.

To sum up, W_mechanical(C, c_R) can be written as

(7)

2.2 Free energy due to mixing

The long polymers and the small solvent molecules form a liquid solution when the long polymers are not cross-linked. Flory^[32] and Huggins^[33] have derived the following expression for the free energy due to mixing,

(8)

where the first term in the bracket is from the entropy of mixing, and the second one is from the enthalpy of mixing. μ⁰ is a reference chemical potential for the fluid; R is the gas constant; and χ is a dimensionless parameter denoting the enthalpy of mixing.

3 Hyperelastic constitutive equations

The tensor field defined in the domain ℜ × [0, +∞) is called the material field, and the one defined in r(t)× [0, +∞) is called the spatial field. A material field Φ in spatial description is Φ_s(x, t)=Φ(X(x, t), t), and at the same time a spatial field Θ has its material description Θ_m =Θ(x(X, t), t). They satisfy the relationship

Table 1 lists the tensor fields used in this paper.

Consider an arbitrary part P in the reference configuration. After the gel is immerged into the solvent, the part P deforms to P_t at time t. According to the first law of thermodynamics, the time-rate of change of the total energy equals the sum of the rate of work done by the external forces

(9)

where n is the outward unit vector normal to ∂P_t.

The nominal stress tensor, also referred to as the first Piola-Kirchhoff stress tensor, gives the current force per unit undeformed area: P=σ_mF^-T det F. The momentum balance law yields the Boussinesq equation of motion

(10)

where ρ₀ =ρ_m det F.

With the aid of the divergence theorem and equation the right term in Eq. (9) can be integrated by parts as

(11)

where N is the unit vector normal to ∂P.

The time derivative of kinetic energy in Eq. (9) is

(12)

Substituting Eqs. (10), (11) and (12) into Eq. (9), we arrive at

(13)

Because the second Piola-Kirchhoff stress tensor S is symmetrical, we can conclude

(14)

which leads to the constitutive equations of the polymer gels

(15)

(16)

3.1 Constitutive equation for chemical potential

We introduce the following polymer volume fraction:

(17)

Thus, the determinants of elastic and swelling distortions can be written as

(18)

According to Eq. (16), the chemical potential μ is given by

(19)

Substituting Eq. (18) into Eq. (19) yields

(20)

3.2 Constitutive equation for stress

According to Eq. (15), the second Piola-Kirchhoff stress tensor S is given by

(21)

where

The first Piola-Kirchhoff stress can be written as

(22)

3.3 Constitutive equation for fluid flux

Let j be the spatial fluid flux, which denotes the number of the small molecules migrating across unit area in deformed state per unit time. We assume that j depends on the spatial gradient of the chemical potential as follows:

(23)

where m is a scalar mobility coefficient which in general is an isotropic function of the stretch and the fluid concentration.

We assume that the mobility at a given temperature T is given by

(24)

where D represents a diffusion coefficient, and c_m=c_R/J with c denoting fluid concentration measured in moles of fluid per unit deformed volume.

4 Finite element implementation 4.1 Governing partial differential equations

Considering the body , the condition of local force balance requires that the inertia effect is negligible and that the viscoelastic process in the body is fully relaxed, so that

(25)

We assume that no chemical reaction occurs, so that the number of the solvent molecules is conserved^[31], namely,

(26)

Substituting (17) into (26) yields

(27)

The stress boundary condition is given by the traction condition

(28)

for all points which lie on the part of the boundary denoted as , where n is the outward unit vector normal to . A quantity with a "bar" denotes a specified function. Similarly, the displacement boundary condition is given by

(29)

for all points which lie on the part of the boundary denoted as . Here, and are complementary surfaces of the boundary of the body in the sense = ∪ and ∪ =∅.

Another pair of boundary conditions in which the chemical potential is specified on and the fluid flux on are

(30)

Similarly, and are also complementary surfaces of the boundary .

To sum up, in the absence of body forces, the strong forms of the coupled partial differential equations can be written as

(31)

(32)

4.2 Weak forms of governing equations

A weak form^[34] for any set of equations can be constructed by multiplying the equation set by an arbitrary function which has the same free indices as in the set of governing equations, and then integrating over the domain of the problem. Firstly, we express the tensor fields in indicial form, namely,

(33)

where e_i(i=1, 2, 3) are the orthogonal bases of tangent space T_XE at the point X ∈ E. Then, multiplying Eqs. (31) and (32) with weighting functions δu_i and δμ, respectively, and integrating over the volume of the gel, we obtain

(34)

(35)

where δu_i and δμ vanish on and , respectively. In this paper, the Einstein summation convention is used, which means that when an index variable appears twice in a term and is not otherwise defined, it implies summation of that term over all the values of the index.

The finite element approximation to the problem starts by dividing the domain of interest, , into a set of subdomains, also called elements, , such that

The approximation for displacements and chemical potentials is given by

(36)

where N_a are the element shape functions, u_i^a and μ^a are time-dependent nodal displacements and chemical potentials, respectively, and the sum ranges over the number of nodes associated with an element. The weighting functions δu_i and δμ are interpolated by the same shape functions

(37)

Substituting Eqs. (36) and (37) into Eqs. (34) and (35), we can obtain the residuals for the displacement and chemical potential

(38)

(39)

To solve this system of coupled equations, a key step is to get the tangent modulus.

4.3 Tangent modulus

The fourth-order tensor field tangent modulus A is given by

(40)

According to Eq. (22), the first Piola-Kirchhoff stress can be written as

(41)

where

Substituting Eq. (41) into Eq. (40) yields

(42)

In Eq. (42), the second-order tensor fields are

(43)

(44)

(45)

where the Kronecker delta δ_ij are used to modify the subscripts in the coefficients of an expression.

The fourth-order tensor fields are

(46)

(47)

(48)

(49)

Substituting Eqs. (43)-(49) into Eq. (42), we can obtain the components A_ijkl of the tangent modulus.

4.4 Newton procedure

The system of coupled equations (38) and (39) is solved using a Newton procedure. So the tangents of the residuals are required for the iterative Newton solver:

(50)

According to the relationship between the first Piola-Kirchhoff stress P and the true stress σ,

(51)

the displacement residual (38) can be written as

(52)

by changing the variable in the integral, where is the deformed element from B^e.

According to the tangent modulus (40) we have

(53)

where (A_s)_{ij1 kn1} is called the spatial tangent modulus,

(54)

Similarly, the remaining tangents are given by

(55)

(56)

(57)

The integrals in the equations are evaluated numerically by Gaussian-quadrature^[35]. The terms and are computed numerically using the finite difference scheme.

We implement our theory in the finite element package ABAQUS by writing the UEL, by which the displacements and chemical potentials at every point of the gel and at every time are solved numerically.

5 Numerical examples

In this section, the reference chemical potential μ⁰ is taken to be 0.0 J/mol, and the bulk modulus of the gel is K=100 N_ckT. Let ν be the volume per solvent molecule, so Ω =νN_A, where N_A is the Avogadro constant. Considering the relationship between the Boltzmann constant and the Avogadro constant R=kN_A, we have Ω K/(RT)=100νN_c. We have normalized the chemical potential μ by RT and the stress by kT/ν.

5.1 Constrained swelling

The kinetic model developed above for the polymeric gels is evaluated by comparing the numerical results with the experimental data reported by Yoon et al.^[36]. In the experiment, thin layers of poly (N-isopropylacrylamide) hydrogels, bounded to a rigid substrate, were immersed in deionized water and swelling kinetics were monitored using epi-fluorescence microscopy, as shown in Fig. 2. The transient change in thickness of the gel normalized by its initial thickness, Δt/L, is reported.

This experiment can be idealized as the case of one-dimensional constrained swelling, as demonstrated in Fig. 3. The dry gel is used as the reference configuration. Let e₁ and e₂ be the material coordinates in the lateral directions, and e₃ be the material coordinate normal to the layer pointing upward.

The diffusion coefficient in the experiment is D_ex=1.5× 10^-11 m²/s. Bouklas and Huang^[37] adjusted the diffusion coefficient to D=6.87× 10^-11 m²/s by comparing the linear poroelastic and nonlinear theory. The thickness of the initial gel is L. The material properties of the gel are νN_s =0.01, νN_c =0.03, α =0.2, η=0.1, and χ =0.46. The gel is constrained in the lateral directions so that the deformation gradient has the matrix representation

(58)

At time t=0 s, we take φ =0.999 to avoid numerical singularity, and the determinant of elastic distortion will be taken as J^e=1. At room temperature, T=298 K. According to Eq. (20), we have μ/(RT)=-5.809. For the mechanical boundary conditions, the traction condition is that the top face is taken to be traction-free, and the displacement condition is u₃ =0 on the bottom face. For the chemical boundary conditions, we prescribe that there is no flux on the bottom face, and on the top face, the chemical potential is μ⁰/(RT)=0.

In the experiment, different thicknesses for the specimens are used, and normalized simulation and experimental data are presented in Figs. 4 and 5, which shows that the agreement is good.

5.2 Influence of entanglements on deformation behavior

Next, we consider the transient swelling of a layer of a gel in one-dimensional constrained swelling. With respect to Fig. 3, the gel is constrained in the lateral e₁ and e₂ directions, so that the motion of the gel, as it absorbs the solvent and swells, points the e₃ direction. An arbitrary material point X can be expressed as X=X₁e₁ +X₂e₂ +X₃e₃. Here, the thickness of the dry gel will be taken as L=0.001 m, and we will take νN_c =10^-3 and χ =0.1. The entanglement parameters will be taken as νN_s =1× 10^-3, α =0.15, and η =0.1.

Figure 6 demonstrates the stretch λ₃ for the gel when it is deforming. There are two limit states during the deformation. Firstly, at short-time limit, the solvent has no time to migrate into the polymer, so the stretch inside the gel remains unperturbed. Secondly, as time progresses, the solvent molecules migrate into the gel gradually, the stretch evolves to the long-time limit, and the change propagates into the depth of the gel.

To investigate the effect of entanglements on the stretches, we fix the time at t=25 000 s. Firstly, Fig. 7 illustrates that the inextensibility parameter α is inversely proportional to the range of the stretch λ₃. Secondly, as can be seen from Fig. 8, the stretch λ₃ rises with increasing values of N_s/(N_c +N_s). Finally, it can be seen from Fig. 9 that the larger the slippage parameter η is, the larger the stretch λ₃ is.

5.3 Plane strain

The initial dry gel is taken to be 0.002 m long and 0.002 m tall, as can be seen in Fig. 10. The gel is allowed to swell freely, and the plane strain conditions are assumed to prevail. According to the symmetry of the problem we only model one fourth of the gel in simulation. The material parameters will be fixed (νN_s=0.041 5, νN_c =0.041 5, α =0.15, χ =0.1, and η =0.1). The diffusion coefficient D takes the value 5× 10^-9 m²/s.

The initial condition for the chemical potential of the dry gel can be computed by Eq. (20) with ϕ =0.999, T=298 K, J^e=1, and μ⁰/(RT)=0, which is μ/(RT)=-5.809. For the mechanical boundary conditions, the traction condition is that the faces AC and CD are taken to be traction-free, and the displacement conditions are u₁ =0 on face AO and u₂ =0 on face DO.

For the chemical boundary conditions, we prescribe that there is no flux on faces AO and DO, and on faces AC and CD, the chemical potential is μ⁰/(RT)=0.

Figure 11 illustrates the contours of the polymer volume fraction ϕ for a gel swelling freely at different time. Right after the dry gel is put into the solvent, the small molecules do not have time to diffuse, as can be seen in Fig. 11(a). Then the corners which contact with the solvent fully swell first, leading to a bowl-like surface of the gel and generating compressive stress to the gel, which can be seen in Figs. 11(b) and 11(c). The bowl-like surfaces in the swelling gels can be seen in the experiment conducted by Achilleos et al.^[38]. After a long time, the swelling process reaches equilibrium, and the fully swelling gel becomes homogeneous with large volumetric change, as shown in Fig. 11(d).

6 Conclusions

In this paper, we have developed a hybrid free energy function of polymeric gels with the Edwards-Vilgis model and the Flory-Huggins solution theory, by which a theory of concurrent large deformation and diffusion in gels is formulated. Based on the theory, a finite element method for the transient processes of gel undergoing large deformations is developed. Firstly, the simulation results are compared with the experimental results of Yoon et al., and the agreement is good. Secondly, we analyze the influence of the entanglements on the stretch during the constrained swelling. Finally, the contours of the polymer volume fraction for a gel swelling freely are given, by which we are able to provide more insight into the transient swelling phenomena of a gel.

The numerical examples show that the entanglements can exert considerable influence upon the mechanical behavior of the gels. The material microstructure parameters can be adjusted to adapt the experiment. The capability of this model is a valuable tool to research into various potential applications of soft matter.