Abstract

A polymeric gel contains a crosslinked polymer network and solvent. Gels can swell or shrink in response to external stimuli. Two typical kinetic processes are involved during the deformation of gels: the viscoelastic and poroelastic responses. Viscoelasticity of gels is generated from local rearrangement of the polymers, while poroelasticity is generated from solvent migration. The coupled time-dependent behaviors of gels can be formulated by coupling a spring-dashpot model with a diffusion–deformation model. Different combinations of spring and dashpot and different ways of dealing with the coupling between solvent migration and rheological models—either through the spring or dashpot—induce significantly different constitutive behaviors and characteristic time-dependent responses of gels. In this work, we quantitatively study how different rheological models coupled with solvent migration affect the transient behavior of gels. We formulate the visco-poroelastic gel theory for the Maxwell model, the Kelvin–Voigt model, and the generalized standard viscoelastic model. In addition, for generalized standard viscoelastic model, we also discuss the different coupling through the secondary spring or the dashpot. The models are implemented into finite element codes, and the transient-state simulations are performed to investigate the time-dependent deformation and frequency-dependent energy dissipation of different rheologically implemented gel models. The result shows that different combinations of spring and dashpot give the gel solid-like properties and liquid-like properties under different time scales; in addition, the coupling of solvent migration with the dashpot in the rheological model results in restrictions of solvent migration under certain length scales.

1 Introduction

Gels are aggregates of crosslinked polymer networks and solvent molecules (Fig. 1). They have good biocompatibility and have been widely used in a variety of biomedical research fields such as drug delivery [13], cell culture scaffold [46], personalized medicine [7], artificial organs [8,9], and transmission media for ultrasound imaging [10]. Gels have also been used as tissue phantoms for various simulation tests [1113]. The discovery of “volume phase transition” phenomena in gels inspired research on utilizing gels as responsive materials [14]. They have been made responsive to diverse environmental cues such as temperature [14], pH [15], light [16], and electric and magnetic fields [17]. Their ability to swell and deswell depending on external environments and stimuli makes them ideal for intelligent devices. Gels have also demonstrated great potential in uses as actuators [18,19], sensors [20,21], soft robotics [22,23], and wearable electronics [24]. A mathematical model that can describe the complex time-dependent behavior of gels is important.

Fig. 1
A gel is composed of a crosslinked polymer network and the solvent molecules. During the evolution of the gel, the polymer network changes the conformation over a short range and the solvent molecules migrate through the network over a long range. These two processes result in the macroscopic behavior of viscoelasticity and poroelasticity, respectively.
Fig. 1
A gel is composed of a crosslinked polymer network and the solvent molecules. During the evolution of the gel, the polymer network changes the conformation over a short range and the solvent molecules migrate through the network over a long range. These two processes result in the macroscopic behavior of viscoelasticity and poroelasticity, respectively.
Close modal

The time-dependent behavior of gels can be originated from two mechanisms: the rearrangement of polymer chains and the migration of solvent through the network, which give rise to the macroscopic viscoelasticity and poroelasticity, respectively (Fig. 1). For most of the covalently crosslinked hydrogels in their highly swollen states, their time-dependent behaviors are primarily dominated by the poroelastic effect [25]. However, for many other gels such as gelatin, collagen, agarose, and other physically and ionically crosslinked gels, their time-dependent behavior is in general coupled visco-poroelastic response. Although the viscoelastic and poroelastic behaviors are intertwined in general, they have distinct macroscopic characteristics. Viscoelastic time is related to intrinsic molecular characters of the polymer chains of the gel and is independent of any macroscopic length, while poroelastic time is related to the diffusion of solvent and thus scales with the square of a macroscopic length that defines the length of diffusion.

Poroelasticity was originally developed by Biot to describe the coupled diffusion and deformation of porous solids filled with fluid [2628]. The theory has been widely applied to study the phenomena of geomaterials [2932]. Later, a fundamentally similar but conceptually different biphasic theory was developed by Mow and coworkers for studying phenomena of biological materials such as cartilage, bones, and soft tissues [33,34]. In recent years, a general framework based on nonequilibrium thermodynamics was developed to formulate the coupled large deformation and diffusion for various gels [3546]. To describe the concurrent reconfiguration of polymer chains and solvent migration of gels, there have been several noticeable attempts to further incorporate the rheological models into the general framework [4750]. A raising question is how the different rheological models and different coupling methods result in different overall time-dependent responses of gels. A systematic study within single general framework is in need.

Following the general framework of the previous gel theory, in this study, we formulate the coupled visco-poroelastic theory for different rheological models, including the Maxwell (MW) model, Kelvin–Voigt (KV) model, and generalized standard viscoelastic (GSV) model. We propose a new phenomenological element that we call as osmotic container and use it to discuss the different coupling of solvent migration with the different parts of the rheological models. The models are implemented into finite element codes to simulate inhomogeneous transient behaviors of gels. The simulation aims to explore how the different coupled visco-poroelastic models give rise to the different time-dependent deformation and frequency-dependent energy dissipation of gels. This work provides a general guideline for future works in choosing different visco-poroelastic models for the different gel systems.

This paper is organized as follows. Section 2 formulates a general framework for the visco-poroelasticity theory based on nonequilibrium thermodynamics; Sec. 3 derives the specific constitutive relations for several different rheological models; Sec. 4 illustrates the nondimensionalization schemes and finite element implementations; finally, Secs. 5 and 6 discuss the different time-dependent characteristics and frequency-dependent energy dissipation of different models through several boundary value problems.

2 Nonlinear Visco-Poroelasticity

We take the dry state of the network as the reference state and mark every material element in this state with a spatial coordinate X. At an arbitrary time t, the material element moves to a new spatial location x(X,t). The deformation of the gel with respect to its reference state is described by the deformation gradient:
(1)

The amount of solvent in the material element can also change over time. It is described by the nominal concentration C(X,t), the number of solvent molecules in the current state per unit reference volume.

During deformation, the polymer chains in the gel undergo local configurational changes, which are characterized by viscoelasticity. For each viscoelastic process, an internal variable needs to be introduced. For simplicity but without losing generality, here, we formulate the theory by considering only one viscoelastic dissipation, i.e., only one dashpot in the rheological model. We define the corresponding internal variable as the deformation gradient of the dashpot ξ=ξ(X,t).

In Secs. 2.1 and 2.2, we first consider the homogeneous deformation of a representative volume of the material and formulate the general constitutive relations for the material with coupled visco-poroelastic behaviors. We then take a body of the material and discuss the field equations that govern the evolution of the inhomogeneous field variables and the mechanical and chemical boundary conditions of the body.

2.1 Homogeneous States.

Take a material element and consider it as an open thermodynamic system (Fig. 2). There are two ways of doing work on the system: through the mechanical and chemical load. Let Ψ be the Helmholtz free energy of the material element per unit reference volume. Thermodynamics requires that
(2)
where s is the nominal stress and μ is the chemical potential of the solvent. The equality of the equation holds if and only if the thermodynamic system experiences reversible processes.
Fig. 2
The material element composed of the polymer network and the solvent molecules is considered as a thermodynamic system. The material element evolves through a series of homogeneous states. (a) In the reference state, the material element is dry and stress free. (b) In the current state, the material element deforms denoted by a deformation gradient FiK and contains C number of solvent molecules per unit reference volume of the element. In addition, the viscoelastic progression of the material element is characterized by the internal variable ξij. The work is done by the external stress siK and the solvent environment where the chemical potential is μ.
Fig. 2
The material element composed of the polymer network and the solvent molecules is considered as a thermodynamic system. The material element evolves through a series of homogeneous states. (a) In the reference state, the material element is dry and stress free. (b) In the current state, the material element deforms denoted by a deformation gradient FiK and contains C number of solvent molecules per unit reference volume of the element. In addition, the viscoelastic progression of the material element is characterized by the internal variable ξij. The work is done by the external stress siK and the solvent environment where the chemical potential is μ.
Close modal
The Helmholtz free energy of the material element is taken to be a function of the independent valuables, Ψ=Ψ(F,C,ξ). The variation of the Helmholtz free energy can be expressed as follows:
(3)
where δ represents the change of the value during a small-time interval.
Combining Eqs. (2) and (3), the inequality equation becomes
(4)
Equation (4) holds for any arbitrary δFiK, δC, and δξij. We assume that the thermodynamic system undergoes a reversible process if the internal variable is fixed. Consequently, we obtain the following constitutive relations of the material element:
(5)
(6)
If the variation of the internal variable is not zero, the free energy of the system is dissipated by viscoelastic deformation. Therefore, Eq. (4) becomes
(7)
The inequality in Eq. (7) holds for any arbitrary δξij. Therefore, we assume that the evolution of the internal variable is governed by
(8)
where Rikjl is the viscosity tensor, which is symmetric and positive definite to guarantee the inequality of Eq. (7). A commonly used specific form for the evolution of the internal variable is given by Pioletti et al. [51]:
(9)
where η is the viscosity and δkl is the Kronecker delta.

2.2 Inhomogeneous Fields.

Take a body of the material that consists of many material elements (Fig. 3). Each material element evolves in time through homogeneous states as governed by Eqs. (5), (6), and (8). Altogether, they develop the inhomogeneous fields of deformation, solvent concentration, stress, and chemical potential in the body. The evolutions of the fields are coordinated through kinematic constraint, force balance, mass conservation, and transport kinetics. Again, there are two ways of doing work to the body: one is through mechanical load and the other is through pumping molecules into the body. Considering the body as an open system, thermodynamic laws require
(10)
where Ti(X,t) is the nominal traction on the boundary of the gel, Q is the number of solvent molecules that are injected into the gel per unit area of the surface in unit time, and μ¯ is the chemical potential of the solvent in the environment. The integration represents the combination of all material elements in the domain.
Fig. 3
A body of the gel consists of many homogeneous material elements and evolves through a series of inhomogeneous states. The inhomogeneous body is subjected to external traction Ti and submerged in the solvent environment where the chemical potential is μ¯. The evolution of the gel is described by the time-dependent fields: the deformation gradient FiK(X,t), the nominal concentration C(X,t), and the internal variable ξij(X,t).
Fig. 3
A body of the gel consists of many homogeneous material elements and evolves through a series of inhomogeneous states. The inhomogeneous body is subjected to external traction Ti and submerged in the solvent environment where the chemical potential is μ¯. The evolution of the gel is described by the time-dependent fields: the deformation gradient FiK(X,t), the nominal concentration C(X,t), and the internal variable ξij(X,t).
Close modal
Substituting Eq. (3) into Eq. (10), the inequality equation becomes
(11)
Let J(X,t) be the nominal flux of the solvent molecule and NK be the unit normal of the boundary of the body in the reference configuration. We add a term of μJKNKdA and meanwhile deduct this term in Eq. (11), so the inequality equation does not change. We then apply divergence theorem to Eq. (11) and meanwhile use the relations of Eqs. (5) and (6), so we get the following mathematical relation:
(12)
This inequality relation (12) holds for any arbitrary δxi/δt, δξij/δt, μ, and JK, and the equal sign is valid when the process is reversible. When the polymer chains rearrange themselves, they slide against each other, and when the solvent molecules transport through the network, they slide against the polymer chains and themselves. These processes dissipate energy. When these processes do not happen, or mathematically when δξij = 0 and JK = 0, the thermodynamic system undergoes the reversible process, and we have
(13)
The mathematical outcome of Eq. (13), in fact, gives the two well-known conservation laws. The first one is the conservation of linear momentum, i.e., the force balance:
(14)
(15)
The second one is the conservation of mass (the number of solvent molecules here):
(16)
(17)
When the polymers undergo local rearrangement, the sliding between polymer chains dissipates energy, which gives rise to the viscoelastic behavior described by the evolution of the internal variable shown in Eqs. (8) and (9). When the solvent flow through the polymer network, the friction between the polymers and the solvent molecules and among the solvent molecules dissipates energy, which is associated with the poroelastic behavior of the material. According to Eq. (12), the following inequality should be satisfied for arbitrary flux JK:
(18)
To ensure this inequality relation, the simplest mathematical form is to assume:
(19)
where MKL is the mobility tensor, which is symmetric and positive definite, so the inequality in Eq. (18) is automatically satisfied. In this work, the mobility tensor is given based on the Einstein relation [52]:
(20)
where D is the diffusivity, k is the Boltzmann constant, T is the temperature, and H is the inverse of the deformation gradient tensor, H=F1.

3 Different Rheological Models and Different Ways of Coupling With Poroelasticity

Rheological models consist of springs and dashpots. They are connected in different ways to describe the different viscoelastic characteristics of different materials. To model the coupling with the solvent, we use an osmotic container that is represented by the shaded box in Fig. 4. In general, there are two ways to couple the osmotic container to the rheological models: coupling through the spring or coupling through the dashpot. They correspond to different physical processes. For coupling through spring, it considers that as the solvent migrates, it causes an immediate stretch or contract to the spring. Then, the force on the spring is held by the presence of the osmotic container and the dashpot does not deform at all. Physically, it corresponds to the case that the volumetric deformation caused by solvent migration is purely elastic, and it does not induce viscous creep of the polymer chains. In this case, the network viscoelasticity is only related to shear deformation and the volumetric change of the pure polymer components due to the mechanical load rather than the chemical load. For coupling through dashpot, it considers that for the solvent to migrate, dashpot must rearrange. Therefore, the solvent can only migrate little by little together with the evolution of the dashpot. This case considers that for the solvent to migrate into or out of the network, the polymers must rearrange viscoelastically to accommodate it, and the network viscoelasticity is fully coupled with the solvent migration.

Fig. 4
Visco-poroelastic models of gels: (a) the maxwell model with the osmotic container coupled with the spring (MW), (b) the Kelvin–Voigt model with the osmotic container coupled with the dashpot (KV), (c) generalized standard viscoelastic models with the osmotic container coupled with the dashpot (GSV-dashpot), and (d) the generalized standard viscoelastic models with the osmotic container coupled with the secondary spring (GSV-spring).
Fig. 4
Visco-poroelastic models of gels: (a) the maxwell model with the osmotic container coupled with the spring (MW), (b) the Kelvin–Voigt model with the osmotic container coupled with the dashpot (KV), (c) generalized standard viscoelastic models with the osmotic container coupled with the dashpot (GSV-dashpot), and (d) the generalized standard viscoelastic models with the osmotic container coupled with the secondary spring (GSV-spring).
Close modal

In this section, we study three widely used rheological models, the Maxwell model, the Kelvin–Voigt model, and the generalized standard viscoelastic model, and discuss how the poroelasticity is coupled with each model. For the generalized standard viscoelastic model, we consider two ways of dealing with the coupling: through dashpot and through spring individually. The total free energy density is specified for each rheological model. In the following, we use “∼” to denote all the variables and parameters of the springs that are connected to the dashpot in series to distinguish them from those that are connected to the dashpot in parallel.

3.1 Maxwell Gel.

The MW model consists of a spring and dashpot connected in series (Fig. 4(a)). The deformation gradient of the spring is F~, and the internal variable ξ is the deformation gradient of the dashpot. The total deformation gradient is given as follows:
(21)
Here, we consider that the coupling of the solvent is through the spring. It means that as the solvent migrates into or out of the network, the spring is stretched and the stress on the spring is held by the osmotic pressure from the osmotic container. In this model, the volumetric deformation due to solvent migration does not invoke viscoelasticity, and the viscoelastic effects are only related to the shear deformation and the volumetric deformation of the polymer components themselves. The deformation of the spring can be decomposed into the volumetric deformation due to solvent migration and the deformation of the polymer components, that is,
(22)
where Ω is the volume occupied by a solvent molecule and I is the identity tensor. Here, the solvent migration is only considered to be associated with the spring but not the dashpot. If the solvent migration is coupled to the deformation of the dashpot, there will be no deterministic amount of solvent in the equilibrium state. This situation might be relevant to simulate dissolving of polymer in solvent but is out of the interest of this work.
The total free energy of the gel is the summation of the elastic energy of the spring and the mixing energy between the polymers and the solvent [53]. For the stretching energy, we adopt the compressible Neo–Hookean model [54]. For the mixing energy, we follow the Flory–Huggins theory [55,56]. Therefore, the total free energy can be expressed as follows:
(23)
where I~1=F~iKF~iK; J~=detF~; Js=detFs=1+ΩC; J~e=J~/Js; G~ and κ~ are the shear modulus and Lamé constant of the spring, respectively; and χ is the Flory–Huggins parameter that characterizes the mixing enthalpy. A smaller χ represents a higher affinity between the polymer and the solvent.
Substituting Eq. (23) into Eqs. (5), (6), and (8), the nominal stress, chemical potential, and the evolution equation of the MW gel are obtained as follows:
(24)
(25)
(26)
where ζ=ξ1. The force balance relation (14), mass conservation relation (16), the solvent migration kinetics (19), the constitutive relations (24) and (25), and the evolution of the dashpot (26) form the complete set of governing equations for the MW gel.

3.2 Kelvin–Voigt Gel.

The KV model consists of a spring and a dashpot connected in parallel (Fig. 4(b)). In this case, the deformation of the dashpot is identical to the deformation gradient of the spring:
(27)

It means that for the solvent to migrate into the network, the network deforms viscoelastically to make room for the solvent.

With the kinematic relation (27), the free energy of the gel is a function of only two variables: the total deformation gradient F and the solvent concentration C, that is, Ψ=Ψ(F,C). The inequality Eq. (4) for the KV gel becomes
(28)
When δFiK = 0, the dashpot is undeformed, and the material element undergoes a reversible thermodynamic process, so μ = ∂Ψ/∂C (Eq. (6)) still applies for the KV model. Thus, the inequality relation (28) becomes
(29)
Similar to the strategy used in deriving Eqs. (8) and (9), to ensure the inequality in Eq. (29), the nominal stress is taken to be
(30)
where the viscosity tensor for the KV gel is taken to be
(31)
The total free energy of the KV gel is the summation of the elastic energy of the spring and the mixing energy of the polymers with the solvent:
(32)
where I1 = FiKFiK; J=detF; Je = J/Js; G and κ are the shear modulus and Lamé constant of the spring, respectively; and χ is the Flory–Huggins parameter.
Substituting the free energy function (32) into Eqs. (30) and (6), we can obtain the nominal stress and chemical potential of the KV gel as follows:
(33)
(34)

The force balance relation (14), the mass conservation relation (16), the kinetics of solvent migration (19), and the constitutive relations (33) and (34) form the complete set of governing equations for the KV gel. The evolution equation for ξ disappears as the evolution of the dashpot is the same as the spring and the evolution of the spring is expressed in Eq. (31). The deformation of the spring is constrained by the dashpot. For the solvent to migrate into or out of the network, it has to push the dashpot, so the solvent can only get into or out of the network little by little, and it is always accompanied by the viscous dissipation, which is different from the MW gel.

3.3 Generalized Standard Viscoelastic Model.

The GSV model consists of a primary spring and a Maxwell segment that are connected in parallel (Figs. 4(c) and 4(d)). The internal variable ξ is defined as the deformation gradient of the dashpot. Both springs contribute to the elastic energy of the gel. As the solvent migrates into or out of the network, the primary spring must deform to accommodate the volume change due to the change of solvent content. However, there are two ways of dealing with the coupling in the Maxwell segment: coupling the osmotic container either to the secondary spring or to the dashpot.

3.3.1 Generalized Standard Viscoelastic Model With the Osmotic Container Coupled to the Dashpot.

In the GSV viscoelastic model with the osmotic container coupled to the dashpot (GSV-dashpot) model, we deal with the coupling of the volumetric change due to the change of solvent content through the dashpot of the Maxwell segment in the GSV model (Fig. 4(c)). This way of coupling considers that as the solvent migrates into or out of the network, initially both springs deform to accommodate the volume change, but over time, the deformation on the secondary spring is released and the deformation related to solvent transport is carried only by the dashpot deformation in the Maxwell segment. In this case, the solvent migration is coupled with the dashpot movement or physically the viscous rearrangement of the polymers (Fig. 4(c)). The visco-poroelastic behavior of the GSV-dashpot model is similar to the KV gel model, but different from the KV gel model that cannot generate deformation in the initial state, the GSV-dashpot model can generate initial deformation due to the presence of the secondary spring. The deformation gradients of each component in the model are related in the following way:
(35)
where F, F~, and ξ are the deformation gradient of the primary spring, the secondary spring, and the dashpot respectively, and Fs is the deformation of the primary spring associated with the change of solvent content.
The free energy of the GSV gels is the summation of the elastic energy of both springs and the mixing energy between polymers and solvent. In the GSV-dashpot gel, the volume change of the secondary spring equals the volume change of the polymeric components of the secondary spring. The free energy function of the GSV-dashpot gel can be expressed as follows:
(36)
where I1 = FiKFiK; J=detF; I~1=F~iKF~iK; J~=detF~; Js=detFs=1+ΩC; Je = J/Js; G and κ are the shear modulus and Lamé constant of the primary spring, respectively; and G~ and κ~ are the shear modulus and Lamé constant of the secondary spring, respectively; and χ is the Flory–Huggins parameter.
Substituting the free energy function (36) into Eqs. (5), (6), and (8), we obtain the constitutive relations and evolution equation for the GSV-dashpot gel as follows:
(37)
(38)
(39)

The force balance relation (14), mass conservation relation (16), the solvent migration kinetics (19), the constitutive relations (37) and (38), and the evolution equation (39) form the complete set of governing equations for the GSV-dashpot gel model.

3.3.2 Generalized Standard Viscoelastic Model With the Osmotic Container Coupled to the Secondary Spring.

In the GSV model with the osmotic container coupled to the secondary spring (GSV-spring) model, we deal with the coupling of the volumetric change due to the change of solvent content through the secondary spring of the Maxwell segment in the GSV model (Fig. 4(d)). This way of coupling considers that as the solvent migrates into or out of the network, it immediately causes the deformation of both springs and the force on the secondary spring is balanced by the osmotic pressure of the osmotic container. Therefore, the volume change due to the change of the solvent content only induce elastic deformation but not viscoelastic deformation. The viscoelasticity is only associated with shear deformation or the volumetric deformation of the pure polymer components. This model corresponds to the case that solvent migration does not induce viscoelastic creep of the polymer chains. The visco-poroelastic behavior of the GSV-spring model is similar to the MW model, but different from the MW model that cannot carry loads in the equilibrium state, the GSV-spring model can carry loads in the equilibrium state due to the presence of the primary spring. The deformation gradients of each component in the model are related in the following way:
(40)
where F, F~, and ξ are the deformation gradient of the primary spring, the secondary spring, and the dashpot respectively, and Fs is the deformation of both springs associated with the change of solvent content.
The free energy of the GSV gels is the summation of the elastic energy of both springs and the mixing energy between polymers and solvent. In the GSV-spring gel, the volumetric strain of the secondary spring consists of both the volume change due to the change of solvent content and the volume change of the polymeric components themselves due to the mechanical load. The free energy function of GSV-spring gel can be expressed as follows:
(41)
where I1 = FiKFiK; J=detF; I~1=F~iKF~iK; J~=detF~; Js=detFs=1+ΩC; Je = J/Js; J~e=detF~e; G and κ are the shear modulus and Lamé constant of the primary spring, respectively; and G~ and κ~ are the shear modulus and Lamé constant of the secondary spring, respectively; and χ is the Flory–Huggins parameter.
Substituting the free energy function (41) into Eqs. (5), (6), and (8), we obtain the constitutive relations and evolution equation for the GSV-spring gel as follows:
(42)
(43)
(44)

The force balance relation (14), mass conservation relation (16), the solvent migration kinetics (19), the constitutive relations (42) and (43), and the evolution equation (44) form the complete set of governing equations for the GSV-spring gel model.

4 Normalized Governing Equations and Finite Element Implementation

The governing equations are normalized based on the constants and the intrinsic parameters of the initial boundary value problem, as presented in Table 1. The normalized variables and parameters are noted by “^.” The time variables and parameters are normalized by an intrinsic time scale τ = ηΩ/kT, and the length variables and parameters are normalized by an intrinsic length scale L0=Dτ.

Table 1

Normalization of the parameters and variables

Parameter/variableNotationSI base unitNormalization
CoordinateX,xmx^,X^=x,X/L0
Timetst^=t/τ
Nominal stresssN/m2s^=s/(kT/Ω)
ModuliG,κ,G~,κ~N/m2G^,κ^,G~^,κ~^=G,κ,G~,κ~/(kT/Ω)
ConcentrationC1/m2C^=ΩC
Chemical potentialμN mμ^=μ/kT
FluxJ1/(s m2)J^=J/(L0/Ωτ)
Frequencyωrad/sω^=ωτ
Parameter/variableNotationSI base unitNormalization
CoordinateX,xmx^,X^=x,X/L0
Timetst^=t/τ
Nominal stresssN/m2s^=s/(kT/Ω)
ModuliG,κ,G~,κ~N/m2G^,κ^,G~^,κ~^=G,κ,G~,κ~/(kT/Ω)
ConcentrationC1/m2C^=ΩC
Chemical potentialμN mμ^=μ/kT
FluxJ1/(s m2)J^=J/(L0/Ωτ)
Frequencyωrad/sω^=ωτ
By using the normalized variables and parameters, the dimensionless governing equations are formulated for the four rheological models. For all the models, the force balance relation and mass conservation relation are the same.
(45)
(46)

For different models, the constitutive relations and kinetic relations are different, which are listed for each model:

  1. MW gel:

    Kinetics of solvent transportation:
    (47)
    Constitutive law for nominal stress:
    (48)
    Constitutive law for chemical potential:
    (49)
    Evolution equation of the internal variable:
    (50)
  • (2)

    KV gel:

    Kinetics of solvent transportation:
    (51)
    Constitutive law for nominal stress:
    (52)
    Constitutive law for chemical potential:
    (53)
  • (3)

    GSV-dashpot gel:

    Kinetics of solvent transportation:
    (54)
    Constitutive law for nominal stress:
    (55)
    Constitutive law for chemical potential:
    (56)
    Evolution equation of the internal variable:
    (57)
  • (4)

    GSV-spring gel:

    Kinetics of solvent transportation:
    (58)
    Constitutive law for nominal stress:
    (59)
    Constitutive law for chemical potential:
    (60)
    Evolution equation of the internal valuable:
    (61)

The dimensionless Eqs. (45)(61) are implemented in the commercial finite element software comsol multiphysics v. 5.4 through the General Form PDE physics. The direct linear solver mumps was used to solve element stiffness equations.

5 Numerical Examples

In this section, we discuss how different visco-poroelastic models give rise to different macroscopic characteristic behaviors under different loading conditions. To characterize the time-dependent behavior of a material, the typical testing methods are creep test, relaxation test, and dynamic oscillation test. Since the creep behavior of a material is reciprocal to its relaxation response, here, in this section, we will only discuss the creep and dynamic oscillatory loading conditions.

As shown in Fig. 5, we formulate a plane strain problem. The in-plane dimension of the gel in its dry state is denoted by the nondimensionalized length L^×L^. As discussed in Sec. 4, all the length scales in the formula are normalized by an intrinsic length scale of the material L0=Dτ. Bigger L^ represents a physically bigger gel block. The gel is first fully swollen in water. Subsequently, the swollen gel is loaded along the X2 direction through either a constant load or a cyclic load. During the process, the gel is kept underwater, and all surfaces are assumed to be permeable. The right half of the configuration is meshed and computed, and the symmetric boundary conditions are applied on the S1 boundary. In addition, the S3 boundary is set to be stress free and the zero vertical displacement constraint is applied on the S4 boundary.

Fig. 5
The configuration of the gel and the meshing of the domain. Symmetric boundary conditions are applied on S1.
Fig. 5
The configuration of the gel and the meshing of the domain. Symmetric boundary conditions are applied on S1.
Close modal

The material parameters for different models are summarized in Table 2. In the computations, all the material parameters are nondimensionalized, so the absolute values of the shear modulus and Lamé constant are not relevant. Since the solid phase is nearly incompressible, we choose the normalized Lamé constant to be three orders of magnitude of the normalized shear modulus. Another nondimensional material parameter is Flory–Huggins interaction parameter χ. Following literature, a representative value of χ = 0.1 is chosen in the following calculation [35]. Another consideration in choosing the material parameters is that for comparison, we choose the parameters for different rheological models to have the same initial swelling ratio and initial solvent concentration.

Table 2

Material parameters for different rheological models

Rheological modelG^,κ^G~^,κ~^χλ0C0
MWG~^=0.01,κ~^=100.12.2229.977
KVG^=0.01,κ^=100.12.2229.977
GSV-dashpotG^=0.01,κ^=10G~^=0.1,κ~^=1000.12.2229.977
GSV-springG^=0.005,κ^=5G~^=0.005,κ~^=50.12.2229.977
Rheological modelG^,κ^G~^,κ~^χλ0C0
MWG~^=0.01,κ~^=100.12.2229.977
KVG^=0.01,κ^=100.12.2229.977
GSV-dashpotG^=0.01,κ^=10G~^=0.1,κ~^=1000.12.2229.977
GSV-springG^=0.005,κ^=5G~^=0.005,κ~^=50.12.2229.977

5.1 Uniaxial Stretching and Creep.

At the time t^=0+, nominal stress s^=0.005 is applied to the S2 surface of the swollen gel and then kept constant. Meanwhile, the displacement of the gel changes over time. We calculate the normalized length of the gel block in the X2 direction at the X1 = 0 position as a function of time for the different rheological models. The calculations are repeated for gels of different initial size L^ to illustrate the coupled or deconvoluted viscoelastic and poroelastic behaviors. Detailed discussions are as follows:

  1. MW gel

    Before the normal stress is applied, the MW gel has been swollen to the equilibrium state. Thus, the deformation gradient of the spring is F~=λ0I, while the dashpot is undeformed. Under the uniaxial normal stress, the spring will be stretched immediately at t^=0+. The initial stretching ratio of the gel is homogeneous everywhere. The values can be calculated from Eqs. (48) and (49) by setting s^22=s^ and s^11=0, that gives λ1 = 2.096 and λ2 = 2.359. Because of the initial stretching, the chemical potential of the gel inside decreases and becomes smaller than the outside, so the solvent starts to migrate into the gel. The outer part of the gel swells first, and the inner part swells last. Therefore, a field of inhomogeneous stress and strain is developed, and the dashpot evolves. The evolution of the dashpot is induced by the mechanical load that causes the volumetric change of the solid polymer components and the shear stress due to the inhomogeneous fields associated with the transient swelling. We used the finite element method through comsol to solve the problem. For the quantitative analysis, we plot the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the normalized time. In Fig. 6(a), the time is normalized as t^=t/τ, and in Fig. 6(b), the time is normalized as t^/L^2. The different curves represent the calculation results for the gels of different dimensionless size L^. Two distinct deformation processes can be observed: the size-dependent poroelastic deformation and the size-independent viscoelastic deformation. As shown in Fig. 6(a), for smaller gels (e.g., the L^ = 0.001 and L^ = 0.01 curves in Fig. 6), the curves first reach a plateau and then increases to infinity. The time scale for the plateau is different for different sizes of gels. If we further normalize the time by L^2, the plateau regions of the L^ = 0.001 and L^ = 0.01 curves collapse when t^/L^2 is close to 1 (Fig. 6(b)), which confirms that the time scale for the plateau region is dominated by the poroelastic effect. For bigger gels (e.g., the L^ = 0.1 curve in Fig. 6), the poroelastic time scale becomes closer to the viscoelastic time scale. In this case, the plateau region disappears and shows only one time scale relevant to the steep increasing region. As shown in Fig. 6(a), the L^ = 0.001, L^ = 0.01, and L^ = 0.1 curves all collapse when the nondimensional time t^ is close to 1, corresponding to the viscoelastic dominant region (Fig. 6(a)). If the size of the gel is even bigger (e.g., the L^ = 1 and L^ = 10 curves in Fig. 6), the viscoelastic and poroelastic time scales are completely intertwined, and the curves show no plateau but only the steep increasing region. But different from the L^ = 0.001, L^ = 0.01, and L^ = 0.1 curves, the time scale for the steep region, in these cases, is not solely determined by the viscoelastic time scale, but also influenced by the poroelastic time scale. Therefore, the time scale for the steep increasing region appears to be bigger than in other cases. Finally, we would like to clarify, the MW gel model should be strictly speaking a visco-poroplastic model. Under constant stress, the dashpot in the MW model can continuously deform to infinity. This model is more suitable for gels with fluid-like behaviors in the long time scale [57].

  2. KV gel

    In the KV model, the overall deformation gradient is identical to the deformation gradient of the dashpot. As the external stress is applied, the gel cannot deform immediately. Therefore, the initial stretching ratio of the gel remains unchanged from λ0 everywhere. Over time, the dashpot deforms and the solvent migrates. The same as before, we plot the normalized height of the gel block l^2(X1=0)/L^ as a function of normalized time t^ and t^/L^2 in Figs. 7(a) and 7(b), respectively. For smaller gels (e.g., the L^ = 0.1 and L^ = 1 curves in Fig. 7), the curves only show one plateau and thus one time scale. This is because for smaller gels, the solvent migrates much faster than the creep of the dashpot, and the time confining factor of the KV model is the viscoelastic time scale related to the dashpot. The magnitude of the deformation is contributed by both solvent migration and the viscoelastic creep, which is 2.429 as shown in Fig. 7. For larger gels (e.g., the L^ = 10 and L^ = 100 curves in Fig. 7), the poroelastic time scale becomes much larger than the viscoelastic time scale, and the curves show two distinct plateaus, corresponding to the viscoelastic and poroelastic time scales, respectively. The first plateau at around 2.360 corresponds to the viscoelastic deformation because the time to reach this plateau is independent of the size of the gel (Fig. 4(a)). The second plateau from 2.360 to 2.429 corresponds to the poroelastic deformation because the time to reach the second plateau scales with the square of the gel size L^2 (Fig. 7(b)). Different from the MW gel that exhibits long-term fluid-like behavior, the KV gel exhibits short-term fluid-like behavior.

  3. GSV-dashpot gel

    Similar to the KV gel, the solvent migration in the GSV-dashpot gel is coupled with the movement of the dashpot. However, as the dashpot is connected to a spring in series, the GSV-dashpot gel can deform at t^=0+. From Eqs. (55) and (56), the initial deformation of the gel is calculated: λ1 = 2.180 and λ2 = 2.265. The same as before, we plot the normalized height of the gel block l^2(X1=0)/L^ as a function of normalized time t^ and t^/L^2 in Figs. 8(a) and 8(b), respectively. As shown, the time-dependent creep behavior of the GSV-dashpot gel is similar to that of the KV gel. The solvent migration of the gel is restricted by the movement of the dashpot. Therefore, the separated poroelastic and viscoelastic regions occur for larger gels. In addition, as time goes to infinity, the secondary spring in the dashpot does not carry any mechanical load. Thus, as the initial swelling ratio of the free swollen gel for both GSV-dashpot gel and the KV gel is identical, the deformation of the GSV-dashpot gel in the equilibrium state is also similar to that of the KV gel l^2(X1=0,t^=)/L^=2.429. The major difference between the GSV-dashpot gel and the KV gel is that at t^=0+, the KV gel is entirely fluid like, while the GSV-dashpot gel can deform as an elastic solid. Therefore, the initial height of the gel in the KV gel model is 2.222, while it is 2.265 in the GSV-dashpot gel model.

  4. GSV-spring gel

    In the GSV-spring gel, the solvent migration is coupled to the motion of the secondary spring. Therefore, the viscoelastic deformation does not constrain the poroelastic deformation in the GSV-spring gel. As the stress is applied, the two springs deform instantaneously. The initial stretching value is given by λ1 = 2.096 and λ2 = 2.359. The same as before, we plot the normalized height of the gel block l^2(X1=0)/L^ as a function of normalized time t^ and t^/L^2 in Figs. 9(a) and 9(b), respectively. Different from all other models, the GSV-spring gel model predicts two distinct time scales for both small gels (the L^ = 0.01 and L^ = 0.1 curves in Fig. 9) and big gels (the L^ = 10 and L^ = 100 curves in Fig. 9). The difference between small and large gel behaviors is that for small gels, the solvent migration time is much shorter than the viscoelastic time, and the curves show two plateaus with the first plateau corresponding to the poroelastic creep and the second plateau corresponding to the viscoelastic creep. It is also confirmed by the scaling relation as the second part of the curves collapse in Fig. 9(a), but the first part of the curves collapse in Fig. 9(b) when the time is normalized by L^2. As the gel reaches the equilibrium state, both springs in the GSV-spring gel are saturated and stretched. Only for the gel size of L^=1, the viscoelastic time is the same as the poroelastic time scale, and the blue curve in Fig. 9 only shows one plateau and one time scale. Comparing with the GSV-dashpot gel where the dashpot is saturated with solvent, the GSV-spring gel generates bigger deformation in the equilibrium state, which is 2.590 compared with 2.429 in the GSV-spring model.

Fig. 6
Creep of the MW gel: the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the (a) dimensionless time t^=t/τ and (b) normalized dimensionless time t^/L^2
Fig. 6
Creep of the MW gel: the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the (a) dimensionless time t^=t/τ and (b) normalized dimensionless time t^/L^2
Close modal
Fig. 7
Creep of the KV gel: the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the (a) dimensionless time t^=t/τ and (b) normalized dimensionless time t^/L^2
Fig. 7
Creep of the KV gel: the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the (a) dimensionless time t^=t/τ and (b) normalized dimensionless time t^/L^2
Close modal
Fig. 8
Creep of the GSV-dashpot gel: the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the (a) dimensionless time t^=t/τ and (b) normalized dimensionless time t^/L^2
Fig. 8
Creep of the GSV-dashpot gel: the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the (a) dimensionless time t^=t/τ and (b) normalized dimensionless time t^/L^2
Close modal
Fig. 9
Creep of the GSV-spring gel: the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the (a) dimensionless time t^=t/τ and (b) normalized dimensionless time t^/L^2
Fig. 9
Creep of the GSV-spring gel: the normalized value of the vertical height of the gel block at X1 = 0 position l^2(X1=0)/L^ as a function of the (a) dimensionless time t^=t/τ and (b) normalized dimensionless time t^/L^2
Close modal

5.2 Dynamic Oscillatory Loading.

As a gel deforms, the elastic energy is stored in the material, but meanwhile, energy is also dissipated because of the viscoelastic and poroelastic processes. When the gel is under cyclic loading, the energy dissipated through one cycle depends on the frequency of deformation. The frequency-dependent energy dissipation is related to the intrinsic viscoelastic and poroelastic properties of the gel. In this section, we quantitatively study the dynamic response of the gel when it is subjected to dynamic cyclic loading for several proposed visco-poroelastic models. The geometry of the gel is shown in Fig. 5. A cyclic load in the form of s^=s^0sinω^t^ is applied to the S2 boundary, where s^0 is a small value. The displacement of the gel in response to the load is also cyclic but with a phase lag δ due to the viscoelastic and poroelastic dissipation, i.e., u^=u^0sin(ω^t^δ). The phase lag amount is positively related to the amount of energy dissipated within one loading cycle in the material. Therefore, to study the dependence of the energy dissipation on the actuation frequency, we quantify the phase lag δ as a function of the normalized frequency ω^ for the different visco-poroelastic models proposed in this work.

  1. MW gel

    The phase lag of MW gel in response to different loading frequencies is shown in Fig. 10. In the low-frequency region, which corresponds to the long-time response of the material, the phase lag of the MW gel is π/2 (Fig. 10(a)). It means that for the MW gel, the work done by the external force is entirely dissipated under one cycle of loading and unloading at extremely low frequencies. Under these frequencies, the MW gel behaves as a fluid that does not store elastic energy. At higher frequencies, local peaks appear in the phase lag curves for the small-size gels. For smaller gels, the solvent migration takes a much shorter time and is separated from the viscoelastic time scale. As these local peaks collapse when the angular frequency is normalized as ω^L^2 (Fig. 10(b)), it means that the energy dissipation corresponding to these local peaks is due to the poroelastic effect of the gels. At extremely high-frequency regions, eventually, the phase lag becomes zero for all the gels of any sizes. It means that the loading is so fast that neither the migration of the solvent nor the movement of the dashpot (i.e., the reconfiguration of the polymers) has enough time to happen during each loading cycle, so the MW gel behaves as an elastic body, and no energy is dissipated during the loading process in the extremely high-frequency region.

  2. KV gel

    The phase lag of the KV gel in response to different loading frequencies is shown in Fig. 10. Different from MW gel in which the dashpot is connected with the spring in series, the spring does not support load in the long time region and the gel’s long-time response is fluid like, in the KV model, the dashpot is connected to the spring in parallel, so the dashpot has to deform to generate any deformation, the time-dependent behavior of the KV gel in a time region that is shorter or close to the viscoelastic time scale is entirely determined by the viscoelastic time scale and the gel is fluid like in the short-time response. As shown in Fig. 11(a), the phase lag in the high-frequency region for different sizes of gels is π/2, which shows that the work done by the external force is entirely dissipated. In the lower frequency region, local peaks appear in the phase lag curves for the larger gels. For larger gels, the solvent migration takes a much longer time and is separated from the viscoelastic time scale. As these local peaks collapse when the angular frequency is normalized as ω^L^2 (Fig. 11(b)), it means that the energy dissipation corresponding to these local peaks is due to the poroelastic effect of the gels. At extremely low-frequency region, eventually, the phase lag becomes zero for all the gels of any sizes. In this long-time limit, the load in the KV gel is entirely supported by the spring and the material behaves as an elastic solid, and no energy is dissipated during the loading and the unloading processes.

  3. GSV-dashpot gel

    The phase lag of the GSV-dashpot gel in response to different loading frequencies is shown in Fig. 12. As discussed previously, the solvent migration in the GSV-dashpot gel is coupled with the deformation of the dashpot. Therefore, solvent transport in a short-time scale is dominated by the viscoelastic time. As shown in Fig. 12(a), for smaller gels, the poroelastic dissipation of the gel is fully coupled with the viscoelastic dissipation and it only shows one peak in the phase lag plot. For larger gels, the poroelastic time scale is much longer than the viscoelastic time scale, and it shows two distinct peaks in the phase lag plot corresponding to the poroelastic and viscoelastic responses. Those peaks that collapse in the δω^ plot are the viscoelastic peaks, while those that collapse in the δω^L^2 are the poroelastic peaks (Fig. 12(b)). The GSV-dashpot gel can be regarded as a generalized KV gel model, as the solvent migration is constrained by the movement of the dashpot, and the poroelastic energy dissipation is restricted by the viscoelastic dissipation rate in both models. However, the GSV-dashpot model does not exhibit a complete fluid-like behavior at any frequency region because of the presence of two springs in parallel that guarantees the pure elastic behavior of the gel in any time scales. For the same reason, instead of a phase lag of π/2 at high-frequency limit as in the KV gel, the phase lag of the GSV-dashpot gel is zero at the high-frequency limit.

  4. GSV-spring gel

    The phase lag of the GSV-spring gel in response to different loading frequencies is shown in Fig. 13. In the GSV-spring gel, the solvent migration is coupled with the deformation of the springs, rather than the dashpot, so the time scale for viscoelasticity and poroelasticity is decoupled. As shown in Fig. 13(a), for both large gels and small gels, there appear two peaks in the phase lag plot, and only for the gel of size L^=1 where the viscoelastic time exactly equals the poroelastic time, there appears one peak in the phase lag plot. Those peaks that collapse in the δω^ plot are the viscoelastic peaks, while those that collapse in the δω^L^2 are the poroelastic peaks (Fig. 13(b)).

Fig. 10
Phase lag of the MW gel with respect to the (a) dimensionless frequency ω^=ωτ and (b) normalized dimensionless frequency ω^L^2
Fig. 10
Phase lag of the MW gel with respect to the (a) dimensionless frequency ω^=ωτ and (b) normalized dimensionless frequency ω^L^2
Close modal
Fig. 11
Phase lag of the KV gel with respect to the (a) dimensionless frequency ω^=ωτ and (b) normalized dimensionless frequency ω^L^2
Fig. 11
Phase lag of the KV gel with respect to the (a) dimensionless frequency ω^=ωτ and (b) normalized dimensionless frequency ω^L^2
Close modal
Fig. 12
Phase lag of the GSV-dashpot gel with respect to the (a) dimensionless frequency ω^=ωτ and (b) normalized dimensionless frequency ω^L^2
Fig. 12
Phase lag of the GSV-dashpot gel with respect to the (a) dimensionless frequency ω^=ωτ and (b) normalized dimensionless frequency ω^L^2
Close modal
Fig. 13
Phase lag of the GSV-spring gel with respect to the (a) dimensionless frequency ω^=ωτ and (b) normalized dimensionless frequency ω^L^2
Fig. 13
Phase lag of the GSV-spring gel with respect to the (a) dimensionless frequency ω^=ωτ and (b) normalized dimensionless frequency ω^L^2
Close modal

The characteristic behavior of the GSV-spring gel is significantly different from that of the GSV-dashpot gel. Unlike the GSV-dashpot gel that can be regarded as a generalized KV gel, the GSV-spring gel is a generalization of the MW gel. In both GSV-spring gel and MW gel, the poroelastic dissipation is decoupled from the dashpot movement, and thus, the poroelastic dissipation can be completely separated from the viscoelastic time scale by controlling the size of the gel. Furthermore, like the GSV-dashpot gel, the existence of the two parallel springs of the GSV-spring gel guarantees the elastic behavior of the gel, and thus, the phase lag of the GSV-spring gel is zero at both the extremely high and extremely low-frequency region as the time scale is far away from both the viscoelastic and poroelastic time scales.

6 Conclusion

In this paper, we formulated the coupled visco-poroelasticy for gels based on different rheological models. We adopted the Maxwell model, the Kelvin–Voigt model, and the generalized standard viscoelastic model for the rheological behaviors of gels and discuss how the coupling of solvent migration with spring or dashpot motions influence the gel’s constitutive behaviors. We used a compressible neo-Hookean hyperelastic solid model for the deformation of the springs and the Flory–Huggins mixing theory for the interaction of the solvent and polymers. All models are implemented in the finite element codes in comsol. The time-dependent behavior of the gel is studied via the simulation of the uniaxial creep of the gel, and the frequency-dependent energy dissipation is studied via the dynamic uniaxial loading of the gel. The solvent migration could be either coupled with the dashpot movement or unrelated to the dashpot. The first case generates a strong couple of poroelastic deformation with the viscoelastic deformation. In this case, the solvent migration could be restricted by the dashpot movement; therefore, the poroelastic behavior might not be separated by changing the size of the gel. On the other hand, if the solvent migration is not coupled with the dashpot movement, the poroelastic deformation and energy dissipation could be separated from the viscoelastic dissipation for both large-size gel and small-size gels as long as the poroelastic time and viscoelastic time are different. Furthermore, different rheological models describe the physical characteristics of the gel in different ways. The Maxwell model and the Kelvin–Voigt model give an intrinsic fluid-like behavior in, respectively, a long time scale (low-frequency region) and a short time scale (high-frequency region), while the generalized standard viscoelastic models (both GSV-dashpot and GSV-spring models) guarantee the elastic behavior of the gel. The theory in this work could be generally applied to gels with different time-dependent characteristics and prospectively inspire future designs of experiments for characterizing gel properties or unraveling the time-dependent and frequency-dependent mechanism of the various practical gel systems.

Acknowledgment

The materials are based upon work supported by the National Science Foundation (NSF) (Grant No. 1554326; Funder ID: 10.13039/501100001809). Y. H. also acknowledges the funding support from Air Force Office of Scientific Research (AFOSR) (Award No. FA9550-19-1-0395; Funder ID: 10.13039/100000181) (Dr. B.-L. “Les” Lee, Program Manager).

References

1.
Hoare
,
T. R.
, and
Kohane
,
D. S.
,
2008
, “
Hydrogels in Drug Delivery: Progress and Challenges
,”
Polymer
,
49
(
8
), pp.
1993
2007
. 10.1016/j.polymer.2008.01.027
2.
Qiu
,
Y.
, and
Park
,
K.
,
2001
, “
Environment-Sensitive Hydrogels for Drug Delivery
,”
Adv. Drug Deliv. Rev.
,
53
(
3
), pp.
321
339
. 10.1016/S0169-409X(01)00203-4
3.
Li
,
J.
, and
Mooney
,
D. J.
,
2016
, “
Designing Hydrogels for Controlled Drug Delivery
,”
Nat. Rev. Mater.
,
1
(
12
), p.
16071
. 10.1038/natrevmats.2016.71
4.
Jen
,
A. C.
,
Wake
,
M. C.
, and
Mikos
,
A. G.
,
1996
, “
Hydrogels for Cell Immobilization
,”
Biotechnol. Bioeng.
,
50
(
4
), pp.
357
364
. 10.1002/(SICI)1097-0290(19960520)50:4<357::AID-BIT2>3.0.CO;2-K
5.
Lee
,
K. Y.
, and
Mooney
,
D. J.
,
2001
, “
Hydrogels for Tissue Engineering
,”
Chem. Rev.
,
101
(
7
), pp.
1869
1880
. 10.1021/cr000108x
6.
Hoffman
,
A. S.
,
2002
, “
Hydrogels for Biomedical Applications
,”
Adv. Drug Deliv. Rev.
,
54
(
1
), pp.
3
12
. 10.1016/S0169-409X(01)00239-3
7.
Aguado
,
B. A.
,
Grim
,
J. C.
,
Rosales
,
A. M.
,
Watson-Capps
,
J. J.
, and
Anseth
,
K. S.
,
2018
, “
Engineering Precision Biomaterials for Personalized Medicine
,”
Sci. Transl. Med.
,
10
(
424
), p.
eaam8645
. 10.1126/scitranslmed.aam8645
8.
Grigoryan
,
B.
,
Paulsen
,
S. J.
,
Corbett
,
D. C.
,
Sazer
,
D. W.
,
Fortin
,
C. L.
,
Zaita
,
A. J.
,
Greenfield
,
P. T.
,
Calafat
,
N. J.
,
Gounley
,
J. P.
,
Ta
,
A. H.
,
Johansson
,
F.
,
Randles
,
A.
,
Rosenkrantz
,
J. E.
,
Louis-Rosenberg
,
J. D.
,
Galie
,
P. A.
,
Stevens
,
K. R.
, and
Miller
,
J. S.
,
2019
, “
Multivascular Networks and Functional Intravascular Topologies Within Biocompatible Hydrogels
,”
Science
,
364
(
6439
), pp.
458
464
. 10.1126/science.aav9750
9.
Lee
,
A.
,
Hudson
,
A. R.
,
Shiwarski
,
D. J.
,
Tashman
,
J. W.
,
Hinton
,
T. J.
,
Yerneni
,
S.
,
Bliley
,
J. M.
,
Campbell
,
P. G.
, and
Feinberg
,
A. W.
,
2019
, “
3D Bioprinting of Collagen to Rebuild Components of the Human Heart
,”
Science
,
365
(
6452
), pp.
482
487
. 10.1126/science.aav9051
10.
Prokop
,
A. F.
,
Vaezy
,
S.
,
Noble
,
M. L.
,
Kaczkowski
,
P. J.
,
Martin
,
R. W.
, and
Crum
,
L. A.
,
2003
, “
Polyacrylamide Gel as an Acoustic Coupling Medium for Focused Ultrasound Therapy
,”
Ultrasound Med. Biol.
,
29
(
9
), pp.
1351
1358
. 10.1016/S0301-5629(03)00979-7
11.
Lafon
,
C.
,
Zderic
,
V.
,
Noble
,
M. L.
,
Yuen
,
J. C.
,
Kaczkowski
,
P. J.
,
Sapozhnikov
,
O. A.
,
Chavrier
,
F.
,
Crum
,
L. A.
, and
Vaezy
,
S.
,
2005
, “
Gel Phantom for Use in High-Intensity Focused Ultrasound Dosimetry
,”
Ultrasound Med. Biol.
,
31
(
10
), pp.
1383
1389
. 10.1016/j.ultrasmedbio.2005.06.004
12.
Forte
,
A. E.
,
Galvan
,
S.
,
Manieri
,
F.
,
Rodriguez y Baena
,
F.
, and
Dini
,
D.
,
2016
, “
A Composite Hydrogel for Brain Tissue Phantoms
,”
Mater. Des.
,
112
, pp.
227
238
. 10.1016/j.matdes.2016.09.063
13.
Jiang
,
S.
,
Liu
,
S.
, and
Feng
,
W.
,
2011
, “
PVA Hydrogel Properties for Biomedical Application
,”
J. Mech. Behav. Biomed. Mater.
,
4
(
7
), pp.
1228
1233
. 10.1016/j.jmbbm.2011.04.005
14.
Tanaka
,
T.
,
1978
, “
Collapse of Gels and the Critical Endpoint
,”
Phys. Rev. Lett.
,
40
(
12
), pp.
820
823
. 10.1103/PhysRevLett.40.820
15.
Tanaka
,
T.
,
Fillmore
,
D.
,
Sun
,
S.-T.
,
Nishio
,
I.
,
Swislow
,
G.
, and
Shah
,
A.
,
1980
, “
Phase Transitions in Ionic Gels
,”
Phys. Rev. Lett.
,
45
(
20
), pp.
1636
1639
. 10.1103/PhysRevLett.45.1636
16.
Tanaka
,
T.
,
Nishio
,
I.
,
Sun
,
S.-T.
, and
Ueno-Nishio
,
S.
,
1982
, “
Collapse of Gels in an Electric Field
,”
Science
,
218
(
4571
), pp.
467
469
. 10.1126/science.218.4571.467
17.
Suzuki
,
A.
, and
Tanaka
,
T.
,
1990
, “
Phase Transition in Polymer Gels Induced by Visible Light
,”
Nature
,
346
(
6282
), pp.
345
347
. 10.1038/346345a0
18.
Richter
,
A.
,
Paschew
,
G.
,
Klatt
,
S.
,
Lienig
,
J.
,
Arndt
,
K.-F.
, and
Adler
,
H.-J. P.
,
2008
, “
Review on Hydrogel-Based PH Sensors and Microsensors
,”
Sensors
,
8
(
1
), pp.
561
581
. 10.3390/s8010561
19.
Buenger
,
D.
,
Topuz
,
F.
, and
Groll
,
J.
,
2012
, “
Hydrogels in Sensing Applications
,”
Prog. Polym. Sci.
,
37
(
12
), pp.
1678
1719
. 10.1016/j.progpolymsci.2012.09.001
20.
Gerlach
,
G.
,
Guenther
,
M.
,
Sorber
,
J.
,
Suchaneck
,
G.
,
Arndt
,
K.-F.
, and
Richter
,
A.
,
2005
, “
Chemical and PH Sensors Based on the Swelling Behavior of Hydrogels
,”
Sens. Actuators B
,
111–112
, pp.
555
561
. 10.1016/j.snb.2005.03.040
21.
Bassil
,
M.
,
Davenas
,
J.
, and
EL Tahchi
,
M.
,
2008
, “
Electrochemical Properties and Actuation Mechanisms of Polyacrylamide Hydrogel for Artificial Muscle Application
,”
Sens. Actuators B
,
134
(
2
), pp.
496
501
. 10.1016/j.snb.2008.05.025
22.
Maeda
,
S.
,
Hara
,
Y.
,
Sakai
,
T.
,
Yoshida
,
R.
, and
Hashimoto
,
S.
,
2007
, “
Self-Walking Gel
,”
Adv. Mater.
,
19
(
21
), pp.
3480
3484
. 10.1002/adma.200700625
23.
Morales
,
D.
,
Palleau
,
E.
,
Dickey
,
M. D.
, and
Velev
,
O. D.
,
2014
, “
Electro-Actuated Hydrogel Walkers With Dual Responsive Legs
,”
Soft Matter
,
10
(
9
), pp.
1337
1348
. 10.1039/C3SM51921J
24.
Kim
,
C.-C.
,
Lee
,
H.-H.
,
Oh
,
K. H.
, and
Sun
,
J.-Y.
,
2016
, “
Highly Stretchable, Transparent Ionic Touch Panel
,”
Science
,
353
(
6300
), pp.
682
687
. 10.1126/science.aaf8810
25.
Hu
,
Y.
,
Zhao
,
X.
,
Vlassak
,
J. J.
, and
Suo
,
Z.
,
2010
, “
Using Indentation to Characterize the Poroelasticity of Gels
,”
Appl. Phys. Lett.
,
96
(
12
), p.
121904
. 10.1063/1.3370354
26.
Biot
,
M. A.
,
1954
, “
Theory of Stress-Strain Relations in Anisotropic Viscoelasticity and Relaxation Phenomena
,”
J. Appl. Phys.
,
25
(
11
), pp.
1385
1391
. 10.1063/1.1721573
27.
Biot
,
M. A.
,
1955
, “
Theory of Elasticity and Consolidation for a Porous Anisotropic Solid
,”
J. Appl. Phys.
,
26
(
2
), pp.
182
185
. 10.1063/1.1721956
28.
Biot
,
M. A.
,
1956
, “
Theory of Deformation of a Porous Viscoelastic Anisotropic Solid
,”
J. Appl. Phys.
,
27
(
5
), pp.
459
467
. 10.1063/1.1722402
29.
Rice
,
J. R.
, and
Cleary
,
M. P.
,
1976
, “
Some Basic Stress Diffusion Solutions for Fluid-Saturated Elastic Porous Media With Compressible Constituents
,”
Rev. Geophys.
,
14
(
2
), pp.
227
241
. 10.1029/RG014i002p00227
30.
Cleary
,
M. P.
,
1978
, “
Elastic and Dynamic Response Regimes of Fluid-Impregnated Solids With Diverse Microstructures
,”
Int. J. Solids Struct.
,
14
(
10
), pp.
795
819
. 10.1016/0020-7683(78)90072-0
31.
McTigue
,
D. F.
,
1986
, “
Thermoelastic Response of Fluid-Saturated Porous Rock
,”
J. Geophys. Res. Solid Earth
,
91
(
B9
), pp.
9533
9542
. 10.1029/JB091iB09p09533
32.
Coussy
,
O.
,
Dormieux
,
L.
, and
Detournay
,
E.
,
1998
, “
From Mixture Theory to Biot’s Approach for Porous Media
,”
Int. J. Solids Struct.
,
35
(
34–35
), pp.
4619
4635
. 10.1016/S0020-7683(98)00087-0
33.
Mow
,
V. C.
,
Kuei
,
S. C.
,
Lai
,
W. M.
, and
Armstrong
,
C. G.
,
1980
, “
Biphasic Creep and Stress Relaxation of Articular Cartilage in Compression: Theory and Experiments
,”
ASME J. Biomech. Eng.
,
102
(
1
), pp.
73
84
. 10.1115/1.3138202
34.
Mak
,
A. F.
,
Lai
,
W. M.
, and
Mow
,
V. C.
,
1987
, “
Biphasic Indentation of Articular Cartilage—I. Theoretical Analysis
,”
J. Biomech.
,
20
(
7
), pp.
703
714
. 10.1016/0021-9290(87)90036-4
35.
Hong
,
W.
,
Zhao
,
X.
,
Zhou
,
J.
, and
Suo
,
Z.
,
2008
, “
A Theory of Coupled Diffusion and Large Deformation in Polymeric Gels
,”
J. Mech. Phys. Solids
,
56
(
5
), pp.
1779
1793
. 10.1016/j.jmps.2007.11.010
36.
Hong
,
W.
,
Liu
,
Z.
, and
Suo
,
Z.
,
2009
, “
Inhomogeneous Swelling of a Gel in Equilibrium With a Solvent and Mechanical Load
,”
Int. J. Solids Struct.
,
46
(
17
), pp.
3282
3289
. 10.1016/j.ijsolstr.2009.04.022
37.
Doi
,
M.
,
2009
, “
Gel Dynamics
,”
J. Phys. Soc. Jpn.
,
78
(
5
), p.
52001
. 10.1143/JPSJ.78.052001
38.
Bouklas
,
N.
, and
Huang
,
R.
,
2012
, “
Swelling Kinetics of Polymer Gels: Comparison of Linear and Nonlinear Theories
,”
Soft Matter
,
8
(
31
), pp.
8194
8203
. 10.1039/c2sm25467k
39.
Bouklas
,
N.
,
Landis
,
C. M.
, and
Huang
,
R.
,
2015
, “
A Nonlinear, Transient Finite Element Method for Coupled Solvent Diffusion and Large Deformation of Hydrogels
,”
J. Mech. Phys. Solids
,
79
, pp.
21
43
. 10.1016/j.jmps.2015.03.004
40.
Cai
,
S.
,
Hu
,
Y.
,
Zhao
,
X.
, and
Suo
,
Z.
,
2010
, “
Poroelasticity of a Covalently Crosslinked Alginate Hydrogel Under Compression
,”
J. Appl. Phys.
,
108
(
11
), p.
113514
. 10.1063/1.3517146
41.
Chester
,
S. A.
, and
Anand
,
L.
,
2010
, “
A Coupled Theory of Fluid Permeation and Large Deformations for Elastomeric Materials
,”
J. Mech. Phys. Solids
,
58
(
11
), pp.
1879
1906
. 10.1016/j.jmps.2010.07.020
42.
Chester
,
S. A.
,
Di Leo
,
C. V.
, and
Anand
,
L.
,
2015
, “
A Finite Element Implementation of a Coupled Diffusion-Deformation Theory for Elastomeric Gels
,”
Int. J. Solids Struct.
,
52
, pp.
1
18
. 10.1016/j.ijsolstr.2014.08.015
43.
Bosnjak
,
N. S.
,
Shuolun
,
W.
, and
Chester
,
S. A.
,
2020
, “
Modeling Deformation-Diffusion in Polymeric Gels
,”
Poromechanics VI
,
Paris, France
,
July 9–13
, pp.
141
148
.
44.
Liu
,
Y.
,
Zhang
,
H.
,
Zhang
,
J.
, and
Zheng
,
Y.
,
2015
, “
Constitutive Modeling for Polymer Hydrogels: A New Perspective and Applications to Anisotropic Hydrogels in Free Swelling
,”
Eur. J. Mech. A. Solids
,
54
, pp.
171
186
. 10.1016/j.euromechsol.2015.07.001
45.
Bosnjak
,
N.
,
Wang
,
S.
,
Han
,
D.
,
Lee
,
H.
, and
Chester
,
S. A.
,
2019
, “
Modeling of Fiber-Reinforced Polymeric Gels
,”
Mech. Res. Commun.
,
96
, pp.
7
18
. 10.1016/j.mechrescom.2019.02.002
46.
Lucantonio
,
A.
,
Nardinocchi
,
P.
, and
Teresi
,
L.
,
2013
, “
Transient Analysis of Swelling-Induced Large Deformations in Polymer Gels
,”
J. Mech. Phys. Solids
,
61
(
1
), pp.
205
218
. 10.1016/j.jmps.2012.07.010
47.
Hu
,
Y.
, and
Suo
,
Z.
,
2012
, “
Viscoelasticity and Poroelasticity in Elastomeric Gels
,”
Acta Mech. Solida Sin.
,
25
(
5
), pp.
441
458
. 10.1016/S0894-9166(12)60039-1
48.
Wang
,
X.
, and
Hong
,
W.
,
2012
, “
A Visco-Poroelastic Theory for Polymeric Gels
,”
Proc. R. Soc. A
,
468
(
2148
), pp.
3824
3841
. 10.1098/rspa.2012.0385
49.
Caccavo
,
D.
, and
Lamberti
,
G.
,
2017
, “
PoroViscoElastic Model to Describe Hydrogels’ Behavior
,”
Mater. Sci. Eng. C
,
76
, pp.
102
113
. 10.1016/j.msec.2017.02.155
50.
Chester
,
S. A.
,
2012
, “
A Constitutive Model for Coupled Fluid Permeation and Large Viscoelastic Deformation in Polymeric Gels
,”
Soft Matter
,
8
(
31
), pp.
8223
8233
. 10.1039/c2sm25372k
51.
Pioletti
,
D. P.
,
Rakotomanana
,
L. R.
,
Benvenuti
,
J.-F.
, and
Leyvraz
,
P.-F.
,
1998
, “
Viscoelastic Constitutive Law in Large Deformations: Application to Human Knee Ligaments and Tendons
,”
J. Biomech.
,
31
(
8
), pp.
753
757
. 10.1016/S0021-9290(98)00077-3
52.
Feynman
,
R. P.
,
Leighton
,
R. B.
, and
Sands
,
M.
,
1965
, “
The Feynman Lectures on Physics; Vol. I
,”
Am. J. Phys.
,
33
(
9
), pp.
750
752
. 10.1119/1.1972241
53.
Flory
,
P. J.
, and
Rehner
,
J.
,
1943
, “
Statistical Mechanics of Crosss; Vol. I Deformations: Application to Human Knee
,
J. Chem. Phys.
,
11
(
11
), pp.
512
520
. 10.1063/1.1723791
54.
Holzapfel
,
G. A.
,
2000
,
Nonlinear Solid Mechanics a Continuum Approach for Engineering
,
John Wiley & Sons Inc.
,
Chichester
,
Chap. 6
.
55.
Huggins
,
M. L.
,
1941
, “
Solutions of Long Chain Compounds
,”
J. Chem. Phys.
,
9
(
5
), p.
440
. 10.1063/1.1750930
56.
Flory
,
P. J.
,
1953
,
Principles of Polymer Chemistry
,
Cornell University Press
,
Ithaca, New York
.
57.
Joanny
,
J. F.
,
Jülicher
,
F.
,
Kruse
,
K.
, and
Prost
,
J.
,
2007
, “
Hydrodynamic Theory for Multi-Component Active Polar Gels
,”
New J. Phys.
,
9
(
11
), p.
422
. 10.1088/1367-2630/9/11/422