A comparative analysis of transient finite-strain coupled diffusion-deformation theories for hydrogels

TL;DR

The paper tackles diffusion–deformation coupling in hydrogels under large strains by unifying several thermodynamically consistent poroelastic theories. It develops a common variational framework, distinguishing compressible and incompressible formulations, and implements monolithic and staggered FEM schemes in FEniCS using Taylor–Hood elements to ensure inf-sup stability. Through 1–3D prototype problems and a reference benchmark, it shows that differences across models primarily arise from the volumetric (energetic) response and the enforcement of incompressibility, with numerical performance (Newton convergence, discretization order) closely linked to the chosen formulation. The work highlights the need for careful model selection and parameter calibration against experiments, and provides open-source tools to enable reproducible comparisons and future extensions to more complex stimuli or reactions. These insights advance the reliable simulation of hydrogel diffusion–deformation and support design of hydrogels with tailored swelling and mechanical responses.

Abstract

This work presents a comparative review and classification between some well-known thermodynamically consistent models of hydrogel behavior in a large deformation setting, specifically focusing on solvent absorption/desorption and its impact on mechanical deformation and network swelling. The proposed discussion addresses formulation aspects, general mathematical classification of the governing equations, and numerical implementation issues based on the finite element method. The theories are presented in a unified framework demonstrating that, despite not being evident in some cases, all of them follow equivalent thermodynamic arguments. A detailed numerical analysis is carried out where Taylor-Hood elements are employed in the spatial discretization to satisfy the inf-sup condition and to prevent spurious numerical oscillations. The resulting discrete problems are solved using the FEniCS platform through consistent variational formulations, employing both monolithic and staggered approaches. We conduct benchmark tests on various hydrogel structures, demonstrating that major differences arise from the chosen volumetric response of the hydrogel. The significance of this choice is frequently underestimated in the state-of-the-art literature but has been shown to have substantial implications on the resulting hydrogel behavior.

Figures (15)

• Figure 1: Representative examples setup. a. One-dimensional transient swelling of a hydrogel bar. The bar is fixed at $Y = 0.0~[m]$, while the opposite end, $Y = 0.01~[m]$, is exposed to a non-reactive solvent. b. Two-dimensional transient swelling of a constrained hydrogel slab. In this example, the hydrogel block is placed in a rigid container with frictionless walls and the deformation in the $X$ direction is constrained. The top surface at $Y = 0.01~[m]$ keeps traction-free and is in contact with the solvent during deformation. At the bottom surface $Y = 0.0~[m]$, the gel is fixed to the container wall and no fluid is allowed to diffuse through it. Due to the solvent absorption, the hydrogel can only swell along the $Y$ direction. c. Two-dimensional hydrogel block is immersed in a non-reactive solvent with a reference chemical potential $\mu^0 = 0$. Only a quarter of the whole model is considered because of the symmetry of the block. For the mechanical boundary conditions, the nodes along edge ab are prescribed to have displacement component $u_y = 0$, while the nodes along edge ad are prescribed to have $u_x = 0$. The edges bc and cd are taken to be traction-free. For the solvent concentration boundary conditions, the edges ab and ad (the symmetry edges) are prescribed a zero fluid flux, and on the edges bc and cd, the chemical potential is prescribed as $\mu = 0$ on $\partial \mathcal{B}_{\mu},~t=\lbrace 0, T \rbrace$. d. Three-dimensional cube immersed in a non-reactive solvent. Only a quarter of the whole model is considered because of the symmetry of the 3D cube. The mechanical boundary conditions are specified such that the $u_y = 0$ in the front face, $u_x = 0$ in the left face, and $u_z = 0$ in the face in the bottom part. For the solvent concentration boundary conditions, the front, left, and bottom faces (the symmetry faces) are prescribed a zero fluid flux, and on the back, right, and top faces, the chemical potential is prescribed as $\mu = 0$ on $\partial \mathcal{B}_{\mu},~t=\lbrace 0, T \rbrace$. Note: the remaining boundary conditions, together with the initial conditions, are defined depending on the specific constitutive theory adopted to study the diffusion-deformation process.
• Figure 2: Constitutive model I:one-dimensional bar (black dots) and two-dimensional hydrogel constrained slab (colored lines) numerical solution comparison for different mesh densities $N_h$ and at different simulation times.a. Deformed two-dimensional constrained slab at $t = 10.0~[s]$. b. Stretch due to swelling $\lambda$. c. Chemical potential $\mu$ normalized by $k_B T$. d. Cauchy compressive stress $\bm{\sigma}_X$. Simulation parameters:$G_0 = 10~[MPa], ~~ \chi = 0.2~[--], ~~ D = 2.0\times10^{-5}~[m^2 s^{-1}]$.
• Figure 3: Constitutive model I:newton iterations along the time steps.a. For different mesh densities ($N_h$). b. For different time step sizes ($N_k$).
• Figure 4: Constitutive model I:convergence analysis of the Newton solver for different mesh densities ($N_h$).a. $N_h = 25$. b. $N_h = 50$. c. $N_h = 100$. d. $N_h = 200$.
• Figure 5: Constitutive model I:convergence analysis of the Newton solver for different time steps size ($N_k$).a. $N_k = 25$. b. $N_k = 50$. c. $N_k = 100$. d. $N_k = 200$.