Nonlinear Visco-Poroelasticity of Gels With Different Rheological Parts

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 [1–3], cell culture scaffold [4–6], 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 [11–13]. 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.

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 [26–28]. The theory has been widely applied to study the phenomena of geomaterials [29–32]. 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 [35–46]. 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 [47–50]. 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.

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.

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.

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.

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

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.

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.

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.

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.

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τ$⁠.

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

1. MW gel:

   Kinetics of solvent transportation:

   $J^K=−C^HKiHLi∂μ^∂X^L$

   (47)
   Constitutive law for nominal stress:

   $s^iK=G~^(F~jKζji−HKi)+κ~^2J^s(J~^e2−1)HKi$

   (48)
   Constitutive law for chemical potential:

   $μ^=−κ~^4(J~^e2−1+2lnJ~^e)+ln(C^1+C^)+χ+1+C^(1+C^)2$

   (49)
   Evolution equation of the internal variable:

   $(ξilξjk+ξimξjmδkl)∂ξjl∂t^=2G~^(F~jNF~iNζjk−ζki)+κ~^J^s(J~^e2−1)ζki(50)$

   (50)

• (2)

  KV gel:

  Kinetics of solvent transportation:

  $J^K=−C^HKiHLi∂μ^∂X^L$

  (51)
  Constitutive law for nominal stress:

  $s^iK=G^(FiK−HKi)+κ^2J^s(J^e2−1)HKi+12(FiLFjK+FiMFjMδKL)∂FjL∂t^(52)$

  (52)
  Constitutive law for chemical potential:

  $μ^=−κ^4(J^e2−1+2lnJ^e)+ln(C^1+C^)+χ+1+C^(1+C^)2$

  (53)

• (3)

  GSV-dashpot gel:

  Kinetics of solvent transportation:

  $J^K=−C^HKiHLi∂μ^∂X^L$

  (54)
  Constitutive law for nominal stress:

  $s^iK=G^(FiK−HKi)+κ^2J^s(J^e2−1)HKi+G~^(F~jKζji−HKi)+κ~^2(J~2−1)HKi(55)$

  (55)
  Constitutive law for chemical potential:

  $μ^=−κ^4(J^e2−1+2lnJ^e)+ln(C^1+C^)+χ+1+C^(1+C^)2$

  (56)
  Evolution equation of the internal variable:

  $(ξilξjk+ξimξjmδkl)∂ξjl∂t^=2G~^(F~jNF~iNζjk−ζki)+κ~^(J~2−1)ζki$

  (57)

• (4)

  GSV-spring gel:

  Kinetics of solvent transportation:

  $J^K=−C^HKiHLi∂μ^∂X^L$

  (58)
  Constitutive law for nominal stress:

  $s^iK=G^(FiK−HKi)+κ^2J^s(J^e2−1)HKi+G~^(F~jKζji−HKi)+κ~^2J^s(J~^e2−1)HKi$

  (59)
  Constitutive law for chemical potential:

  $μ^=−κ^4(J^e2−1+2lnJ^e)−κ~^4(J~^e2−1+2lnJ~^e)+ln(C^1+C^)+χ+1+C^(1+C^)2(60)$

  (60)
  Evolution equation of the internal valuable:

  $(ξilξjk+ξimξjmδkl)∂ξjl∂t^=2G~^(F~jNF~iNζjk−ζki)+κ~^J^s(J~^e2−1)ζki(61)$

  (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.

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 X₂ 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 S₁ boundary. In addition, the S₃ boundary is set to be stress free and the zero vertical displacement constraint is applied on the S₄ boundary.

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.

At the time $t^=0+$⁠, nominal stress $s^=0.005$ is applied to the S₂ 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 X₂ direction at the X₁ = 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 λ₁ = 2.096 and λ₂ = 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 X₁ = 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 λ₀ 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: λ₁ = 2.180 and λ₂ = 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 λ₁ = 2.096 and λ₂ = 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.