# A finite element implementation of finite deformation surface and bulk
poroelasticity

Jaemin Kim
Sibley School of Mechanical and Aerospace Engineering
Cornell University
Ithaca, NY 14853, USA
&Ida Ang
Sibley School of Mechanical and Aerospace Engineering
Cornell University
Ithaca, NY 14853, USA
&Francesco Ballarin
Dipartimento di Matematica e Fisica “N. Tartaglia"
Università Cattolica del Sacro Cuore
via della Garzetta 48, Brescia 25133, Italy
&Chung-Yuen Hui
Sibley School of Mechanical and Aerospace
Engineering, Cornell University
Ithaca, NY 14853, USA
Field of Theoretical and Applied Mechanics
Cornell University
Ithaca, NY 14853, USA
Global Station for Soft Matter, GI-CoRE
Hokkaido University
Sapporo, Japan
&Nikolaos Bouklas
Sibley School of Mechanical and Aerospace
Engineering, Cornell University
Ithaca, NY 14853, USA
Field of Theoretical and Applied Mechanics
Cornell University
Ithaca, NY 14853, USA Corresponding author<EMAIL_ADDRESS>
###### Abstract

We present a theoretical and computational model for the behavior of a porous
solid undergoing two interdependent processes, the finite deformation of a
solid and species migration through the solid, which are distinct in bulk and
on surface. Nonlinear theories allow us to systematically study porous solids
in a wide range of applications, such as drug delivery, biomaterial design,
fundamental study of biomechanics and mechanobiology, and the design of
sensors and actuators. As we aim to understand the physical phenomena at a
smaller length scale towards comprehending the fundamental biological
processes and the miniaturization of devices, the surface effect becomes more
pertinent. Although existing methodologies provide the necessary tools to
study coupled bulk effects for deformation and diffusion; however, very little
is known about fully coupled bulk and surface poroelasticity at finite strain.
Here we develop a thermodynamically consistent formulation for multiphysics
processes of surface and bulk poroelasticity, specialized for soft hydrated
solids, along with a corresponding finite element implementation. Our
multiphysical approach captures the interplay between competing processes of
finite deformation and species diffusion through the bulk and surface, and
provides invaluable insight when surface effects are important.

_Keywords_ Hydrogels $\cdot$ Surface energy $\cdot$ Surface diffusion $\cdot$
Swelling

##  1 Introduction

Poroelasticity is the theory that describes the coupling of deformations of a
continuum body with diffusive flow of a solvent through the body itself due to
its permeable/porous nature. The theory of linear poroelasticity as pioneered
by Biot (Biot, 1955, 1956, 1962) has been applied to describe hydrogels and
hydrated biological tissues (Hong et al., 2008, 2015; Li and Mooney, 2016),
but in many cases nonlinear variants of the theory of poroelasticity are
necessary due to bulk deformation of finite deformations and nonlinear complex
constitutive responses (Chester and Anand, 2010; Yoon et al., 2010; Li et al.,
2012; Bouklas and Huang, 2012; Bouklas et al., 2015; Ang et al., 2020). These
nonlinear theories have been highly relevant in the past few decades as the
study of hydrated biomaterials, tissues and hydrogels at finite strains.
Applications in drug delivery, tissue regeneration, biomaterials design,
fundamental studies in biomechanics and mechanobiology, design of morphing
structures as well as sensors and actuators, nonlinear theories of
poroelasticity have provided a deep understanding of the rich observed
material and structural responses. However, the majority of these studies have
been in the macroscale, where surface effect can often be neglected as bulk
energetic contributions dominate the response of the system. But, as we aim
towards miniaturization of devices, and also aim to probe and fundamentally
understand the effects that dominate the material response at smaller length
scales surface effect become more relevant. When the length scales that are
probed are sufficiently small, one can expect the coupling of elasticity and
solvent diffusion to be significantly influenced by surface contributions
toward all aspects of these multiphysical processes.

Following free energy-based phenomenological approaches for the bulk response
of materials, similar approaches can can be developed to study the surface and
interface response of materials and materials systems. The pioneering works of
Gurtin (Gurtin and Murdoch, 1975, 1978) established a systematic framework to
derive kinematics and balance principles for surfaces that focused on the
thermodynamics of purely elastic systems. Further, this framework was extended
by Gurtin to account for thermal effects and mass transport particularized for
surfaces (Gurtin and Jabbour, 2002; Cermelli et al., 2005). Many studies have
followed since exploring the mechanics of surfaces and interfaces for continua
involving purely mechanical as well as coupled processes (Dingreville et al.,
2005; Steinmann, 2008; Henann and Bertoldi, 2014; Javili et al., 2014; Liu et
al., 2020).

The competition of bulk and surface energetics is length scale dependent. For
elastic solids with shear modulus $\mathnormal{G}_{0}$ and surface energy
$\widetilde{\gamma}$, a relevant length scale is the elastocapillary length
$\ell=\widetilde{\gamma}/\mathnormal{G}_{0}$ (Style et al., 2017). Generally
speaking, if the characteristic dimension is much larger than elastocapillary
length ($\mathnormal{H}\gg\ell$), surface effects are negligible. On the other
hand, if the characteristic dimension is much smaller than elastocapillary
length ($\mathnormal{H}\ll\ell$), surface effects are dominant. Focusing on
soft solids, work from Style et al. (2017) and Bico et al. (2018) focused on
the implication of surface energetics on the mechanical behavior of this class
of materials, looking at the subtle relationships between surface stresses and
the bulk response. Surface and interface energetics have been further explored
in the context of soft solids to study their implications to adhesion,
fracture, instabilities as well as the design of composites.

Specialized boundary value problems allow for the development of analytical
solutions in this context (Bouklas and Huang, 2012; Ang et al., 2020).
However, finite element implementation that account for surface effects is
necessary to tackle general problems. Notably, Javili and Steinmann (Javili
and Steinmann, 2009, 2010) developed the finite element framework for finite
deformation accounting for surface energetics. Steinmann and co-authors
(McBride et al., 2011; Javili et al., 2013) also developed finite element
schemes to study several classes of coupled problems, including mass transport
and thermomechanics. These types of frameworks allowed the study of
morphogenesis and wound healing in microtissue systems, where “active” surface
and bulk effects were shown to emanate from the cellular activity of the
extracellular matrix (Mailand et al., 2019; Kim et al., 2020, 2023). Focused
on drug delivery applications, Bouklas and co-authors (Ang et al., 2020)
investigated the mechanical surface effect in the context of finite
deformation poroelasticity, uncovering the importance of such effects on the
transient response of swelling hydrogel microspheres of the porous solids
coupled with solvent diffusion using a mixed finite element method. More
recently, Dortdivanlioglu and co-authors (Rastogi and Dortdivanlioglu, 2022)
proposed a computational formulation that allows the consideration of
curvature-dependent surface energies using isogeometric analysis.

There is an extensive body of work on the development of mixed finite element
frameworks for linear poroelasticity, including stabilization approaches for
the treatment of numerical instabilities that are known to arise due to the
violation of the Ladyzhenskaya–Babuška–Brezzi (LBB) condition (Babuška, 1971;
Brezzi, 1974; Murad and Loula, 1994; Bathe, 2001). For different variants of
nonlinear poroelasticity, an array of approaches has been proposed (Hong et
al., 2008; Zhang et al., 2009; Chester and Anand, 2010; Duda et al., 2010;
Bouklas and Huang, 2012; Lucantonio et al., 2013; Chester et al., 2015;
Bouklas et al., 2015; MacMinn et al., 2016; Ang et al., 2020; Leronni and
Bardella, 2021). Although the existing methodologies provide the necessary
tools to study coupled bulk effects for deformation and diffusion, to the best
of the authors’ knowledge an implementation of the finite element methods that
fully couples bulk and surface poroelasticity at finite strain has not been
previously presented. Motivated by the theoretical work of McBride et al.
(2011) that considered coupled thermomechanics with species diffusion
accounting for individual surface energetic contributions from elasticity,
heat transfer and mass diffusion, here we aim to develop a robust and open-
source finite element framework for the treatment of surface and bulk
poroelasticity. In our previous work (Ang et al., 2020), we considered the
surface elasticity coupled with poroelasticity in bulk, not accounting for
surface diffusion.

The ultimate goal of this paper is the development of a thermodynamically-
consistent formulation for surface and bulk poroelasticity specialized for
soft hydrated solids, along with the development of a corresponding finite
element implementation. Focusing on soft solids, and more specifically
hydrogels, we discover a set of numerical complications that are not evident
when the general theoretical framework is presented. This is because of the
nature of the coupling of surface and bulk poroelasticity, as well as the
nearly incompressibility assumption that is utilized to model the elastomeric
component of the hydrogel. The manuscript is organized as follows: In Section
2, we describe the coupled theory, including kinematics, balance laws and the
general form of the constitutive equations considering both bulk and surface
energetics. In Section 3, the free energies are specialized for a hydrogel
material. Section 4 presents the finite element implementation including the
weak form, normalization and solution procedure. Section 5 elaborates on the
numerical simulations for two three-dimensional (3D) boundary value problems
corresponding to free-standing response of a cube-like structure as surface
effects dominate, as well as the response under external loading highlighting
unique contributions of surface effects, and the bulk and surface diffusion
pathways. Section 6 provides the conclusions and outlook for future work.

##  2 A continuum multiphysics approach for bulk and surface poroelasticity

In this section, we present a brief overview of key concept of differential
geometry towards describing the kinematics, mechanical equilibrium and mass
conservation, and obtain the general relationships for the constitutive
relations for the bulk and surface responses. It is noteworthy that a surface
quantity to be introduced is not necessarily the same as the bulk quantity
evaluated on the surface. Let $\mathnormal{V}$ be a fixed reference
configuration of a continuum body $\mathcal{B}$, where material points in the
reference configuration are tracked through the vector
$\textbf{X}\in\mathnormal{V}$.

######  Remark 1 (Notation).

The dyad ($\otimes$) and dot ($\cdot$) denote the tensor and dot products,
respectively. The double dot ($:$) denotes the double contraction over two
second-order tensors; e.g.
$\mathbf{A}:\mathbf{B}=\mathnormal{A}_{ij}\mathnormal{B}_{ij}$. For
convenience, we adopt the convention that other operations involving repeated
indices will not be indicated by a specific symbol, specifically meaning for
$\mathbf{v}=\mathbf{A}\mathbf{u}$ ($v_{i}=A_{ij}u_{j}$) and
$\mathbf{A}=\mathbf{B}\mathbf{C}$ ($A_{ij}=B_{jk}C_{kj}$), we do not write
$\mathbf{v}=\mathbf{A}\cdot\mathbf{u}$ and
$\mathbf{A}=\mathbf{B}\cdot\mathbf{C}$, respectively. Throughout the text, we
use the Einstein summation convention: when an index variable appears twice in
an expression, it implies summation of that variable over all values of the
index. Latin indices take the value 1, 2, 3, and are used in variables
describing the ambient, three-dimensional, Euclidean space. Greek indices take
the value 1, 2, and are used in variables in the embedded, two-dimensional,
surface. $\\{\bullet\\}$ and $\\{\widetilde{\bullet}\\}$ denote bulk and
surface quantities, respectively, for a body occupying volume $\mathnormal{V}$
bounded by outer surface denoted as $\mathnormal{S}$. The overhat
$\\{\widehat{\bullet}\\}$ denotes normalized quantities, and the superposed
dot $\\{\dot{\bullet}\\}$ denotes the material time derivative. We write
$\mathbf{A}^{-1}=\left(\mathbf{A}\right)^{-1}$ and
$\mathrm{tr}\left(\mathbf{A}\right)$ for the inverse and trace of a tensor
$\mathbf{A}$.

###  2.1 Key concepts of differential geometry

We briefly review the key concepts of differential geometry required to
describe the kinematics of the bulk and surface. Further details on
differential geometry can be found in books and papers, for example, Gurtin
(Gurtin and Murdoch, 1975), Green (Green and Zerna, 1992), Steinmann
(Steinmann, 2008; Javili et al., 2014) and Carmo (Do Carmo, 2016).

We first focus on the bulk description. Let $\mathbf{X}$ be an arbitrary
vector expressed in terms of the individual components of the curvilinear
coordinate system $\xi^{i}$, i.e.,
$\mathbf{X}=\mathbf{X}(\xi^{1},\xi^{2},\xi^{3})$. We can always associate two
new triplets of vectors ($\textbf{G}_{i}$ and $\textbf{G}^{i}$) with the
general curvilinear coordinate $\xi^{i}$, by Green and Zerna (1992).

$\displaystyle\textbf{G}_{i}=\frac{\partial\mathbf{X}}{\partial\xi^{i}}\quad$
$\displaystyle\text{and}\quad\textbf{G}^{i}=\frac{\partial\xi^{i}}{\partial\mathbf{X}}$
(1a)
$\displaystyle\textbf{I}=\delta^{i}_{j}\;\textbf{G}_{i}\otimes\textbf{G}^{j}\quad$
$\displaystyle\text{where}\quad\delta^{i}_{j}=\textbf{G}^{i}\cdot\textbf{G}_{j}$
(1b)

where $\textbf{G}_{i}$ and $\textbf{G}^{i}$ are referred to as covariant and
contravariant basis vectors, respectively. I and $\delta^{i}_{j}$ are the
mixed-variant identity tensor and the mixed Kronecker delta, respectively.
Note that $\delta^{i}_{j}=1$ for $i=j$, and $\delta^{i}_{j}=0$ for $i\neq j$
($i,j=1,2,3$).

The covariant and contravariant basis vectors are not necessarily orthogonal
to each other but linearly independent, requiring
$(\textbf{G}_{1}\times\textbf{G}_{2})\cdot\textbf{G}_{3}\neq 0$ (the vectors
do not lie in a plane) (Holzapfel, 2000). The vectors $\textbf{G}_{i}$ and
$\textbf{G}^{i}$ are connected to metric tensors $\textit{G}_{ij}$,
$\textit{G}^{ij}$ by

$\displaystyle\textbf{G}_{i}=\textit{G}_{ij}\textbf{G}^{j}\quad\text{with}\quad\textit{G}_{ij}=\textit{G}_{ji}=\textbf{G}_{i}\cdot\textbf{G}_{j}$
(2a)
$\displaystyle\textbf{G}^{i}=\textit{G}^{ij}\textbf{G}_{j}\quad\text{with}\quad\textit{G}^{ij}=\textit{G}^{ji}=\textbf{G}^{i}\cdot\textbf{G}^{j}$
(2b)

where two tensors are inverse to each other, i.e.,
$\textit{G}^{ik}\textit{G}_{kj}=\delta^{i}_{j}$. Note that metric tensors with
the identical indices represent the squares of the lengths; for example,
$\textit{G}_{11}=\lvert\textbf{G}_{1}\rvert^{2}$, whereas the metric tensors
with different indices represent the product of the lengths and the cosine of
the angle $\theta$; for example,
$\textit{G}_{12}=\lvert\textbf{G}_{1}\rvert\lvert\textbf{G}_{2}\rvert\cos\left(\theta(\textbf{G}_{1},\textbf{G}_{2})\right)$.

If a quantity $\\{\bullet\\}$ is a vector field defined throughout a volume V
which is bounded by a closed surface S, then bulk gradient and divergence
operators in 3D general curvilinear coordinates are defined as follows:

$\displaystyle\boldsymbol{\nabla}_{\textbf{X}}\\{\bullet\\}$
$\displaystyle=\frac{\partial\\{\bullet\\}}{\partial\xi^{i}}\otimes\textbf{G}^{i}$
(3a) $\displaystyle\boldsymbol{\nabla}_{\textbf{X}}\cdot\\{\bullet\\}$
$\displaystyle=\frac{\partial\\{\bullet\\}}{\partial\xi^{i}}\cdot\textbf{G}^{i}=\boldsymbol{\nabla}_{\textbf{X}}\\{\bullet\\}:\textbf{I}$
(3b)

With the above equations, the divergence theorem is (Green and Zerna, 1992)

$\int_{V}\boldsymbol{\nabla}_{\textbf{X}}\cdot\,\\{\bullet\\}\,\mathrm{d}V=\int_{S}\\{\bullet\\}\cdot\textbf{N}\,\mathrm{d}S$
(4)

where N is a unit outward normal vector on the surface $\mathnormal{S}$.

We now move our attention to the surface description. Let
$\widetilde{\mathbf{X}}$ denote an arbitrary vector expressed in terms of the
individual components of the general surface curvilinear coordinate system
$\widetilde{\xi}{}^{\alpha}$, i.e.,
$\widetilde{\mathbf{X}}=\widetilde{\mathbf{X}}(\widetilde{\xi}{}^{1},\widetilde{\xi}{}^{2}$).
We can always associate two new doublets of vectors
($\widetilde{\textbf{G}}_{\alpha}$ and $\widetilde{\textbf{G}}{}^{\alpha}$)
with the general curvilinear coordinate $\widetilde{\xi}{}^{\alpha}$, by
coordinate transformations (Green and Zerna, 1992).

$\displaystyle\widetilde{\textbf{G}}_{\alpha}=\frac{\partial\widetilde{\mathbf{X}}}{\partial\widetilde{\xi}{}^{\alpha}}\quad$
$\displaystyle\text{and}\quad\widetilde{\textbf{G}}{}^{\alpha}=\frac{\partial\widetilde{\xi}{}^{\alpha}}{\partial\widetilde{\mathbf{X}}}$
(5a)
$\displaystyle\widetilde{\textbf{I}}=\delta{}^{\alpha}_{\beta}\;\widetilde{\textbf{G}}_{\alpha}\otimes\widetilde{\textbf{G}}{}{}^{\beta}\quad$
$\displaystyle\text{where}\quad\widetilde{\delta}{}^{\alpha}_{\beta}=\widetilde{\textbf{G}}{}^{\alpha}\cdot\widetilde{\textbf{G}}_{\beta}$
(5b)

where $\widetilde{\textbf{G}}_{\alpha}$ and
$\widetilde{\textbf{G}}{}^{\alpha}$ are referred to as covariant and
contravariant surface basis vectors, respectively. $\widetilde{\textbf{I}}$
and $\delta{}^{\alpha}_{\beta}$ are the mixed-variant surface identity tensor
and the mixed Kronecker delta, respectively ($\alpha,\beta=1,2$). Note that
the surface identity tensor also can be expressed by
$\widetilde{\textbf{I}}=\textbf{I}-\textbf{N}\otimes\textbf{N}$.

The vectors $\widetilde{\textbf{G}}_{\alpha}$ and
$\widetilde{\textbf{G}}{}^{\alpha}$ are connected with the geometrical
characteristics

$\displaystyle\widetilde{\textbf{G}}_{\alpha}=\widetilde{\textit{G}}_{\,\alpha\beta}\widetilde{\textbf{G}}{}^{\beta}\quad\text{with}\quad\widetilde{\textit{G}}_{\,\alpha\beta}=\widetilde{\textit{G}}_{\,\beta\alpha}=\widetilde{\textbf{G}}_{\alpha}\cdot\widetilde{\textbf{G}}_{\beta}$
(6a)
$\displaystyle\widetilde{\textbf{G}}{}^{\alpha}=\widetilde{\textit{G}}{}^{\,\alpha\beta}\widetilde{\textbf{G}}_{\beta}\quad\text{with}\quad\widetilde{\textit{G}}{}^{\,\alpha\beta}=\widetilde{\textit{G}}{}^{\,\beta\alpha}=\widetilde{\textbf{G}}{}^{\alpha}\cdot\widetilde{\textbf{G}}{}^{\beta}$
(6b)

where $\widetilde{\textit{G}}_{\alpha\beta}$ and
$\widetilde{\textit{G}}{}^{\,\alpha\beta}$ are surface metric tensors, and
these two mappings are inverse to each other, i.e.
$\widetilde{\textit{G}}{}^{\,\alpha\gamma}\widetilde{\textit{G}}_{\gamma\beta}=\widetilde{\delta}{}^{\alpha}_{\beta}$.

If a quantity $\\{\widetilde{\bullet}\\}$ is a surface vector field defined
throughout a surface S, which is bounded by a closed curve L, then the surface
gradient and divergence operators in a 2D general curvilinear coordinates are
defined as follows:

$\displaystyle\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\\{\widetilde{\bullet}\\}$
$\displaystyle=\frac{\partial\\{\widetilde{\bullet}\\}}{\partial\widetilde{\xi}{}^{\alpha}}\otimes\widetilde{\textbf{G}}{}^{\alpha}$
(7a)
$\displaystyle\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\cdot\\{\widetilde{\bullet}\\}$
$\displaystyle=\frac{\partial\\{\widetilde{\bullet}\\}}{\partial\widetilde{\xi}{}^{\alpha}}\cdot\widetilde{\textbf{G}}{}^{\alpha}=\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\\{\widetilde{\bullet}\\}:\widetilde{\textbf{I}}$
(7b)

With the above equations, we introduce the surface divergence theorem as
follows (Steinmann, 2008):

$\int_{\mathnormal{S}}\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\cdot\\{\widetilde{\bullet}\\}\,\text{d}S=\int_{\mathnormal{L}}\\{\widetilde{\bullet}\\}\cdot\widetilde{\textbf{N}}\,\text{d}\mathnormal{L}-\int_{\mathnormal{S}}\widetilde{\kappa}\\{\widetilde{\bullet}\\}\cdot\textbf{N}\,\text{d}S$
(8)

where $\widetilde{\textbf{N}}$ is the unit outward bi-normal vector to the
boundary curve $\mathnormal{L}$ (see Figure 1), and
$\widetilde{\kappa}=-\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\cdot\textbf{N}$
is the total curvature (twice the mean surface curvature) (Gurtin and Murdoch,
1975).

Figure 1: Schematic illustration of the reference, initial, and current state
of a continuum body. The initial state is assumed to be isotropically scaled
from the reference state. The reference volume and surface, and boundary are
denoted by $\mathnormal{V}$ and $\mathnormal{S}$, and $\mathnormal{L}$
respectively. The normal vector to the surface in the reference and current
configuration ($\mathbf{N}$ and $\mathbf{n}$) and the bi-normal vector to the
boundary ($\widetilde{\mathbf{N}}$) are shown, where the over-tilde indicates
surface quantities.

###  2.2 Kinematics

We use the notation $\chi:\mathnormal{V}\rightarrow\mathbb{R}^{3}$ for the
deformation of body $\mathcal{B}$. A motion $\chi$ is the vector field of the
mapping $\textit{{x}}=\chi(\textbf{X},\mathnormal{t})$, of a material point in
the reference configuration $\textbf{X}\in\mathnormal{V}$ to a position in the
deformed configuration $\textit{{x}}\in\mathnormal{v}$. The kinematics of a
typical particle are described by the displacement vector field in the spatial
description,
$\textbf{u}(\textbf{X},t)=\textit{{x}}\,(\textbf{X},t)-\textbf{X}$. The
kinematics of an infinitesimal bulk element are described by

$\displaystyle\textbf{F}(\textbf{X},\textit{t})$
$\displaystyle=\frac{\partial\chi(\textbf{X},\textit{t})}{\partial\textbf{X}}=\boldsymbol{\nabla}_{\textbf{X}}\,\textit{{x}}\,(\mathbf{X},t)$
(9a) $\displaystyle\textbf{F}^{-1}(\textit{{x}},\textit{t})$
$\displaystyle=\frac{\partial\chi^{-1}(\textit{{x}},\textit{t})}{\partial\textit{{x}}}=\boldsymbol{\nabla}_{\textit{{x}}}\,\textbf{X}(\textit{{x}},\textit{t})$
(9b)

where $\textbf{F}(\textbf{X},\textit{t})$ and
$\textbf{F}^{-1}(\textit{{x}},\textit{t})$ are the deformation gradient and
inverse deformation gradient, respectively. Note that
$\mathnormal{J}(\textbf{X},\textit{t})=\textrm{d}\textit{v}/\textrm{d}\textit{V}=\text{det}\,\textbf{F}(\textbf{X},\textit{t})>0$
is the Jacobian determinant of the deformation gradient defining the ratio of
a volume element between the material and spatial configuration.

Assuming conformity, the surface displacement
$\widetilde{\textbf{u}}(\widetilde{\textbf{X}},t)$ can be determined by
$\mathbf{u}(\mathbf{X},t)|_{\mathnormal{S}}=\widetilde{\textbf{u}}(\widetilde{\textbf{X}},t)$.
The motion of an arbitrary differential vector element $\mathrm{d}\mathbf{X}$
can be mapped by the deformation gradient F to a vector $\mathrm{d}\mathbf{x}$
in the deformed configuration. However, a unit normal vector N in the material
configuration cannot be transformed into a unit normal vector n in the spatial
configuration using solely the deformation gradient (Holzapfel, 2000;
Steinmann, 2008). This motivates us to follow the kinematics of an
infinitesimal surface element (Kim et al., 2020, 2023).

$\displaystyle\widetilde{\textbf{F}}(\widetilde{\textbf{X}},\textit{t})$
$\displaystyle=\frac{\partial\chi(\widetilde{\textbf{X}},\textit{t})}{\partial\widetilde{\textbf{X}}}\cdot\widetilde{\mathbf{I}}=\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\,\widetilde{\textit{{x}}}\,(\widetilde{\mathbf{X}},t)$
(10a)
$\displaystyle\widetilde{\textbf{F}}^{-1}(\widetilde{\textit{{x}}},\textit{t})$
$\displaystyle=\frac{\partial\chi^{-1}(\widetilde{\textit{{x}}},\textit{t})}{\partial\widetilde{\textit{{x}}}}\cdot\widetilde{\textit{{i}}}=\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textit{{x}}}}\,\widetilde{\textbf{X}}(\widetilde{\textit{{x}}},\textit{t})$
(10b)

where $\widetilde{\textbf{F}}(\widetilde{\textbf{X}},\textit{t})$ and
$\widetilde{\textbf{F}}^{-1}(\widetilde{\textit{{x}}},\textit{t})$ are the
surface deformation gradient and inverse surface deformation gradient. Note
$\widetilde{\textbf{F}}$ is a mapping from the surface (2D) to the bulk (3D),
so it is rank-deficient (it has rank 2, whereas full-rank would be 3). The
$\widetilde{\textbf{I}}=\textbf{I}-\textbf{N}\otimes\textbf{N}$ and
$\widetilde{\textit{{i}}}=\textit{{i}}-\textbf{n}\otimes\textbf{n}$ are the
mixed surface unit tensors with the outward unit normal vectors N and n, where
$\widetilde{\mathbf{I}}$ and $\widetilde{\textit{{i}}}$ act as a surface
(idempotent) projection tensors in material and spatial configurations,
respectively. Note that
$\widetilde{\mathnormal{J}}(\widetilde{\textbf{X}},\textit{t})=\textrm{d}\textit{a}/\textrm{d}\textit{A}=\text{det}\,\widetilde{\textbf{F}}(\widetilde{\textbf{X}},\textit{t})>0$
is the Jacobian determinant of the surface deformation gradient defining the
area ratio of a surface element between material and spatial configuration.

We introduce $\mathbf{C}$ and $\widetilde{\mathbf{C}}$, the right Cauchy-Green
tensors in the bulk and on the surface, respectively, as

$\displaystyle\mathbf{C}$ $\displaystyle=\mathbf{F}^{\mathrm{T}}\mathbf{F}$
(11a) $\displaystyle\widetilde{\textbf{C}}$
$\displaystyle=\widetilde{\textbf{F}}^{\text{T}}\,\widetilde{\textbf{F}}$
(11b)

and $\mathnormal{I}_{1}=\text{tr}(\mathbf{C})$ is the first principal
invariant. The detailed derivation for the surface kinematics can be found in
Green and Zerna (1992); Steinmann (2008); Do Carmo (2016). Note that we cannot
perform the inverse of the surface right Cauchy-Green tensor due to its rank
deficiency. Nevertheless, we can still obtain its inverse form in the
generalized sense,

$\widetilde{\mathbf{C}}^{-1}=\widetilde{\mathbf{I}}\mathbf{C}^{-1}\widetilde{\mathbf{I}}$
(12)

which will be utilized in the forthcoming developments for defining the
surface kinetic law.

###  2.3 Mechanical equilibrium

Mechanical equilibrium is assumed to be maintained at all times during the
motion. The strong form of the corresponding governing equation is

$\displaystyle\boldsymbol{\nabla}_{\textbf{X}}\cdot\mathbf{P}+\mathbf{B}$
$\displaystyle=0\quad$ $\displaystyle\text{in}\quad\mathnormal{V}$ (13a)
$\displaystyle\mathbf{P}\mathbf{N}-\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\cdot\widetilde{\mathbf{P}}$
$\displaystyle=\mathbf{T}\quad$
$\displaystyle\text{on}\quad\mathnormal{S}_{\mathrm{T}}$ (13b)
$\displaystyle\mathbf{u}$ $\displaystyle=\mathbf{u}_{p}\quad$
$\displaystyle\text{on}\quad\mathnormal{S}_{\mathrm{u}}$ (13c)
$\displaystyle[[\widetilde{\mathbf{P}}\widetilde{\mathbf{N}}]]$
$\displaystyle=0\quad$ $\displaystyle\text{on}\quad\mathnormal{L}$ (13d)

where $\mathbf{P}$ and $\widetilde{\mathbf{P}}$ are the first Piola-Kirchoff
stresses in bulk and on surface, $\mathbf{B}$ and $\mathbf{T}$ are the body
force and the traction vector, and $\mathbf{u}_{p}$ is the prescribed
displacement. Note that a Neumann-type boundary condition is also defined on
boundary curves that $[[\bullet]]$ indicates summation over surfaces
intersecting on boundary curves (Steinmann, 2008).

###  2.4 Species conservation

Through species (or mass) conservation, the strong form for the corresponding
governing equation is (McBride et al., 2011):

$\displaystyle\dot{\mathnormal{C}}+\boldsymbol{\nabla}_{\textbf{X}}\cdot\mathbf{J}$
$\displaystyle=\mathnormal{r}\quad$
$\displaystyle\text{in}\quad\mathnormal{V}$ (14a)
$\displaystyle\dot{\widetilde{\mathnormal{C}}}+\widetilde{\boldsymbol{\nabla}}_{\textbf{X}}\cdot\widetilde{\mathbf{J}}-\mathbf{J}\cdot\mathbf{N}$
$\displaystyle=\mathnormal{i}\quad$
$\displaystyle\text{on}\quad\mathnormal{S}_{\widetilde{\mathnormal{C}}}$ (14b)
$\displaystyle\widetilde{\mathbf{J}}$
$\displaystyle=\widetilde{\mathbf{J}}_{p}\quad$
$\displaystyle\text{on}\quad\mathnormal{S}_{\widetilde{\mathbf{J}}}$ (14c)
$\displaystyle[[\widetilde{\mathbf{J}}\cdot\widetilde{\mathbf{N}}]]$
$\displaystyle=0\quad$ $\displaystyle\text{on}\quad\mathnormal{L}$ (14d)

where $\mathnormal{C}$ is the bulk nominal concentration (the number of
species per unit reference volume), $\widetilde{\mathnormal{C}}$ is the
surface nominal concentration (the number of species per unit reference area),
$\mathbf{J}$ is the bulk nominal flux (the number of species per unit time per
unit area), and $\widetilde{\mathbf{J}}$ is the surface nominal flux (the
number of species per unit time per unit length). The $\mathnormal{r}$ and
$\mathnormal{i}$ are the source/sink terms for the number of species injected
into the reference volume and area per unit time, and
$\widetilde{\mathbf{J}}_{p}$ is the prescribed surface flux. Equation 14a is
the standard species conservation equation prescribed for the bulk in the
reference configuration. Equation 14b describes the species balance on the
surface, where the third term on the left hand is similar to the source
$\mathnormal{i}$, but describes the outward normal flux from the bulk to the
surface; note that species can get to the surface either from the exterior
(environment) or from the interior (bulk). Additionally, Equation 14 are
supplemented with initial conditions,

$\displaystyle\mathnormal{C}(\mathbf{X},\mathnormal{t}=0)$
$\displaystyle=\mathnormal{C}_{0}$ (15a)
$\displaystyle\widetilde{\mathnormal{C}}(\widetilde{\mathbf{X}},\mathnormal{t}=0)$
$\displaystyle=\widetilde{\mathnormal{C}}_{0}$ (15b)

where $\mathnormal{C}_{0}$ and $\widetilde{\mathnormal{C}}_{0}$ are the
initial species concentration in the bulk and on the surface at time $t=0$.
Note that we assume homogeneous initial conditions; this choice will become
apparent as we discuss the general solution process for the corresponding
boundary value problems in Section 5.

###  2.5 Thermodynamic considerations

In addition to the free energy density in the bulk, we consider the free
energy density on the surface. We will denote the bulk and surface free energy
densities as

$\Psi(\mathbf{F},\mathnormal{C})\quad\text{and}\quad\widetilde{\Psi}(\widetilde{\mathbf{F}},\widetilde{\mathnormal{C}})$
(16)

where we assume that they are functions of the deformation gradient and
nominal concentration in the bulk and on the surface, respectively.

Considering a system that includes an elastic and porous solid component
coupled with species that are free to migrate in the porous network, the rate
of change of the system’s free energy $\mathcal{G}$ has to account for several
effects (Holzapfel, 2000; Gurtin et al., 2010; Hong et al., 2008; Bouklas et
al., 2015; Ang et al., 2020). This can be expressed as

$\dot{\mathcal{G}}=\int_{\mathnormal{V}}\dot{\Psi}\,\mathrm{d}V+\int_{\mathnormal{S}}\dot{\widetilde{\Psi}}\,\mathrm{d}S-\int_{\mathnormal{V}}\mathbf{B}\dot{\mathbf{x}}\,\mathrm{d}V-\int_{\mathnormal{S}}\mathbf{T}\dot{\mathbf{x}}\,\mathrm{d}S-\int_{\mathnormal{V}}\mu\mathnormal{r}\,\mathrm{d}V-\int_{\mathnormal{S}}\widetilde{\mu}\mathnormal{i}\,\mathrm{d}S$
(17)

where the third and fourth terms are the rate of mechanical work by the body
force $\mathbf{B}$ and traction vector $\mathbf{T}$, and the fifth and sixth
terms are the rate of chemical work by the bulk chemical potential $\mu$ and
the surface chemical potential $\widetilde{\mu}$. Note that thermodynamics
dictate that the free energy of the system should not increase, i.e.,
$\dot{\mathcal{G}}\leq 0$.

Substituting Equation 13 and Equation 14 into Equation 17, the rate of change
of the free energy of the system can be expressed as follows,

$\displaystyle\dot{\mathcal{G}}$
$\displaystyle=\int_{\mathnormal{V}}\dot{\Psi}\,\mathrm{d}V+\int_{\mathnormal{S}}\dot{\widetilde{\Psi}}\,\mathrm{d}S-\int_{\mathnormal{V}}\mathbf{P}:\dot{\mathbf{F}}\,\mathrm{d}V-\int_{\mathnormal{S}}\widetilde{\mathbf{P}}:\dot{\widetilde{\mathbf{F}}}\,\mathrm{d}S$
$\displaystyle-\int_{\mathnormal{V}}\mu\dot{\mathnormal{C}}\,\mathrm{d}V-\int_{\mathnormal{S}}\widetilde{\mu}\dot{\widetilde{\mathnormal{C}}}\,\mathrm{d}S+\int_{\mathnormal{S}}\left(\widetilde{\mu}-\mu\right)\mathbf{J}\cdot\mathbf{N}\,\mathrm{d}S$
$\displaystyle+\int_{\mathnormal{V}}\mathbf{J}\cdot\boldsymbol{\nabla}_{\mathbf{X}}\mu\,\mathrm{d}V+\int_{\mathnormal{S}}\widetilde{\mathbf{J}}\cdot\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\widetilde{\mu}\,\mathrm{d}S\leq
0$ (18)

Using the chain-rule, the rate of bulk and surface free energy densities can
be expressed as

$\displaystyle\dot{\Psi}$
$\displaystyle=\frac{\partial\Psi}{\partial\mathbf{F}}:\dot{\mathbf{F}}+\frac{\partial\Psi}{\partial\mathnormal{C}}\dot{\mathnormal{C}}$
(19a) $\displaystyle\dot{\widetilde{\Psi}}$
$\displaystyle=\frac{\partial\widetilde{\Psi}}{\partial\widetilde{\mathbf{F}}}:\dot{\widetilde{\mathbf{F}}}+\frac{\partial\widetilde{\Psi}}{\partial\widetilde{\mathnormal{C}}}\dot{\widetilde{\mathnormal{C}}}$
(19b)

By substituting Equation 19 into Equation 18, and rearranging terms yields

$\displaystyle\dot{\mathcal{G}}$
$\displaystyle=\int_{\mathnormal{V}}\left(\frac{\partial\Psi}{\partial\mathbf{F}}-\mathbf{P}\right):\dot{\mathbf{F}}\,\mathrm{d}V+\int_{\mathnormal{S}}\left(\frac{\partial\widetilde{\Psi}}{\partial\widetilde{\mathbf{F}}}-\widetilde{\mathbf{P}}\right):\dot{\widetilde{\mathbf{F}}}\,\mathrm{d}S$
$\displaystyle+\int_{\mathnormal{V}}\left(\frac{\partial\Psi}{\partial\mathnormal{C}}-\mu\right)\dot{\mathnormal{C}}\,\mathrm{d}V+\int_{\mathnormal{S}}\left(\frac{\partial\widetilde{\Psi}}{\partial\widetilde{\mathnormal{C}}}-\widetilde{\mu}\right)\dot{\widetilde{\mathnormal{C}}}\,\mathrm{d}S$
$\displaystyle+\int_{\mathnormal{S}}\left(\widetilde{\mu}-\mu\right)\mathbf{J}\cdot\mathbf{N}\,\mathrm{d}S+\int_{\mathnormal{V}}\mathbf{J}\cdot\boldsymbol{\nabla}_{\mathbf{X}}\mu\,\mathrm{d}V+\int_{\mathnormal{S}}\widetilde{\mathbf{J}}\cdot\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\widetilde{\mu}\,\mathrm{d}S\leq
0$ (20)

where each integral represents a distinct mechanism of energy dissipation,
associated with mechanical and chemical works. The inequality must hold at
every point of the continuum body and at all times during a thermodynamic
process. To satisfy the constraint, the Coleman-Noll procedure (Holzapfel,
2000) states that each integrand in Equation 20 to be either negative or equal
to zero.

From the first four terms in Equation 20 we obtain the following constitutive
relations

$\displaystyle\mathbf{P}$
$\displaystyle=\frac{\partial\Psi(\mathbf{F},\mathnormal{C})}{\partial\mathbf{F}}\quad$
$\displaystyle\text{and}\quad\mu$
$\displaystyle=\frac{\partial\Psi(\mathbf{F},\mathnormal{C})}{\partial\mathnormal{C}}$
(21a) $\displaystyle\widetilde{\mathbf{P}}$
$\displaystyle=\frac{\partial\widetilde{\Psi}(\widetilde{\mathbf{F}},\widetilde{\mathnormal{C}})}{\partial\widetilde{\mathbf{F}}}\quad$
$\displaystyle\text{and}\quad\widetilde{\mu}$
$\displaystyle=\frac{\partial\widetilde{\Psi}(\widetilde{\mathbf{F}},\widetilde{\mathnormal{C}})}{\partial\widetilde{\mathnormal{C}}}$
(21b)

for the bulk and surface first Piola-Kirchhoff stresses tensors, and the bulk
and surface chemical potentials, respectively.

For the fifth term in Equation 20, the most obvious way to guarantee the
dissipation inequality on surface is to impose the condition (McBride et al.,
2011),

$\mu(\mathbf{X},t)=\widetilde{\mu}(\widetilde{\textbf{X}},t)\quad\text{on}\quad\mathnormal{S}$
(22)

which prescribes local chemical equilibrium between the surface and the bulk.
This point is discussed in McBride et al. (2011), where also the possibility
of an alternate strategy to satisfy that constraint is highlighted.

Finally, to maintain that the last two terms in Equation 20 remain to be
negative or zero, we adopt kinetic laws for diffusion (Hong et al., 2008) in
bulk and on surface. This allows us to maintain negative semi-definiteness and
describes the consistent species diffusion that is driven by gradients of
chemical potential:

$\displaystyle\mathbf{J}$
$\displaystyle=-\mathbf{M}\boldsymbol{\nabla}_{\mathbf{X}}\mu$ (23a)
$\displaystyle\widetilde{\mathbf{J}}$
$\displaystyle=-\widetilde{\mathbf{M}}\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\widetilde{\mu}$
(23b)

where $\mathbf{M}$ and $\widetilde{\mathbf{M}}$ are the bulk and surface
mobility tensors. In the upcoming section, we will specialize the symmetric
and positive definite mobility tensors $\mathbf{M}$ and
$\widetilde{\mathbf{M}}$ to fully define the constitutive laws for the fluxes.

##  3 Specific considerations for hydrogels

In this work, we will focus on a coupled bulk and surface poroelastic
framework for hydrogels. We have to specialize our choices for the surface and
bulk free energy densities, corresponding constitutive laws, and definition of
mobility tensors for our theory to be complete and to be able to proceed to
the development of the numerical solution scheme.

###  3.1 Particularizing the surface and bulk free energy densities

For the free energy densities of polymer in the bulk, we adopt the Flory-
Huggins model Flory (1942, 1961); Treloar (1975); Hong et al. (2008), assuming
the following additive decomposition into elastic and mixing contributions

$\Psi(\mathbf{F},\mathnormal{C})=\Psi_{e}(\mathbf{F})+\Psi_{m}(\mathnormal{C})\,.$
(24)

These can be individually specialized to

$\displaystyle\Psi_{e}(\mathbf{F})$
$\displaystyle=\frac{\mathnormal{1}}{2}\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}(\mathnormal{I}_{1}-3-2\ln{\mathnormal{J}})$
(25a) $\displaystyle\Psi_{m}(\mathnormal{C})$
$\displaystyle=-\frac{\mathnormal{k}_{B}\mathnormal{T}}{\Omega}\left[\Omega\mathnormal{C}\,\ln\left(\frac{1+\Omega\mathnormal{C}}{\Omega\mathnormal{C}}\right)+\frac{\chi}{1+\Omega\mathnormal{C}}\right]$
(25b)

where $\Omega$ is the volume of a solvent molecule, $\chi$ is a dimensionless
parameter of polymer-solvent mixing, $\mathnormal{N}$ is the number of polymer
chains per unit reference volume, $\mathnormal{k}_{B}$ is Boltzmann’s
constant, and $\mathnormal{T}$ is the absolute temperature. Note that the
first and second terms in Equation 25b represent the entropy and enthalpy of
mixing, respectively Hong et al. (2008).

Following Hong et al. (2008, 2009); Bouklas and Huang (2012); Bouklas et al.
(2015); Ang et al. (2020), we assume that the polymer chains and the diffusion
species are individually incompressible. Furthermore, the gel is a condensed
matter with negligible void space, so any volume change of the hydrogel is due
to species diffusion, thus

$1+\Omega\mathnormal{C}=\mathnormal{J}\Rightarrow\mathnormal{C}=\frac{\mathnormal{J}-1}{\Omega}\,.$
(26)

Similar to Equation 24, we assume that the surface free energy density of the
hydrogel is also decomposed into elastic and mixing contributions as (see
Remark 2))

$\widetilde{\Psi}(\widetilde{\mathbf{F}},\widetilde{\mathnormal{C}})=\widetilde{\Psi}_{e}(\widetilde{\mathbf{F}},\widetilde{\mathnormal{C}})+\widetilde{\beta}\widetilde{\Psi}_{m}(\widetilde{\mathnormal{C}})$
(27)

which can be individually expressed as

$\displaystyle\widetilde{\Psi}_{e}(\widetilde{\mathbf{F}},\widetilde{\mathnormal{C}})$
$\displaystyle=\frac{\widetilde{\kappa}}{2}\left(\widetilde{\mathnormal{J}}-1-\widetilde{\Omega}\widetilde{\mathnormal{C}}\right)^{2}+\widetilde{\gamma}\widetilde{\mathnormal{J}}$
(28a) $\displaystyle\widetilde{\Psi}_{m}(\widetilde{\mathnormal{C}})$
$\displaystyle=-\frac{\mathnormal{k}_{B}\mathnormal{T}}{\widetilde{\Omega}}\left[\widetilde{\Omega}\widetilde{\mathnormal{C}}\,\ln\left({\frac{1+\widetilde{\Omega}\widetilde{\mathnormal{C}}}{\widetilde{\Omega}\widetilde{\mathnormal{C}}}}\right)+\frac{\widetilde{\chi}}{1+\widetilde{\Omega}\widetilde{\mathnormal{C}}}\right]$
(28b)

In a microscopic description, the length scales of polymer chains and solvent
molecules are very different. Focusing on the surface, where the interaction
of polymer chains and solvent molecules is crucial, and more specifically the
finite thickness of the surface zone, here we consider the the idealization of
a constant number of layers $\widetilde{\beta}$ of potential sites that
solvent molecules can occupy in the surface zone. Here, $\widetilde{\Omega}$
is the area that a solvent molecule occupies on the surface,
$\widetilde{\chi}$ is a dimensionless parameter of polymer-solvent mixing. In
the second term of the right hand side of Equation 28a, we consider a constant
surface energy per unit current area
$\widetilde{\gamma}\,[\mathrm{J}/\mathrm{m}^{2}]$, which leads to a fluid-like
response for the aggregate. Additionally, it is not realistic to enforce the
equivalent of exact incompressibility for the surface, as one cannot constrain
the number of layers of species on the surface, even though the constituents
are individually incompressible. Thus, we introduce a penalty term to loosely
introduce a constraint for species concentration in the surface connected with
the area change of the boundary, with a penalty coefficient
$\widetilde{\kappa}\,[\mathrm{J}/\mathrm{mol}/\mathrm{m}^{2}]$.

######  Remark 2 (The application of Flory-Huggins model to surface
poroelasticity).

In the Flory-Huggins model Huggins (1941); Flory (1942), the solvent and
polymer molecules are considered to be arranged in the 3D lattice sites such
that each site may be occupied either by a solvent molecule or by a segment of
the polymer chain, and the equations of entropy and enthalpy of mixing are
derived in the statistical thermodynamics approach (see Chapter 7 of Treloar
(1975) for the details). Without loss of generality, the Flory-Huggins model
can be also considered on the 2D lattice model, where the probability of
encountering a site occupied by the solvent and polymer molecules is the
function of the area. Although the number of adjacent sites would be different
in 2D and 3D models, the Flory parameters $\chi$ and $\widetilde{\chi}$ are
introduced to eliminate the number of adjacent sites, and their dependence on
the dimensionality is to be fitted by experimental data. We introduce the
parameter $\widetilde{\beta}$ in the mixing part of surface free energy
$\widetilde{\Psi}_{m}$ to account for the number of 2D lattice layers that
arise due to the size difference of solvent molecules and polymer chains
considering the surface zone, but we do not consider the interaction energy
between the layers. Note that the elastic part of surface free energy
$\widetilde{\Psi}_{e}$ does not include the number of layer parameter
$\widetilde{\beta}$ which we assume to be a constant.

###  3.2 Constitutive relations

Using Equations 21, 25 and 28, the specific constitutive relations are
obtained as follows:

$\displaystyle\mathbf{P}$
$\displaystyle=\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}\left(\mathbf{F}+\alpha_{1}\mathnormal{J}\mathbf{F}^{-\mathrm{T}}\right)\quad\text{with}\quad\alpha_{1}=-\frac{1}{\mathnormal{J}}+\frac{1}{\mathnormal{N}\Omega}\left[\frac{1}{\mathnormal{J}}+\ln\left(\frac{\mathnormal{J}-1}{\mathnormal{J}}\right)+\frac{\chi}{\mathnormal{J}^{2}}\right]$
(29a) $\displaystyle\mu$
$\displaystyle=\mathnormal{k}_{B}\mathnormal{T}\left[\ln\left(\frac{\Omega\mathnormal{C}}{1+\Omega\mathnormal{C}}\right)+\frac{1}{1+\Omega\mathnormal{C}}+\frac{\chi}{\left(1+\Omega\mathnormal{C}\right)^{2}}\right]$
(29b) $\displaystyle\widetilde{\mathbf{P}}$
$\displaystyle=\left[\widetilde{\kappa}\left(\widetilde{\mathnormal{J}}-1-\widetilde{\Omega}\widetilde{\mathnormal{C}}\right)+\widetilde{\gamma}\right]\widetilde{\mathnormal{J}}\widetilde{\mathbf{F}}^{-\mathrm{T}}$
(29c) $\displaystyle\widetilde{\mu}$
$\displaystyle=\mathnormal{k}_{B}\mathnormal{T}\left[\ln\left(\frac{\widetilde{\Omega}\widetilde{\mathnormal{C}}}{1+\widetilde{\Omega}\widetilde{\mathnormal{C}}}\right)+\frac{1}{1+\widetilde{\Omega}\widetilde{\mathnormal{C}}}+\frac{\widetilde{\chi}}{\left(1+\widetilde{\Omega}\widetilde{\mathnormal{C}}\right)^{2}}\right]-\widetilde{\Omega}\,\widetilde{\kappa}\left(\widetilde{J}-1-\widetilde{\Omega}\widetilde{\mathnormal{C}}\right)\
\Rightarrow\widetilde{\mathnormal{C}}=\widetilde{\mathnormal{C}}\,(\widetilde{\mu},\widetilde{\mathnormal{J}})$
(29d)

which are the specific forms of the constitutive relations for the first
Piola-Kirchhoff stresses and the chemical potential in the bulk and on the
surface. Note that we do not use the constitutive relation of Equation 29b,
but the incompressibility condition of Equation 26 to determine the bulk
concentration $\mathnormal{C}$. On the other hand, we do not have an
incompressibility condition for solvent molecules on surface, and we should
use the constitutive relation Equation 29d to determine the surface
concentration $\widetilde{\mathnormal{C}}$. However, there is no closed-form
solution for the surface concentration in Equation 29d, so we need to solve
the constitutive relation numerically.

###  3.3 Bulk and surface diffusion

Particularizing the expression for the bulk mobility, again following Hong et
al. (2008), and prescribing an equivalent definition for surface mobility

$\displaystyle\mathbf{M}$
$\displaystyle=\frac{\mathnormal{C}\mathnormal{D}}{\mathnormal{k}_{B}\mathnormal{T}}\mathbf{C}^{-1}$
(30a) $\displaystyle\widetilde{\mathbf{M}}$
$\displaystyle=\frac{\widetilde{\mathnormal{C}}\widetilde{\mathnormal{D}}}{\mathnormal{k}_{B}\mathnormal{T}}\widetilde{\mathbf{C}}^{-1}$
(30b)

where $\mathnormal{k}_{B}$ is the Boltzmann’s constant, $\mathnormal{T}$ is an
absolute temperature. The coefficient of diffusion of the solvent molecules
$\mathnormal{D}$ in the bulk and $\widetilde{\mathnormal{D}}$ on the surface
are assumed to be isotropic and independent of the deformation and
concentration as the simplest approximation Hong et al. (2008). Note that the
current concentrations are related to the nominal concentrations as
$\mathnormal{c}=\mathnormal{C}/\mathnormal{J}$ and
$\widetilde{\mathnormal{c}}=\widetilde{\mathnormal{C}}/\widetilde{\mathnormal{J}}$.

##  4 Mixed finite element formulation

This section presents a finite element formulation based on the nonlinear
theory in Sections 2 and 3. The main aim here is to provide an accessible
open-source implementation that will be utilized by the research community,
and thus, we choose FEniCS Logg et al. (2012a); Alnæs et al. (2015) for the
implementation; a choice which also affects the set-up of the mixed finite
element formulation. The formulation starts with the strong form of the
governing equations and initial and boundary conditions. We first introduce
the weak form of the problem and subsequently describe the normalization,
discretization, solution steps, and some specific implementation details for
FEniCS. As the main focus here is the coupled mechanics between diffusion and
deformation in bulk and on surface, we here-on neglect the source terms, i.e.,
$\mathnormal{r}=\mathnormal{i}=0$.

###  4.1 Three-field weak form

Since the chemical boundary conditions for hydrogels are often specified in
terms of chemical potential or diffusion flux (proportional to the gradient of
chemical potential), it is convenient to use the chemical potential (instead
of solvent concentration) as an independent variable in the finite element
formulation. Additionally, this allows us to avoid $\mathcal{C}^{1}$
continuity requirements Bouklas et al. (2015). For this purpose, we rewrite
the free energy densities as a function of the deformation gradient and
chemical potential through a Legendre transform Hong et al. (2008); McBride et
al. (2011); Bouklas et al. (2015), which replaces a variable with its
thermodynamic conjugate:

$\displaystyle\Phi(\mathbf{F},\mu)$
$\displaystyle=\Psi(\mathbf{F},\mathnormal{C})-\mu\mathnormal{C}$ (31a)
$\displaystyle\widetilde{\Phi}(\widetilde{\mathbf{F}},\widetilde{\mu})$
$\displaystyle=\widetilde{\Psi}(\widetilde{\mathbf{F}},\widetilde{\mathnormal{C}})-\widetilde{\mu}\widetilde{\mathnormal{C}}$
(31b)

After the Legendre transform, similar to Equation 21, the constitutive
relations can be rewritten as follows:

$\displaystyle\mathbf{P}$
$\displaystyle=\frac{\partial\Phi(\mathbf{F},\mu)}{\partial\mathbf{F}}\quad$
$\displaystyle\text{and}\quad\mathnormal{C}$
$\displaystyle=-\frac{\partial\Phi(\mathbf{F},\mu)}{\partial\mu}$ (32a)
$\displaystyle\widetilde{\mathbf{P}}$
$\displaystyle=\frac{\partial\widetilde{\Phi}(\widetilde{\mathbf{F}},\widetilde{\mu})}{\partial\widetilde{\mathbf{F}}}\quad$
$\displaystyle\text{and}\quad\widetilde{\mathnormal{C}}$
$\displaystyle=-\frac{\partial\widetilde{\Phi}(\widetilde{\mathbf{F}},\widetilde{\mu})}{\partial\widetilde{\mu}}$
(32b)

which yields the same constitutive relation for the surface first Piola-
Kirchhoff $\widetilde{\mathbf{P}}$, but the bulk first Piola-Kirchhoff stress
$\mathbf{P}$ should be modified because the incompressibility condition was
enforced by substituting the Equation 26 into Equation 25 to eliminate the
bulk concentration $\mathnormal{C}$ Hong et al. (2009):

$\displaystyle\mathbf{P}$
$\displaystyle=\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}\left(\mathbf{F}+\alpha_{2}\mathnormal{J}\mathbf{F}^{-\mathrm{T}}\right)\quad\text{with}\quad\alpha_{2}=-\frac{1}{\mathnormal{J}}+\frac{1}{\mathnormal{N}\Omega}\left[\frac{1}{\mathnormal{J}}+\ln\left(\frac{\mathnormal{J}-1}{\mathnormal{J}}\right)+\frac{\chi}{\mathnormal{J}^{2}}-\frac{\mu}{\mathnormal{k}_{B}\mathnormal{T}}\right]$
(33)

Note that the bulk concentration $\mathnormal{C}$ can still be obtained by
incompressibility condition of Equation 26, but the surface concentration
$\widetilde{\mathnormal{C}}$ is now given implicitly by solving a nonlinear
algebraic equation in Equation 29d.

Keeping into account the inability of the FEniCS framework to solve nonlinear
equations at the Gauss-point level, we proceed to solve Equation 29d directly
with the mixed finite element method. The weak form of the problem is obtained
by using a set of test functions, which satisfy the necessary integrability
conditions. By multiplying Equation 13a, Equation 14a and Equation 29d with
the test functions $\delta\mathbf{u}$, $\delta\mu$ and
$\delta\widetilde{\mathnormal{C}}$, and integrating over the domain,
respectively, then we obtain that

$\displaystyle\int_{\mathnormal{V}}\mathbf{P}:\boldsymbol{\nabla}_{\mathbf{X}}\delta\mathbf{u}\,\mathrm{d}V+\int_{\mathnormal{S}}\widetilde{\mathbf{P}}:\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\delta\mathbf{u}\,\mathrm{d}S=0$
(34a)
$\displaystyle\int_{\mathnormal{V}}\dot{\mathnormal{C}}\,\delta\mu\,\mathrm{d}V-\int_{\mathnormal{V}}\mathbf{J}\cdot\boldsymbol{\nabla}_{\mathbf{X}}\delta\mu\,\mathrm{d}V+\int_{\mathnormal{S}}\dot{\widetilde{\mathnormal{C}}}\,\delta\mu\,\mathrm{d}S-\int_{\mathnormal{S}}\widetilde{\mathbf{J}}\cdot\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\delta\mu\,\mathrm{d}S=0$
(34b)
$\displaystyle\int_{\mathnormal{S}}\Bigg{[}\widetilde{\mu}-\mathnormal{k}_{B}\mathnormal{T}\Big{[}\ln\left(\frac{\widetilde{\Omega}\widetilde{\mathnormal{C}}}{1+\widetilde{\Omega}\widetilde{\mathnormal{C}}}\right)+\frac{1}{1+\widetilde{\Omega}\widetilde{\mathnormal{C}}}+\frac{\widetilde{\chi}}{\left(1+\widetilde{\Omega}\widetilde{\mathnormal{C}}\right)^{2}}\Big{]}-\widetilde{\Omega}\,\widetilde{\kappa}\left(\widetilde{J}-1-\widetilde{\Omega}\widetilde{\mathnormal{C}}\right)\Bigg{]}\delta\widetilde{\mathnormal{C}}\,\mathrm{d}S=0$
(34c)

where the derivation of Equation 34b can be found in Appendix A. The statement
of the weak form is to find the trial functions, $\mathbf{u}$, $\mu$ and
$\widetilde{\mathnormal{C}}$, such that the integrals in Equation 34 are
satisfied for any permissible test functions, $\delta\mathbf{u}$, $\delta\mu$
and $\delta\widetilde{\mathnormal{C}}$.

###  4.2 Normalization

For the finite element simulations in the following section, all variables and
parameters are normalized, as denoted by $\\{\widehat{\bullet}\\}$. All
lengths are normalized by a characteristic dimension, $\mathnormal{H}$ (e.g.
the length of the edge of a cube in the reference configuration). The chemical
potential and stresses are normalized as follows,

$\widehat{\mu}=\frac{\mu}{\mathnormal{k}_{B}\mathnormal{T}},\quad\widehat{\mathbf{P}}=\frac{\mathbf{P}}{\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}},\quad\widehat{\widetilde{\mathbf{P}}}=\frac{\widetilde{\mathbf{P}}}{\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}\mathnormal{H}}$
(35)

Following the normalization of the surface stress, surface energy and penalty
coefficient should be normalized in the same way,

$\widehat{\widetilde{\gamma}}=\frac{\widetilde{\gamma}}{\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}\mathnormal{H}},\quad\widehat{\widetilde{\kappa}}=\frac{\widetilde{\kappa}}{\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}\mathnormal{H}}$
(36)

Solvent concentration and time are normalized as follows,

$\widehat{\mathnormal{C}}=\Omega\mathnormal{C},\quad\widehat{\widetilde{\mathnormal{C}}}=\widetilde{\Omega}\widetilde{\mathnormal{C}},\quad\hat{\mathnormal{t}}=\frac{\mathnormal{t}}{\tau}$
(37)

where $\tau=\mathnormal{H}^{2}/\mathnormal{D}$ is the characteristic time
scale of diffusion.

Recall that the shear modulus linearized about the undeformed dry state
($\mathbf{F}=\mathbf{I}$) is
$\mathnormal{G}=\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}$ Bouklas and
Huang (2012). In a homogeneously swollen stress-free state, the deformation
gradient takes the form $\mathbf{F}=\lambda_{0}\mathbf{I}$, with $\lambda_{0}$
the swelling stretch, the shear modulus about that state is defined as follows
Bouklas and Huang (2012):

$\mathnormal{G}_{0}=\frac{\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}}{\lambda_{0}}$
(38)

Considering again a homogeneously swollen state, the elastocapillary length
scale is defined by $\ell=\widetilde{\gamma}/\mathnormal{G}_{0}$ for constant
surface free energy per unit current surface area. Taking into account the
characteristic dimension $\mathnormal{H}$ the normalized elastocapillary
length scale is given as

$\widehat{\ell}=\frac{\ell}{\mathnormal{H}}=\widehat{\widetilde{\gamma}}\lambda_{0}$
(39)

By substituting the normalized variables into the weak forms in Equation 34,
we can obtain the normalized weak forms, for the bulk and the surface.

$\displaystyle\int_{\mathnormal{V}}\widehat{\mathbf{P}}:\boldsymbol{\nabla}_{\widehat{\textbf{X}}}\delta\widehat{\mathbf{u}}\,\mathrm{d}\widehat{V}+\int_{\mathnormal{S}}\widehat{\widetilde{\mathbf{P}}}:\widetilde{\boldsymbol{\nabla}}_{\widehat{\widetilde{\textbf{X}}}}\delta\widehat{\mathbf{u}}\,\mathrm{d}\widehat{S}=0$
(40a)
$\displaystyle\frac{\widetilde{\Omega}\mathnormal{H}}{\Omega}\int_{\mathnormal{V}}\left(\frac{\mathrm{d}\widehat{\mathnormal{C}}}{\mathrm{d}\hat{\mathnormal{t}}}\,\delta\widehat{\mu}-\int_{\mathnormal{V}}\widehat{\mathbf{J}}\cdot\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\delta\widehat{\mu}\right)\mathrm{d}\widehat{V}+\int_{\mathnormal{S}}\left(\frac{\mathrm{d}\widehat{\widetilde{\mathnormal{C}}}}{\mathrm{d}\hat{\mathnormal{t}}}\,\delta\widehat{\widetilde{\mu}}-\int_{\mathnormal{S}}\widehat{\widetilde{\mathbf{J}}}\cdot\widetilde{\boldsymbol{\nabla}}_{\widehat{\widetilde{\textbf{X}}}}\delta\widehat{\widetilde{\mu}}\right)\mathrm{d}\widehat{S}=0$
(40b)
$\displaystyle\int_{\mathnormal{S}}\Bigg{[}\widehat{\widetilde{\mu}}-\ln\left(\frac{\widehat{\widetilde{\mathnormal{C}}}}{1+\widehat{\widetilde{\mathnormal{C}}}}\right)-\frac{1}{1+\widehat{\widetilde{\mathnormal{C}}}}-\frac{\widetilde{\chi}}{\left(1+\widehat{\widetilde{\mathnormal{C}}}\right)^{2}}+\mathnormal{N}\widetilde{\Omega}\mathnormal{H}\,\widehat{\widetilde{\kappa}}\left(\widetilde{J}-1-\widehat{\widetilde{\mathnormal{C}}}\right)\Bigg{]}\delta\widehat{\widetilde{\mathnormal{C}}}\,\mathrm{d}\widehat{S}=0$
(40c)

where $\mathnormal{N}\Omega$ and
$\mathnormal{N}\widetilde{\Omega}\mathnormal{H}$ are the dimensionless
parameters, and the normalized fluxes in bulk and on surface are obtained by

$\displaystyle\widehat{\mathbf{J}}$
$\displaystyle=-\widehat{\mathnormal{C}}\mathbf{C}^{-1}\cdot\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\widehat{\mu}$
(41a) $\displaystyle\widehat{\widetilde{\mathbf{J}}}$
$\displaystyle=-\widehat{\widetilde{\mathnormal{C}}}\frac{\widetilde{\mathnormal{D}}}{\mathnormal{D}}\widetilde{\mathbf{C}}^{-1}\cdot\widetilde{\boldsymbol{\nabla}}_{\widehat{\widetilde{\textbf{X}}}}\widehat{\widetilde{\mu}}$
(41b)

where $\widetilde{\mathnormal{D}}/\mathnormal{D}$ is the ratio of diffusivity
between in bulk and on surface.

###  4.3 Temporal discretization

The backward Euler scheme is used to integrate Equation 40b over time:

$\displaystyle\frac{\widetilde{\Omega}\mathnormal{H}}{\Omega}\int_{\mathnormal{V}}\left[\frac{1}{\Delta\hat{\mathnormal{t}}}\left(\widehat{\mathnormal{C}}{}^{\hat{\mathnormal{t}}+\Delta\hat{\mathnormal{t}}}-\widehat{\mathnormal{C}}{}^{\hat{\mathnormal{t}}}\right)\delta\widehat{\mu}-\int_{\mathnormal{V}}\widehat{\mathbf{J}}{}^{\hat{\mathnormal{t}}+\Delta\hat{\mathnormal{t}}}\cdot\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\delta\widehat{\mu}\right]\mathrm{d}\widehat{V}$
$\displaystyle+\int_{\mathnormal{S}}\left[\frac{1}{\Delta\hat{\mathnormal{t}}}\left(\widehat{\widetilde{\mathnormal{C}}}{}^{\hat{\mathnormal{t}}+\Delta\hat{\mathnormal{t}}}-\widehat{\widetilde{\mathnormal{C}}}{}^{\hat{\mathnormal{t}}}\right)\delta\widehat{\widetilde{\mu}}-\int_{\mathnormal{S}}\widehat{\widetilde{\mathbf{J}}}{}^{\hat{\mathnormal{t}}+\Delta\hat{\mathnormal{t}}}\cdot\widetilde{\boldsymbol{\nabla}}_{\widehat{\widetilde{\textbf{X}}}}\delta\widehat{\widetilde{\mu}}\right]\mathrm{d}\widehat{S}=0$
(42)

where the superscripts indicate the time step, at the current time step
($\hat{\mathnormal{t}}+\Delta\hat{\mathnormal{t}}$) or the previous step
$\hat{\mathnormal{t}}$. We can combine Equation 40a, Equation 40c and Equation
42 as

$\displaystyle\int_{\mathnormal{V}}\widehat{\mathbf{P}}:\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\delta\widehat{\mathbf{u}}\,\mathrm{d}\widehat{V}+\int_{\mathnormal{S}}\widehat{\widetilde{\mathbf{P}}}:\widetilde{\boldsymbol{\nabla}}_{\widehat{\widetilde{\textbf{X}}}}\delta\widehat{\mathbf{u}}\,\mathrm{d}\widehat{S}$
$\displaystyle+\frac{\widetilde{\Omega}\mathnormal{H}}{\Omega}\int_{\mathnormal{V}}\left[\left(\widehat{\mathnormal{C}}-\widehat{\mathnormal{C}}{}^{\hat{\mathnormal{t}}}\right)\delta\widehat{\mu}-\Delta\hat{\mathnormal{t}}\,\widehat{\mathbf{J}}\cdot\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\delta\widehat{\mu}\right]\,\mathrm{d}\widehat{V}+\int_{\mathnormal{S}}\left[\left(\widehat{\widetilde{\mathnormal{C}}}-\widehat{\widetilde{\mathnormal{C}}}{}^{\hat{\mathnormal{t}}}\right)\delta\widehat{\widetilde{\mu}}-\Delta\hat{\mathnormal{t}}\,\widehat{\widetilde{\mathbf{J}}}\cdot\widetilde{\boldsymbol{\nabla}}_{\widehat{\widetilde{\textbf{X}}}}\delta\widehat{\widetilde{\mu}}\right]\,\mathrm{d}\widehat{S}$
$\displaystyle+\int_{\mathnormal{S}}\Bigg{[}\widehat{\widetilde{\mu}}-\ln\left(\frac{\widehat{\widetilde{\mathnormal{C}}}}{1+\widehat{\widetilde{\mathnormal{C}}}}\right)-\frac{1}{1+\widehat{\widetilde{\mathnormal{C}}}}-\frac{\widetilde{\chi}}{\left(1+\widehat{\widetilde{\mathnormal{C}}}\right)^{2}}+\mathnormal{N}\widetilde{\Omega}\mathnormal{H}\,\widehat{\widetilde{\kappa}}\left(\widetilde{J}-1-\widehat{\widetilde{\mathnormal{C}}}\right)\Bigg{]}\delta\widehat{\widetilde{\mathnormal{C}}}\,\mathrm{d}\widehat{S}=0$
(43)

where the superscript ($\hat{\mathnormal{t}}+\Delta\hat{\mathnormal{t}}$) is
omitted for all the terms at the current time step and
$\widehat{\mathnormal{C}}{}^{\hat{\mathnormal{t}}}$ and
$\widehat{\widetilde{\mathnormal{C}}}{}^{\hat{\mathnormal{t}}}$ are the
species concentration at the previous time step in the bulk and on the
surface.

###  4.4 Spatial discretization

A mixed finite element method is utilized to solve for the normalized
displacement, chemical potential and surface concentration fields
concurrently. To avoid the numerical instability with the mixed method, we
should employ proper spatial discretization techniques Bouklas et al. (2015).
We utilize a Taylor-Hood element Taylor and Hood (1973) for the bulk
quantities (displacement and chemical potential), where the interpolation
order for chemical potential (linear) is one-order lower than for the
displacement (quadratic). We note that a linear interpolation is utilized for
the and surface concentration field.

The normalized displacement and chemical potential are interpolated through
the domain of interest as

$\widehat{\mathbf{u}}=\mathbf{H}{}^{\widehat{\mathbf{u}}}\widehat{\mathbf{u}}{}^{n},\quad\widehat{\mu}=\mathbf{H}{}^{\widehat{\mu}}\widehat{\boldsymbol{\mu}}{}^{n},\quad\widehat{\widetilde{\mathnormal{C}}}=\widetilde{\mathbf{H}}{}^{\widehat{\widetilde{\mathnormal{C}}}}\widehat{\widetilde{\boldsymbol{\mathnormal{C}}}}{}^{n}$
(44)

where $\mathbf{H}{}^{\widehat{\mathbf{u}}}$, $\mathbf{H}{}^{\widehat{\mu}}$
and $\widetilde{\mathbf{H}}{}^{\widehat{\widetilde{\mathnormal{C}}}}$ are the
shape functions, $\widehat{\mathbf{u}}{}^{n}$,
$\widehat{\boldsymbol{\mu}}{}^{n}$ and
$\widehat{\widetilde{\boldsymbol{\mathnormal{C}}}}{}^{n}$ are the nodal values
of the normalized displacement, chemical potential and surface concentration,
respectively. Note that the shape function
$\widetilde{\mathbf{H}}{}^{\widehat{\widetilde{\mathnormal{C}}}}$ is only
defined on surface elements. The test functions are discretized in the same
way

$\delta\widehat{\mathbf{u}}=\mathbf{H}{}^{\widehat{\mathbf{u}}}\delta\widehat{\mathbf{u}}{}^{n},\quad\delta\widehat{\mu}=\mathbf{H}{}^{\widehat{\mu}}\delta\widehat{\boldsymbol{\mu}}{}^{n},\quad\delta\widehat{\widetilde{\mathnormal{C}}}=\widetilde{\mathbf{H}}{}^{\widehat{\widetilde{\mathnormal{C}}}}\delta\widehat{\widetilde{\boldsymbol{\mathnormal{C}}}}{}^{n}$
(45)

The stresses, concentrations, and fluxes are evaluated at integration points,
depending on the gradients of the displacement and chemical potential via the
constitutive relations. Taking the gradient of Equation 44, we obtain that

$\displaystyle\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\widehat{\mathbf{u}}$
$\displaystyle=\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\mathbf{H}{}^{\widehat{\mathbf{u}}}\widehat{\mathbf{u}}{}^{n}$
$\displaystyle=\mathbf{B}{}^{\widehat{\mathbf{u}}}\widehat{\mathbf{u}}{}^{n}$
(46a) $\displaystyle\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\widehat{\mu}$
$\displaystyle=\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\mathbf{H}{}^{\widehat{\mu}}\delta\widehat{\boldsymbol{\mu}}{}^{n}$
$\displaystyle=\mathbf{B}{}^{\widehat{\mu}}\widehat{\boldsymbol{\mu}}{}^{n}$
(46b)
$\displaystyle\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\widehat{\widetilde{\mathnormal{C}}}$
$\displaystyle=\boldsymbol{\nabla}_{\widehat{\mathbf{X}}}\widetilde{\mathbf{H}}{}^{\widehat{\widetilde{\mathnormal{C}}}}\widehat{\widetilde{\boldsymbol{\mathnormal{C}}}}{}^{n}$
$\displaystyle=\widetilde{\mathbf{B}}{}^{\widehat{\widetilde{\mathnormal{C}}}}\widehat{\widetilde{\boldsymbol{\mathnormal{C}}}}{}^{n}$
(46c)

where $\mathbf{B}{}^{\widehat{\mathbf{u}}}$ and $\mathbf{B}{}^{\widehat{\mu}}$
are the gradients of the shape functions in the bulk, and
$\widetilde{\mathbf{B}}{}^{\widehat{\widetilde{\mathnormal{C}}}}$ is the one
on the surface.

###  4.5 Nonlinear solution

The weak form in Section 4.3 can be expressed as a system of nonlinear
equations,

$\mathcal{N}(\mathbf{d})=\mathbf{f}\quad\text{with}\quad\mathbf{d}=\left[\widehat{\mathbf{u}}{}^{n}\;\widehat{\boldsymbol{\mu}}{}^{n}\;\widehat{\widetilde{\boldsymbol{\mathnormal{C}}}}{}^{n}\right]^{\mathrm{T}}$
(47)

Note that $\mathcal{N}(\mathbf{d})$ denotes the part of the weak form that is
not known at current time step, and we take all the known quantities to the
right-hand side and denote as $\mathbf{f}$ (which is known from previous time
step). The residual of nonlinear equations in iteration step $\mathnormal{i}$
is given by
$\mathbf{R}_{\mathnormal{i}}=\mathbf{f}-\mathcal{N}(\mathbf{d}_{\mathnormal{i}})$,
which can be solved using the Newton–Raphson method. In particular, the
procedure requires the calculation of the tangent Jacobian matrix at each
iteration, namely,

$\frac{\partial\mathcal{N}}{\partial\mathbf{d}}\Bigg{|}_{\mathbf{d}_{i}}=\begin{bmatrix}\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\mathbf{u}}}&\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\mu}}&\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\widetilde{\mathnormal{C}}}}\\\
\mathbf{K}^{\widehat{\mu}\widehat{\mathbf{u}}}&\mathbf{K}^{\widehat{\mu}\widehat{\mu}}&\mathbf{K}^{\widehat{\mu}\widehat{\widetilde{\mathnormal{C}}}}\\\
\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\mathbf{u}}}&\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\mu}}&\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\widetilde{\mathnormal{C}}}}\end{bmatrix}$
(48)

FEniCS version 2019.2.0 Logg et al. (2012a); Alnæs et al. (2015) is used to
numerically solve the coupled non-linear equations via the Portable Extensible
Toolkit for Scientific Computations (PETSc) Scalable Nonlinear Equations
Solvers (SNES) interface Balay et al. (2019). This process repeats until a
level of convergence specified within the SNES solver. At each iteration, the
block Jacobian matrices are set up using multiphenics Ballarin (2019), a
python library that also facilitates the definition of boundary restricted
variables within FEniCS, a feature necessary for dealing in our case with the
surface concentration. We note that the surface chemical potential is tied
with the definition of the corresponding bulk quantity, but the surface
concentration is not.

###  4.6 FEniCS implementation details

Although many documents Logg et al. (2012b); Langtangen and Logg (2017);
Bleyer (2018) provide FEniCS implementation details, it is still not trivial
to deal with surface kinematics using FEniCS as several surface quantities
have a rank deficiency. For educational purposes, we review a few important
mathematical definitions of surface quantities, and then we provide their
corresponding FEniCS definitions.

The surface unit tensor is defined by
$\widetilde{\textbf{I}}=\textbf{I}-\textbf{N}\otimes\textbf{N}$, which can be
interpreted as surface (idempotent) projection tensor, and surface deformation
gradient can be obtained from surface projection of bulk deformation gradient
$\widetilde{\textbf{F}}=\mathbf{F}\widetilde{\mathbf{I}}$. Similarly, the
surface gradient is obtained by surface projection of the bulk gradient
$\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\mu=\boldsymbol{\nabla}_{\textbf{X}}\mu\widetilde{\textbf{I}}$.

⬇

N = FacetNormal(mesh)

I_bulk = Identity(3)

I_surf = I_bulk - outer(N,N)

F_surf = dot(F_bulk,I_surf)

grad_mu_surf = dot(grad(mu),I_surf)

Note that the surface unit tensor $\widetilde{\mathbf{I}}$ is absent from the
(principal) normal component $\textbf{N}\otimes\textbf{N}$; as a result, the
surface projection always leads to a rank deficiency, and we cannot perform
the inverse operations for the surface quantities. Nevertheless, they possess
an inverse in the generalized sense, i.e.,
$\widetilde{\textbf{F}}^{-1}:=\widetilde{\mathbf{I}}\mathbf{F}^{-1}$ and
$\widetilde{\textbf{C}}^{-1}:=\widetilde{\mathbf{I}}\mathbf{C}^{-1}\widetilde{\mathbf{I}}$.

⬇

F_surf_inv = dot(I_surf,inv(F_bulk))

C_surf_inv = dot(I_surf,dot(C_bulk,I_surf))

In addition, we cannot perform the determinant operations for the surface
quantities due to rank deficiency. We have the alternative way of using
Nanson’s formula to obtain the determinant of surface quantities, i.e.,
$\widetilde{\mathnormal{J}}=\text{det}\,\widetilde{\textbf{F}}:=\lvert\mathrm{cof}\mathbf{F}\mathbf{N}\rvert$.

⬇

A_surf = dot(cofac(F_bulk),N)

J_surf = sqrt(dot(A_surf,A_surf))

##  5 Numerical examples

To study the performance of the suggested framework, we analyze the transient
responses for two initial boundary value problems involving bulk and surface
poroelasticity of hydrogels: (1) free contraction of a cube with smooth edges
and (2) uniaxial tension of a cube with sharp edges. In both cases the edge
length of the cubes is $10^{-3}\,\mathrm{m}$. In the following numerical
study, we take the initial swelling ratio at $\lambda_{0}=3.215$, and the
Flory interaction parameter at $\chi=\widetilde{\chi}=0.2$. At room
temperature $\mathnormal{k}_{B}\mathnormal{T}=4\times
10^{-21}\,\mathrm{J}/\mathrm{mol}$, and the representative value of the volume
and area per molecule are $\Omega=10^{-28}\,\mathrm{m}^{3}$ and
$\widetilde{\Omega}=10^{-19}\,\mathrm{m}^{2}$. In the absence of solvent
molecules, the dry network has shear moduli of
$\mathnormal{N}\mathnormal{k}_{B}\mathnormal{T}=4\times
10^{4}\,\mathrm{N}/\mathrm{m}^{2}$, which gives $\mathnormal{N}\Omega=10^{-3}$
and $\mathnormal{N}\widetilde{\Omega}\mathnormal{H}=10^{3}$ Ang et al. (2020).
The characteristic dimension and time scales are set as
$\mathnormal{H}=10^{-3}\,\mathrm{m}$ and $\tau=1.0\,\mathrm{sec}$. Without
loss of generality, we prescribe the following values for the parameters
$\widetilde{\beta}=1.0$, $\widehat{\widetilde{\gamma}}=1.0$,
$\widehat{\widetilde{\kappa}}=10^{-3}$ and
$\widetilde{\mathnormal{D}}/\mathnormal{D}=1.0$.

For robustness of the numerical procedure, we follow a two-stage process. We
initialize from a homogeneously swollen state, with initial homogeneous
swelling stretch $\lambda_{0}$, neglecting elastic surface energy
contributions. In the first stage, we ramp the surface energy
$\widehat{\widetilde{\gamma}}$ linearly, from zero to its prescribed value,
for the time interval $\mathnormal{t}/\tau\in[0.0,1.0]$. This time interval is
significantly smaller compared to the time the system needs for equilibration.
Then in the second stage, we exponentially increase the time steps
$\Delta\mathnormal{t}/\tau$ until equilibrium is attained, while holding the
value of the surface energy $\widehat{\widetilde{\gamma}}$ fixed.

We denote the free-swelling stretch by $\mathbf{F}=\lambda_{0}\mathbf{I}$,
which relates to the chemical potential $\widehat{\mu}_{0}$ by setting stress
Equation 33 to be zero Hong et al. (2008, 2009),

$\mathnormal{N}\Omega\left(\frac{1}{\lambda_{0}}-\frac{1}{\lambda_{0}^{3}}\right)+\ln\left(1-\frac{1}{\lambda_{0}^{3}}\right)+\frac{1}{\lambda_{0}^{3}}+\frac{\chi}{\lambda_{0}^{6}}=\widehat{\mu}_{0}$
(49)

We can obtain the initial free-swelling stretch $\lambda_{0}=3.215$ by setting
the immersed condition, $\widehat{\mu}_{0}=0.0$ Hong et al. (2008). We can
also obtain the initial normalized surface concentration by numerically
solving the normalized constitutive relation of Equation 29d, which yields the
initial value, $\widehat{\widetilde{\mathnormal{C}}}_{0}\approx 9.3673$. It is
important to note that the boundary of swollen hydrogel is assumed to be
impermeable for the cases studies here, which is formulated by the flux
boundary condition, $\widetilde{\mathbf{J}}_{p}=0$ (see Equation 14). That is,
the species on the surface are not allowed to migrate into the exterior
(environment). We note that for immersed conditions, we would need to assume
local equilibrium instantaneously on the surface which would override the
surface diffusion mechanism.

Figure 2: Temporal sequence of the chemical potential during the free
contraction of cubes. Figure 3: Temporal sequence of the bulk contraction
during the free contraction of cubes. .

To investigate the numerical stability when time step is small
($\Delta\mathnormal{t}\rightarrow 0$), the block Jacobian matrix in Equation
48 can be reduced to

$\frac{\partial\mathcal{N}}{\partial\mathbf{d}}\Bigg{|}_{\mathbf{d}_{i}}=\begin{bmatrix}\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\mathbf{u}}}&\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\mu}}&\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\widetilde{\mathnormal{C}}}}\\\
\mathbf{K}^{\widehat{\mu}\widehat{\mathbf{u}}}&\mathbf{K}^{\widehat{\mu}\widehat{\mu}}&\mathbf{K}^{\widehat{\mu}\widehat{\widetilde{\mathnormal{C}}}}\\\
\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\mathbf{u}}}&\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\mu}}&\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\widetilde{\mathnormal{C}}}}\end{bmatrix}\approx\begin{bmatrix}\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\mathbf{u}}}&\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\mu}}&\mathbf{K}^{\widehat{\mathbf{u}}\widehat{\widetilde{\mathnormal{C}}}}\\\
\mathbf{K}^{\widehat{\mu}\widehat{\mathbf{u}}}&\mathbf{0}&\mathbf{K}^{\widehat{\mu}\widehat{\widetilde{\mathnormal{C}}}}\\\
\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\mathbf{u}}}&\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\mu}}&\mathbf{K}^{\widehat{\widetilde{\mathnormal{C}}}\widehat{\widetilde{\mathnormal{C}}}}\end{bmatrix}$
(50)

where $\mathbf{K}^{\widehat{\mu}\widehat{\mu}}$ is proportional to time step
($\Delta\mathnormal{t}$), and it approaches zero at the short time limit. This
indicates the saddle point problem structure inherently, which is known to
lead to numerical oscillations; however, as the time progresses in the
transient process, the parabolic nature of the diffusion equations regularizes
the problem (see Wan (2003) Chapter 3 for the details of the saddle point
problem of $3\times 3$ block matrices). It is worth noting, that this type of
oscillations for the bulk poroelasticity problem, stemming from the inf-sup
problem are known to be aleviated by several approaches, including the choice
of Taylor-Hood spaces (as shown for the nonlinear problem in Bouklas et al.
(2015)). Here we have also chosen a Taylor-Hood space between the displacement
and chemical potential, but also between the displacement and surface
concentration.

Figure 4: Temporal sequence of the surface concentration during the free
contraction of cubes. Figure 5: Temporal sequence of the surface flux during
the free contraction of cubes. Figure 6: Surface energy drives the species
migration between in bulk and on surface. (a) The numbers of species in bulk
and on surface are tracked with respect to the normalized time. (b) The change
of volume and surface area are tracked with respect to the normalized time.
During the surface energy ramping ($\mathnormal{t}/\tau<10^{0}$), the bulk is
gaining the species, and the volume increases, while the surface is losing the
species and the surface area decreases. While surface energy is being fixed
($\mathnormal{t}/\tau\geq 10^{0}$), the species still migrate from the surface
into the bulk due to chemical potential. At equilibrium
($\mathnormal{t}/\tau\sim 1.6\times 10^{5}$), there is no migration between
the bulk and the surface. We note that the total number of species is constant
at all times, and the volume increase of bulk is identical to the sum of the
volume of species molecules moving into the bulk, i.e.,
$\int_{V}\left(\mathnormal{J}-1\right)\mathrm{d}\widehat{\mathnormal{V}}=\int_{V}\left(\widehat{\mathnormal{C}}-\widehat{\mathnormal{C}}_{0}\right)\mathrm{d}\widehat{\mathnormal{V}}$
due to the incompressibility condition of Equation 26.

###  5.1 Free contraction of a cube with smooth edges

In this example, we first investigate the surface and bulk poroelastic effects
on a hydrogel cube. The cube is considered to be cast in a a fully swollen
state, and then released in an environment in which it is free to deform, but
with which it cannot exchange solvent ($\widetilde{\mathbf{J}}_{p}=0$). The
two-stage solution procedure detailed above is followed here.

We plot the temporal sequence of the finite element simulation for the
chemical potential, bulk concentration, surface concentration, and surface
flux in Figures 2, 3, 4 and 5, respectively, where all the values are
normalized. We remove a quarter of the domain from the images of chemical
potential in Figure 2 and bulk concentration in Figure 3 to show the contour
plot on two sections of the interior along with the field on the surface. The
surface flux in Figure 5 is shown on only one face of the cube for
visualization purposes. Images are taken at normalized time
$\mathnormal{t}/\tau=0.0$, $1.0$, $8.4$, $7.8\times 10^{4}$ in the clockwise,
which are denoted by Step 1, 2, 3 and 4 in the following discussion. In Step 1
of Figures 2, 3, 4 and 5, the cubes are at the initial state by a homogeneous
swelling stretch, $\lambda_{0}=3.215$. It is noteworthy that the cube is taken
to have smooth edges (a case with sharp edges will be discussed in the
following example).

Step 2 of Figures 2, 3, 4 and 5 corresponding to the time steps that the
linearly ramping of the surface energy $\widehat{\widetilde{\gamma}}$ are
included. This is for time $\mathnormal{t}/\tau=1.0$ which is significantly
smaller (4 orders of magnitude) compared to the equilibration time. We note
that the current interpolation scheme leads to unstable results if the ramping
time is significantly shorter, pointing to the need for the development of a
more advanced stabilization scheme in the future. This pathology might be a
consequence of the three-field formulation. As the surface energy increases,
the elastocapillary length $\widehat{\ell}$ becomes larger than the
characteristic dimension $\mathnormal{H}$; as a result, the surface effects
become dominant that the geometry becomes almost spherical. In Step 2 of
Figure 2, the chemical potential is heterogeneous due to the deformation
driven by surface energy. In Step 2 of Figures 3 and 4, the distribution of
the bulk concentration on boundary is distinct from the surface concentration,
and has not evolved as rapidly. This is owing to the fact that the species
migration in bulk is slower than the one on surface, a specific modeling
choice to highlight the lack of volumetric constraints for surface diffusion.
However, in the equilibrium state, as shown in Step 4 of Figures 3 and 4,
their distributions are similar on the surface.

On top of surface and bulk diffusion, in this system there is an exchange of
species between the surface and the bulk, as captured in the strong form of
the balance law for species conservation Equation 14. The species exchange
between bulk and surface is easily tracked by accounting for the molecular
incompressiblity condition of Equation 26, indicating that the change of
volume of the bulk is the same as the total volume of the species migrating
from bulk to surface. In Figure 6, we track the normalized total species
population in bulk and on surface at all time steps by integrating the
normalized concentrations over the bulk and the surface accordingly. We can
observe the species migration between $\mathnormal{t}/\tau=0.0$ and $7.8\times
10^{4}$ while the total number of species is conserved. In Step 2 of Figure 4,
the surface concentration is high near the smoothened vertex of the cube
because the deformation at that location is high compared to the facet, which
subsequently leads to higher chemical potential and species migration. In Step
2 of Figure 5, the surface diffusion is mainly directed from the edges to the
interior of the facet. The direction of the surface diffusion is governed by
the surface deformation and chemical potential (see Equation 41).

Beyond Step 2, and for the rest of the transient response the surface energy
($\widehat{\widetilde{\gamma}}=1$) is maintained at the same level, and we
exponentially ramp the time steps towards equilibrium. The deformation of the
cubes is small beyond the initial surface energy ramping, but the
concentration (especially in the bulk) significantly changes, driven by
gradients of the chemical potential. In Step 3 of Figure 3 and Figure 4, we
can see the diffusion in bulk and on surface, and the species migrate toward
the interior of the facet of the cubes. In Step 4 of Figure 2, the simulation
reaches the equilibrium state, which is confirmed by the homogeneous chemical
potential. the distribution of the bulk concentration on the boundary is
similar to the surface concentration, as shown in Step 4 of Figure 3 and
Figure 4. In Figure 6, we can observe the significant migration between bulk
and surface during ramping the surface energy ($\mathnormal{t}/\tau<1$), but
the migration between bulk and surface is negligible with surface energy being
fixed ($1<\mathnormal{t}/\tau\sim 7.8\times 10^{4}$).

Figure 7: Temporal sequence of the chemical potential during the tension of
cubes. Figure 8: Temporal sequence of the bulk contraction during the tension
of cubes. Figure 9: Temporal sequence of the surface concentration during the
tension of cubes. Figure 10: Temporal sequence of the surface flux during the
tension of cubes.

###  5.2 Uniaxial tension of a cube with sharp edges

In this example, we first investigate the surface and bulk poroelastic effects
on a hydrogel cube subject to external load. Similar to the previous example
the cube is considered to be cast in a a fully swollen state, and released in
an environment with which it cannot exchange solvent, but in this case it is
subject to displacement controlled loading. The two-stage solution procedure
detailed above is followed here as well. The lateral surfaces at
$\mathnormal{x}(\mathnormal{t}=0)=\pm\mathnormal{H}/2$ are the clamped
boundaries, and the others are free. In this case, the displacement boundary
conditions at lateral surfaces, are ramped linearly from zero in the same
timeframe that ramping of the surface energy takes place.

In this case external mechanical loading along with surface and bulk energy
contributions together dictate the transient response of the system. The other
conditions are the same as in the previous example of Section 5.1. We note
that for this example, we have taken a cube geometry with sharp corners. We
plot the temporal sequence of finite element simulations for the chemical
potential, bulk concentration, surface concentration, and surface flux in
Figures 7, 8, 9 and 10, respectively, where all the values are normalized. The
surface flux in Figure 10 is shown on only one face of the cube for
visualization purposes. Images are taken at normalized time
$\mathnormal{t}/\tau=0.0$, $1.0$, $5.4\times 10^{2}$, $4.1\times 10^{6}$ in
the clockwise, which are denoted by Step 1, 2, 3 and 4 in the following
discussion.

In Step 1 of Figures 7, 8, 9 and 10, the cubes are at the initial state by a
homogeneous swelling stretch, $\lambda_{0}=3.215$. In Step 2 of Figures 7, 8,
9 and 10, the cubes are deformed by the surface energy and stretch, after the
initial load and surface energy ramping. The cubes are initially clamped at
lateral sides, and we linearly ramp the surface energy
$\widehat{\widetilde{\gamma}}=1$ and stretch $\varepsilon=68.9\%$ between
$\mathnormal{t}/\tau=0.0$ and $1.0$. In Step 3 of Figures 7, 8, 9 and 10, we
maintain the surface energy and the stretch, and we exponentially ramp the
time-step size towards equilibrium. In Step 4 of Figures 7, 8, 9 and 10, the
simulation approaches the equilibrium state. Through the second example, we
can observe that surface concentration is distinct from the bulk concentration
evaluated on the boundary.

Unlike the previous example in Section 5.1, the cube at initial state has the
sharp corners to investigate the impact of sharp geometric features to the
numerical scheme. We can observe an oscillation of the surface concentration
near the edges of the cube in Figure 9, which is partially regularized with
time, reminiscent of the enforcement of an initial condition in a heat
conduction problem. This originates from the weakly enforced continuity
condition in Equation 14d, which states that the net surface flux on the
boundary curve (between the adjacent surface elements) is always balanced. Of
course, this can easily be resolved (as in the previous example), where sharp
features can be smoothened at the expense of computational cost owing to a
finer non-structured mesh. Alternative, one could embark on deriving an
appropriate stabilization scheme for this. We do not consider the impact of
these initial oscillations to be significant to the transient response,
similar to imposing discontinuous chemical potential boundary (with respect to
the initial condition) for swelling of hydrogels Bouklas et al. (2015). Again,
this might be specific to the three-field formulation, and a two-field
formulation could also alleviate the issue.

##  6 Conclusion

A continuum multiphysics formulation and corresponding finite element
implementation has been presented to account for bulk and surface nonlinear
poroelasticity in hydrogels. The governing equations for the response of a
continuum body have been presented, and the general constitutive relations are
derived in a thermodynamically-consistent manner, which is subsequently
specialized for hydrogels. The proposed nonlinear theory is numerically solved
through the mixed finite element method by employing the open-source finite
element framework FEniCS. In addition, we provide FEniCS implementation
details that could be important for problems where surface kinematics are of
interest.

Two numerical examples are investigated to understand the effect of coupling
bulk and surface poroelasticity. In the first example, where a cube with
smoothened geometric features is investigated (see Figures 2, 3, 4, 5 and 6),
we probe the variation of concentration, surface flux, and change of volume
and surface area along with the chemical potential, which follow accounting
for a fluid-like surface energy for the soft solid. As a result of the
transient procedure and species exchange between bulk and surface, a smooth
cube-like object is gradually transformed into an almost spherical one. From
this observation, one may expect that the volume of the cube is decreased.
Interestingly, the result is the opposite, as shown in Figure 6. It implies
that the species collectively migrate from the surface into the bulk because
the volume increase of bulk is identical to the volume of total species moving
into the bulk; consequently, the volume of the cube rather increases in spite
of the surface area decreasing. This result showcases that the multiphysical
complication could play an important role in understanding the response of
soft solids in length scales that surface effects dominate.

In the second example (see Figures 7, 8, 9 and 10), we further investigate the
response of a cube with sharp features, under a external load, while still
accounting for a fluid-like surface energy. Although the distributions of
species concentration in bulk and on the surface are similar on the boundary
in the first example (see Figures 3 and 4), the second example showcases that
the pathways of species migrations can be significantly different (see Figures
8 and 9). In this example, we choose the same mobility constant in bulk and on
the surface as the simplest case, but this result still provides an important
insight into hydrogel-based applications at small scales. For example, one may
use the distinct migration mechanism for the advanced design of hydrogel-based
sensors and actuators.

Even though we focus on hydrogels in this work, the general framework and
finite element implementation developed here can allow studying several
problems for solids as well in mechanobiology, where diffusion processes in
the bulk and on the surface are coupled with the elastic response. One may
adopt our numerical framework, but different constitutive relations must be
specialized for studying other materials. One specific example that the
authors plan to pursue is the modeling of morphogenic processes in tissue
mechanics and, more specifically modeling of contractile microtissues Mailand
et al. (2019); Kim et al. (2020, 2023), where cells apply forces and contract
the extracellular matrix (ECM), but at the same time can diffuse within the
bulk of the ECM but also on its periphery, having significant implications in
problems like wound healing. In this paper, we report numerical issues that
can arise due to several reasons and propose approaches to alleviate them. For
the authors’ future work, it would also be interesting to extend our model to
consider growth mechanisms to account for surface morphology in evolving
natural systems, such as encountered during tumor growth.

## Acknowledgments

JK and NB acknowledge the support by the National Science Foundation under
grant no. CMMI-2129776. CYH acknowledges the support by the National Science
Foundation under grant no. CMMI-1903308. FB thanks the project “Numerical
modeling of flows in porous media” funded by Università Cattolica del Sacro
Cuore, and the INDAM-GNCS project “Metodi numerici per lo studio di strutture
geometriche parametriche complesse” (CUP_E53C22001930001, PI Dr. Maria
Strazzullo).

## Declarations

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##  Appendix A Derivation of weak form of species conservation

Staring form Equation 14a,

$\int_{\mathnormal{V}}\dot{\mathnormal{C}}\,\delta\mu\,\mathrm{d}V+\int_{\mathnormal{V}}\left(\boldsymbol{\nabla}_{\mathbf{X}}\cdot\mathbf{J}\right)\delta\mu\,\mathrm{d}V=0$
(51)

By applying the product rule and divergence theorem Holzapfel (2000), and then
substituting Equation 14b into Equation 51,

$\int_{\mathnormal{V}}\dot{\mathnormal{C}}\,\delta\mu\,\mathrm{d}V-\int_{\mathnormal{V}}\mathbf{J}\cdot\boldsymbol{\nabla}_{\mathbf{X}}\delta\mu\,\mathrm{d}V+\int_{\mathnormal{S}}\dot{\widetilde{\mathnormal{C}}}\,\delta\mu\,\mathrm{d}S+\int_{\mathnormal{S}}\left(\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\cdot\widetilde{\mathbf{J}}\right)\delta\mu\,\mathrm{d}S=0$
(52)

where the last term in Equation 52 can be rewritten by product rule and
surface divergence theorem,

$\int_{\mathnormal{S}}\left(\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\cdot\widetilde{\mathbf{J}}\right)\delta\mu\,\mathrm{d}S=-\int_{\mathnormal{S}}\widetilde{\mathbf{J}}\cdot\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\delta\mu\,\mathrm{d}S+\int_{\mathnormal{L}}\left(\widetilde{\mathbf{J}}\cdot\widetilde{\mathbf{N}}\right)\delta\mu\,\mathrm{d}L$
(53)

where the second term in the right-hand side of Equation 53 is assumed to be
zero. By substituting Equation 53 into Equation 52,

$\int_{\mathnormal{V}}\dot{\mathnormal{C}}\,\delta\mu\,\mathrm{d}V-\int_{\mathnormal{V}}\mathbf{J}\cdot\boldsymbol{\nabla}_{\mathbf{X}}\delta\mu\,\mathrm{d}V+\int_{\mathnormal{S}}\dot{\widetilde{\mathnormal{C}}}\,\delta\mu\,\mathrm{d}S-\int_{\mathnormal{S}}\widetilde{\mathbf{J}}\cdot\widetilde{\boldsymbol{\nabla}}_{\widetilde{\textbf{X}}}\delta\mu\,\mathrm{d}S=0$
(54)

This equation corresponds to Equation 34b in the manuscript.

## References

  * Biot [1955] Maurice A Biot.  Theory of elasticity and consolidation for a porous anisotropic solid.  _Journal of applied physics_ , 26(2):182–185, 1955.
  * Biot [1956] Maurice A Biot.  Theory of propagation of elastic waves in a fluid-saturated porous solid. ii. higher frequency range.  _The Journal of the acoustical Society of america_ , 28(2):179–191, 1956.
  * Biot [1962] Maurice A Biot.  Mechanics of deformation and acoustic propagation in porous media.  _Journal of applied physics_ , 33(4):1482–1498, 1962.
  * Hong et al. [2008] Wei Hong, Xuanhe Zhao, Jinxiong Zhou, and Zhigang Suo.  A theory of coupled diffusion and large deformation in polymeric gels.  _Journal of the Mechanics and Physics of Solids_ , 56(5):1779–1793, 2008.
  * Hong et al. [2015] Sungmin Hong, Dalton Sycks, Hon Fai Chan, Shaoting Lin, Gabriel P Lopez, Farshid Guilak, Kam W Leong, and Xuanhe Zhao.  3d printing of highly stretchable and tough hydrogels into complex, cellularized structures.  _Advanced materials_ , 27(27):4035–4040, 2015\.
  * Li and Mooney [2016] Jianyu Li and David J Mooney.  Designing hydrogels for controlled drug delivery.  _Nature Reviews Materials_ , 1(12):1–17, 2016\.
  * Chester and Anand [2010] Shawn A Chester and Lallit Anand.  A coupled theory of fluid permeation and large deformations for elastomeric materials.  _Journal of the Mechanics and Physics of Solids_ , 58(11):1879–1906, 2010.
  * Yoon et al. [2010] Jinhwan Yoon, Shengqiang Cai, Zhigang Suo, and Ryan C Hayward.  Poroelastic swelling kinetics of thin hydrogel layers: comparison of theory and experiment.  _Soft Matter_ , 6(23):6004–6012, 2010.
  * Li et al. [2012] Bo Li, Yan-Ping Cao, Xi-Qiao Feng, and Huajian Gao.  Mechanics of morphological instabilities and surface wrinkling in soft materials: a review.  _Soft Matter_ , 8(21):5728–5745, 2012.
  * Bouklas and Huang [2012] Nikolaos Bouklas and Rui Huang.  Swelling kinetics of polymer gels: comparison of linear and nonlinear theories.  _Soft Matter_ , 8(31):8194–8203, 2012.
  * Bouklas et al. [2015] Nikolaos Bouklas, Chad M Landis, and Rui Huang.  A nonlinear, transient finite element method for coupled solvent diffusion and large deformation of hydrogels.  _Journal of the Mechanics and Physics of Solids_ , 79:21–43, 2015.
  * Ang et al. [2020] Ida Ang, Zezhou Liu, Jaemin Kim, Chung-Yuen Hui, and Nikolaos Bouklas.  Effect of elastocapillarity on the swelling kinetics of hydrogels.  _Journal of the Mechanics and Physics of Solids_ , page 104132, 2020\.
  * Gurtin and Murdoch [1975] Morton E Gurtin and A Ian Murdoch.  A continuum theory of elastic material surfaces.  _Archive for rational mechanics and analysis_ , 57(4):291–323, 1975.
  * Gurtin and Murdoch [1978] Morton E. Gurtin and Ian A. Murdoch.  Surface stress in solids.  _International Journal of Solids and Structures_ , 14(6):431 – 440, 1978.  ISSN 0020-7683.
  * Gurtin and Jabbour [2002] Morton E Gurtin and Michel E Jabbour.  Interface evolution in three dimensions with curvature-dependent energy and surface diffusion: Interface-controlled evolution, phase transitions, epitaxial growth of elastic films.  _Archive for rational mechanics and analysis_ , 163:171–208, 2002.
  * Cermelli et al. [2005] Paolo Cermelli, Eliot Fried, and Morton E Gurtin.  Transport relations for surface integrals arising in the formulation of balance laws for evolving fluid interfaces.  _Journal of Fluid Mechanics_ , 544:339–351, 2005.
  * Dingreville et al. [2005] Rémi Dingreville, Jianmin Qu, and Mohammed Cherkaoui.  Surface free energy and its effect on the elastic behavior of nano-sized particles, wires and films.  _Journal of the Mechanics and Physics of Solids_ , 53(8):1827–1854, 2005.
  * Steinmann [2008] Paul Steinmann.  On boundary potential energies in deformational and configurational mechanics.  _Journal of the Mechanics and Physics of Solids_ , 56(3):772–800, 2008.
  * Henann and Bertoldi [2014] David L Henann and Katia Bertoldi.  Modeling of elasto-capillary phenomena.  _Soft Matter_ , 10(5):709–717, 2014.
  * Javili et al. [2014] A Javili, A McBride, P Steinmann, and BD Reddy.  A unified computational framework for bulk and surface elasticity theory: a curvilinear-coordinate-based finite element methodology.  _Computational Mechanics_ , 54(3):745–762, 2014\.
  * Liu et al. [2020] Zezhou Liu, Anand Jagota, and Chung-Yuen Hui.  Modeling of surface mechanical behaviors of soft elastic solids: theory and examples.  _Soft Matter_ , 16(29):6875–6889, 2020.
  * Style et al. [2017] Robert W Style, Anand Jagota, Chung-Yuen Hui, and Eric R Dufresne.  Elastocapillarity: Surface tension and the mechanics of soft solids.  _Annual Review of Condensed Matter Physics_ , 8:99–118, 2017.
  * Bico et al. [2018] José Bico, Étienne Reyssat, and Benoît Roman.  Elastocapillarity: When surface tension deforms elastic solids.  _Annual Review of Fluid Mechanics_ , 50:629–659, 2018.
  * Javili and Steinmann [2009] Ali Javili and Paul Steinmann.  A finite element framework for continua with boundary energies. part i: the two-dimensional case.  _Computer Methods in Applied Mechanics and Engineering_ , 198(27-29):2198–2208, 2009.
  * Javili and Steinmann [2010] A Javili and Paul Steinmann.  A finite element framework for continua with boundary energies. part ii: The three-dimensional case.  _Computer Methods in Applied Mechanics and Engineering_ , 199(9-12):755–765, 2010.
  * McBride et al. [2011] AT McBride, A Javili, P Steinmann, and Swantje Bargmann.  Geometrically nonlinear continuum thermomechanics with surface energies coupled to diffusion.  _Journal of the Mechanics and Physics of Solids_ , 59(10):2116–2133, 2011.
  * Javili et al. [2013] Ali Javili, Andrew McBride, and Paul Steinmann.  Thermomechanics of solids with lower-dimensional energetics: on the importance of surface, interface, and curve structures at the nanoscale. a unifying review.  _Applied Mechanics Reviews_ , 65(1), 2013.
  * Mailand et al. [2019] Erik Mailand, Bin Li, Jeroen Eyckmans, Nikolaos Bouklas, and Mahmut Selman Sakar.  Surface and bulk stresses drive morphological changes in fibrous microtissues.  _Biophysical journal_ , 117(5):975–986, 2019\.
  * Kim et al. [2020] Jaemin Kim, Erik Mailand, Ida Ang, Mahmut Selman Sakar, and Nikolaos Bouklas.  A model for 3d deformation and reconstruction of contractile microtissues.  _Soft Matter_ , 2020.
  * Kim et al. [2023] Jaemin Kim, Erik Mailand, Mahmut Selman Sakar, and Nikolaos Bouklas.  A model for mechanosensitive cell migration in dynamically morphing soft tissues.  _Extreme Mechanics Letters_ , 58:101926, 2023.
  * Rastogi and Dortdivanlioglu [2022] Animesh Rastogi and Berkin Dortdivanlioglu.  Modeling curvature-resisting material surfaces with isogeometric analysis.  _Computer Methods in Applied Mechanics and Engineering_ , 401:115649, 2022.
  * Babuška [1971] Ivo Babuška.  Error-bounds for finite element method.  _Numerische Mathematik_ , 16(4):322–333, 1971\.
  * Brezzi [1974] Franco Brezzi.  On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers.  _Publications mathématiques et informatique de Rennes_ , (S4):1–26, 1974.
  * Murad and Loula [1994] Márcio A Murad and Abimael FD Loula.  On stability and convergence of finite element approximations of biot’s consolidation problem.  _International Journal for Numerical Methods in Engineering_ , 37(4):645–667, 1994.
  * Bathe [2001] Klaus-Jürgen Bathe.  The inf–sup condition and its evaluation for mixed finite element methods.  _Computers & structures_, 79(2):243–252, 2001\.
  * Zhang et al. [2009] Jiaping Zhang, Xuanhe Zhao, Zhigang Suo, and Hanqing Jiang.  A finite element method for transient analysis of concurrent large deformation and mass transport in gels.  _Journal of Applied Physics_ , 105(9):093522, 2009\.
  * Duda et al. [2010] Fernando P Duda, Angela C Souza, and Eliot Fried.  A theory for species migration in a finitely strained solid with application to polymer network swelling.  _Journal of the Mechanics and Physics of Solids_ , 58(4):515–529, 2010.
  * Lucantonio et al. [2013] Alessandro Lucantonio, Paola Nardinocchi, and Luciano Teresi.  Transient analysis of swelling-induced large deformations in polymer gels.  _Journal of the Mechanics and Physics of Solids_ , 61(1):205–218, 2013.
  * Chester et al. [2015] Shawn A Chester, Claudio V Di Leo, and Lallit Anand.  A finite element implementation of a coupled diffusion-deformation theory for elastomeric gels.  _International Journal of Solids and Structures_ , 52:1–18, 2015.
  * MacMinn et al. [2016] Christopher W MacMinn, Eric R Dufresne, and John S Wettlaufer.  Large deformations of a soft porous material.  _Physical Review Applied_ , 5(4):044020, 2016\.
  * Leronni and Bardella [2021] Alessandro Leronni and Lorenzo Bardella.  Modeling actuation and sensing in ionic polymer metal composites by electrochemo-poromechanics.  _Journal of the Mechanics and Physics of Solids_ , 148:104292, 2021.
  * Green and Zerna [1992] Albert Edward Green and Wolfgang Zerna.  _Theoretical elasticity_.  Courier Corporation, 1992.
  * Do Carmo [2016] Manfredo P Do Carmo.  _Differential geometry of curves and surfaces: revised and updated second edition_.  Courier Dover Publications, 2016.
  * Holzapfel [2000] Gerhard A Holzapfel.  _Nonlinear solid mechanics II_.  John Wiley & Sons, Inc., 2000.
  * Gurtin et al. [2010] Morton E Gurtin, Eliot Fried, and Lallit Anand.  _The mechanics and thermodynamics of continua_.  Cambridge University Press, 2010.
  * Flory [1942] Paul J Flory.  Thermodynamics of high polymer solutions.  _The Journal of chemical physics_ , 10(1):51–61, 1942.
  * Flory [1961] PJ Flory.  Thermodynamic relations for high elastic materials.  _Transactions of the Faraday Society_ , 57:829–838, 1961\.
  * Treloar [1975] LR G Treloar.  The physics of rubber elasticity.  1975\.
  * Hong et al. [2009] Wei Hong, Zishun Liu, and Zhigang Suo.  Inhomogeneous swelling of a gel in equilibrium with a solvent and mechanical load.  _International Journal of Solids and Structures_ , 46(17):3282–3289, 2009.
  * Huggins [1941] Maurice L Huggins.  Solutions of long chain compounds.  _The Journal of chemical physics_ , 9(5):440–440, 1941.
  * Logg et al. [2012a] Anders Logg, Kent-Andre Mardal, and Garth N. Wells.  _Automated Solution of Differential Equations by the Finite Element Method_.  Springer, 2012a.
  * Alnæs et al. [2015] Martin S. Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells.  The FEniCS project version 1.5.  _Archive of Numerical Software_ , 3(100), 2015.
  * Taylor and Hood [1973] Cedric Taylor and Paul Hood.  A numerical solution of the navier-stokes equations using the finite element technique.  _Computers & Fluids_, 1(1):73–100, 1973.
  * Balay et al. [2019] Satish Balay, Shrirang Abhyankar, Mark Adams, Jed Brown, Peter Brune, Kris Buschelman, Lisandro Dalcin, Alp Dener, Victor Eijkhout, W Gropp, et al.  PETSc users manual.  2019\.
  * Ballarin [2019] F Ballarin.  multiphenics – easy prototyping of multiphysics problems in FEniCS.  https://github.com/multiphenics/multiphenics, 2019.
  * Logg et al. [2012b] Anders Logg, Kent-Andre Mardal, and Garth Wells.  _Automated solution of differential equations by the finite element method: The FEniCS book_ , volume 84.  Springer Science & Business Media, 2012b.
  * Langtangen and Logg [2017] Hans Petter Langtangen and Anders Logg.  _Solving PDEs in python: the FEniCS tutorial I_.  Springer Nature, 2017.
  * Bleyer [2018] Jeremy Bleyer.  Numerical tours of computational mechanics with fenics.  _Zenodo_ , 2018.
  * Wan [2003] Jing Wan.  _Stabilized finite element methods for coupled geomechanics and multiphase flow_.  Stanford university, 2003.