[1,2]Jorge-Humberto Urrea-Quintero

[3]Michele Marino

[4]Thomas Wick

[1]Udo Nackenhorst

1]Leibniz University Hannover, IBNM - Institute of Mechanics and Computational
Mechanics, Appelstraße 9a, Hannover, 30167, Lower Saxony, Germany

2]Technische Universtität Braunschweig, iRMB - Institute for Computational
Modeling in Civil Engineering, Pockelsstr. 3, Braunschweig, 38106, Lower
Saxony, Germany

3]Department of Civil Engineering and Computer Science Engineering, University
of Rome Tor Vergata, Via del Politecnico 1, Rome, 00133, Italy

4]Leibniz University Hannover, Institut für Angewandte Mathematik,
Welfengarten 1, Hannover, 30167, Lower Saxony, Germany

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

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [

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

###### keywords:

hydrogels, finite element method, thermodynamic consistency, mathematical
classification, large deformations, nonlinear coupled problems, FEniCS.

##  1 Introduction

Hydrogels, a class of versatile soft materials with the unique ability to
absorb and retain fluid within their three-dimensional network structures,
have gained widespread attention in recent years across various industrial
applications. Their diverse uses include serving as drug carriers in
biomedical applications [43, 27], absorbents of pollutants in agriculture
[27], and smart sensors or actuators in engineering [10, 32]. These properties
can be fine-tuned by manipulating chemical composition, crosslinking density,
and fluid content [9]. In this work, we aim to provide an in-depth examination
of diffusion-deformation hydrogel theories and the implications of selecting
underlying constitutive theory in the large deformation setting.

Several comprehensive reviews have emerged in the study of hydrogels,
discussing topics ranging from hydrogel deformation theories to the
microstructural impact on mechanical behaviors. [36] provided an exhaustive
review of hydrogel deformation theories, blending both theoretical analyses
and practical applications, emphasizing mechanics’ role. [26] delve into
advances in constitutive models for hydrogels and shape memory polymers
(SMPs), categorizing six primary hydrogel types and highlighting potential
hyperelastic model adaptations. Meanwhile, [30] offered insights into hydrogel
network models, discussing the microstructural impact on their mechanical
behaviors like swelling, elasticity, and fracture. This research underscores
the potential synergy between network modeling and continuum mechanics in
capturing hydrogel dynamics.

Modeling the diffusion-deformation process in hydrogels necessitates a
comprehensive understanding of the material’s behavior during swelling and
mechanical deformation. This understanding hinges on a mathematical
formulation that simultaneously accounts for the diffusion of fluid molecules
within the hydrogel’s polymer network and the corresponding mechanical
deformations induced by swelling. Two primary approaches can be employed to
describe these phenomena: mixture theories and macro-scale theories. Mixture
theories, representing the porous medium as spatially superimposed interacting
layers [18], often involve a complex array of model parameters and
constitutive choices, making them challenging to calibrate in practical
applications. Consequently, this paper focuses on macro-scale poroelastic
theories, which view the medium as a homogenous material characterized by a
coupled deformation-diffusion response, as developed in the seminal works by
Biot [e.g., 3]. Within this framework, the key equations, such as the mass
balance and mechanical equilibrium equations, govern both the conservation of
mass and the balance of forces within the hydrogel. An appropriate chemo-
mechanical constitutive description for the macroscopic continuum completes
the modeling framework. Constitutive laws couple the mechanical response of
the polymer network (accounting, for example, for hyperelastic or viscoelastic
effects) with changes in fluid distribution within the hydrogel mesh,
describing mixing effects and the interaction between diffusion and
deformation.

Over the years, various models have been developed to describe coupled
diffusion-deformation effects in elastomeric gels. These models couple a
common general poroelastic framework with the effects of different physical
fields, resulting in a wide range of stimuli-based responses encountered in
various applications. Regarding the general poroelastic framework, for
instance, [25] and [49] provided a continuum mechanics framework, as well as
analytical and finite element solutions, of the coupled problem in the large
deformation setting. [12] provided a comprehensive thermodynamics-consistent
formulation of the diffusion-deformation theory under isothermal conditions.
In contrast, [38] benchmarked the diffusion-deformation theory against some
experiments involving localized exposure of the gel boundary to a solvent,
where large bending deformations appear during solvent absorption. With
respect to more complex stimuli-based responses, [13] extended the theory and
accounted for temperature effects as well. [14] summarized the main
developments of the thermo-mechanically coupled theory for fluid permeation in
elastomeric materials and provided an open-access Abaqus implementation to
simulate hydrogel’s response in 3D.

A shared mathematical characteristic of these models is the saddle point
structure of the underlying variational formulation. As a result of the finite
element implementations, the Ladyzhenskaya-Babuska-Brezzi (LBB), i.e.,
discrete inf-sup, condition might be violated, requiring equilibrated inf-sup
stable finite element spaces. Otherwise, oscillatory distributions of chemical
potential as the primary variable controlling species diffusion arise. These
challenges can be addressed in several ways. To name a few of these
approaches, in [7], an in-depth analysis of Taylor-Hood elements (see, e.g.,
[20]) for gel modeling is provided, while [29] proposed an enhanced assumed
strain (EAS) method for the coupled problem. [4] proposed a minimization
formulation for the coupled diffusion-deformation problem of polymeric
hydrogels at large strains and compared this variational framework to
classical saddle point-based structures. Large volume changes with instability
patterns in the presence of geometrical constraints were successfully modeled.
[44] developed a numerical platform to simulate the dynamic behaviors of
responsive gels. The authors of this study addressed some of the challenges
that were not previously resolved, such as how to handle time-dependent and
coupled mass diffusion and deformation fields in a very short time, as well as
numerical instability issues.

Overall, there is no consensus on formulating a constitutive model for
hydrogels and its numerical implementation. This study emphasizes the
essential building blocks found in common across various hydrogel theories,
mainly focusing on the standard poroelastic description. It is unclear whether
the underlying constitutive choices at the free energy level are equivalent.
The significance of fundamental choices in coupling diffusion and deformation
is frequently overlooked. For example, such choices may relate to whether
volume changes in the hydrogel are solely attributed to fluid content or
involve elastic mechanisms as well. In the latter scenario, the detailed
description of the volumetric response plays a crucial role. We aim to analyze
the impact of such assumptions on the constitutive model.

Not only do these fundamental constitutive choices lead to differences in the
observed physical response, but they also influence the feasibility and
appropriateness of different numerical solutions. Some theories lend
themselves to a monolithic implementation of the coupled problem, while others
favor a staggered approach. Despite their importance, these considerations
have been largely neglected in existing literature.

###  1.1 General goals

This paper reviews some of the most notable models that describe the
diffusion-deformation process of elastomeric gels. The models were selected
based on their consistency with thermodynamic principles, their ability to
represent different behaviors, and the easiness to reproduce the results using
the open-source computing platform FEniCS [2], serving as a benchmark for
future studies. We provide a common mathematical framework and classifications
from which both staggered and monolithic numerical algorithms are derived for
studying the diffusion-deformation behavior of hydrogels despite their various
derivation methods. We carefully examine the distinctions and similarities
between different constitutive equations and analyze their impact on the
deformation of hydrogels and the evolution of associated fields through
simulations. This lays the foundation for future theoretical extensions,
including additional mechanisms such as chemical reactions, degradation, or
damage.

###  1.2 Open science context

Traditionally, researchers have relied on custom numerical solvers or
commercial software like Abaqus or COMSOL. While [14] made their Abaqus code
publicly available, a recent push has been towards open-source platforms for
addressing coupled multiphysics models. Emerging software libraries such as
deal.II111https://www.dealii.org/, OpenFOAM222https://www.openfoam.com/,
MOOSE333https://mooseframework.inl.gov/, and
FEniCS444https://fenicsproject.org/ exemplify this trend, signaling a shift in
the way researchers approach the complex challenges of modeling hydrogels.

Unlike commercial software packages, open-source projects, like FEniCS,
generally allow the user to have a direct hand in implementing the problem
statement, the discretization, the numerical solution, and specific
manipulations since all instances can be accessed. To this end, more
information about the mathematical-numerical classification and design of
algorithms is required. While this provides an advantage, granting users more
liberty to develop and test novel numerical algorithms and discretizations, it
also constitutes a challenge. The need for more hands-on implementation and
debugging is time-consuming.

To include readerships interested in open-source software and as part of the
broader scope of our paper, we establish mathematical classifications to
derive numerical algorithms for the coupled problem, enabling the precise
study of well-posedness and facilitating rigorous numerical analysis. This
includes introducing inf-sup stable Taylor-Hood finite elements for spatial
discretization and using a strongly $A$-stable, first-order, implicit Euler
scheme for temporal discretization. These technical developments are
summarized into compact mathematical formulations that serve as the starting
point for our FEniCS implementation, allowing for meticulous examination of
the numerical algorithms.

###  1.3 Main contributions and paper outline

The main contributions of this paper are summarized as follows:

  1. 1.

Unification in the formulation and notation of some representative diffusion-
deformation theories applied to the large deformation of hydrogels.

  2. 2.

Mathematical classification and subsequent derivation and computational
comparison of monolithic and staggered numerical solutions accuracy in a
variational setting using the finite element method (FEM).

  3. 3.

One- to three-dimensional numerical simulations of well-known prototype
problems to study the theories’ capabilities and robustness of the FEM
implementation.

  4. 4.

Comparison of the different constitutive models in a single benchmark problem
highlighting their main characteristics.

The paper is organized as follows. Section 2 provides a common background
knowledge of gels and the basic formulations of the chosen theories using a
unified notation. Section 3 is devoted to the mathematical classification and
numerical approximation using the FEM in variational settings. Section 4
explains the general setting of the numerical simulations campaign. Section 5
presents the simulation results for some well-known prototype problems of the
diffusion-deformation of hydrogels. Section 6 presents a unified benchmark
problem for the diffusion-deformation of hydrogels. Concluding remarks are
given in Section 7.

##  2 Nonlinear theory for the diffusion-deformation of elastomeric gels

In this section, we summarize some of the well-known theories describing the
diffusion-deformation mechanisms in elastomeric gels that undergo large
deformations under isothermal conditions.

As notation rules, we denote gradient in the reference and current
configuration by $\text{Grad}(\bullet)$ and $\text{grad}(\bullet)$,
respectively, whereas the divergence in the reference and current
configuration is denoted by $\text{Div}(\bullet)$ and $\text{div}(\bullet)$,
respectively. The time derivative of any field is denoted by
$\partial_{t}(\bullet)$. The operator $\text{tr}(\mathbf{A})$ refers to the
trace of the second-order tensor $\mathbf{A}$. We denote the spatial dimension
with $d$, and in this paper, we exclusively work with $d=3$. Finally, let
$I:=(0,T)$ be the time interval with end time value $T>0$ and $\bar{I}$ its
closure.

###  2.1 Kinematics of the deformation

Consider a continuum homogeneous elastomeric body $\mathcal{B}$ living in the
Euclidean space $\mathbb{R}^{3}$ and its boundary
$\partial\mathcal{B}=\partial\mathcal{B}_{\mathbf{u}}\cup\partial\mathcal{B}_{\bar{\mathbf{t}}}=\partial\mathcal{B}_{\varphi}\cup\partial\mathcal{B}_{\bar{J}_{R}}$.
Here, $\partial\mathcal{B}_{\mathbf{u}}$ denotes the displacement boundary
(Dirichlet), $\partial\mathcal{B}_{\bar{\mathbf{t}}}$ the traction boundary
(Neumann), $\partial\mathcal{B}_{\varphi}$ the fluid concentration-related
boundary (Dirichlet), and $\partial\mathcal{B}_{\bar{J}_{R}}$ the fluid flux
boundary (Neumann). The outward normal vector to the domain boundaries in the
reference configuration is denoted by $\mathbf{n}_{R}\in\mathbb{R}^{d}$. The
current configuration of this homogeneous body at any instant of time can be
described by a one-to-one transformation mapping
$\bm{\varphi}_{t}:\mathcal{B}_{R}\rightarrow\mathbb{R}^{d}$, where
$\mathcal{B}_{R}$ refers to the reference configuration. In principle,
selecting any state as the reference state $\mathcal{B}_{R}$ (e.g., the
stress-free dry gel) is possible. For convenience, authors have also chosen to
define the initial state $\mathcal{B}_{R}$ as an isotropically swollen
configuration from the dry state, leading to better reproduction of
experimental conditions [6, 24].

The position vector in the reference configuration
$\mathbf{X}\in\mathcal{B}_{R}$ is related to the one in the current
configuration $\mathbf{x}\in\mathcal{B}$ by
$\mathbf{x}=\bm{\varphi}_{t}(\mathbf{X},t)$. The displacement field reads
$\mathbf{u}(\mathbf{X},t)=\mathbf{x}-\mathbf{X}$. The transformation map
$\bm{\varphi}_{t}(\mathbf{X},t)$ can be described in terms of the deformation
gradient,

$\mathbf{F}=\text{Grad}(\mathbf{x})=\mathbf{I}+\text{Grad}(\mathbf{u}),$ (1)

with $J=J(\mathbf{X},t)=\text{det}\mathbf{F}>0$, the determinant, representing
the volume change of a volume element $dv=JdV$ from the reference ($dV$) to
the current ($dv$) configuration, and $\mathbf{I}$ being the identity tensor.

As is standard,

$\displaystyle\mathbf{C}$ $\displaystyle=\mathbf{F}^{T}\mathbf{F},$ (2)
$\displaystyle\mathbf{b}$ $\displaystyle=\mathbf{F}\mathbf{F}^{T},$ (3)

denote the right and left Cauchy-Green tensors, respectively. Additionally,
the first invariant of $\mathbf{C}$ is given by

$I_{1}({\bf F})=\text{tr}(\mathbf{C})=\text{tr}(\mathbf{F}^{T}\mathbf{F})\,.$
(4)

####  2.1.1 Chemical potential and swelling deformation

The solvent component of the hydrogel is described by introducing its chemical
potential $\mu$, that is, the energy absorbed or released due to a change in
its content. Fluid absorption/desorption is kinematically described through an
inelastic part of the deformation ${\bf F}_{f}$. The volume change
$J_{f}=\text{det}{\bf F}_{f}$ associated with fluid absorption/desorption is
linked with the referential fluid concentration variable $c_{R}$, i.e., the
number of absorbed fluid molecules per unit volume of the reference
configuration, by enforcing:

$J_{f}=1+\Omega c_{R}\,,$ (5)

where $\Omega$ denotes the volume of a mole of fluid molecules. This
relationship can also be equivalently described by introducing the polymer
volume fraction variable $\phi$, defined as:

$\phi=(1+\Omega c_{R})^{-1}=J_{f}^{-1}\,,$ (6)

resulting $0\leq\phi\leq 1$. The dry state of the gel corresponds to $\phi=1$,
and $\phi<1$ represents a swollen state.

A generally adopted choice is to assume an isotropic swelling deformation
${\bf F}_{f}$, hence reading as:

${\bf
F}_{f}=\lambda^{s}\mathbf{I}\quad\text{with}\quad\lambda^{s}=J_{f}^{1/3}=\phi^{-1/3}=(1+\Omega
c_{R})^{1/3}\,,$ (7)

where $\lambda^{s}$ represents the polymer network stretch due to swelling. It
is noteworthy that, since $c_{R}>0$ by definition, it results in
$J_{f},\phi,\lambda^{s}>0$.

####  2.1.2 Elastic deformation

The total deformation of a hydrogel is obtained from the superimposition of
the fluid-related deformation gradient ${\bf F}_{f}$ and the elastic one ${\bf
F}_{e}$. The latter originates from the effects of mechanical actions to
restore compatibility. Based on the previously introduced choices, we obtain:

${\bf F}=(1+\Omega c_{R})^{1/3}{\bf F}_{e}.\,$ (8)

At this standpoint, depending on the volume change associated with the elastic
deformation $J_{e}=\text{det}{\bf F}_{e}=J_{f}^{-1}J$, a general
classification between two classes of models is introduced:

  1. 1.

_elastic compressible models (non-constrained formulations)_. In this case,
the elastic part of the deformation is assumed to be compressible. It is then
allowed ${J}_{e}\neq 1$;

  2. 2.

_elastic incompressible models (constrained formulations)_. In this case, the
elastic part of the deformation is assumed to be perfectly incompressible, and
the total volume change of the hydrogel is related only to fluid volume
changes. In other words, it results in $J_{e}=1$ and the kinematic constraint,

$J\stackrel{{\scriptstyle!}}{{=}}J_{f}\,,$ (9)

has to be enforced within the theoretical formulation.

###  2.2 Governing partial differential equations

The two governing partial differential equations (PDEs) for the vector-valued
displacements $\mathbf{u}:\bar{\mathcal{B}}_{R}\to\mathbb{R}^{d}$ and scalar-
valued, time-dependent, fluid content in terms of the concentration
$c_{R}:\bar{\mathcal{B}}_{R}\times\bar{I}\to\mathbb{R}$ when expressed in the
reference configuration, consist of:

  1. 1.

The local form of the balance of linear momentum reads:

$\begin{cases}\text{Div}(\mathbf{P})+\mathbf{b}_{R}=\bm{0},&~{}\text{in}~{}\mathcal{B}_{R},\\\
\mathbf{u}=\bar{\mathbf{u}},&~{}\text{on}~{}\partial\mathcal{B}_{\mathbf{u}},\\\
\mathbf{P}\mathbf{n}_{R}=\bar{\mathbf{t}},&~{}\text{on}~{}\partial\mathcal{B}_{\bar{\mathbf{t}}},\\\
\mathbf{P}|_{t=0}=\mathbf{P}_{0},&~{}\text{in}~{}\mathcal{B}_{R}\,.\end{cases}$
(10)

Here, $\mathbf{P}:\mathcal{B}_{R}\to\mathbb{R}^{d\times d}$ denotes the first
Piola-Kirchhoff stress tensor, and
$\mathbf{P}_{0}:\mathcal{B}_{R}\to\mathbb{R}^{d\times d}$ its initial value at
$t=0$. An alternative stress measure also commonly employed is the Cauchy
stress tensor $\bm{\sigma}=J^{-1}\mathbf{P}\mathbf{F}^{T}$. External actions
consist of body forces per unit deformed volume in the reference configuration
$\mathbf{b}_{R}:\mathcal{B}_{R}\to\mathbb{R}^{d}$. Moreover, boundary
conditions prescribe displacement
$\bar{\mathbf{u}}:\mathcal{B}_{R}\to\mathbb{R}^{d}$ and traction
$\bar{\mathbf{t}}:\partial\mathcal{B}_{\bar{\mathbf{t}}}\to\mathbb{R}^{d}$ on
separate portions of the boundary. Notice that inertial effects have been
neglected due to the considerably slow dynamics of the fluid diffusion
evolution w.r.t. the time scale of the wave propagation.

  2. 2.

The local form of the mass balance of fluid content inside the hydrogel reads:

$\begin{cases}\partial_{t}c_{R}+\text{Div}(\mathbf{J}_{R})=0,&~{}\text{in}~{}\mathcal{B}_{R}\times
I,\\\ \mu=\bar{\mu},&~{}\text{on}~{}\partial\mathcal{B}_{\varphi}\times I,\\\
-\mathbf{J}_{R}\cdot\mathbf{n}_{R}=\bar{J}_{R},&~{}\text{on}~{}\partial\mathcal{B}_{\bar{J}_{R}}\times
I,\\\
\mu|_{t=0}=\mu_{0},&~{}\text{in}~{}\mathcal{B}_{R}\times\\{t=0\\}\,.\end{cases}$
(11)

Here, $\mathbf{J}_{R}$ is the fluid flux vector in the reference
configuration, related to the one ${\bf j}$ in the current configuration via
$\mathbf{J}_{R}=J\mathbf{F}^{-1}{\bf j}$. Boundary conditions prescribe values
of chemical potential $\bar{\mu}:\partial\mathcal{B}_{\varphi}\times
I\to\mathbb{R}$ and fluid flux
$\bar{J}_{R}:\partial\mathcal{B}_{\bar{J}_{R}}\times I\to\mathbb{R}^{d}$ in
the reference configuration on the boundaries. Moreover,
$\mu_{0}:\partial\mathcal{B}_{\varphi}\to\mathbb{R}$ refers to the initial
value of the chemical potential inside the hydrogel. Notice that the mass
balance of fluid content is written in terms of $c_{R}$ and ${\bf J}_{R}$, but
its corresponding boundary and initial conditions involve a different
variable, the chemical potential $\mu$. The connection of equation (11) with
$\mu$ becomes clear by introducing a constitutive relation for ${\bf J}_{R}$,
e.g., through Fick’s laws of diffusion.

###  2.3 Constitutive equations: stress and chemical potential

This section introduces constitutive equations for the stress and chemical
potential, addressing both compressible or perfectly incompressible
formulations.

####  2.3.1 Compressible formulations

Following standard thermodynamic arguments [48, 21, 12, 13], the local form of
the second law of thermodynamics reads:

${\bf P}:\dot{\bf F}+\mu\dot{c}_{R}-{\bf
J}_{R}\cdot\text{Grad}(\mu)-\dot{\psi}_{R}\geq 0,$ (12)

where $\psi_{R}:\mathcal{B}_{R}\times I\to\mathbb{R}$ is the free energy
density function (per unit reference volume). Guided by equation (12), the
free energy density function can be regarded as a function of the total
deformation ${\bf F}$ and fluid concentration $c_{R}$, that is:

$\psi_{R}=\Psi_{R}({\bf F},c_{R})\,.$ (13)

Consequently, equation (12) can be reformulated as:

${\bf P}:\dot{\bf F}+\mu\dot{c}_{R}-{\bf
J}_{R}\cdot\text{Grad}(\mu)-\frac{\partial\Psi_{R}}{\partial{\bf F}}:\dot{\bf
F}-\frac{\partial\Psi_{R}}{\partial c_{R}}:\dot{c}_{R}\geq 0,$ (14)

and the following thermodynamically-consistent constitutive relations can be
established for the first Piola-Kirchhoff stress tensor (PK1):

$\mathbf{P}=\dfrac{\partial\Psi_{R}}{\partial\mathbf{F}},$ (15)

and for the chemical potential:

$\mu=\dfrac{\partial\Psi_{R}}{\partial c_{R}}\,.$ (16)

Alternative formulations can be found in the state-of-the-art for defining
stresses and chemical potential, built upon elastic stress and active chemical
potential concepts. These are reviewed and discussed in appendix A, showing
that both approaches lead to identical results.

####  2.3.2 Incompressible formulations

Incompressible models require a special treatment of the kinematic constraint
in equation (9) that has to be satisfied _a priori_ within the formulation. A
first possibility is to extend the free energy $\Psi_{R}$ in equation (13) in
the context of Lagrangian formulations by introducing a constrained free-
energy $\Psi_{R}^{c}$ that reads:

$\Psi_{R}^{c}({\bf F},c_{R},P)=\Psi_{R}({\bf F},c_{R})+p(J-J_{f}(c_{R})),$
(17)

where $p$ represents a pressure-like Lagrange multiplier to enforce the
kinematic constraint in equation (9) through a variational framework. By
replacing $\Psi_{R}$ with $\Psi_{R}^{c}$, equation (14) reads:

${\bf P}:\dot{\bf F}+\mu\dot{c}_{R}-{\bf
J}_{R}\cdot\text{Grad}(\mu)-\left(\frac{\partial\Psi_{R}}{\partial{\bf
F}}-pJ{\bf F}^{-T}\right):\dot{\bf F}-\left(\frac{\partial\Psi_{R}}{\partial
c_{R}}+\Omega p\right)\dot{c}_{R}\geq 0,$ (18)

since $\partial J/\partial{\bf F}=J{\bf F}^{-T}$ and $\partial J_{f}/\partial
c_{R}=\Omega$. From equation (18), the following thermodynamically consistent
choices can be introduced for the first Piola-Kirchhoff stress tensor:

$\mathbf{P}=\dfrac{\partial\Psi_{R}}{\partial\mathbf{F}}-pJ{\bf F}^{-T},$ (19)

and for the chemical potential:

$\mu=\dfrac{\partial\Psi_{R}}{\partial c_{R}}+\Omega p.$ (20)

Such an approach has been followed, for instance, by [6, 12].

Alternatively, the kinematic constraint in equation (9) can be embedded
directly within the free-energy function. For instance, this rationale is
described and adopted by [33, 34]. In this case, the fluid concentration
$c_{R}$ is regarded as a dependent variable since it is related to the
deformation gradient ${\bf F}$ (or displacements ${\bf u}$) through equation
(9). Hence, the free energy can be reformulated in the form:

$\psi_{R}=\tilde{\Psi}_{R}({\bf F})\,.$ (21)

Then, equation (14) reads:

${\bf P}:\dot{\bf F}+\frac{\mu}{\Omega}\dot{J}-{\bf
J}_{R}\cdot\text{Grad}(\mu)-\left(\frac{\partial\tilde{\Psi}_{R}}{\partial{\bf
F}}\right):\dot{\bf F}\geq 0,$ (22)

since $\dot{c}_{R}=\dot{J}_{f}/\Omega=\dot{J}/\Omega$ from the
incompressibility constraint. From equation (22) and noting that
$\dot{J}=J{\bf F}^{-T}:\dot{\bf F}$, the first Piola-Kirchhoff stress tensor
can be expressed as:

$\mathbf{P}=\frac{\partial\tilde{\Psi}_{R}}{\partial\mathbf{F}}-\frac{\mu
J}{\Omega}{\bf F}^{-T}.$ (23)

As previously noted, the fluid concentration $c_{R}$ can no longer be
considered as an independent variable. Therefore, the chemical potential $\mu$
shall now be considered as the problem’s primary variable. Consequently, a
special treatment is required for the $\partial_{t}c_{R}$ term in the mass
balance equation (28). As introduced by [33, 34], the following relationship
holds true from equation (9) between the time derivative of the fluid
concentration and the determinant of the deformation gradient:

$\partial_{t}c_{R}=\dfrac{1}{\Omega}\partial_{t}J_{f}=\dfrac{J}{\Omega}\text{div}(\partial_{t}\mathbf{u})=\dfrac{J}{\Omega}\text{Grad}(\partial_{t}\mathbf{u}):\mathbf{F}^{-T},$
(24)

which yields

$\dfrac{J}{\Omega}\text{Grad}(\partial_{t}\mathbf{u}):\mathbf{F}^{-T}+\text{Div}(\mathbf{J}_{R})=0,~{}\text{in}~{}\mathcal{B}_{R}\times
I\,.$ (25)

###  2.4 Constitutive equations: fluid flux

From equation (14), a thermodynamically motivated choice for the fluid flux
$\mathbf{j}$ is to assume that it is proportional to the gradient of the
chemical potential in the current configuration, namely,

$\mathbf{j}=-\dfrac{cD}{k_{B}T}\text{grad}(\mu),$ (26)

with $c$ the solvent concentration in the current configuration, related to
the nominal concentration by $c_{R}=Jc$, and $D$ is the solvent diffusivity
assumed to be a constant. Notice that $\text{grad}(\mu)$ can be pulled back to
the reference configuration by
$\text{grad}(\mu)=\mathbf{F}^{-T}\text{Grad}(\mu)$. Hence, the flux in the
reference configuration reads:

$\mathbf{J}_{R}=-J\mathbf{F}^{-1}\left(\dfrac{cD}{k_{B}T}\text{grad}(\mu)\right)=-\dfrac{c_{R}D}{k_{B}T}\mathbf{b}^{-1}\text{Grad}(\mu)\,.$
(27)

Inserting equation (27) into (11) yields for the mass balance equation:

$\partial_{t}c_{R}-\text{Div}\left(\dfrac{c_{R}D}{k_{B}T}\mathbf{b}^{-1}\text{Grad}(\mu)\right)=0,~{}\text{in}~{}\mathcal{B}_{R}\times
I\,,$ (28)

with $\mu=\mu({\bf F},c_{R})$ from equation (16).

Notice that to make the units of equation (28) consistent, $c_{R}$ must have
units of [mol], $k_{B}$ units of [J K-1], $T$ units of [K], $\mu$ units of [J
mol-1], $D$ units of [m2 s-1], and $\bm{b}$ is dimensionless. Some authors
have defined constitutive equation (26) in terms of $RT$, where $R$ is the
gasses constant, (see, e.g., [13] or [14]). In this case, $\mu$ units of [J].

###  2.5 Specialization of constitutive theories

Specific choices for the free energy function characterize different state-of-
the-art models linking stress and chemical potential variations with
diffusion-deformation mechanisms. We introduce constitutive models in
accordance with the framework outlined in equation (13) for compressible
formulations and equations (17) or (21) for incompressible formulations.

Unless explicitly specified otherwise, we choose to characterize fluid content
using the concentration variable $c_{R}$. As a result, we express fluid-
related volume changes and polymer volume fractions as functions of $c_{R}$,
denoting them as $J_{f}=J_{f}(c_{R})$ and $\phi=\phi(c_{R})$, respectively,
based on equations (5) and (6), respectively. By inverting these
relationships, we can reformulate the proposed theories using different
primary variables to describe fluid content whenever needed.

The free energy function is, in general, written in a separable additive form:

$\psi_{R}=\psi_{R}^{\text{mix}}+\psi_{R}^{\text{mech}}\,.$ (29)

Here, $\psi_{R}^{\text{mix}}$ describes the mixing of the solvent with the
polymer network. Overall, there is a rather consensus agreement that this is
well described by the Flory-Huggins/Flory-Rehner theory (see, e.g., [12, 6,
34]), reading:

$\psi_{R}^{\text{mix}}=\Psi_{R}^{\text{mix}}(c_{R})=\mu^{0}c_{R}+\dfrac{k_{B}T}{\Omega}\left[\Omega
c_{R}\ln\left(\dfrac{\Omega c_{R}}{1+\Omega
c_{R}}\right)+\chi\left(\dfrac{\Omega c_{R}}{1+\Omega c_{R}}\right)\right],$
(30)

where $\mu^{0}$ is the chemical potential of the unmixed pure solvent, $k_{B}$
refers to Boltzmann’s constant, $T$ is the absolute temperature, and $\chi$ is
a dimensionless parameter named Flory’s interaction parameter. The latter
represents the disaffinity between the polymer and the fluid. In particular,
if $\chi$ is increased, the fluid molecules are expelled from the gel, and it
shrinks, while if $\chi$ is decreased, the gel swells.

Furthermore, $\psi_{R}^{\text{mech}}$ is the contribution to the change in the
free energy due to the deformation of the polymer network, for which elastic
incompressible or compressible formulations differ by considering:

$\psi_{R}^{\text{mech}}=\begin{cases}\Psi_{R}^{\text{s}}(\mathbf{F},J)&\text{incompressible},\\\
\Psi_{R}^{\text{s}}(\mathbf{F},J)+\Psi_{R}^{\text{en}}(J,c_{R})&\text{compressible},\end{cases}$
(31)

where $\Psi_{R}^{\text{s}}$ is an entropic component and
$\Psi_{R}^{\text{en}}$ an energetic contribution.

The entropic component is usually defined following the arguments of classical
statistical mechanics models for rubber elasticity, [13]. For small to
moderate values of stretching, Gaussian statistics provide an estimate of the
entropy change due to mechanical stretching of the polymer network resulting
in the form of a Neo-Hooke material that takes the form [26]:

$\Psi_{R}^{\text{s}}(\mathbf{F},J)=\dfrac{G_{0}}{2}\left[I_{1}({\bf
F})-3-2\ln\left(J\right)\right]\,,$ (32)

where $I_{1}$ is given in equation (4) and $G_{0}\approx Nk_{B}T$ represents
the shear modulus, with $N$ being the number of polymer chains per unit
reference volume, i.e., crosslink polymer network density.

In contrast, different choices have been made by authors when it comes to
defining the energetic component $\Psi_{R}^{\text{en}}$ of the free energy due
to the deformation of the polymer network. Some available solutions will be
discussed in Section 2.5.2, and after that, some state-of-the-art perfectly
incompressible models will be presented.

####  2.5.1 Incompressible constitutive models

Two incompressible constitutive models are presented here.

Constitutive model I: Following the ideas by [25, 49, 34], this model assumes
a perfectly incompressible elastic material response, by introducing a
constrained material response within the rationale presented in equation (21).
Therefore, the constraint in equation (9) is directly embedded in the free-
energy function, leading to:

$\tilde{\Psi}_{R}({\bf F})=\left.\left(\Psi_{R}^{\text{s}}({\bf
F},J)+\tilde{\Psi}_{R}^{\text{mix}}(J)\right)\right|_{J=\text{det}({\bf
F})}\,,$ (33)

where the mixing part of the energy $\tilde{\Psi}_{R}^{\text{mix}}$ respects
$\tilde{\Psi}_{R}^{\text{mix}}(J)=\Psi_{R}^{\text{mix}}(c_{R}(J))$ with
$c_{R}(J)=(J-1)/\Omega$ from equation (5) under the condition of equation (9).
Considering the entropic component in equation (32), the first Piola-Kirchhoff
stress tensor is derived from equation (23) as:

$\mathbf{P}({\bf F},\mu)=\frac{\partial\Psi_{R}^{\text{s}}}{\partial{\bf
F}}+\left(\frac{\partial\Psi_{R}^{\text{s}}}{\partial
J}+\frac{\partial\tilde{\Psi}_{R}^{\text{mix}}}{\partial J}\right)J{\bf
F}^{-T}-\frac{J\mu}{\Omega}{\bf
F}^{-T}=\Big{[}G_{0}\left(\mathbf{b}-\mathbf{I}\right)+Jp_{\mu}(J,\mu)\mathbf{I}\Big{]}\mathbf{F}^{-T}\,.$
(34)

Here, $p_{\mu}$ is the equivalent volumetric stress:

$p_{\mu}(J,\mu)=-\dfrac{\mu}{\Omega}+\dfrac{k_{B}T}{\Omega}\left[\ln\left(1-\dfrac{1}{J}\right)+\dfrac{1}{J}+\dfrac{\chi}{J^{2}}\right].$
(35)

Constitutive model II: This model has been introduced by [12]. Also, in this
case, a perfectly incompressible elastic material response is assumed,
introducing a constrained material response by following a Lagrangian approach
as described in equation (17). By exploiting the kinematic constraints in
equations (5) and (9), the entropic component in equation (32) is re-
formulated as $\hat{\Psi}_{R}^{\text{e}}({\bf
F},c_{R})=\Psi_{R}^{\text{s}}({\bf F},J_{f}(c_{R}))$ and the total free-energy
reads:

$\Psi_{R}=\Psi_{R}^{\text{mix}}(c_{R})+\hat{\Psi}_{R}^{\text{e}}({\bf
F},c_{R})\,.$ (36)

The first Piola-Kirchhoff stress tensor is derived from equation (19)
yielding:

$\mathbf{P}({\bf F},P)=\frac{\partial\hat{\Psi}_{R}^{\text{s}}}{\partial{\bf
F}}-pJ{\bf F}^{-T}=G_{0}\left(\mathbf{b}-P\mathbf{I}\right)\mathbf{F}^{-T}\,,$
(37)

where $P=pJ$ is the pressure term in the reference configuration. To be
consistent with the original model, the fluid concentration $c_{R}$ is
replaced in the formulation with the polymer volume fraction $\phi$ by means
of equation (6). Accordingly, the chemical potential can be obtained from
equations (20), (30), (32) and (36) as:

$\mu(\phi,P)=\left.\left(\frac{\partial\hat{\Psi}_{R}^{\text{mix}}}{\partial
c_{R}}+\frac{\partial\hat{\Psi}_{R}^{\text{s}}}{\partial c_{R}}+\Omega
p\right)\right|_{c_{R}(\phi)}=\mu^{0}+k_{B}T\left[\ln\left(1-\phi\right)+\phi+\chi\phi^{2}\right]-\Omega
G_{0}\phi+\Omega P\phi\,,$ (38)

where $P=pJ=pJ_{f}=p\phi^{-1}$ follows directly from equations (5), (6) and
(9).

####  2.5.2 Compressible models

Compressible models are characterized by the superposition of an entropic
energetic component (given in equation (32)) and an energetic part of the free
energy. The latter reflects changes in the internal energy associated with the
volumetric mechanical deformation of the swollen elastomer. Three well-
established constitutive models for $\Psi_{R}^{\text{en}}$ are here discussed:

$\Psi_{R}^{\text{en}}(J,c_{R})=\begin{cases}\Psi_{R,1}^{\text{en}}(J,c_{R})=\dfrac{K}{2}\left[J-J_{f}(c_{R})\right]^{2}&\text{\cite[citep]{[\@@bibref{Number}{Bouklas2015Nonlinear}{}{}]}}\\\
\Psi_{R,2}^{\text{en}}(J,c_{R})=\dfrac{K}{2}\ln(J\phi(c_{R}))^{2}&\text{\cite[citep]{[\@@bibref{Number}{Chester2011thermo}{}{}]}}\\\
\Psi_{R,3}^{\text{en}}(J,c_{R})=\phi(c_{R})^{-1}\left[\dfrac{K}{2}\ln(J\phi(c_{R}))^{2}\right]&\text{\cite[citep]{[\@@bibref{Number}{Chester2015Abaqus}{}{}]}}\end{cases}$
(39)

where $K>0$ represents the bulk modulus of the gel.

Constitutive model III: The following formulation is based on the ideas by
[7]. To be consistent with the original model, the swelling volume change
$J_{f}$ is considered in place of $c_{R}$ within the formulation, but we
highlight that this is straight linked through equation (5). Recalling that
$J=\text{det}({\bf F})$, the first Piola-Kirchhoff stress tensor considers the
entropic and energetic components in equations (32) and (39)1, reading from
equation (15):

$\mathbf{P}({\bf F},J_{f})=\frac{\partial\Psi_{R}^{\text{s}}}{\partial{\bf
F}}+\left(\frac{\partial\Psi_{R}^{\text{s}}}{\partial
J}+\frac{\partial\Psi_{R,1}^{\text{en}}}{\partial J}\right)J{\bf
F}^{-T}=G_{0}\left(\mathbf{F}+\alpha({\bf F},J_{f})J{\bf F}^{-T}\right),$ (40)

with:

$\alpha({\bf F},J_{f})=-\dfrac{1}{J}+\dfrac{K}{G_{0}}\left(J-J_{f}\right).$
(41)

The chemical potential is obtained from equations (16), (29), (30) and (39)1
as:

$\mu({\bf F},J_{f})=\left.\frac{\partial\Psi_{R}^{\text{mix}}}{\partial
c_{R}}\right|_{c_{R}(J_{f})}=\mu^{0}+k_{B}T\left[\ln\left(1-\dfrac{1}{J_{f}}\right)+\dfrac{1}{J_{f}}+\dfrac{\chi}{J_{f}^{2}}\right]-\Omega
K\left(J-J_{f}\right).$ (42)

It is worth noting that the mass balance equation (28) reads in this context
as:

$\dfrac{1}{\Omega}\partial_{t}J_{f}-\text{Div}\Bigl{(}{\bf M}({\bf
F},J_{f})\text{Grad}(\mu({\bf
F},J_{f}))\Bigr{)}=0\quad\text{in}~{}\mathcal{B}_{R}\times I,$ (43)

where ${\bf M}={\bf M}({\bf F},J_{f})$ denotes the species mobility tensor:

${\bf M}({\bf
F},J_{f})=\dfrac{D}{k_{B}T}\dfrac{J_{f}-1}{\Omega}\mathbf{b}^{-1}\,.$ (44)

Constitutive models IV and V: Starting from the original incompressible model
introduced in [12] (_constitutive model II_), the same group of authors
presented alternative compressible formulations by adding energetic components
as given in equations (39)2 ([13], _constitutive model IV_) and (39)3 ([14],
_constitutive model V_). To be consistent with the original models, the
polymer volume fraction $\phi$ is considered in place of $c_{R}$ within the
formulation, but we highlight that this is straight linked through equation
(6). Hence, recalling that $J=\text{det}({\bf F})$, the first Piola-Kirchhoff
stress tensor follows, for _constitutive model IV_ , from equation (15) with
equations (32) and (39)2 as:

$\mathbf{P}({\bf F},\phi)=\dfrac{\partial\Psi_{R}^{\text{s}}}{\partial{\bf
F}}+\left(\dfrac{\partial\Psi_{R}^{\text{s}}}{\partial
J}+\dfrac{\partial\Psi_{R,2}^{\text{en}}}{\partial J}\right)J{\bf
F}^{-T}=\left[G_{0}\left(\mathbf{b}-\mathbf{I}\right)+K\left(\ln(J\phi)\right)\mathbf{I}\right]\mathbf{F}^{-T}\,,$
(45a) and, for _constitutive model V_ with (39)3, as: $\mathbf{P}({\bf
F},\phi)=\dfrac{\partial\Psi_{R}^{\text{s}}}{\partial{\bf
F}}+\left(\dfrac{\partial\Psi_{R}^{\text{s}}}{\partial
J}+\dfrac{\partial\Psi_{R,3}^{\text{en}}}{\partial J}\right)J{\bf
F}^{-T}=\left[G_{0}\left(\mathbf{b}-\mathbf{I}\right)+\phi^{-1}K\left(\ln(J\phi)\right)\mathbf{I}\right]\mathbf{F}^{-T}\,.$
(45b)

Equation (16), together with equation (30), (39)2 and (39)3, yields the
chemical potential:

$\mu({\bf
F},\phi)=\begin{cases}\left.\left(\dfrac{\partial\Psi_{R}^{\text{mix}}}{\partial
c_{R}}+\dfrac{\partial\Psi_{R,2}^{\text{en}}}{\partial
c_{R}}\right)\right|_{c_{R}(\phi)}=\mu^{0}+g({\bf F},\phi),&\text{\emph{model
IV} }\\\ \left.\left(\dfrac{\partial\Psi_{R}^{\text{mix}}}{\partial
c_{R}}+\dfrac{\partial\Psi_{R,3}^{\text{en}}}{\partial
c_{R}}\right)\right|_{c_{R}(\phi)}=\mu^{0}+g({\bf F},\phi)+\dfrac{1}{2}\Omega
K\left(\ln(J\phi)\right)^{2}&\text{\emph{model V}},\end{cases}$ (46)

with:

$g({\bf
F},\phi)=RT\left[\ln\left(1-\phi\right)+\phi+\chi\phi^{2}\right]+\begin{cases}-\Omega
K\left(\ln(J\phi)\right)\phi&\text{\emph{model IV}}\\\ -\Omega
K\left(\ln(J\phi)\right)&\text{\emph{model V}}\end{cases}\,.$ (47)

It is noteworthy that, in the original papers, authors formulate the theories
based on the elastic PK1 and the active chemical potential, leading, however,
to identical results as proved in appendix A. Moreover, in these contexts, the
mass balance equation (28) is conveniently reformulated in terms of $\phi$,
instead of $c_{R}$, as:

$-\dfrac{1}{\Omega\phi^{2}}\partial_{t}\phi-\text{Div}\Bigl{(}{\bf M}({\bf
F},\phi^{-1})\text{Grad}(\mu({\bf
F},\phi))\Bigr{)}=0\quad\text{in}~{}\mathcal{B}_{R}\times I\,,$ (48)

where ${\bf M}$ is given in equation (44) and $J_{f}=\phi^{-1}$ from equation
(6).

####  2.5.3 Preliminary comparisons between models and final considerations

In the context of compressible models, we can gain valuable insights by
examining how different choices of the energy component $\Psi_{R}^{\text{en}}$
impact the stress constitutive relationships. To illustrate this, let us focus
on equation (40) within _constitutive model III_ , which primarily penalizes
substantial elastic deformations, that is, when $|J-J_{f}|=|J_{e}-1|J_{f}\gg
0$. In contrast, _constitutive models IV and V_ in equations (45) incorporate
penalties for both significant elastic deformations and the fully swollen
state, as evidenced by the behavior of $|\ln(J\phi)|=|\ln(J_{e})|$, which
approaches $+\infty$ when $|J_{e}-1|\gg 0$, and inversely, it approaches
$-\infty$ as the polymer volume fraction $\phi$ approaches $0$. However, it is
worth noting that _constitutive models IV and V_ diverge from each other in
treating large swelling deformations. Specifically, _constitutive model V_
penalizes these deformations, occurring when
$\phi^{-1}=J_{f}\rightarrow+\infty$, whereas _constitutive model IV_ does not.

Furthermore, as previously highlighted, _constitutive models IV and V_ serve
as the compressible counterparts to _constitutive model II_. This connection
becomes evident when we observe that the Lagrange multiplier $P$ in equation
(37) is effectively replaced by a term related to bulk modulus in equations
(45). However, from a numerical implementation standpoint, this substitution
carries significant implications. In general, $P$ represents an additional
primary variable that must be determined through the stationary conditions of
the constrained functional (17) with respect to the Lagrange multiplier. Only
under specific circumstances, such as a traction-free condition in the
presence of plane stresses, can $P$ be directly defined from equilibrium
conditions. In these instances, constitutive model II simplifies to having
only two primary variables, namely, u and $\phi$ (as demonstrated in, for
example, equation (78) in the subsequent Section 5.2). From this point
forward, we will exclusively focus on these special cases.

###  2.6 Weak formulations

The solution of the coupled PDE system consists of a vector-valued field of
displacements ($\mathbf{u}$) and, depending on the formulation, a scalar-
valued field ($\varphi$) given either by the concentration ($c_{R}$), polymer
volume fraction ($\phi$), or chemical potential ($\mu$). Hence, it results
$\varphi\in\\{c_{R},\phi,\mu\\}$.

Here we adopt standard notation for the usual Lebesgue and Sobolev spaces,
e.g., [47]. The functional space $H^{1}(\mathcal{B}_{R})^{d}$ is a Sobolev
space that consists of functions defined on a bounded domain
$X\subset\mathbb{R}^{d}$, with square integrable partial derivatives up to the
first order.

Here, $X:=\mathcal{B}_{R}$. Specifically, a function $w\in
H^{1}(\mathcal{B}_{R})^{d}$, if it satisfies the following conditions, namely
$w$ is square integrable: $\int_{\mathcal{B}_{R}}|w(x)|^{2}dx<\infty$ and the
first-order partial derivatives of $w$ exist and are square-integrable such
that $\int_{\mathcal{B}_{R}}|\nabla w(x)|^{2}dx<\infty$.

The norm associated with this space is given by

$||w||_{H^{1}(\mathcal{B}_{R})^{d}}:=\left(\int_{\mathcal{B}_{R}}|w(x)|^{2}dx+\int_{\mathcal{B}_{R}}|\nabla
w(x)|^{2}dx\right)^{1/2}.$ (49)

This norm induces a complete metric space with respect to which the functions
in $H^{1}(\mathcal{B}_{R})^{d}$ can be well-defined and approximated.

The coupled system of equations is formulated in terms of a variational
coupled system. To this end, we define the trial and test spaces as follows:

$\displaystyle Q$
$\displaystyle:=H^{1}\left(\mathcal{\mathcal{B}_{R}}\right),\quad
Q_{\bar{\varphi}}:=\Big{\\{}\varphi\in
Q\Big{|}~{}h(\varphi)=\bar{h}~{}~{}\text{on}~{}~{}\partial\mathcal{B}_{\varphi}\Big{\\}},\quad
Q_{0}:=\Big{\\{}\varphi\in
Q\Big{|}~{}\varphi=0~{}~{}\text{on}~{}~{}\partial\mathcal{B}_{\varphi}\Big{\\}},$
$\displaystyle V$
$\displaystyle:=H^{1}\left(\mathcal{\mathcal{B}_{R}}\right)^{d},\quad
V_{\bar{\mathbf{u}}}:=\Big{\\{}\mathbf{u}\in
V\Big{|}~{}\mathbf{u}=\bar{\mathbf{u}}~{}~{}\text{on}~{}~{}\partial\mathcal{B}_{\mathbf{u}}\Big{\\}},\quad
V_{0}:=\Big{\\{}\mathbf{u}\in
V\Big{|}~{}\mathbf{u}=0~{}~{}\text{on}~{}~{}\partial\mathcal{B}_{\mathbf{u}}\Big{\\}}.$

We notice that in $Q_{\bar{\varphi}}$, the boundary conditions may be given
explicitly or implicitly through the relation $h(\varphi)=\bar{h}$. This is
seen in the set of equations (11). To elucidate this, let us consider two
distinct scenarios: i. for constitute models II, IV, and V, we seek $\phi$,
and the boundary condition is $\phi=\bar{\phi}$, i.e., $h(\varphi):=\phi$ and
$\bar{h}:=\bar{\phi}$. On the other hand, ii. in constitutive model III, we
have $\varphi:=J_{f}$ and $\bar{h}:=\bar{\mu}$. Then, we solve the coupled
problem for $J_{f}$ and $\mu$, where $J_{f}$ is implicitly obtained from
equation (42). Once we know $J_{f}$ or $\phi$, then $c_{R}$ can be recovered
from (5).

To formulate both problem statements in an abstract fashion, we introduce for
the displacement system the semi-linear form
$a((\mathbf{u},\varphi))(\mathbf{v})$, which is nonlinear in the first
argument (trial function) and linear in the second argument (test function).
Furthermore, let $b(\mathbf{v})$ be the given right-hand side data. Next, for
the balance of fluid concentration, we use $c((\mathbf{u},\varphi))(q)$ and
$d(q)$. Then, the weak formulation can be written as:

######  Formulation 1.

(Diffusion-deformation of gels in $\mathcal{B}_{R}$). Find
$(\mathbf{u},\varphi)\in V_{\bar{\mathbf{u}}}\times Q_{\bar{\varphi}}$, with
$\varphi(0)=\varphi_{0}$, such that for $t\in I$ it holds

$\displaystyle a((\mathbf{u},\varphi))(\mathbf{v})+b(\mathbf{v})$
$\displaystyle=\mathbf{0},~{}~{}\forall\mathbf{v}\in V_{0},$ (50)
$\displaystyle c((\mathbf{u},\varphi))(q)+d(q)$ $\displaystyle=0,~{}~{}\forall
q\in Q_{0},$

where

$\displaystyle a((\mathbf{u},\varphi))(\mathbf{v})$
$\displaystyle=\int_{\mathcal{B}_{R}}\mathbf{P}(\mathbf{u},\varphi):\nabla\mathbf{v}dV,$
(51) $\displaystyle b(\mathbf{v})$
$\displaystyle=-\int_{\mathcal{B}_{R}}\mathbf{b}_{R}\cdot\mathbf{v}dV,$ (52)
$\displaystyle c((\mathbf{u},\varphi))(q)$
$\displaystyle=\int_{\mathcal{B}_{R}}\partial_{t}c_{R}(\varphi)\cdot
qdV+\int_{\mathcal{B}_{R}}{\bf M}({\bf
u},\varphi)\nabla\mu\left(\mathbf{u},\varphi\right)\cdot\nabla qdV,$ (53)
$\displaystyle d(q)$ $\displaystyle=0,$ (54)

with $``~{}:~{}"$ denoting the double contraction of the second-order tensors
$\mathbf{P}$ and $\nabla\mathbf{v}$, where $\mathbf{P}(\mathbf{u},\varphi)$ is
defined by either equation (34), (37), (40), or (45), depending of the
constitutive model adopted. Whereas ${\bf M}({\bf u},\varphi)$ is given by
equation (44).

Notice that the two balance equations are fully coupled through the
constitutive equations of the stress $\mathbf{P}(\mathbf{u},\varphi)$ and the
species mobility tensor ${\bf M}({\bf u},\varphi)$.

##  3 Classifications, discretization, and numerical solution

In this section, based on Formulation 1, we explain numerical coupling
strategies, provide mathematical classifications, and introduce spatial and
temporal discretizations. These derivations serve as a starting point for the
implementation in FEniCS. The reader is referred to the introduction to
understand the importance of this section, particularly when testing novel
algorithms, comparing them, and pursuing our own numerical developments,
including code debugging. This section provides the link between the strong
form problem statements in Section 2 and the numerical simulations carried out
in Section 5 by following the road map outlined in [46][Section 12.3].

###  3.1 Coupling strategies

There exist several ways for realizing numerically the coupling of several
PDEs. Here, we discuss two fundamental strategies to be implemented later,
namely, monolithic and partitioned approaches. In the former, Formulation 1,
it treated all-at-once, while in a partitioned approach, the two subproblems
are decoupled and solved in an iterative fashion. Here, we closely follow the
concepts and notation introduced by [45][Chapter 3].

Variational-monolithic coupling: In the variational-monolithic setting, the
coupling conditions are realized in an exact fashion in the weak (i.e.,
variational) formulation. Formulation 1 is given in such a variational-
monolithic fashion and, more specifically, the coupling conditions are of
volume-coupling type [45][Section 3.3.3].

In the monolithic approach, the entire PDE system can be either solved all-at-
once, which usually requires physics-based preconditioners. Either the system
is decoupled on the solver level within an outer monolithic iteration, e.g.,
GMRES (generalized minimal residual) or multigrid, and the preconditioner is
constructed based on decoupled subproblems. In general, monolithic solutions
can be computationally demanding depending on the complexity of the problem at
hand.

The monolithic scheme can be naturally extended to account for time-dependent
problems. In this case, we need to introduce a suitable time discretization
scheme and solve the previous problem at each time step. As an alternative, a
global space-time formulation can be formulated, discretized, and solved
accordingly; see e.g., [19] for a specific space-time multigrid realization
and analysis on the linear level.

Partitioned (staggered) approach: Conversely, in the partitioned approach, the
system of PDEs is broken down into smaller subsystems, and each subsystem is
solved independently using its own numerical method. The solutions of these
subsystems are then coupled to obtain the solution of the entire system. The
partitioned approach typically involves the following steps:

  1. 1.

Initialization: provide initial guesses $\hat{\bf u}^{0}$ and
$\hat{\varphi}^{0}$ for the unknown fields $\mathbf{u}$ and $\varphi$.

  2. 2.

Iteration: Let $(\hat{\bf u}^{j},\hat{\varphi}^{j})$ and $(\hat{\bf
u}^{j-1},\hat{\varphi}^{j-1})$ be the trial values of $\mathbf{u}$ and
$\varphi$ at the current and previous iterations, respectively.

######  Algorithm 1.

For $j=1,2,3,\dots$, iterate:

     * •

Given $\hat{\varphi}^{j-1}$, find $\hat{\bf u}^{j}$:

$a((\hat{\bf
u}^{j},\hat{\varphi}^{j-1}))(\mathbf{v})+b(\mathbf{v})=\mathbf{0}$ (55)

     * •

Given $\hat{\bf u}^{j}$, find $\hat{\varphi}^{j}$:

$c((\hat{\bf u}^{j},\hat{\varphi}^{j}))(q)+d(q)=0$ (56)

  3. 3.

Check for convergence: compare the updated and previous trial values. The
iteration is considered converged if the difference is below a specified
tolerance. That is, check whether

$\max\left(||\hat{\bf u}^{j}-\hat{\bf
u}^{j-1}||,||\hat{\varphi}^{j}-\hat{\varphi}^{j-1}||\right)<TOL.$ (57)

If the criterion is fulfilled, stop and assign ${\bf u}=\hat{\bf u}^{j}$ and
$\varphi=\hat{\varphi}^{j}$. If not, increment $j\rightarrow j+1$ and return
to step 2.

The previous procedure extends to discrete formulations of time-dependent
problems, for which a sequence of unknown fields $\mathbf{u}^{n}$ and
$\varphi^{n}$ at time points $t_{n}$ with $n=1,2,3,...,N$ is sought for. Here,
the previously introduced algorithm remains identical and performed at each
time point $t_{n}$. In this context, trial quantities at iteration step $j$
and time $t_{n}$ can be denoted as $\hat{\bf u}^{n,j}$ and
$\hat{\varphi}^{n,j}$.

###  3.2 Mathematical classification

We follow the ideas presented by [45] and classify the diffusion-deformation
prototype problem in Formulation 1. This is the starting point to design
appropriate algorithms, which are of interest in this work and have been
introduced before. Moreover, such classifications are required for
mathematical and numerical analysis, which both exceed the focus of this work.
To this end, we can formally analyze the problem statement as follows:

  1. 1.

Orders in time and space: Equation (10) represents a quasi-static problem with
no time derivatives and second order in space. In the thermodynamic context,
this quasi-static problem is in a thermodynamic equilibrium at each instance.
Equation (11) is a nonlinear, time-dependent, and of advection-diffusion type
when solved for $c_{R}$ or $\lambda^{s}$. Equation (11) is first order in time
and second order in space, i.e., a nonlinear parabolic PDE. On the other hand,
when it is solved for $\phi$, a minus sign appears in front of the time
derivative term (see equation (48)), which is rather unusual.

  2. 2.

Nonlinearities: They appear due to two reasons in Formulation 1. First, the
constitutive equations for $\mathbf{P}$ and $\mu$ are nonlinear. The
constitutive theories are formulated in a large deformation setting with
compressibility constraints. Second, the coupling terms enter in a nonlinear
fashion in the respective other problem.

  3. 3.

Type of coupling: The coupling is of domain-type and occurs via coefficients
and solution variables. In a partitioned approach (see Section 3.1), assuming
that one variable is given, the displacement PDE displays only the geometric
nonlinearity coming from the Neo-Hookean type of constitutive equation for the
deformation. On the other hand, the fluid balance concentration PDE can become
quasi-linear if the problem is solved for $c_{R}$ or $\lambda^{s}$ because a
lower-order term of the solution variable is multiplied with the highest
derivative. It can also remain fully nonlinear if the problem is solved for
$\phi$.

###  3.3 Temporal and spatial discretization

In this section, we discuss the discretization in time and space. A classical
finite difference scheme is employed for the temporal discretization of the
fluid balance concentration PDE, resulting in a quasi-stationary solution in
space at each time point. The spatial discretization is based on a Galerkin
FEM formulation [15]. Here, due to the structure of Formulation 1 of saddle-
point type, the inf-sup stability, i.e., LBB (Ladyzhenskaya-Babuska-Brezzi),
must be guaranteed, which requires on the discrete level using Taylor-Hood
elements, also known as $\mathcal{P}_{2}/\mathcal{P}_{1}$ elements. The
Taylor-Hood elements consist of quadratic basis functions ($\mathcal{P}_{2}$)
for the displacement field and linear basis functions ($\mathcal{P}_{1}$) for
the fluid concentration-related field. As just mentioned, the main reason for
using Taylor-Hood elements lies in their ability to satisfy the LBB stability
condition, which helps avoid numerical instabilities like spurious chemical
potential oscillations [7]. The notation in this section mainly matches the
notation adopted by [37] and [45].

Assume the computational domain $\mathcal{B}_{R}$ is partitioned into open
elements $K$ that depend on the spatial dimension $d$. A mesh consists of
quadrilateral, triangular, or hexahedron cells $K$, all of them available in
FEniCS. Here, we employ hexahedron cells. They perform a non-overlapping cover
of the computational domain $\mathcal{B}_{R}\subset\mathbb{R}^{d}$. Let
$\mathcal{T}_{h}=\\{K\\}$ be a conforming mesh of the bounded domain
$\mathcal{B}_{R}\subset\mathbb{R}^{d}$, with mesh size $h$.

We employ Taylor-Hood elements, i.e., a pair of finite element spaces
$(V_{h}^{TH},Q_{h}^{TH})$, where:

  * •

$Q_{h}^{TH}\subseteq Q$ is the space of continuous, piecewise linear functions
on $\mathcal{T}_{h}$, i.e.,

$Q_{h}^{TH}(\mathcal{T}_{h}):=\Big{\\{}q_{h}\in\left[C\left(\mathcal{B}_{h}\right)\right]\Big{|}~{}q_{h}|_{K}\in\left[\mathcal{P}_{1}\left(K\right)\right]~{}~{}\forall
K\in\mathcal{T}_{h},~{}q_{h}|_{\partial\mathcal{B}_{q_{h}}}=0\Big{\\}}.$ (58)

  * •

$V_{h}^{TH}\subseteq V$ is the space of continuous, piecewise quadratic
functions on $\mathcal{T}_{h}$, i.e.,

$V_{h}^{TH}(\mathcal{T}_{h}):=\Big{\\{}\mathbf{v}_{h}\in\left[C\left(\mathcal{B}_{h}\right)\right]^{d}\Big{|}~{}\mathbf{v}_{h}|_{K}\in\left[\mathcal{P}_{2}\left(K\right)\right]^{d}~{}~{}\forall
K\in\mathcal{T}_{h},~{}\mathbf{v}_{h}|_{\partial\mathcal{B}_{\mathbf{v}_{h}}}=0\Big{\\}}.$
(59)

Here, $\mathcal{P}_{1}(K)$ denotes the space of linear polynomials over the
element $K$ in $\mathbb{R}$ and $\mathcal{P}_{2}(K)^{d}$ denotes the space of
quadratic polynomials over the element $K$ in $\mathbb{R}^{d}$.

Thus, the discrete variational monolithic formulation for the diffusion-
deformation model reads:

######  Formulation 2.

(Semi-discrete in space variational monolithic diffusion-deformation of gels
in $\mathcal{T}_{h}$). Find
$(\mathbf{u}_{h},\varphi_{h})\in\\{\bar{\mathbf{u}}_{h}|_{\partial\mathcal{B}_{\bar{u}}}+V_{h}^{TH}\\}\times\\{\bar{\varphi}_{h}|_{\partial\mathcal{B}_{\bar{\varphi}}}+Q_{h}^{TH}\\}$
, with $\varphi_{h}(0)=\varphi_{0}$, such that for $t\in I$ it holds

$\displaystyle
a((\mathbf{u}_{h},\varphi_{h}))(\mathbf{v}_{h})+b(\mathbf{v}_{h})$
$\displaystyle=\mathbf{0},~{}~{}\forall\mathbf{v}_{h}\in V_{h}^{TH},$ (60)
$\displaystyle c((\mathbf{u}_{h},\varphi_{h}))(q_{h})+d(q_{h})$
$\displaystyle=0,~{}~{}\forall q_{h}\in Q_{h}^{TH}.$

Moreover, the discrete variational formulation using the staggered approach
reads:

######  Formulation 3.

(Semi-discrete in space variational diffusion-deformation of gels in
$\mathcal{T}_{h}$). Find
$\mathbf{u}_{h}\in\\{\bar{\mathbf{u}}_{h}|_{\partial\mathcal{B}_{\bar{u}}}+V_{h}^{TH}\\}$
and
$\varphi_{h}\in\\{\bar{\varphi}_{h}|_{\partial\mathcal{B}_{\bar{\varphi}}}+Q_{h}^{TH}\\}$,
with $\varphi_{h}(0)=\varphi_{0}$, such that, for the nonlinear iterations
$j=1,2,\dots$ and $t\in I$ it holds

$\displaystyle a((\hat{\bf
u}_{h}^{j},\hat{\varphi}_{h}^{j-1}))(\mathbf{v}_{h})+b(\mathbf{v}_{h})$
$\displaystyle=\mathbf{0},~{}~{}\forall\mathbf{v}_{h}\in V_{h}^{TH},$ (61)
$\displaystyle c((\hat{\bf u}_{h}^{j},\hat{\varphi}_{h}^{j}))(q_{h})+d(q_{h})$
$\displaystyle=0,~{}~{}\forall q_{h}\in Q_{h}^{TH},$

until the iteration converges, i.e., equation (57) is fulfilled.

The time-dependent term in the balance of fluid concentration is approximated
using the first-order implicit Euler discretization for $n=1,2,\dots,N_{f}$,
with $N_{f}$ the final simulation time index at the final time $T$, as

$\partial_{t}c_{R}(\varphi)\approx\frac{c_{R}(\varphi^{n})-c_{R}(\varphi^{n-1})}{\Delta
t},$ (62)

where $\varphi^{n}$ and $\varphi^{n-1}$ are the value of $\varphi$ at the
current and previous time step, respectively, and $\Delta t$ is the time step
increment. Following [45][Chapter 5, Definition 52, p. 90], we split the semi-
linear form $c((\mathbf{u},\varphi))(q)$ into time derivative and non-
derivative terms. To this end, we have

$\displaystyle c((\mathbf{u},\varphi))(q)$
$\displaystyle:=c_{T}((\varphi))(q)+c_{E}((\mathbf{u},\varphi))(q)$
$\displaystyle:=\int_{\mathcal{B}_{R}}\partial_{t}c_{R}(\varphi)\cdot
q\,dV+\int_{\mathcal{B}_{R}}{\bf M}({\bf
u},\varphi)\nabla\mu\left(\mathbf{u},\varphi\right)\cdot\nabla qdV\,.$ (63)

Then, the difference approximation of the time derivative with time step
increment $\Delta t$ yields:

$c_{T}((\varphi))(q)\approx c_{T}^{\Delta
t}((\varphi^{n}))(q):=\int_{\mathcal{B}_{R}}\frac{c_{R}(\varphi^{n})-c_{R}(\varphi^{n-1})}{\Delta
t}\cdot qdV\,.$ (64)

Then, the fully discrete variational monolithic formulation for the diffusion-
deformation model reads:

######  Formulation 4.

(Fully-discrete in space variational monolithic diffusion-deformation of gels
in $\mathcal{T}_{h}$). Let $h$ be the spatial discretization parameter and $n$
the current time point index. Find
$(\mathbf{u}_{h}^{n},\varphi_{h}^{n})\in\\{\bar{\mathbf{u}}_{h}|_{\partial\mathcal{B}_{\bar{u}}}+V_{h}^{TH}\\}\times\\{\bar{\varphi}_{h}|_{\partial\mathcal{B}_{\bar{\varphi}}}+Q_{h}^{TH}\\}$,
with $\varphi(0)=\varphi_{0}$, such that for $n=1,\ldots,N_{f}$ it holds

$\displaystyle
a((\mathbf{u}_{h}^{n},\varphi_{h}^{n}))(\mathbf{v}_{h})+b(\mathbf{v}_{h})$
$\displaystyle=\mathbf{0},~{}~{}\forall\mathbf{v}_{h}\in V_{h}^{TH},$ (65)
$\displaystyle c_{T}^{\Delta
t}((\varphi_{h}^{n}))(q_{h})+c_{E}((\mathbf{u}_{h}^{n},\varphi_{h}^{n}))(q_{h})+d(q_{h})$
$\displaystyle=0,~{}~{}\forall q_{h}\in Q_{h}^{TH}.$

We notice that the abstract cycle from monolithic problem statements until the
final linear solution is outlined in [45][Section 7.8.4]. Using the same
abstract concept, but replacing the monolithic nonlinear solution with some
iteration (Formulation 3), we obtain the fully discrete variational
formulation using the staggered approach:

######  Formulation 5.

(Fully-discrete in space variational diffusion-deformation of gels in
$\mathcal{T}_{h}$). Let $h$ be the spatial discretization parameter and $n$
the current time point index. Find
$\mathbf{u}_{h}^{n}\in\\{\bar{\mathbf{u}}_{h}|_{\partial\mathcal{B}_{\bar{u}}}+V_{h}^{TH}\\}$
and
$\varphi_{h}^{n}\in\\{\bar{\varphi}_{h}|_{\partial\mathcal{B}_{\bar{\varphi}}}+Q_{h}^{TH}\\}$,
with $\varphi(0)=\varphi_{0}$, such that, for the nonlinear iterations
$j=1,2,\dots$ and $n=1,\dots,N_{f}$ it holds

$\displaystyle a((\hat{\bf
u}_{h}^{n,j},\hat{\varphi}_{h}^{n,j-1}))(\mathbf{v}_{h})+b(\mathbf{v}_{h})$
$\displaystyle=\mathbf{0},~{}~{}\forall\mathbf{v}_{h}\in V_{h}^{TH},$ (66)
$\displaystyle c_{T}^{\Delta
t}((\hat{\varphi}_{h}^{n}))(q_{h})+c_{E}((\mathbf{u}_{h}^{n},\varphi_{h}^{n}))(q_{h})+d(q_{h})$
$\displaystyle=0,~{}~{}\forall q_{h}\in Q_{h}^{TH},$

until the iteration converges at time point $t_{n}$, i.e., equation (57) is
fulfilled, and then proceeds to $t_{n+1}$. Then, we set as initial guesses of
the nonlinear iterative scheme at the next time step $(\hat{\bf
u}_{h}^{n+1,0},\hat{\varphi}_{h}^{n+1,0}):=(\mathbf{u}_{h}^{n},\varphi_{h}^{n})$.

###  3.4 Numerical solution

The numerical consequences of our previous classification (specifically the
type of nonlinearities) are that we will always have to solve at least one
fully nonlinear PDE, independently of whether the problem is solved in a
monolithic, Formulation 4, or staggered way, Formulation 5. Here, we utilize a
Newton-type solver [17] and the consistent linearization of the system of PDEs
at each nonlinear iteration step. Specifically, for Formulation 4, a
convenient way (see [45][Section 3.3.3.2 and Section 3.3.3.3]) is to formulate
a common semi-linear form

$A(\mathbf{U}_{h}^{n})(\mathbf{\Psi}_{h}):=a((\mathbf{u}_{h}^{n},\varphi_{h}^{n}))(\mathbf{v}_{h})+c_{T}^{\Delta
t}((\varphi_{h}^{n}))(q_{h})+c_{E}((\mathbf{u}_{h}^{n},\varphi_{h}^{n}))(q_{h})\,,$
(67a) and corresponding right-hand side:
$F(\mathbf{\Psi}_{h}):=b(\mathbf{v}_{h})+d(q_{h})$ (67b)

with the joint unknown, trial function,
$\mathbf{U}_{h}^{n}:=(\mathbf{u}_{h}^{n},\varphi_{h}^{n})$ and the joint test
function $\mathbf{\Psi}_{h}:=(\mathbf{v}_{h},q_{h})$. This corresponds to Step
4 [45][Section 7.8.4.1]. Then, we can proceed with Step 5 (Newton’s method) in
[45][Section 7.8.4.1] for the nonlinear solution. The resulting solving system
is analogous to equations (67) for Formulation 5, but with
$\mathbf{U}_{h}^{n}$ replaced by $\hat{\mathbf{U}}_{h}^{n,j}:=(\hat{\bf
u}_{h}^{n,j},\hat{\varphi}_{h}^{n,j})$ that collects trial values of unknown
fields at step $j$ within the iterative solution scheme.

In this work, automatic differentiation offered by FEniCS was employed rather
than calculating the Jacobian by hand. Sparse LU decomposition (Gaussian
elimination) is used inside each Newton step to solve the arising linear
equation systems. As we have the specific derivations in Formulation 4 and
Formulation 5 at hand, a future extension is to employ iterative methods, like
GMRES - generalized minimal residuals [41], or multigrid solvers [23], for
which however, preconditioners need to be developed.

Finally, it is well known that Dirichlet boundary conditions on the chemical
potential, see equation (11)4, might be the source of spurious numerical
oscillations due to large sudden pressures within the hydrogel at the start of
the simulation. Hence, whenever needed, such boundary condition is
incrementally applied during the simulation, a strategy known as time-ramping
boundary condition. In each case, the boundary condition is increased as fast
as possible to meet a good compromise between numerical stability and physical
reality. This is achieved by introducing a time-dependent exponential term
that multiplies the Dirichlet boundary condition, i.e.,
$h(\varphi,t)=\bar{h}(\varphi)\left(1-\text{exp}(-\alpha_{r}t)\right)$, where
$\alpha_{r}$ is a positive constant that determines the rate of ramping.

##  4 General settings for the numerical simulation campaign

The coupled diffusion-deformation problem is faced by solving equations (10)
with null body forces ${\bf b}_{R}={\bf 0}$ and (11) and discretized either as
in equation (65) (Formulation 4, monolithic) or equation (66) (Formulation 5,
staggered) together with the specific choices of Constitutive models I to V
introduced in Section 2. Table 1 summarizes the main features of the models
considered in this work and refers to their respective equations.

Two campaigns of numerical simulations will be presented. Section 5 addresses
different representative prototype problems by adopting parameter settings
presented in the original papers where each constitutive model has been
originally presented. The aim is to explore the robustness of the developed
numerical implementations. Section LABEL:res:benchmark presents a unified
benchmark problem to have a unique reference simulation example that shows the
differences in the response for each constitutive equation.

The solution strategy outlined in Section 3.1 is adopted in all simulation
cases. The algorithms are implemented in FEniCS, and the code is provided so
the reader can reproduce and verify the results (online repository link:
https://doi.org/10.25835/5v49yfk0). The following material parameters remain
constant in all cases: $k_{B}=1.38065\times 10^{-23}~{}[Jmol^{-1}]$,
$T=298~{}[K]$, $\Omega=1.7\times 10^{-28}~{}[mol^{-1}]$ and
$N\Omega=0.001~{}[m^{3}mol^{-1}]$. The remaining material parameters vary
depending on the specific constitutive model and the adopted simulation setup.
They can be found in the caption of the figures associated with each numerical
result. For the sake of notation, let ${\bf E}_{1},{\bf E}_{2},{\bf E}_{3}$ be
introduced as a Cartesian coordinate system in the reference configuration
(resp., ${\bf e}_{1},{\bf e}_{2},{\bf e}_{3}$ in the current one),
parametrized in $X$, $Y$ and $Z$ (resp., $x$, $y$ and $z$). Moreover, let the
following stress components be introduced:

$P_{X}={\bf P}:({\bf e}_{1}\otimes{\bf E}_{1})\,,\quad P_{Y}={\bf P}:({\bf
e}_{2}\otimes{\bf E}_{2})\,,\quad\sigma_{x}=\bm{\sigma}:({\bf
e}_{1}\otimes{\bf e}_{1})\,.$ (68)

Numerical results will be reported in terms of displacements $\mathbf{u}$,
stretch $\lambda$, and the evolution of $\mu$ since they are the most relevant
to analyze from a physical viewpoint.

Table 1: Models’ main features summary. Constitutive model | Reference | Equations |  Primary variables | Solution strategy |  Incompressible
---|---|---|---|---|---
monolithic | staggered | yes | no
I | [34] | (10), (25), and (34) |  $\mu$ and $\mathbf{u}$ | x |  | x |
II |  [12]111Constitutive model II originally has $P$ as a primary variable but, as discussed in Section 2.5.3, a functional dependency $P=P(\phi)$ can be analytically determined in some special cases, like those addressed in this work.
Note: The rationale behind the selected solution strategies is detailed in the text of Section 5 for each constitutive model. | (10), (37), (38), and (48) |  $\phi$ and $\mathbf{u}$ | x | x | x |
III | [7] | (10), (40), (42), and (43) |  $J_{f}$, $\mu$ and $\mathbf{u}$ |  | x |  | x
IV | [13] | (10), (45a), (46)1, and (48) |  $\phi$ and $\mathbf{u}$ | x | x |  | x
V | [14] | (10), (45b), (46)2, and (48) |  $\phi$ and $\mathbf{u}$ | x | x |  | x

##  5 Representative prototype problems

Four representative prototype problems are introduced. First, we consider a
one-dimensional transient swelling of a hydrogel bar along the $Y$-direction.
The bar is fixed at $Y=0$ and free at $Y=0.01~{}[m]$ where $P_{Y}=0$. At this
latter end, the bar is exposed to a non-reactive solvent (see Figure 1a). The
deformation gradient takes the form

$\text{1D:}\quad\mathbf{F}=\text{diag}(1,\lambda,1)\,,$ (69)

occurring along the $Y$ direction. For perfectly incompressible models (i.e.,
_constitutive models I_ and _II_), it then results in
$\text{det}(\mathbf{F})=\lambda=J_{f}=\phi^{-1}$, that is the total stretch is
equal to the swelling volume change and inversely proportional to polymer
volume fraction. Furthermore, from the equilibrium condition (10), it follows
that $\partial\mathbf{P}/\partial Y=0$, which, after considering traction free
boundary conditions, leads to $P_{Y}=0$ and, in consequence, $\sigma_{y}=0$.
This problem has been studied previously by both linear and nonlinear
theories, e.g., [6, 7]. Due to the simplicity of its numerical settings, this
case study will serve as a reference, providing an estimate of the correct
order of magnitude of the quantities of interest, i.e., $\lambda$, $\mu$,
$P_{X}$ or $\sigma_{x}$, in more complex case studies.

As a second example, we investigate the transient swelling of a constrained
hydrogel slab in a two-dimensional setting in-plane strain (in the $(X,Y)$
plane). In this example, the hydrogel block is placed in a rigid container
with frictionless walls, and the deformation in the $X$ direction is
constrained. Only the upper part of the hydrogel is exposed to a non-reactive
solvent (see Figure 1b). This example has been previously considered by, e.g.,
[12] and [34], and represents the 2D counterpart of the previously introduced
1D example. However, it is noteworthy that the one-dimensional problem is
numerically solved for a single scalar field (representing either $\lambda$ or
$\phi$ depending on the constitutive model), and the other quantities of
interest are computed in the post-processing stage. On the other hand, the
numerical solution of the two-dimensional problem is obtained by considering
the complete sets of unknowns, that is, a scalar field describing the fluid
content and the vector displacement field $\mathbf{u}$. Then, quantities of
interest (e.g., the chemical potential $\mu$) are computed in the post-
processing stage.

As a third example, we consider the transient free-swelling in a two-
dimensional setting (in the $(X,Y)$ plane) of a polymer gel with an initially
square cross-section. A free hydrogel block is immersed in a non-reactive
solvent. Due to the symmetry of the deformation, only a quarter of the whole
model needs to be considered (see Figure 1c). A similar setup to this example
can be found in, e,g., [13] and [34]. At the steady state, the deformation
gradient is uniform within the domain, reading in the Cartesian representation
as:

$\text{Free
swelling:}\qquad\mathbf{F}_{\infty}^{2D}=\text{diag}\left(\lambda_{\infty}^{2D},\lambda_{\infty}^{2D},1\right)\,,$
(70)

with $\lambda_{\infty}^{2D}$ referring to the (constant) 2D final stretch at
the steady state.

The fourth and last example corresponds to the extension of the two-
dimensional block example into three-dimensions, namely, a cube is immersed in
a non-reactive solvent to swell due to the solvent absorption freely (see
Figure 1d). In the steady state, the deformation gradient results:

$\text{Free
swelling:}\qquad\mathbf{F}_{\infty}^{3D}=\lambda_{\infty}^{3D}\mathbf{I}$ (71)

with $\lambda_{\infty}^{3D}$ referring to the (constant) 3D final stretch at
the steady state. This simulation resembles that presented, e.g., by [38] and
will be only performed considering the constitutive model I.

We estimate the convergence order with a well-known heuristic formula. Let us
denote the errors by $E_{h}$, $E_{h/2}$, and $E_{h/4}$, where $h$ is the mesh
size parameter as before. Under the assumption that our discretized problem
has a convergence order of $p$, then each error should be roughly $(1/2)^{p}$
times the previous error. Therefore, we can estimate $p$ by taking the
logarithm base $2$ of the error ratios as:

$p=\dfrac{1}{\log(2)}\log\left(\left|\dfrac{E_{h}-E_{h/2}}{E_{h/2}-E_{h/4}}\right|\right).$
(72)

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=\\{0,T\\}$. 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=\\{0,T\\}$. 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.

###  5.1 Constitutive model I

One-dimensional transient swelling. Here, we follow the ideas presented by
[34]. Recalling equations (34) and (35) with equation (69), the balance of
linear momentum (10) reduces to:

$P_{Y}=-\dfrac{\mu}{\Omega}+\dfrac{k_{B}T}{\Omega}\left[\ln\left(1-\dfrac{1}{\lambda}\right)+\dfrac{1}{\lambda}+\dfrac{\chi}{(\lambda)^{2}}\right]+G_{0}\left(\lambda-\dfrac{1}{\lambda}\right)=0.$
(73)

Additionally, the fluid balance equation (11) in a one-dimensional setting
yields

$\begin{cases}\dfrac{\partial_{t}\lambda}{\Omega}=\dfrac{\partial}{\partial
Y}\left(\dfrac{D}{\Omega k_{B}T}(\lambda)^{-1}\dfrac{\partial\mu}{\partial
Y}\right),&~{}\text{in}~{}\mathcal{B}_{R},\\\
\lambda=\bar{\lambda},&~{}\text{at}~{}Y=0.01,\\\
\dfrac{\partial\lambda}{\partial Y}=0,&~{}\text{at}~{}Y=0.0,\\\
\lambda(t=0)=\lambda_{0},&~{}\text{in}~{}\mathcal{B}_{R}.\end{cases}$ (74)

Differentiating equation (73) with respect to $Y$ yields

$\dfrac{\partial\mu}{\partial
Y}=k_{B}T\left[\dfrac{1}{\lambda(\lambda-1)}-\dfrac{1}{(\lambda)^{2}}-\dfrac{2\chi}{(\lambda)^{3}}\right]\dfrac{\partial\lambda}{\partial
Y}+\Omega
G_{0}\left(1+\dfrac{1}{(\lambda)^{2}}\right)\dfrac{\partial\lambda}{\partial
Y}\,.$ (75)

Substituting equation (75) into equation (74) leads to a nonlinear partial
differential equation with respect to $\lambda(Y,t)$.

The weak discretized form of equation (74) reads

$\int_{\mathcal{B}_{R}}\left(\dfrac{\lambda_{h}^{n}-\lambda_{h}^{n-1}}{\Delta
t}q_{h}+\dfrac{D}{\Omega
k_{B}T}(\lambda_{h}^{n})^{-1}\dfrac{\partial\mu_{h}^{n}}{\partial
Y}\dfrac{\partial q_{h}}{\partial Y}\right)dY=0,~{}\forall q_{h}\in Q_{h}\,,$
(76)

where $\partial\mu_{h}^{n}/\partial Y$ is given by equation (75) evaluated at
$\lambda_{h}^{n}$. After solving equation (76), the chemical potential can be
obtained from (73) and the stress component $\sigma_{x}$ in the direction
orthogonal to swelling as $\sigma_{x}=p_{\mu}(\lambda,\mu)$, with $p_{\mu}$
from equation (35).

The numerical results after solving equation (76) with the FEM are presented
in Figure 2 (black dots). A one-dimensional mesh is created to discretize the
hydrogel bar, with the number of elements and time steps defined from a
convergence study on $\lambda(Y,t)$ (results not shown). In the end, the
reference solution is obtained with $200$ elements and $50$ time steps. The
boundary condition at $Y=0.01~{}[m]$ is obtained after solving equation (73)
for $\lambda$ with $\mu=0$, which yields
$\lambda(Y=0.01,t)=\bar{\lambda}=1.498$. The initial condition is defined as
$\lambda_{0}=1.0$. It is worth mentioning that the obtained results accurately
capture the ones presented by [34] in Figure 3 therein for the same problem
setup.

Two-dimensional constrained hydrogel slab. Constitutive model I is now solved
for the constrained hydrogel slab considered in Figure 1b.

Figure 2 presents the comparison between the solution previously obtained for
the one-dimensional bar (black dots) and the two-dimensional constrained slab
examples (different colored lines). For the two-dimensional case, we show
results from different simulations with increasing number of elements and time
steps. The deformed hydrogel at $t=10.0~{}[s]$ is displayed in Figure 2a, and
a time step equal to $0.2$ ($50$ time steps) was used to get Figure 2a. It is
observed from Figures 2b \- d that differences in numerical solution are
rather small for mesh densities higher than $50$. A zoom-in is included in
Figure 2b to better distinguish between the different mesh densities. The
discrepancy between the one- and two-dimensional cases can be assessed from
Figure 2, which shows the effect of approximating the time derivative in the
fluid concentration balance through the displacement and increasing the
problem’s dimension.

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\times
10^{-5}~{}[m^{2}s^{-1}]$.

From Figure 2c, it is observed that the chemical potential’s rate increases
fast in the beginning when the solvent starts entering the gel and becomes
slower when it tends to the steady state. Stress $\sigma_{x}$ follows a
similar pattern as observed in Figure 2d. This behavior can be understood as
follows. The gradient of $\mu$ is relatively large close to the top surface,
at $Y=0.01~{}[m]$, and the compressive stress in the interior part of the
hydrogel is the smallest, which are both helpful for the diffusion of the
solvent content. The hydrogel’s network is relaxed and allows easy solvent
absorption. The gradient of $\mu$ becomes smaller and $\sigma_{x}$ larger as
the solvent’s concentration increases in the hydrogel. This prevents the
solvent from penetrating the hydrogel further and slows the diffusion process.
There is less space within the hydrogel network, and it starts to saturate.
Therefore, each quantity approaches the corresponding steady solution at a
decreasing rate.

Next, a computational convergence analysis is performed to investigate the
robustness and computational cost of the monolithic approach. We focus on
investigating the effect of mesh density and the time step size on the
behavior of the implemented numerical algorithm. We aim to understand how
these parameters influence the performance of the algorithm in solving the
nonlinear system of equations associated with the FEM discretization.

First, we conduct performance studies of the Newton solver for each time step
to assess its efficiency in solving the nonlinear equations. We observe
whether the mesh density and time step size influence the number of Newton
iterations, determining if it remains consistent throughout the simulation.
Figure 3 depicts time versus the number of Newton iterations, where Figure 3a
showcases variations in mesh density and Figure 3b displays variations in time
step size. From these findings, it is clear that at the beginning of the
simulation, the number of Newton iterations is highest but subsequently
decreases until it reaches a steady value of three iterations.

Additionally, we examine the Newton algorithm’s convergence behavior at five
different time points during the simulation. We aim to confirm whether the
Newton iteration has quadratic convergence in error. We also investigate
whether changes in mesh density or time steps impact the Newton iteration’s
convergence. Figures 4 and 5 demonstrate the convergence behavior, considering
both absolute and relative errors. Figure 4 presents the results for various
mesh density values, while Figure 5 displays the results for different time
step sizes. Figures 4 and 5 reveal that the error decays very quickly for all
the cases, except for the first iteration. But, the Newton iteration also
converges in only $8$ iterations in this case. These results support the
efficiency and robustness of the Newton algorithm. These findings are relevant
as they assure the reliability of the numerical solution for the coupled
problem solved using a monolithic approach.

After establishing that the Newton solver is reliable, we test whether using
Taylor-Hood elements and the Euler method leads to the expected convergence in
space and time. Specifically, we measure the $L_{2}$ error for various mesh
densities and time step sizes with respect to the highest fidelity solution
and the values of both displacement and chemical potential at the center of
the top face of the two-dimensional hydrogel slab.

The convergence analysis is detailed in Tables 2 and 3, where the results
demonstrate a second-order convergence for spatial discretization and a first-
order convergence for temporal discretization. These findings are in
accordance with the theoretically predicted orders from the FEM and Euler
discretization scheme.

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

Table 2: Spatial discretization convergence analysis for the two-dimensional constrained hydrogel slab example considering constitutive model I. Level |  | Mesh density
---
$N_{h}$
| Time steps
---
$N_{k}$
$h$ | Elements | DoFs |  $\mu$(0.5.1.0) | L2 error for $\mu$ |  $\mathbf{u}$(1.0.1.0) | L2 error for $\mathbf{u}$
1 | 25 | 50 | 0.04 | 1250 | 5878 | -0.44017195 | 1.3896e-4 | 0.33459041 | 3.5884e-5
2 | 50 | 50 | 0.02 | 5000 | 23003 | -0.43999694 | 3.4425e-5 | 0.33461516 | 9.0131e-6
3 | 100 | 50 | 0.01 | 20000 | 91003 | -0.43994541 | 7.1092e-6 | 0.33462156 | 1.9465e-6
4 | 200 | 50 | 0.005 | 80000 | 362003 | -0.43993062 | – | 0.33462318 | –
conv. order |  |  |  |  |  | 1.80 | 1.93 | 1.98 | 1.93

Table 3: Time discretization convergence analysis for the two-dimensional constrained hydrogel slab example considering constitutive model I. Level |  | Mesh density
---
$N_{h}$
| Time steps
---
$N_{k}$
$k$ |  $\mu$(0.5.1.0) | L2 error for $\mu$ |  $\mathbf{u}$(1.0.1.0) | L2 error for $\mathbf{u}$
1 | 25 | 25 | 0.4 | -0.44915384 | 1.2353e-2 | 0.33195688 | 3.1447e-3
2 | 25 | 50 | 0.2 | -0.44017195 | 5.5713e-3 | 0.33459041 | 1.4267e-3
3 | 25 | 100 | 0.1 | -0.43534375 | 1.9333e-3 | 0.33601343 | 4.9668e-4
4 | 25 | 200 | 0.05 | -0.43277477 | – | 0.33677258 | –
conv. order |  |  |  | 0.91 | 0.90 | 0.91 | 0.88

Finally, it is noteworthy that we tested both monolithic and partitioned
approaches for the solution of the _constitutive model I_ by [34]. However,
only the monolithic implementation produced physically consistent results. The
reason seems to be due to the reformulation of the concentration-time
derivative. The time-dependent concentration term is expressed in terms of
$\text{Grad}(\partial_{t}\mathbf{u}):\mathbf{F}^{-T}$ in equation (25). Hence,
such a term is highly nonlinear in terms of displacements, and its accuracy is
highly affected by the adopted spatial discretization. Moreover, if the
deformation process is very slow, i.e., $\partial_{t}\mathbf{u}\approx 0$,
which is the case here, the time-dependent term in the balance of fluid
concentration may become negligible, and this equation becomes quasi-static.
These issues cause some numerical difficulties in the staggered approach to
capture the transient behavior of the coupled problem.

Two-dimensional free swelling of a square block. A square block is immersed in
a solvent with a reference chemical potential $\mu^{0}=0$, as illustrated in
Figure 1c. Recalling equations (34) and (35), and following the same procedure
as in [34], the theoretical value of the steady state stretching results
$\lambda_{\infty}^{2D}=1.35$ (see equation (70)).

Figure 6: _Constitutive model I_ : two-dimensional hydrogel block of an
initially square cross-section immersed in a non-reactive solvent at different
simulation times. a. Deformed two-dimensional block at three different time
steps $t=1.0~{}[s]$, $t=5.0~{}[s]$, and $t=10.0~{}[s]$. b. Stretch due to
swelling $\lambda$ at different times across $X=0.005$. c. Transient evolution
of $\lambda$ measured at different corners of the two-dimensional block.
Simulation parameters:
$G_{0}=10~{}[MPa],~{}~{}\chi=0.2~{}[--],~{}~{}D=5.0\times
10^{-5}~{}[m^{2}s^{-1}]$.

Figure 6 displays the simulation results of the transient diffusion-
deformation process for the two-dimensional square block. The thick gray line
in each subfigure in Figure 6a indicates the reference body. From Figure 6a,
it is evidenced that the initially square block gets distorted at the
beginning of the deformation process while swelling. The origin is the
pronounced $\lambda$ gradient in the early stage of the transient behavior as
evidenced in Figure 6b. This distortion vanishes as time progresses and all
corners reach a similar stretching value as Figure 6c shows. The two-
dimensional blocks exhibit a diffusion-deformation process that aligns with
the one observed in the two-dimensional slab. The final stretch reaches the
steady state at the previously computed theoretical value
$\lambda_{\infty}^{2D}=1.35$.

Three-dimensional free swelling cubic block. The last example corresponds to
the extension of the previous example from two to three dimensions as
illustrated in Figure 1d. The steady-state stretching is estimated to be
$\lambda_{\infty}^{3D}=1.45$, which can be used to verify the simulation
results at the steady state.

Figure 7: _Constitutive model I_ : three-dimensional hydrogel cubic block
immersed in a non-reactive solvent at different simulation times. a. Deformed
three-dimensional cube projected to planes $XY$ and $YZ$ at two different time
steps $t=0.1$ and $t=1.0$. b. Stretch due to swelling $\lambda$ along $Y$ at
different times across the plane $XZ=0.01~{}[m]$. c. Transient evolution of
$\lambda$ measured at different corners of the three-dimensional cubic block.
Simulation parameters:
$G_{0}=10~{}[MPa],~{}~{}\chi=0.2~{}[--],~{}~{}D=7.5\times
10^{-5}~{}[m^{2}s^{-1}]$.

Figure 7 presents the simulation results of the transient diffusion-
deformation process for the three-dimensional block. Compared to the two-
dimensional block, the cube is less distorted, as seen in the $XY$ and $YZ$
plane projections in Figure 7a. This can be explained by the increment in the
diffusivity coefficient $D$ from $2.5\times 10^{-5}$ to $5.0\times
10^{-5}~{}[m^{2}s^{-1}]$, where the latter is the minimum value that leads to
the Newton algorithm iteration convergence. The diffusivity effect is observed
in Figure 7b when compared to Figure 6b. The former displays a faster
evolution of the stretch over time. Figure 7c shows the stretch evolution at
different corners. It is observed that each observed corner presents a
different stretch at the beginning of the simulation, but it fades over time.
Consequently, as illustrated in Figure 7a, the three-dimensional block
recovers its cube shape.

It is noted that the two-dimensional and three-dimensional problems are
addressed without altering the diffusion-deformation model or its constitutive
equations. Nonetheless, adjusting the diffusivity coefficient value was
necessary. This is because a pronounced stretch gradient results in high-
stress values and excessive distortion of the elements used for domain
discretization as more degrees of freedom are added to the problem. This
stability issue is commonly encountered in diffusion-deformation studies of
hydrogels with non-reactive solvent absorption (refer to [7] for an in-depth
discussion). Although beyond the scope of this study, stabilization methods
can be employed to overcome this issue (see [29, 4] for more information). Our
results provide a foundation for testing some of the approaches.

###  5.2 Constitutive model II

The second model under consideration is composed of the balance equations (10)
and (11) together with constitutive equations (37) and (38).

One-dimensional transient swelling. Here, we follow the ideas presented by
[12]. By considering the 1D deformation gradient in equation (69) under the
incompressibility condition, the stress equation (37) yields:

$P_{Y}=G_{0}\phi^{-1}-P\phi=0\,.$ (77)

Hence, it results:

$P=G_{0}\phi^{-2}\,.$ (78)

Notice that this expression for $P$ is case specific and equation (78) holds
true exactly only in the present 1D case. After some manipulations, the
balance of fluid concentration given in equation (48) can be rewritten as

$\begin{cases}\partial_{t}\phi=\dfrac{D}{k_{B}T}\left[\dfrac{\partial\phi}{\partial
Y}\dfrac{\partial\mu}{\partial
Y}-\phi^{2}\left(\dfrac{1-\phi}{\phi}\right)\dfrac{\partial^{2}\mu}{\partial
Y^{2}}\right],&~{}\text{in}~{}\mathcal{B}_{R},\\\
\phi=\bar{\phi},&~{}\text{at}~{}Y=0.01,\\\ \dfrac{\partial\phi}{\partial
Y}=0,&~{}\text{at}~{}Y=0.0,\\\
\phi(t=0)=\phi_{0},&~{}\text{in}~{}\mathcal{B}_{R},\end{cases}$ (79)

with

$\dfrac{\partial\mu}{\partial
Y}=k_{B}T\left[-\dfrac{1}{1-\phi}+1+2\chi\phi\right]\dfrac{\partial\phi}{\partial
Y}-\Omega G_{0}\left(1+\dfrac{1}{\phi^{2}}\right)\dfrac{\partial\phi}{\partial
Y},$ (80)

obtained after differentiating equation (38) with respect to $Y$.

The discretized weak form of equation (79) reads

$\int_{\mathcal{B}_{R}}\left(\dfrac{\phi_{h}^{n}-\phi_{h}^{n-1}}{\Delta
t}q_{h}-\dfrac{D}{k_{B}T}\left[\dfrac{\partial\phi_{h}^{n}}{\partial
Y}\dfrac{\partial\mu_{h}^{n}}{\partial
Y}q_{h}-(\phi_{h}^{n})^{2}\left(\dfrac{1-\phi_{h}^{n}}{\phi_{h}^{n}}\right)\dfrac{\partial\mu_{h}^{n}}{\partial
Y}\dfrac{\partial q_{h}}{\partial Y}\right]\right)dY=0,~{}\forall q_{h}\in
Q_{h},$ (81)

where $\partial\mu_{h}^{n}/\partial Y$ is given by equation (80) evaluated at
$\phi_{h}^{n}$. Again, equation (81) can be directly implemented in FEniCs
with the corresponding boundary and initial conditions and solved using the
FEM. Solution of equation (81) serves as a reference solution for the
diffusion-deformation problem adopting constitutive model II. The stress can
be computed replacing $P$ into equation (37), yielding:

$\sigma_{x}=J^{-1}P_{X}(\phi,\mathbf{u})=G_{0}\phi\left[1-\phi^{-2}\right],$
(82)

which for $\phi\in[0,1]$ results in a compressive stress state. Whereas the
chemical potential in equation (38) becomes

$\mu(\phi)=\mu^{0}+k_{B}T\left[\ln\left(1-\phi\right)+\phi+\chi\phi^{2}\right]+\Omega
G_{0}\left[\phi^{-1}-\phi\right],$ (83)

after the definition of $P$.

The numerical solution of equation (81) is presented in Figure 8 (black dots).
A one-dimensional mesh is created to discretize the hydrogel bar. As a result
of a convergence study, the final number of elements equals $200$, and $50$
time steps are chosen ($\Delta t=0.02$). The boundary condition at
$Y=0.01~{}[m]$ is determined from equation (38), which gives
$\phi(Y=0.01,t)=0.2$. The initial condition is defined as $\phi_{0}=0.75$,
which corresponds to a pre-swollen state conveniently chosen to alleviate the
numerics and allow for a comparison between the one and two-dimensional
models. In the original work by [12], this one-dimensional problem is
numerically approximated using a finite difference method in space. The
simulation results here obtained with the FEM can be quantitatively compared
to those reported by [12] for a similar example (see Figures 3, 5, and 7 in
[12]).

Two-dimensional constrained hydrogel slab. Constitutive model II is solved for
the constrained hydrogel slab considered in Figure 1b. Initial and boundary
conditions for $\phi$ are kept analogous to the one-dimensional case. Hence,
since we have a dominant diffusion along the $Y$-direction also in this case
and an analogous stress state, we approximate the Lagrange multiplier $P$ by
employing the one obtained in the 1D case, i.e. with equation (78). Under the
limitations of such approximation, the stress for the two-dimensional case
hence reads:

$\sigma_{x}=J^{-1}P_{X}(\phi,\mathbf{u})=J^{-1}G_{0}\left[\mathbf{b}-\phi^{-2}\mathbf{I}\right],$
(84)

after replacing equation (78) into equation (37). Since the same approximation
for $P$ would be inaccurate in a free swelling condition, we will not face
these case studies for constitutive model II.

Figure 8 presents the comparison between the one-dimensional bar (black dots)
and two-dimensional constrained slab examples (different colored lines). The
one-dimensional simulation is taken as the reference solution to the problem.
The two-dimensional slab problem is solved for different mesh densities and
different number of time steps. It was concluded that a mesh density of
$N_{h}=30$ and $25$ time steps are enough to produce accurate results, so the
following simulation results are produced using this simulation setup. The
deformed hydrogel at $t=10.0[s]$ is displayed in Figure 8a. Figures 8b \- d
show the time evolution of $\phi$, $\mu$ and $\sigma_{x}$. It is observed that
the numerical solutions for the one and two-dimensional problems coincide.
This is not surprising since the large deformation displayed by the hydrogel
slab makes it behave like a one-dimensional bar.

Figure 8: _Constitutive model II_ : 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=1.0$. b. Polymer
volume fraction $\phi$. c. Chemical potential $\mu$ normalized by $k_{B}T$. d.
Cauchy compressive stress $\bm{\sigma}_{X}$. Simulation parameters:
$G_{0}=1~{}[MPa],~{}~{}\chi=0.2~{}[--],~{}~{}D=2.5\times
10^{-5}~{}[m^{2}s^{-1}]$.

In Figure 8c, it is observed that the chemical potential’s rate increases fast
in the beginning when the solvent starts entering the gel and becomes slower
when it tends to the steady state. The pronounced gradient of $\mu$ along the
spatial dimension in the beginning and a relaxed polymer network facilitates
solvent absorption. As the gradient of $\mu$ flattens out and the compressive
stress increases, solvent diffusion slows down.

[6] and [14] proposed to solve an additional nonlinear equation either for
$J_{f}$ or $\phi$ at the Gauss integration point to fully determine the time
derivative term and then solve the balance of fluid concentration for the
chemical potential. If equations (43) and (48) are solved for $\mu$, either
equation (42), (38), or (46) must be solved as an implicit nonlinear equation
at each Gauss integration point for either $J_{f}$ or $\phi$, depending of the
adopted constitutive model, as suggested by [7] and [12, 14].

In [14], the coupled problem was solved using the commercial software Abaqus,
via an UMAT routine. The coupled problem was solved for $\phi$, $\mu$, and
$\mathbf{u}$. In particular, $\phi$ was defined as a local variable, and the
nonlinear equation (38) was solved at each Gauss integration point. This
approach allowed to fully determine the time derivative in equation (48)2 at
each time instance. Thus, a static problem for $\mu$ and $\mathbf{u}$ is
solved also at each time step. A similar idea was adopted by [7].

We decided to test an alternative approach. That is, to find a suitable
expression for $\text{Grad}(\mu)$ in equation (48) from the constitutive
equation (38). Then, solve the couple problem only for $\phi$ and
$\mathbf{u}$. Such an approach allows us to avoid the computation of $\phi$ at
each integration point. By following this approach, we noticed that the
solution becomes less prompt to numerical instabilities, and monolithic and
staggered formulations can be adopted to solve the resulting coupled problem.
It is worth remarking that without the automatic differentiation capabilities
of FEniCs, the computation of $\text{Grad}(\mu)$ can be a very tedious task,
and it is not difficult to imagine that was the reason why [14] opted for
their approach.

We conduct an investigation to assess the staggered approach’s behavior. We
verify the number of staggered iterations between the displacement and polymer
volume fraction sub-systems and, as for the monolithic case, we examine the
effects of mesh density and time step size on the algorithm’s performance,
focusing particularly on the Newton solver’s efficiency and convergence.

For the staggered solution, we modify the strategy slightly. Instead of
computing the total number of Newton iterations required for convergence, we
compute the average number of Newton iterations over the inner iteration loop
for each time increment for the displacement and polymer volume fraction sub-
systems. This approach provides us with insights into the staggered solution’s
convergence behavior and whether the number of staggered and Newton iterations
remains stable throughout the simulation time. It allows us to gauge the
efficiency and stability of the staggered approach, and compare it to the
monolithic approach.

Figure 9 illustrates the time versus the staggered and Newton iterations for
different mesh densities and time step sizes. Figures 9a and 9d show how many
staggered iterations are necessary at each time step to reach an $L_{2}$ error
lower than $10^{-10}$ both for $\mathbf{u}$ and $\phi$, respectively. It is
observed that the number of iterations is never higher than three, independent
of the mesh density and time step size. This highlights the effectiveness of
the staggered approach.

The convergence patterns for the Newton iterations and the impact of
variations in mesh density and time step size were consistent with those
observed in the monolithic approach. Figures 9b, e focuses on the Newton
iteration for the displacement sub-system, whereas Figure 9c,f reports the
Newton iterations along time for the polymer volume fraction sub-system, for
different mesh densities and time step sizes, respectively. Results in Figure
9 reveal a similar trend to the monolithic case regarding the Newton
iterations at each time increment.

Additionally, we explored the Newton algorithm’s convergence behavior at five
distinct simulation moments, examining the quadratic convergence in error and
the influence of mesh density and time step size. Figure 10 shows the error
decay with respect to the Newton iterations. Sub-figures 10a, b, c, d
correspond to the displacement, while sub-figures 10e, f, g, h refer to the
polymer volume fraction, both for different mesh densities. The results are
analogous to the monolithic approach, where the absolute and relative errors
decay fast as the Newton iteration increases, demonstrating the staggered
approach’s robustness.

The observed similarities in the convergence behaviors between the monolithic
and staggered approaches further underline the robustness and reliability of
the numerical solutions, confirming the applicability of both methods to the
coupled problem at hand.

Concerning convergence analyses for the discretization itself, for the Taylor-
Hood elements and Euler method, our findings for spatial and temporal
discretization yielded a second-order and first-order convergence,
respectively. These outcomes are summarized in Tables 4 and 5, reinforcing the
staggered approach’s validity in solving the coupled diffusion-deformation
problem.

Figure 9: _Constitutive model II_ : staggered algorithm and Newton iterations
along the time steps. For different mesh densities ($N_{h}$): a. for the
staggered scheme, b. for displacement, and c. for the polymer volume fraction.
For different time step sizes ($N_{k}$): d. for the staggered scheme, e. for
displacement, and f. for the polymer volume fraction. Figure 10: _Constitutive
model II_ : convergence analysis of the Newton solver for different mesh
densities ($N_{h}$). For the displacement: a. $N_{h}=30$. b. $N_{h}=60$. c.
$N_{h}=120$. d. $N_{h}=240$. For the polymer volume fraction: e. $N_{k}=25$.
f. $N_{h}=50$. g. $N_{h}=100$. h. $N_{h}=200$.

Table 4: Spatial discretization convergence analysis for the two-dimensional constrained hydrogel slab example considering constitutive model II. Level |  | Mesh density
---
$N_{h}$
| Time steps
---
$N_{k}$
$h$ | Elements | DoFs |  $\phi$(0.5.1.0) | L2 error for $\phi$ |  $\mathbf{u}$(1.0.1.0) | L2 error for $\mathbf{u}$
1 | 30 | 40 | 0.033 | 1800 | 8403 | 0.31318141 | 2.2590e-5 | 2.83246074 | 1.4800e-4
2 | 60 | 40 | 0.016 | 7200 | 33003 | 0.31309735 | 5.4560e-6 | 2.83288465 | 3.4709e-5
3 | 120 | 40 | 0.008 | 28800 | 130803 | 0.31307316 | 1.2288e-6 | 2.83299197 | 9.0328e-6
4 | 240 | 40 | 0.004 | 115200 | 520803 | 0.31306639 | – | 2.83301882 | –
conv. order |  |  |  |  |  | 1.83 | 2.02 | 2.00 | 2.14

Table 5: Time discretization convergence analysis for the two-dimensional constrained hydrogel slab example considering constitutive model II. Level |  | Mesh density
---
$N_{h}$
| Time steps
---
$N_{k}$
$k$ |  $\phi$(0.5.1.0) | L2 error for $\phi$ |  $\mathbf{u}$(1.0.1.0) | L2 error for $\mathbf{u}$
1 | 30 | 20 | 0.5 | 0.31996593 | 8.2017e-3 | 2.72595105 | 9.6975e-02
2 | 30 | 40 | 0.25 | 0.31309911 | 3.7086e-3 | 2.80367321 | 4.3200e-02
3 | 30 | 80 | 0.125 | 0.30939811 | 1.2936e-3 | 2.84445820 | 1.4871e-02
4 | 30 | 160 | 0.0625 | 0.30741375 | – | 2.86582984 | –
conv. order |  |  |  | 0.90 | 0.90 | 0.93 | 0.92

###  5.3 Constitutive model III

The third model under consideration comprises the balance equations (10) and
(11) together with constitutive equations (40) and (42).

One-dimensional transient swelling. The basic ingredients to describe the
hydrogel deformation considering an energetic constraint were presented by [6,
7]. By considering the 1D deformation gradient in equation (69), the balance
of fluid concentration reads

$\begin{cases}\gamma(\lambda)\partial_{t}\lambda=D\dfrac{\partial}{\partial
Y}\left[\xi(\lambda)\dfrac{\partial\lambda}{\partial
Y}\right],&~{}\text{in}~{}\mathcal{B}_{R},\\\
\lambda=\bar{\lambda},&~{}\text{at}~{}Y=0.01,\\\
\dfrac{\partial\lambda}{\partial Y}=0,&~{}\text{at}~{}Y=0.0,\\\
\lambda(t=0)=\lambda_{0},&~{}\text{in}~{}\mathcal{B}_{R},\end{cases}$ (85)

with

$\gamma(\lambda)=\lambda_{0}+\dfrac{Nk_{B}T}{K\lambda_{0}}\left(1+\dfrac{1}{\lambda^{2}}\right),$
(86)

and

$\xi(\lambda)=\left(\dfrac{1}{(\lambda_{0}^{2}\lambda)^{2}}-\dfrac{2\chi(\lambda_{0}^{2}\lambda-1)}{(\lambda_{0}^{2}\lambda)^{3}}\right)\left(\frac{Nk_{B}T(1+\lambda^{2})}{K(\lambda_{0}\lambda^{2})^{2}}+\frac{\lambda_{0}^{2}}{\lambda^{2}}\right)+\dfrac{N\Omega(\lambda_{0}^{2}\lambda-1)(1+\lambda^{2})}{(\lambda_{0}\lambda^{2})^{2}},$
(87)

where $\lambda_{0}$ refers to the initial swelling ratio of the gel.

The discretized weak form of equation (85) reads

$\int_{\mathcal{B}_{R}}\left(\gamma(\lambda_{h}^{n})\dfrac{\lambda_{h}^{n}-\lambda_{h}^{n-1}}{\Delta
t}q_{h}+D\xi(\lambda_{h}^{n})\dfrac{\partial\lambda_{h}^{n}}{\partial
Y}\dfrac{\partial q_{h}}{\partial Y}\right)dY=0,~{}\forall q_{h}\in Q_{h},$
(88)

where $\gamma(\lambda_{h}^{n})$ and $\xi(\lambda_{h}^{n})$ are given by
equations (86) and (87) evaluated at $\lambda_{h}^{n}$, respectively. Equation
(88) can be solved using the FEM with FEniCS. The postprocessing quantities
$\mu$ and $\sigma_{x}$ results:

$\mu(\lambda_{h})=k_{B}T\left[\ln\left(\dfrac{\lambda_{0}^{2}\lambda_{h}-1}{\lambda_{0}^{2}\lambda_{h}}\right)+\dfrac{\lambda_{0}^{2}\lambda_{h}+\chi}{(\lambda_{0}^{2}\lambda_{h})^{2}}+\dfrac{N\Omega}{\lambda_{0}^{2}}\left(\lambda_{h}-\dfrac{1}{\lambda_{h}}\right)\right],$
(89)

and

$\sigma_{x}(\lambda_{h})=Nk_{B}T\left[\left(\lambda_{0}-\dfrac{1}{\lambda_{0}}\right)-K\left(\dfrac{\lambda_{h}}{\lambda_{0}}\right)\left(\lambda_{h}-\dfrac{1}{\lambda_{h}}\right)\right],$
(90)

respectively. Notice that $J_{f}$ in equation (42) was replaced in by
$\lambda_{0}^{2}\lambda_{h}$ to get equation (89). More details about the
derivation of equations (89) and (90) can be found in appendix A of [7].

Two-dimensional constrained hydrogel slab. Example’s geometry corresponds to
that illustrated in Figure 1b. Initial and boundary conditions correspond to
those defined for the one-dimensional case.

Figure 11 presents the comparison between the one-dimensional bar (black dots)
and two-dimensional constrained slab examples (different colored lines). The
two-dimensional slab problem is solved for a mesh density $N_{h}=25$ and $10$
time steps ($\Delta t=1.0~{}[s]$). This setup yields an acceptable numerical
accuracy for the two-dimensional configuration. The deformed hydrogel at
$t=10.0~{}[s]$ is displayed in Figure 11a. Figures 11b \- d show the time
evolution of $\lambda=J_{f}$, $\mu$ and $\sigma_{x}$. It is observed that the
numerical solutions for the one and two-dimensional problems are on the same
order of magnitude and get closer as the two-dimensional slab becomes larger
and progressively resembles a one-dimensional domain. The same pattern was
observed for the constitutive model I. However, the hydrogel experiences a
lower level of stretch when considering the two-dimensional setup. This is the
effect of determining first $\lambda$ through the nonlinear equation and then
$\mu$ from the static PDE. The boundary condition can only be imposed on
$\mu$, not $\lambda$ as in the one-dimensional case.

Figure 11: _Constitutive model III_ : one-dimensional bar (black dots) and
two-dimensional hydrogel constrained slab (colored lines) numerical solution
comparison at different simulation times. a. Deformed two-dimensional
constrained slab at $t=10.0~{}[m]$. 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}=Nk_{B}T\approx
41~{}[Pa],K=100~{}Nk_{B}T~{}[Pa],~{}~{}\chi=0.4~{}[--],~{}~{}D=2.5\times
10^{-6}~{}[m^{2}s^{-1}]$.

The transient behavior of the hydrogel follows the already observed evolution
for constitutive models I and II. At the beginning of the deformation, the
solvent enters the gel faster because of the high $\mu$ gradient and low
compressible stress. From Figure 11, it is noticed that the hydrogel slab
undergoes a rather large deformation, i.e., the final size is about $3$ times
the original one, despite the presence of the energetic constrain in the free
energy function. The small shear modulus $G_{0}$ value can justify this large
deformation.

Constitutive model III as considered in this work was formally presented by
[6] and solved using the FEM by [7]. The coupled problem was solved for
$c_{R}$, $\mu$, and $\mathbf{u}$. In particular, $c_{R}$ was defined as a
local variable, and the nonlinear equation (42) was solved at each Gauss
integration point to determine the time derivative in equation (43) at each
time instant. Next, $\mu$ and $\mathbf{u}$ are computed by the solution of a
static problem defined by equations (10) and (11) at each time step.

Here we tried to find a suitable expression for $\text{Grad}(\mu)$ in equation
(48) from the constitutive equation (42) using FEniCS automatic
differentiation tools and avoid the solution of the nonlinear equation at the
Gauss points for the two-dimensional example. Although we got some numerical
results, these appear highly inaccurate compared to the reference 1D outcomes.
Consequently, we only report the results obtained following [7] original
formulation.

It should be highlighted that among the constitutive models accounted for in
this work, constitutive model III numerical solution is the most fragile.
There is not much room to explore different values for the model parameters.
Small changes in their values would lead to the divergence of the Newton-
Raphson algorithm. It is also not possible to arbitrarily increase the mesh
density or make the time step smaller. This is not surprising since [7]
presented a detailed analysis of the induced instabilities in hydrogels in the
presence of geometrical constraints. These instabilities arise because the
exposure of a gel to a fluid not only leads to a large increase in volume but
also to a wave-like buckling pattern. Free surfaces expand due to the species
influx and are bonded to unswollen inner parts of the gel simultaneously. For
high osmotic pressure, this mechanism leads to buckling patterns, that have
extensively been analyzed within experimental setups (see, e.g., [22, 16, 31]
for more details on hydrogels instabilities).

Due to this numerical fragility, it was not possible to perform a convergence
analysis in the same manner as for the two previous constitutive models.
However, for the chosen simulation setup presented in this sub-section, both
the absolute and relative errors decay faster as the Newton iterations
increase, in a manner similar to that observed for constitutive model II (see
Figures 4 and 5).

Two-dimensional free swelling of a square block. Figure 1c illustrates the
considered example. 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 set equal to zero. The initial
condition for $J_{f}$ is defined as $J_{f_{0}}=1.4$. In this case, it was not
necessary to apply a time-ramping strategy.

Figure 12: _Constitutive model III_ : two-dimensional hydrogel block of an
initially square cross-section immersed in a non-reactive solvent at different
simulation times. a. Deformed two-dimensional block at three different time
steps $t=1.0~{}[s]$, $t=5.0~{}[s]$, and $t=10.0~{}[s]$. b. Polymer volume
fraction $\phi$ at different times across $X=0.005~{}[m]$. c. Transient
evolution of $\phi$ measured at different corners of the two-dimensional
block. Simulation parameters: $G_{0}=Nk_{B}T\approx
0.41~{}[MPa],K=100~{}Nk_{B}T~{}[MPa],~{}~{}\chi=0.4~{}[--],~{}~{}D=2.5\times
10^{-6}~{}[m^{2}s^{-1}]$.

The simulation results for the two-dimensional square block are presented in
Figure 12. The box plotted with a thick gray line in Figure 12a corresponds to
the reference body before undergoing deformation. In contrast to the
constitutive model I, the block tends to keep its square shape along the
simulation time. This can be attributed to the value of the diffusivity that
prevents high $\mu$ gradients as observed in Figure 12b. However, we chose
this simulation setup because it leads to numerical convergence. For
simulation setup producing more pronounced $\mu$ gradients, bigger distortion,
or larger deformations, suitable stabilization techniques are required (see,
e.g., [29, 4] for some approaches), which are out of the scope of the current
research. An important difference between constitutive model III and the
previous constitutive models is that $\lambda$ does not tend to a steady state
value as time simulation progresses. This can be appreciated from Figure 12c.

###  5.4 Constitutive models IV and V

The fourth and fifth models under consideration are composed of the balance
equations (10) and (11) together with constitutive equations (45) and (46) for
the stress and chemical potential, respectively.

To avoid redundancies, we only present results for the two-dimensional free
swelling of a square block with constitutive models IV and V. It is worth
recalling that the main difference between constitutive model II and IV - V is
the addition of an energetic component in the free energy function, thus
eliminating the need to introduce a Lagrange multiplier as in constitutive
model II.

Figure 13: _Constitutive model IV_ : Two-dimensional hydrogel block of an
initially square cross-section immersed in a non-reactive solvent at different
simulation times. a. Deformed two-dimensional block at three different time
steps $t=1.0~{}[s]$, $t=5.0~{}[s]$, and $t=10.0~{}[s]$. b. Polymer volume
fraction $\phi$ at different times across $X=0.005~{}[m]$. c. Transient
evolution of $\phi$ measured at different corners of the two-dimensional
block. Simulation parameters:
$G_{0}=1~{}[MPa],~{}~{}K=100~{}[MPa],~{}~{}\chi=0.2~{}[--],~{}~{}D=7.5\times
10^{-9}~{}[m^{2}s^{-1}],~{}~{}\alpha_{r}=10.0~{}[--]$. Figure 14:
_Constitutive model V_ : Two-dimensional hydrogel block of an initially square
cross-section immersed in a non-reactive solvent at different simulation
times. a. Deformed two-dimensional block at three different time steps
$t=1.0~{}[s]$, $t=5.0~{}[s]$, and $t=10.0~{}[s]$. b. Polymer volume fraction
$\phi$ at different times across $X=0.005~{}[m]$. c. Transient evolution of
$\phi$ measured at different corners of the two-dimensional block. Simulation
parameters:
$G_{0}=1~{}[MPa],~{}~{}K=100~{}[MPa],~{}~{}\chi=0.2~{}[--],~{}~{}D=7.5\times
10^{-9}~{}[m^{2}s^{-1}],~{}~{}\alpha_{r}=10.0~{}[--]$.

The simulation results obtained for constitutive models IV and V are displayed
in Figures 13 and 14, respectively. From these results, it becomes evident
that both constitutive models lead to a similar level of deformation and time
evolution of $\phi$.

##  6 A reference benchmark problem

In this section, the two-dimensional square block problem is studied again as
a unified benchmark problem for the diffusion-deformation of hydrogels. As a
matter of fact, previous results cannot be directly compared due to the
disparate value of the material parameters selected by authors in the original
papers (and adopted for the prototype problems presented in the previous
section). Our goal here is to have a unique reference simulation example and
show the differences in the response for each constitutive equation. We are
interested in testing how much the diffusion-deformation behavior is affected
by the different constitutive choices under the same simulation setup.

Hydrogels are a type of material that can vary significantly in their physical
properties due to factors such as the particular polymer used, the degree of
crosslinking, and the presence of any additives [9]. As such, the shear
modulus ($G_{0}$) and bulk modulus ($K$) can have a wide range of values.
Typically, the shear modulus of hydrogels can range from $1~{}[Pa]$ to
$1~{}[MPa]$, with softer, more water-rich gels being at the lower end of the
scale and more crosslinked or polymer-rich gels being at the higher end. Among
the studies considered in this research, we found very different values for
the shear modulus ranging from $0.1~{}[MPa]$ in [12, 13] to $10~{}[MPa]$ as
reported by [35].

Similarly, the bulk modulus, which is a measure of a material’s resistance to
uniform compression, can also vary widely for hydrogels. The bulk modulus is
typically larger than the shear modulus. For many hydrogels, the bulk modulus
can range from around $10~{}[kPa]$ to several $[MPa]$. However, these values
should only be used as a rough guide, as the specific values for any
particular hydrogel can vary significantly depending on its formulation and
preparation details. One needs to refer to specific experimental data for the
particular hydrogel of interest for a precise value. The most common value for
$K$ reported in the reviewed papers was $K=100~{}G_{0}$.

In the present study, motivated by state-of-the-art values, we choose
$G_{0}=1~{}[MPa]$ and $K=100~{}[MPa]$. The diffusion coefficient is used as
the free parameter to tune the simulation such that the results become
comparable and the numerical stability is guaranteed. All other parameters
remain the same as for the previously reported simulation results. Notice that
this simulation setup is similar to the one used to solve constitutive models
II, and IV - V and it was inspired by [14].

H

Figure 15: _Reference benchmark_ : two-dimensional hydrogel block of an
initially square cross-section immersed in a non-reactive solvent at different
simulation times. a. constitutive model I, b. constitutive model III, and c.
constitutive model IV. Common simulation parameters:
$G_{0}=1~{}[MPa],~{}~{}K=100~{}[MPa],~{}~{}\chi=0.2~{}[--]$.

Figure 15 shows the simulation results for the two-dimensional square block
problem. Results are reported for constitutive models I, III, and IV. The
reason for this selection is twofold: i. constitutive model II introduces an
additional field represented by the Lagrange multiplier, thus requiring ad hoc
numerical treatments outside of the scopes of present work; ii. the same
example was solved for constitutive model V with the same parameters and
results in Figure 15 can be directly compared with those plotted in Figure 14.

From Figure 15, first column, it is possible to observe that the level of
deformation achieved in each case by the end of the simulation is different
despite the similarity of the parameters. However, it is not on the
deformation or stretch where the main difference between the different models
becomes evident. It is on the transient evolution of the solvent as
appreciated in the second and third columns of the figure. Constitutive models
I and IV reach their steady state in about $4.0~{}[s]$ (Figure 15a, third
column) and $6.0~{}[s]$ (Figure 15c, third column), respectively, whereas
constitutive model III is far from reaching the steady state (Figure 15b,
third column).

It is important to note that the diffusivity values used in the benchmark
problem for constitutive models I and III are the lowest possible before
encountering numerical issues. For constitutive models III - IV, it is
feasible to set the diffusivity to $8.5\times 10^{-3}[m^{2}s^{-1}]$, matching
constitutive model I. However, given our simulation setup, this is an
exceptionally high value. Consequently, the system will reach steady in both
scenarios almost immediately, eliminating any transient behavior. In summary,
the mentioned diffusivity values differ significantly from each other. They
belong to different time scales. This difference can be attributed to the
distinct solution strategies adopted for solving the diffusion equation in
each case. We delved into the specifics of this when introducing each
constitutive model.

One takeaway drawn from this comparison study is that further refinement may
be required for the constitutive models such that closer predictions are
achieved. In particular, each model predicts a different level of deformation
and requires a fine-tuning of some parameters to ensure numerical convergence.
This can complicate the model’s experimental validation because it could lead
to different values of the material parameters and, for instance, ambiguity
regarding which values are correct and their interpretability. Nonetheless,
there is still room for improvement in selecting the energetic component of
the constitutive model, which can fix the issue. In the end, all constitutive
models capture the diffusion-deformation process in a reasonable way and offer
valuable insights into the coupled problem.

##  7 Conclusion

This study presented a detailed classification and analysis of various
nonlinear models that depict the diffusion-deformation process in hydrogels
caused by non-reactive solvent absorption. We have consolidated these theories
into a unified framework, demonstrating that, despite not being evident, all
theories follow equivalent thermodynamic arguments. For instance, while having
a common set of governing equations, each model showcases differences in the
enforcement of incompressibility and formulation of the constitutive model —
particularly in terms of the free energy function’s components, mainly at the
energetic level. At present, further research appears necessary to
conveniently account for the energetic component, ultimately aiming for a
cohesive and unified constitutive model.

Different numerical approaches adopted by various researchers to solve these
models were analyzed. Our implemented numerical methods presented reasonable
predictions for the diffusion-deformation process. However, our findings
suggest that the superiority of one strategy over another remains
inconclusive. The selection of an appropriate model should be grounded in a
comprehensive understanding of the hydrogel’s composition and must be
validated experimentally. Notice that some of the theories here introduced
have been extended to account for temperature variations, chemical reactions,
and damage (refer to, e.g., [42, 39, 28, 24] for more details). So, they
represent a solid foundation for the study of elastomeric materials.

We emphasized on discerning the differences among the leading models in the
literature and verifying whether the results presented by different authors
are still valid in light of open-source general-purpose software available
nowadays, such as FEniCS. To this end, an important part of this work is
Section 3 with the mathematical classification and resulting numerical
discretization and algorithms. These allowed us to carefully design
mathematical formulations, which were then implemented and used to study the
previously mentioned models.

Advances in automated solution techniques for the finite element method (FEM),
such as FEniCS, provide the user with a streamlined approach for solving
systems of partial differential equations compared to traditional model
development. Unlike models implemented on commercial software (e.g., Abaqus),
the implementation of our models grants users considerable control over
several components.

The simulation results indicate that there is not a one-size-fits-all model
compatible with the parameters across all scenarios. Depending on the specific
problem configuration being simulated, adjustments to one or multiple
parameters are necessary to guarantee the numerical stability of the solution.
Our observations in Section 5 revealed a strong dependency of the numerical
solution’s convergence on the diffusion coefficient and bulk modulus value. In
fact, very different diffusion coefficients for each example featured in this
study were required. Nonetheless, once a suitable value for the material
parameters is identified, the numerical solution achieves reliable accuracy
and robustness, as evidenced by the convergence analysis. Here, a consistent
number of Newton iterations was displayed along the simulation time to reach
absolute and relative errors smaller than $10^{-10}$. Moreover, concerning the
spatial and temporal discretization, quadratic and linear convergence orders,
respectively, for the FEM and Euler schemes for both the monolithic and
staggered approaches were observed.

Furthermore, constitutive model III stands to be the more fragile in terms of
numerical stability. Because constitutive model III directly penalizes the
material compressibility, it calls for a numerical solution using a mixed
formulation where the stress becomes a primary variable as it is standard in
the literature dealing with incompressible materials (see, e.g., [8] or [40]).

Another observation was regarding the application of boundary conditions.
Depending on the specific simulation scenario, it might be necessary to adopt
a time-ramping strategy to prevent numerical instabilities, which is well
known in continuum mechanics and due to the mathematical functional framework
in order to have compatible conditions. However, these mathematical
assumptions might not be met in many engineering applications. Therefore, one
must be careful with any assumption made while numerically solving the problem
at hand to avoid unphysical numerical results.

As an overall outcome, in Section 5, we found that each model presents diverse
deformation states and solute concentrations within the hydrogel, highlighting
the complexity of the investigated problem. It is evident that pinpointing a
single appropriate model to describe the diffusion-deformation of hydrogels
remains a challenge, given the different calibration mechanisms each offers.
These differences underscore the need for more comprehensive experimental data
to reconcile these theoretical distinctions with actual observations. With the
models’ validation, we can transition from merely describing the process to
predicting hydrogel behavior, thereby using the model for designing new
materials or optimizing the mechanical properties of existing ones.

Some efforts have been made to validate the diffusion-deformation process of
hydrogels as predicted by the coupled model. For example, [11] performed a
validation of a very similar theory as presented by [6] with a major focus on
the linear theory. Several measurements of the strain experienced by a
gelatin-glycerol-water hydrogel under free-swelling conditions are reported,
and the best-fitting parameters are identified. Neither the time evolution of
the diffusion process nor the internal stress were investigated. [5, 1]
performed an experimental work to validate an extension of the original models
presented by [14]. The extended models account for the viscoelastic response
of elastomeric gels under isothermal and non-isothermal conditions. Again, the
transient behavior of the diffusion process is neglected, and only the
stresses are computed for the steady state under external loading conditions.
Consequently, no study has been carried out to validate the models as
presented in this study.

Supplementary information No supplementary information was produced from this
research.