## 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 [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.

## 2 Nonlinear Visco-Poroelasticity

*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:

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 $\xi =\xi (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.

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

*δ*represents the change of the value during a small-time interval.

*δF*

_{iK},

*δ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:

*δξ*

_{ij}. Therefore, we assume that the evolution of the internal variable is governed by

*R*

_{ikjl}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]:

*η*is the viscosity and

*δ*

_{kl}is the Kronecker delta.

### 2.2 Inhomogeneous Fields.

*Q*is the number of solvent molecules that are injected into the gel per unit area of the surface in unit time, and $\mu \xaf$ is the chemical potential of the solvent in the environment. The integration represents the combination of all material elements in the domain.

*N*

_{K}be the unit normal of the boundary of the body in the reference configuration. We add a term of $\u222b\mu 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:

*δx*

_{i}/

*δt*,

*δξ*

_{ij}/

*δt*,

*μ*, and

*J*

_{K}, 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

*J*

_{K}= 0, the thermodynamic system undergoes the reversible process, and we have

*J*

_{K}:

*M*

_{KL}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]:

*D*is the diffusivity,

*k*is the Boltzmann constant,

*T*is the temperature, and $H$ is the inverse of the deformation gradient tensor, $H=F\u22121$.

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

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.

*χ*is the Flory–Huggins parameter that characterizes the mixing enthalpy. A smaller

*χ*represents a higher affinity between the polymer and the solvent.

### 3.2 Kelvin–Voigt Gel.

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

*C*, that is, $\Psi =\Psi (F,C)$. The inequality Eq. (4) for the KV gel becomes

*I*

_{1}=

*F*

_{iK}

*F*

_{iK}; $J=detF$;

*J*

_{e}=

*J*/

*J*

_{s};

*G*and

*κ*are the shear modulus and Lamé constant of the spring, respectively; and

*χ*is the Flory–Huggins parameter.

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 $\xi $ 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 $\xi $ 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.

*I*

_{1}=

*F*

_{iK}

*F*

_{iK}; $J=detF$; $I~1=F~iKF~iK$; $J~=detF~$; $Js=detFs=1+\Omega C$;

*J*

_{e}=

*J*/

*J*

_{s};

*G*and

*κ*are the shear modulus and Lamé constant of the primary spring, respectively; and $G~$ and $\kappa ~$ are the shear modulus and Lamé constant of the secondary spring, respectively; and

*χ*is the Flory–Huggins parameter.

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

*I*

_{1}=

*F*

_{iK}

*F*

_{iK}; $J=detF$; $I~1=F~iKF~iK$; $J~=detF~$; $Js=detFs=1+\Omega C$;

*J*

_{e}=

*J*/

*J*

_{s}; $J~e=detF~e$;

*G*and

*κ*are the shear modulus and Lamé constant of the primary spring, respectively; and $G~$ and $\kappa ~$ are the shear modulus and Lamé constant of the secondary spring, respectively; and

*χ*is the Flory–Huggins parameter.

## 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\tau $.

Parameter/variable | Notation | SI base unit | Normalization |
---|---|---|---|

Coordinate | $X,x$ | m | $x^,X^=x,X/L0$ |

Time | t | s | $t^=t/\tau $ |

Nominal stress | $s$ | N/m^{2} | $s^=s/(kT/\Omega )$ |

Moduli | $G,\kappa ,G~,\kappa ~$ | N/m^{2} | $G^,\kappa ^,G~^,\kappa ~^=G,\kappa ,G~,\kappa ~/(kT/\Omega )$ |

Concentration | C | 1/m^{2} | $C^=\Omega C$ |

Chemical potential | μ | N m | $\mu ^=\mu /kT$ |

Flux | $J$ | 1/(s m^{2}) | $J^=J/(L0/\Omega \tau )$ |

Frequency | ω | rad/s | $\omega ^=\omega \tau $ |

Parameter/variable | Notation | SI base unit | Normalization |
---|---|---|---|

Coordinate | $X,x$ | m | $x^,X^=x,X/L0$ |

Time | t | s | $t^=t/\tau $ |

Nominal stress | $s$ | N/m^{2} | $s^=s/(kT/\Omega )$ |

Moduli | $G,\kappa ,G~,\kappa ~$ | N/m^{2} | $G^,\kappa ^,G~^,\kappa ~^=G,\kappa ,G~,\kappa ~/(kT/\Omega )$ |

Concentration | C | 1/m^{2} | $C^=\Omega C$ |

Chemical potential | μ | N m | $\mu ^=\mu /kT$ |

Flux | $J$ | 1/(s m^{2}) | $J^=J/(L0/\Omega \tau )$ |

Frequency | ω | rad/s | $\omega ^=\omega \tau $ |

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

MW gel:

Kinetics of solvent transportation:(47)$J^K=\u2212C^HKiHLi\u2202\mu ^\u2202X^L$Constitutive law for nominal stress:(48)$s^iK=G~^(F~jK\zeta ji\u2212HKi)+\kappa ~^2J^s(J~^e2\u22121)HKi$Constitutive law for chemical potential:(49)$\mu ^=\u2212\kappa ~^4(J~^e2\u22121+2lnJ~^e)+ln(C^1+C^)+\chi +1+C^(1+C^)2$Evolution equation of the internal variable:(50)$(\xi il\xi jk+\xi im\xi jm\delta kl)\u2202\xi jl\u2202t^=2G~^(F~jNF~iN\zeta jk\u2212\zeta ki)+\kappa ~^J^s(J~^e2\u22121)\zeta ki(50)$

- (2)
KV gel:

Kinetics of solvent transportation:(51)$J^K=\u2212C^HKiHLi\u2202\mu ^\u2202X^L$Constitutive law for nominal stress:(52)$s^iK=G^(FiK\u2212HKi)+\kappa ^2J^s(J^e2\u22121)HKi+12(FiLFjK+FiMFjM\delta KL)\u2202FjL\u2202t^(52)$Constitutive law for chemical potential:(53)$\mu ^=\u2212\kappa ^4(J^e2\u22121+2lnJ^e)+ln(C^1+C^)+\chi +1+C^(1+C^)2$

- (3)
GSV-dashpot gel:

Kinetics of solvent transportation:(54)$J^K=\u2212C^HKiHLi\u2202\mu ^\u2202X^L$Constitutive law for nominal stress:(55)$s^iK=G^(FiK\u2212HKi)+\kappa ^2J^s(J^e2\u22121)HKi+G~^(F~jK\zeta ji\u2212HKi)+\kappa ~^2(J~2\u22121)HKi(55)$Constitutive law for chemical potential:(56)$\mu ^=\u2212\kappa ^4(J^e2\u22121+2lnJ^e)+ln(C^1+C^)+\chi +1+C^(1+C^)2$Evolution equation of the internal variable:(57)$(\xi il\xi jk+\xi im\xi jm\delta kl)\u2202\xi jl\u2202t^=2G~^(F~jNF~iN\zeta jk\u2212\zeta ki)+\kappa ~^(J~2\u22121)\zeta ki$

- (4)
GSV-spring gel:

Kinetics of solvent transportation:(58)$J^K=\u2212C^HKiHLi\u2202\mu ^\u2202X^L$Constitutive law for nominal stress:(59)$s^iK=G^(FiK\u2212HKi)+\kappa ^2J^s(J^e2\u22121)HKi+G~^(F~jK\zeta ji\u2212HKi)+\kappa ~^2J^s(J~^e2\u22121)HKi$Constitutive law for chemical potential:(60)$\mu ^=\u2212\kappa ^4(J^e2\u22121+2lnJ^e)\u2212\kappa ~^4(J~^e2\u22121+2lnJ~^e)+ln(C^1+C^)+\chi +1+C^(1+C^)2(60)$Evolution equation of the internal valuable:(61)$(\xi il\xi jk+\xi im\xi jm\delta kl)\u2202\xi jl\u2202t^=2G~^(F~jNF~iN\zeta jk\u2212\zeta ki)+\kappa ~^J^s(J~^e2\u22121)\zeta ki(61)$

## 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^\xd7L^$. As discussed in Sec. 4, all the length scales in the formula are normalized by an intrinsic length scale of the material $L0=D\tau $. 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*_{2} 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*_{1} boundary. In addition, the *S*_{3} boundary is set to be stress free and the zero vertical displacement constraint is applied on the *S*_{4} 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.

Rheological model | $G^,\kappa ^$ | $G~^,\kappa ~^$ | χ | λ_{0} | C_{0} |
---|---|---|---|---|---|

MW | – | $G~^=0.01,\kappa ~^=10$ | 0.1 | 2.222 | 9.977 |

KV | $G^=0.01,\kappa ^=10$ | – | 0.1 | 2.222 | 9.977 |

GSV-dashpot | $G^=0.01,\kappa ^=10$ | $G~^=0.1,\kappa ~^=100$ | 0.1 | 2.222 | 9.977 |

GSV-spring | $G^=0.005,\kappa ^=5$ | $G~^=0.005,\kappa ~^=5$ | 0.1 | 2.222 | 9.977 |

Rheological model | $G^,\kappa ^$ | $G~^,\kappa ~^$ | χ | λ_{0} | C_{0} |
---|---|---|---|---|---|

MW | – | $G~^=0.01,\kappa ~^=10$ | 0.1 | 2.222 | 9.977 |

KV | $G^=0.01,\kappa ^=10$ | – | 0.1 | 2.222 | 9.977 |

GSV-dashpot | $G^=0.01,\kappa ^=10$ | $G~^=0.1,\kappa ~^=100$ | 0.1 | 2.222 | 9.977 |

GSV-spring | $G^=0.005,\kappa ^=5$ | $G~^=0.005,\kappa ~^=5$ | 0.1 | 2.222 | 9.977 |

### 5.1 Uniaxial Stretching and Creep.

At the time $t^=0+$, nominal stress $s^=0.005$ is applied to the *S*_{2} 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*_{2} direction at the *X*_{1} = 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:

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~=\lambda 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*X*_{1}= 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/\tau $, 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].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.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^=\u221e)/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.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.

### 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\omega ^t^$ is applied to the *S*_{2} 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(\omega ^t^\u2212\delta )$. 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 $\omega ^$ for the different visco-poroelastic models proposed in this work.

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 $\omega ^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.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 $\omega ^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.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 $\delta \u223c\omega ^$ plot are the viscoelastic peaks, while those that collapse in the $\delta \u223c\omega ^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.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 $\delta \u223c\omega ^$ plot are the viscoelastic peaks, while those that collapse in the $\delta \u223c\omega ^L^2$ are the poroelastic peaks (Fig. 13(b)).

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