Swelling-induced pattern transformations of periodic hydrogels – from the
wrinkling of internal surfaces to the buckling of thin films

E. Polukhov, L. Pytel & M.-A. Keip∗∗∗Corresponding author: marc-
<EMAIL_ADDRESS>Phone: +49 711 685 66233, Fax: +49 711
685 66347

Institute of Applied Mechanics

Department of Civil and Environmental Engineering

University of Stuttgart, Stuttgart, Germany

Abstract. We investigate pattern transformations of periodic hydrogel systems
that are triggered by swelling-induced structural instabilities. The types of
microstructures considered in the present work include single-phase and two-
phase voided hydrogel structures as well as reinforced hydrogel thin films.
While the observed transformations of the single-phase structures show good
agreement with experimental findings, the two-phase materials provide novel
patterns associated with wrinkling of internal surfaces. Furthermore, an
extensive parametric study on the reinforced hydrogel thin films reveals new
opportunities for the design of complex out-of-plane surface modes caused by
swelling-induced instabilities. Next to the mentioned buckling-type
instabilities, we encountered the development of micro-creases at the internal
surfaces of periodic media before the loss of strong ellipticity of effective
moduli.

Keywords. periodic hydrogel structures, pattern transformations, buckling,
Bloch-Floquet analysis, swelling-induced instabilities, strong ellipticity

##  1 Introduction

Hydrogels are polymeric, soft and biocompatible materials that undergo large
deformations under swelling [33]. At swollen state, their mechanical response
depends on both the properties of the polymer and the solvent phase. While the
elastic stiffness correlates with the density of cross-links of the polymer
network, viscoelasticity is driven by the viscosity of the solvent [1].
Depending on how the polymer network is held together, hydrogels can be
categorized as physical (reversible) or chemical (permanent), see [68] and
[33] for specific examples related to each category. In the former group, the
network is interconnected by polymer-chain entanglements as well as possible
additional forces resulting from ionic or hydrophobic interactions and
hydrogen bonding [83]. The latter group of hydrogels is based on covalently
cross-linked networks. Both physical and chemical hydrogels can be
inhomogeneous due to clustered regions of cross-links and swelling, as well as
the existence of defects in the polymer network. Furthermore, hydrogels can be
classified as conventional or (multi-)stimuli responsive [33]. While both are
hydrophilic and undergo swelling in a solvent, (multi-)stimuli responsive
hydrogels may also respond to changes in pH, temperature or electric field
[32]. In the present work, we will focus on conventional hydrogels and
describe their chemo-mechanical response related to the diffusion of a solvent
through their polymer network.

###  1.1 Engineering and biomedical applications of hydrogels

Since hydrogels have favorable properties and show tunable, versatile
response, they have a wide range of applications, such as fertilizer
encapsulaters in agriculture [69], soft actuators [40] in engineering as well
as artificial muscles [45], tissue-engineered auricles [4], contact lenses and
wound dressings [8] in biomedicine, see also [68, 5, 33] and the references
therein. By virtue of intricate fabrication techniques, further applications
are being realized not only at macroscopic scale, but also at micro- and
nanoscopic scales [53, 70]. Here, we mention in particular the fabrication by
means of photolithography, which has attracted increased attention for the
manufacturing of patterned and composite hydrogels exhibiting controllable
deformations [53, 51]. We refer to [77, 52] for the use of photolithography to
fabricate composite periodic hydrogel structures composed of high-swelling and
non-swelling regions. By means of the latter, the authors were able to induce
out-of-plane buckling in water yielding various periodic buckling patterns.
Consequently, by designing associated hydrogel microstructures and by tuning
material parameters, a rich set of deformation modes can be activated which
increases their potential application as soft devices exploiting chemo-
mechanical interactions.

Besides pattern transforming instabilities in periodic systems, surface
wrinkling and creasing are observed in hydrogel structures. These
instabilities can be exploited to further advance engineering and biomedical
applications. In particular, wrinkling of hydrogels can be utilized for
controlled formation of microgears [81], generation of multicellular spheroids
[86], controlled cell spreading [41] among many other possibilities, see also
[49, 26] and the detailed review [15] on instabilities in gels. We further
refer to the recent experimental realization of nano-wrinkled architectures by
means of laser direct assembly [20].

###  1.2 Continuum-mechanical models for hydrogels

Not least because of the above mentioned favorable properties and
applications, hydrogels have gained increased attention in theoretical and
computational research. Here, one main goal is to develop versatile and
reliable models to simulate their response. Related models can give further
insights into underlying features, provide access to the design of pertinent
microstructures for specific applications and allow the investigation of the
influence of various parameters on the respective behavior. In particular
continuum-mechanical models are established and widely used to capture the
response of hydrogels at macroscopic as well as microscopic scales. These
models take into account chemo-mechanical interactions within the polymer
network of a hydrogel, as well as between the polymer network and the
diffusing solvent. We refer to [85, 36] for continuum-mechanical models of
gels based on Flory–Huggins theory [21, 38], see also the Flory–Rehner type
models [22, 23]. Finite-element implementations of the models [85, 36] have
been applied to the simulation of transient diffusion processes in [82] and
drying-induced instabilities in [35]. We further refer to the continuum
approaches addressed in [11, 12].

Note that the finite-element formulation of the above mentioned models has
been based on the chemical potential as a primary unknown. In combination with
the mechanical deformation map, this corresponds to a saddle-point
formulation. Because of the a priori indefiniteness of the underlying system
of equations, related implementations could result in numerical challenges
associated with _inf-sup_ _instabilities_ [74]. Such challenges can be
overcome by employing a minimization-based variational formulation as proposed
by [7] for hydrogels; see [55] for a basis rooted in poromechanics.

A further similarity of the above mentioned approaches is the use of Biot’s
theory [6, 13] relying on a Flory–Rehner-type energy function and Fick’s law
of diffusion. Alternative continuum-mechanical models have been established on
the basis of the so-called Theory of Porous Media [17, 18, 19]. We further
refer to recent multiscale modeling approaches to describe the diffusion
process in heterogeneous microstructures [50, 44, 65, 64, 43].

###  1.3 Instabilities of hydrogels and their prediction

Figure 1: Examples of swelling-induced instabilities in constrained hydrogel
structures.

Continuum-mechanical models of hydrogels have further been exploited to
investigate complex instability phenomena such as surface wrinkling, creasing
and pattern transformations as a result of deformation-diffusion processes
[15] (see Fig. 1 for an illustration). Experimental investigations of surface
instabilities have been documented, for example, in [73, 76, 28]. In [28],
wrinkling and creasing-type instabilities have been studied in pHEMA
gels111The acronym pHEMA stands for poly(2-hydroxyethyl methacrylate) [28]. ,
where it has been shown that the resulting surface morphology strongly depends
on cross-linker concentration and solvent-polymer interaction. In particular,
the authors demonstrated that while a specific hydrogel layer develops
wrinkling in water, the same layer shows a transition from wrinkling to
creasing when exposed to alcohol and alcohol-water mixture. In [26], surface
wrinkling was exploited to control stem-cell morphology and differentiation.
Further studies on wrinkling and creasing of hydrogel bilayers have been
performed analytically in [47, 46, 79] and computationally in [16, 39, 72]. In
contrast to analytical investigations, computational models usually allow for
the examination of instabilities within a transient framework. Furthermore,
they are not limited with respect to the complexity of underlying geometries.

In addition to surface instabilities, pattern-transforming instabilities in
periodic hydrogel structures have been investigated experimentally in [84, 87,
78]. In [84], the authors have shown that when PDMS gels222The acronym PDMS
stands for poly(dimethyl siloxane) [84].  with periodically embedded circular
pores of 1 $\mu\text{m}$ diameter are swollen in toluene, an organic solvent,
elastic buckling of circular pores into elliptic shapes is observed. In the
buckled state, the major axes of the ellipses are mutually orthogonal to one
another developing a so-called diamond-plate pattern (see Fig. 1). The authors
have utilized this instability for the one-step assembly of patterned
functional surfaces using superparamagnetic nanoparticles. In [87], similar
buckling patterns are exploited to manipulate the optical properties of pHEMA-
based periodic hydrogel membranes. In this case, the authors observed diamond-
plate patterns during drying of the membrane below a critical temperature
after immersing it in deionized water. Further experimental studies on
periodic hydrogel films are documented in [77, 52], where microstructures with
high-swelling and non-swelling regions developed out-of-plane buckling
patterns. It can be shown that the observed patterns depend significantly on
the morphological parameters of the microstructure.

Similar to the mentioned experimental studies, [62, 61] deal with the 2D
numerical investigation of the effect of prestrain and imperfections on the
instability of periodically perforated hydrogels under plane-strain
conditions. Their simulations are performed on a larger domain of a periodic
structure contained from $2\times 2$ and $10\times 10$ unit-cell
representative volume elements. In line with experiments, they have found that
the diamond-plate patterns are a dominant instability mode, which are obtained
considering random imperfections of the voids without eigenvalue analysis.
Following the theoretical and numerical modeling approaches of these works
supplemented further with eigenvalue buckling analysis in Abaqus, in [59, 48]
3D pattern transformations on a gel film bonded to a soft substrate are
investigated. The authors have studied sequential bifurcations on bifurcated
paths considering that the first bifurcation is associated with a hexagonal
dimple mode, which led to complex patterns such as herringbone and labyrinth
and gave further insight into the nature of the formation of these patterns.
We note that in these contributions, the chemical potential is assumed
spatially constant within the hydrogel domain, so that transient effects (such
as time-dependent diffusion) was not considered. Furthermore, the instability
analyses are usually considered on larger domains, which can be reduced to
unit-cell computations based on Bloch-Floquet representation theorem as
discussed in this work.

All mentioned experimental and numerical studies highlight the wide-ranging
opportunities to tune the response of periodic systems by exploiting
instabilities and hence to obtain materials with enhanced functionalities that
could be considered, for example, in the development of soft and smart
devices.

###  1.4 Novel features and objectives of the present contribution

Motivated by the above mentioned works, we study numerically structural as
well as material instabilities in periodic hydrogel systems. All considered
numerical setups are based on corresponding experimental investigations and
the prescription of experimentally motivated boundary conditions. In
particular, we study instabilities of perforated composite hydrogel
microstructures in two dimensions (cf. [84, 87, 78]) as well as in composite
hydrogel thin films in three spatial dimensions (cf. [77, 52]). Here, we do
not only observe short-wavelength structural instabilities, which alter the
periodicity of the microstructure, but also long-wavelength structural
instabilities considering a transient diffusion process.

To model the transient diffusion process in hydrogels, we employ a rate-type
variational minimization formulation following [55, 7], see also [65] for an
embedding into a multiscale framework. For spatial discretization in two
dimensions, we exploit Raviart–Thomas-type finite elements coupled with an
implicit Euler scheme for time discretization. In three spatial dimensions, we
discretize the domain with (non-conforming) hexahedral finite elements. Within
both space-time discrete formulations, we implement Bloch–Floquet wave
analysis to determine pattern-transforming instabilities together with the
altered periodicity of the microstructure on a unit-cell based computations
without considering any imperfections apriori. Furthermore, we derive a
condition for the loss of strong ellipticity based on effective moduli of
perforated microstructures. The latter is based on an extension of the
surface-averaging approach of [54] to the coupled problem of transient chemo-
mechanics.

Our studies reveal various experimentally confirmed as well as novel buckling
patterns of periodic hydrogel systems depending on their morphology and
associated material properties. For example, it will be shown that micro-
coated hydrogel microstructures develop complex combinations of the well-known
diamond-plate pattern coupled with intricate micro-wrinkling along internal
surfaces. To our best knowledge, such pattern transformations have not been
recorded yet and could serve for additional applications of periodic
hydrogels. Furthermore, our numerical studies in three spatial dimensions
indicate new opportunities for the design of complex instability-induced
pattern transformations driven by the morphology of the underlying
microstructure.

###  1.5 Outline

The present work is structured as follows. As a theoretical basis for the
computational analysis of hydrogels, we discuss a variational formulation and
its numerical implementation in Section 2. Based on that, we address
fundamental aspects of structural and material stability analysis as well as
the associated numerical implementation in Section 3. The developed framework
will then be applied to the computational investigation of instabilities of
two- and three-dimensional hydrogel microstructures as well as corresponding
parametric studies in Section 4. Here, we also document a set of newly
observed pattern transformations. We close with a discussion of the main
contributions and findings in Section 5.

##  2 Variational formulation of deformation-diffusion processes

In this section, we discuss the continuum modeling and computational
implementation of the coupled chemo-mechanical problem. We adopt a variational
minimization formulation as it has proven convenient for the investigation of
instabilities [72]. In the underlying variational principle, the deformation
map and the solvent-volume flux are chosen as the primary unknowns [55, 7].
The main ingredients of this formulation are an energy-storage functional, a
dissipation-potential functional and an external-power functional, where the
latter arises due to applied traction and chemical-potential boundary
conditions. For the latter, we propose additional boundary conditions in terms
of the jumps of primary fields, which proves useful for the investigation of
instabilities in periodic hydrogel structures.

###  2.1 Primary state variables and the second law of thermodynamics

We choose the deformation map ${\boldsymbol{\varphi}}$ and the solvent-volume
flux ${\IKbb H}$ as the primary state variables of the considered minimization
principle. If we denote by ${\mathcal{B}}_{0}$ and ${\mathcal{B}}_{t}$ the
reference and the current configuration of the body, respectively, we can
describe the mapping of these fields as

${\boldsymbol{\varphi}}:\begin{cases}{\mathcal{B}}_{0}\times{\mathcal{T}}\to{\mathcal{}B}_{t}\subset{\mathcal{R}}^{3}\\\\[2.15277pt]
({\boldsymbol{\mathnormal{X}}},t)\mapsto{\boldsymbol{\varphi}}({\boldsymbol{\mathnormal{X}}},t)\end{cases}\quad\mbox{and}\quad{\IKbb
H}:\begin{cases}{\mathcal{B}}_{0}\times{\mathcal{T}}\to{\mathcal{R}}^{3}\\\\[2.15277pt]
({\boldsymbol{\mathnormal{X}}},t)\mapsto{\IKbb
H}({\boldsymbol{\mathnormal{X}}},t)\end{cases}\,,$ (1)

where the deformation map characterizes the deformed position
${\boldsymbol{\mathnormal{x}}}={\boldsymbol{\varphi}}({\boldsymbol{\mathnormal{X}}},t)$
of a material point ${\boldsymbol{\mathnormal{X}}}\in{\mathcal{B}}_{0}$ at a
given time $t\in{\mathcal{T}}$. The solvent-volume flux is a relative
Lagrangian field characterizing local transport of the solvent volume relative
to the motion of the polymer333We refer to [14] for a micromechanically
motivated description of poromechanics. .

An additional state variable associated with the solvent-volume flux is given
by the solvent-volume content
$s:{\mathcal{B}}_{0}\times{\mathcal{T}}\to{\mathcal{R}}^{+}$. For an arbitrary
subdomain ${\mathcal{P}}_{0}\subseteq{\mathcal{B}}_{0}$, the solvent-volume
content is determined from the balance of solvent volume [25]

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{{\mathcal{P}}_{0}}s\,\mathrm{d}V=-\int_{\partial{\mathcal{P}}_{0}}{\IKbb
H}\cdot{\boldsymbol{\mathnormal{N}}}\,\mathrm{d}A\quad\Rightarrow\quad\dot{s}=-\mathop{\operator@font
Div}{\IKbb H}\quad\mbox{in}\ {\mathcal{P}}_{0}.$ (2)

Note that the flux field is defined in the Lagrangian configuration, so that
no convective terms related to the spatial change of the solvent-volume
content arise.

Based on the second law of thermodynamics, we can write the following
inequality for the evolution of the stored-energy density of an arbitrary
subdomain ${\mathcal{P}}_{0}\subseteq{\mathcal{B}}_{0}$

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{{\mathcal{P}}_{0}}\psi\,\mathrm{d}V\leq\underbrace{\int_{\partial{\mathcal{P}}_{0}}{\boldsymbol{\mathnormal{T}}}\cdot\dot{{\boldsymbol{\varphi}}}\,\mathrm{d}A-\int_{\partial{\mathcal{P}}_{0}}\mu
h\,\mathrm{d}A}_{P_{ext}(\dot{{\boldsymbol{\varphi}}},{\IKbb H})}\,,$ (3)

where we have introduced the traction vector ${\boldsymbol{\mathnormal{T}}}$
acting on the surface of the deformed subdomain, defined per unit undeformed
area $\partial{\mathcal{P}}_{0}$, as well as the surface flux $h={\IKbb
H}\cdot{\boldsymbol{\mathnormal{N}}}$; $\mu$ is the chemical potential of the
solvent and characterizes the supplied energy to the system as a result of the
influx of solvent volume. We refer to [14] for an analogous interpretation in
poromechanics. Localizing the above inequality and distinguishing local and
diffusive parts of the dissipation, we obtain

${\mathcal{D}}_{loc}:={\boldsymbol{\mathnormal{P}}}:\dot{{\boldsymbol{\mathnormal{F}}}}+\mu\dot{s}-\dot{\psi}\geq
0\quad\mbox{and}\quad{\mathcal{D}}_{dif}:={\IKbb H}\cdot{\IKbb M}\geq 0\,,$
(4)

where ${\boldsymbol{\mathnormal{F}}}:=\mathop{\operator@font Grad}\varphi$ is
the deformation gradient with $J:=\det{\boldsymbol{\mathnormal{F}}}>0$ and
${\IKbb M}:=-\mathop{\operator@font Grad}\mu$ denotes the negative gradient of
the chemical potential. Since we assume diffusion to take place in an
otherwise elastic solid, the local part of the dissipation must vanish, i.e.,
${\mathcal{D}}_{loc}=0$. Thus, (4)1 can be exploited to obtain the
constitutive relations for the first Piola-Kirchhoff stress tensor
${\boldsymbol{\mathnormal{P}}}=\partial_{\boldsymbol{\mathnormal{F}}}\hat{\psi}({\boldsymbol{\mathnormal{F}}},s)$
and for the chemical potential
$\mu=\partial_{s}\hat{\psi}({\boldsymbol{\mathnormal{F}}},s)$. The inequality
(4)2 characterizes the dissipation due to the diffusion process, where ${\IKbb
M}$ is the driving force for the solvent-volume transport. The latter can be
determined from a dissipation-potential density $\hat{\phi}({\IKbb
H};{\boldsymbol{\mathnormal{F}}},s)$ for a given constitutive state, i.e.,
${\IKbb M}=\partial_{\IKbb H}\hat{\phi}$.

###  2.2 Constitutive functions: The free-energy and the dissipation-potential
density

As we have seen above, the stress tensor, the chemical potential and the
negative gradient of the chemical potential are local constitutive fields.
They are determined from a free-energy density
$\hat{\psi}({\boldsymbol{\mathnormal{F}}},s)$ and a dissipation-potential
density $\hat{\phi}({\IKbb H};{\boldsymbol{\mathnormal{F}}},s)$.

In what follows, we consider an additive free-energy function of the form [37,
34, 7, 11]

$\hat{\psi}({\boldsymbol{\mathnormal{F}}},s)=\hat{\psi}_{mech}({\boldsymbol{\mathnormal{F}}})+\hat{\psi}_{chem}(s)+\hat{\psi}(J,s)\,.$
(5)

While the mechanical contribution of the free-energy function is assumed to be
of Neo-Hookean-type, the chemical contribution is specified as a Flory–Rehner-
type function [23]. These contributions together with a coupling term yield
the overall free-energy function

$\hat{\psi}({\boldsymbol{\mathnormal{F}}},s)={\displaystyle\frac{\gamma}{2}\big{[}{\boldsymbol{\mathnormal{F}}}:{\boldsymbol{\mathnormal{F}}}-3-2\ln(\det{\boldsymbol{\mathnormal{F}}})\big{]}}+\alpha\big{[}{\displaystyle
s\ln\big{(}\frac{s}{1+s}\big{)}+\frac{\chi
s}{1+s}}\big{]}+{\displaystyle\frac{\epsilon}{2}\big{(}\det{\boldsymbol{\mathnormal{F}}}-1-s\big{)}^{2}}\,,$
(6)

where $\gamma$ is the shear modulus and $\alpha:=kT/\nu$ is a mixing modulus
with Boltzmann constant $k$, absolute temperature $T$ and volume of a solvent
molecule $\nu$; $\chi$ is the Flory–Huggins interaction parameter, $\epsilon$
is a penalty parameter to constrain $\det{\boldsymbol{\mathnormal{F}}}=1+s$
for incompressible hydrogels. We refer to Table 1 for further information on
the given parameters.

The dissipation-potential function is chosen such that (4)2 is satisfied a
priori. Therefore, $\hat{\phi}({\IKbb H};{\boldsymbol{\mathnormal{F}}},s)$ is
assumed to be a convex and normalized function of the solvent-volume flux
${\IKbb H}$ according to $\hat{\phi}({\boldsymbol{\mathit{0}}})=0$ and
$\partial_{\IKbb
H}\hat{\phi}({\boldsymbol{\mathit{0}}})={\boldsymbol{\mathit{0}}}$ at the
reference configuration. In the present work, we will consider a dissipation
function related to Fickian diffusion given by

${\displaystyle\hat{\phi}({\IKbb
H};{\boldsymbol{\mathnormal{C}}},s)=\frac{1}{2Ms}{\boldsymbol{\mathnormal{C}}}:({\IKbb
H}\otimes{\IKbb H})}\,,$ (7)

where $M$ is referred to as a mobility parameter. The dissipation-potential
function is formulated in terms of given values
$\\{{\boldsymbol{\mathnormal{F}}},s\\}$ and therefore evaluated numerically at
the previous time step $(\cdot)_{n}$ within an incremental variational
formulation [55, 7], see Section 2.4 for details.

The above introduced constitutive functions are given with respect to a dry
reference configuration. Since (6) is singular for the dry state $s(t=0)=0$, a
pre-swollen state is usually considered as the reference configuration [37,
34, 7]. According to that, the constitutive functions (6) and (7) defined per
unit volume of the pre-swollen reference configuration and appear as

$\hat{\psi}=\frac{\gamma}{2J_{0}}[J_{0}^{2/3}{\boldsymbol{\mathnormal{F}}}:{\boldsymbol{\mathnormal{F}}}-3-2\ln(JJ_{0})]+\frac{\alpha}{J_{0}}[s\ln(\frac{s}{1+s})+\frac{\chi
s}{1+s}]+\frac{\epsilon}{2J_{0}}(JJ_{0}-1-s)^{2}$ (8)

and

$\hat{\phi}=\frac{1}{2J_{0}^{1/3}Ms_{n}}{\boldsymbol{\mathnormal{C}}}_{n}:({\IKbb
H}\otimes{\IKbb H})\,,$ (9)

where the pre-swollen state is characterized by a Jacobian $J_{0}$. The
initial condition for the solvent-volume content associated with the pre-
swollen reference configuration reads

$s(t=0)=s_{0}=\frac{\gamma}{\epsilon}\big{(}J_{0}^{-1/3}-\frac{1}{J_{0}}\big{)}+J_{0}-1\,.$
(10)

The chemical potential at the pre-swollen reference configuration is
determined as [7]

$\mu_{0}=-\frac{\epsilon}{J_{0}}(J_{0}-1-s_{0})+\frac{\alpha}{J_{0}}\big{[}\ln{\big{(}\frac{s_{0}}{1+s_{0}}\big{)}}+\frac{1}{1+s_{0}}+\frac{\chi}{(1+s_{0})^{2}}\big{]}\,.$
(11)

The latter relation is particularly required when applying chemical-potential
boundary conditions in the minimization formulation, see Section 4.

###  2.3 Boundary conditions for periodic hydrogel structures

Figure 2: Boundary conditions for a two-dimensional periodic hydrogel
microstructure. The response of periodic materials can be determined by means
of unit-cell representative volume elements (unit-cell RVEs) of their
microstructure. In the present approach, the unit-cell problem is solved in
consideration of the jump conditions for the deformation map
$\left[\kern-1.33334pt\left[{\boldsymbol{\varphi}}\right]\kern-1.33334pt\right]=\overset{\rule[-0.17223pt]{4.47778pt}{0.25832pt}}{{\boldsymbol{\mathnormal{F}}}}{\boldsymbol{\mathnormal{X}}}$
and the solvent-volume flux $\left[\kern-1.33334pt\left[{\IKbb
H}\right]\kern-1.33334pt\right]\cdot{\boldsymbol{\mathnormal{N}}}=0$ across
the external boundary according to (12). As can be seen in the figure, the
solvent is assumed to enter the RVE through the boundary of the voids. In our
computations, this is realized by applying the chemical potential at the
void’s boundary according to (47).

Since our main aim is to investigate instabilities in periodic hydrogels, we
will now provide a suitable setup for associated boundary conditions. As we
will see, these are partly inspired from homogenization theory such that they
allow for the definition of effective quantities, which could be used in the
material stability analysis to be discussed later.

For illustration, let us consider the two-dimensional representation of a
periodic hydrogel as depicted in Fig. 2. Here we assume that any out-of-plane
deformation could be neglected and therefore typical plane-strain conditions
apply. Since the structure is periodic, any stability analysis could be
reduced to the smallest periodic domain, here given in terms of a square-
shaped matrix with a circular void such that
${\mathcal{B}}_{0}={\mathcal{D}}_{0}\cup{\mathcal{H}}_{0}$ with
${\mathcal{H}}_{0}$ denoting the domain of the void. In the following, we will
refer to such a smallest possible periodic domain as a unit-cell
representative volume element (unit-cell RVE).

For the chemo-mechanical loading of the unit-cell RVE we prescribe boundary
conditions of the generalized form

$\left[\kern-1.66669pt\left[{\boldsymbol{\varphi}}\right]\kern-1.66669pt\right]\,=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}{\boldsymbol{\mathnormal{X}}}\quad\text{and}\quad\left[\kern-1.66669pt\left[{\IKbb
H}\right]\kern-1.66669pt\right]\cdot{\boldsymbol{\mathnormal{N}}}=0\quad\text{on}\quad\partial{\mathcal{B}}_{0}\,,$
(12)

where we have introduced an _effective_ or _macroscopic_ deformation gradient
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}$
[56, 65]. The dual variable to the effective deformation gradient is the
effective first Piola-Kirchhoff stress tensor
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}$.
Both quantities can be determined as the volume averages of their microscopic
counterparts as444We note that when we prescribe the effective deformation
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}$
via the RVE boundary, the effective stresses can be determined via (13)2. Vice
versa, it is also possible to apply effective stresses
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}$
on the RVE. In that case, the effective deformation gradient can be determined
via (13)1, see [54, 57].

$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}=\frac{1}{|{\mathcal{B}}_{0}|}\int_{\partial{\mathcal{B}}_{0}}{\boldsymbol{\varphi}}\otimes{\boldsymbol{\mathnormal{N}}}\,\mathrm{d}A\quad\mbox{and}\quad\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}=\frac{1}{|{\mathcal{B}}_{0}|}\int_{\partial{\mathcal{B}}_{0}}({\boldsymbol{\mathnormal{P}}}{\boldsymbol{\mathnormal{N}}})\otimes{\boldsymbol{\mathnormal{X}}}\,\mathrm{d}A\,,$
(13)

where $|{\mathcal{B}}_{0}|$ denotes the volume of a unit-cell RVE. Note that
according to (12)2, we have defined vanishing jump conditions for the solvent-
volume flux such that no effective flux is applied on the RVE. This setting is
motivated by the assumption that the diffusion process is driven exclusively
by suitable chemical-potential boundary conditions acting at the internal
boundary of unit-cell RVEs. Please refer to Section 4 for more details.

###  2.4 Rate-type variational formulation of the deformation-diffusion
problem

Above, we have introduced the energy-storage function
$\hat{\psi}({\boldsymbol{\mathnormal{F}}},s)$, the dissipation-potential
function $\hat{\phi}({\IKbb H})$ and the external-loading functional
$P_{ext}(\dot{{\boldsymbol{\varphi}}},{\IKbb H})$. Having these functions at
hand, we define a rate-type potential to characterize the deformation-
diffusion response [55]

$\displaystyle\Pi(\dot{{\boldsymbol{\varphi}}},{\IKbb H})$
$\displaystyle:=\int_{{\mathcal{B}}_{0}}\pi(\dot{{\boldsymbol{\varphi}}},{\IKbb
H})\,\mathrm{d}V-P_{ext}(\dot{{\boldsymbol{\varphi}}},{\IKbb
H})\quad\mbox{with}\quad\pi(\dot{{\boldsymbol{\varphi}}},{\IKbb
H}):=\frac{\mbox{d}}{\mbox{d}t}\hat{\psi}({\boldsymbol{\mathnormal{F}}},s)+\hat{\phi}({\IKbb
H};{\boldsymbol{\mathnormal{F}}},s),$ (14)

where $\pi(\dot{{\boldsymbol{\varphi}}},{\IKbb H})$ is the internal-power
density. As a result, the variational principle for the considered problem in
the continuous setting can be introduced as

$\\{\dot{{\boldsymbol{\varphi}}}^{\star},{\IKbb
H}^{\star}\\}=\mathop{\operator@font
Arg}\big{\\{}\inf_{\dot{{\boldsymbol{\varphi}}}\in{\mathcal{W}}_{\dot{{\boldsymbol{\varphi}}}}}\inf_{{\IKbb
H}\in{\mathcal{W}}_{{\IKbb H}}}\Pi(\dot{{\boldsymbol{\varphi}}},{\IKbb
H})\big{\\}}\,,$ (15)

where ${\mathcal{W}}_{\dot{{\boldsymbol{\varphi}}}}$ and
${\mathcal{W}}_{{\IKbb H}}$ are the admissible spaces for the rate of the
deformation map and the solvent-volume flux

$\displaystyle{\mathcal{W}}_{\dot{{\boldsymbol{\varphi}}}}$
$\displaystyle:=\\{\dot{{\boldsymbol{\varphi}}}\in
H^{1}({\mathcal{B}}_{0})\,\big{|}\,\dot{{\boldsymbol{\varphi}}}=\dot{{\boldsymbol{\varphi}}}_{D}\,\,\text{on}\,\,\partial{\mathcal{B}}_{0}^{\boldsymbol{\varphi}}\\}\quad\mbox{and}\quad$
(16) $\displaystyle{\mathcal{W}}_{{\IKbb H}}$ $\displaystyle:=\\{{\IKbb H}\in
H(\mathop{\operator@font Div},{\mathcal{B}}_{0})\,\big{|}\,{\IKbb H}={\IKbb
H}_{D}\,\,\text{on}\,\,\partial{\mathcal{B}}_{0}^{h}\\}\,.$

The Euler-Lagrange equations follow as the necessary conditions of the
variational formulation as

$\boxed{\begin{array}[]{lrclcl}\mbox{Balance of linear
momentum:}&\mathop{\operator@font
Div}[\partial_{{\boldsymbol{\mathnormal{F}}}}\hat{\psi}]&=&{\boldsymbol{\mathit{0}}}&\text{in}&{\mathcal{B}}_{0}\\\\[4.2679pt]
\mbox{Inverse Fickian law:}&\nabla[\partial_{s}\hat{\psi}]+\partial_{\IKbb
H}\hat{\phi}&=&{\boldsymbol{\mathit{0}}}&\text{in}&{\mathcal{B}}_{0}\\\\[4.2679pt]
\mbox{Mechanical
traction:}&\partial_{{\boldsymbol{\mathnormal{F}}}}\hat{\psi}\cdot{\boldsymbol{\mathnormal{N}}}-{{\boldsymbol{\mathnormal{T}}}}&=&{\boldsymbol{\mathit{0}}}&\text{on}&\partial{\mathcal{B}}_{0}^{T}\\\\[4.2679pt]
\mbox{Chemical
potential:}&-\partial_{s}\hat{\psi}+{\mu}&=&0&\text{on}&\partial{\mathcal{B}}_{0}^{\mu}\end{array}}$
(17)

###  2.5 Space-time discretization of the variational formulation

We solve (15) using a conforming finite-element formulation. Before
considering the discretization in space, we apply a modified version of the
implicit Euler scheme in a time interval $[t_{n},t]$ with $t$ being the
current time step with $0<t\leq T$ in a considered overall simulation time
$T$. Note that for the variables evaluated at the current time step, we do not
use subscripts. Consequently, the time-discrete version of (15) reads555We
remark that the variational formulation (18) should also contain the energy-
storage function as well as the external-power functional evaluated at the
previous time step. However, since these terms are constant they vanish under
variation and are therefore omitted from the equations to arrive at a compact
representation. [55]

$\Pi^{\tau}({\boldsymbol{\varphi}},{\IKbb
H}):=\int_{{\mathcal{B}}_{0}}\hat{\psi}({\boldsymbol{\mathnormal{F}}},s)+\tau\hat{\phi}({\IKbb
H};{\boldsymbol{\mathnormal{F}}}_{n},s_{n})\>\mathrm{d}V-P{}^{\tau}_{ext}({{\boldsymbol{\varphi}}},{\IKbb
H})\,,$ (18)

where the solvent-volume content at the current time step is determined by
$s=s_{n}-\tau\mathop{\operator@font Div}[{\IKbb H}]$. The dissipation-
potential function $\hat{\phi}({\IKbb
H};{\boldsymbol{\mathnormal{F}}}_{n},s_{n})$ is evaluated using the
deformation gradient and the solvent-volume content from the previous time
step $t_{n}$ in order to stay consistent with the continuous variational
setting [55, 7]. For the explicit form of $\hat{\phi}$ we refer to (9). The
time-discrete external-power functional reads

$P{}^{\tau}_{ext}({{\boldsymbol{\varphi}}},{\IKbb
H})=\int_{\partial{\mathcal{B}}_{0}}{\boldsymbol{\mathnormal{T}}}\cdot{\boldsymbol{\varphi}}\>\mathrm{d}A-\int_{\partial{\mathcal{B}}_{0}}\tau\mu
h\>\mathrm{d}A.$ (19)

As a result, we obtain the time-discrete variational formulation with an
initial condition for the solvent-volume content as

$\\{{{\boldsymbol{\varphi}}^{\star}},{\IKbb
H}^{\star}\\}=\mathop{\operator@font
Arg}\\{\inf_{{{\boldsymbol{\varphi}}}\in{\mathcal{W}}_{{{\boldsymbol{\varphi}}}}}\inf_{{\IKbb
H}\in{\mathcal{W}}_{{\IKbb H}}}\Pi^{\tau}({{\boldsymbol{\varphi}}},{\IKbb
H})\\}\quad\text{with}\quad s(t_{0})=s_{0}.$ (20)

The above formulation is implemented into a conforming finite-element
formulation by using a combination of Lagrange- and Raviart–Thomas-type shape
functions [7, 65]. This then results in a minimization principle in the space-
time discrete form as

${\underline{\boldsymbol{\mathsf{d}}}^{\star}}=\mathop{\operator@font
Arg}\\{\min_{\underline{\boldsymbol{\mathsf{d}}}}\Pi^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}})\\},$ (21)

where $\underline{\boldsymbol{\mathsf{d}}}$ contains all degrees of freedom
comprising the node-based displacements and the edge-based normal-projected
solvent-volume fluxes.

##  3 Stability analysis of periodic hydrogel systems

In this section, we describe methods to investigate the structural and the
material stability of periodic hydrogel systems. For the structural stability
analysis, we utilize the minimization structure of the underlying incremental
system of equations [30, 31, 56]. In doing so, we investigate in particular
_short-wavelength_ and _long-wavelength_ instabilities. The former are usually
associated with a change in periodicity of an underlying microstructure, which
could arise, for example, due to wrinkling of internal surfaces in perforated
systems. We will detect these types of instabilities via Bloch–Floquet
analysis. Long-wavelength instabilities are related to the loss of strong
ellipticity of a corresponding homogenized continuum. These kinds of
instabilities can be detected by checking the positive definiteness of an
underlying acoustic tensor. In what follows, we will also present a
corresponding numerical method for the analysis of related instabilities.

###  3.1 Minimization-based structural stability analysis

Within a minimization-based approach, hydrogel systems are considered to be
structurally stable at an equilibrium state
$\underline{\boldsymbol{\mathsf{d}}}^{\star}$, if all variations of the
incremental potential from this equilibrium state under infinitesimally small
perturbations are positive, that is [56]

$\Delta\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star},\delta\underline{\boldsymbol{\mathsf{d}}}):=\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star}+\epsilon\delta\underline{\boldsymbol{\mathsf{d}}})-\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star})>0\,,$ (22)

where $\epsilon$ is a small parameter. Considering a Taylor approximation of
(22) and truncating higher-order terms, we obtain

$\Delta\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star},\delta\underline{\boldsymbol{\mathsf{d}}})\approx\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{\epsilon=0}\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star}+\epsilon\delta\underline{\boldsymbol{\mathsf{d}}})+\frac{1}{2}\frac{\mathrm{d}^{2}}{\mathrm{d}\epsilon^{2}}\bigg{|}_{\epsilon=0}\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star}+\epsilon\delta\underline{\boldsymbol{\mathsf{d}}})>0\,.$
(23)

Since the first variation vanishes at equilibrium, we arrive at

$\Delta\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star},\delta\underline{\boldsymbol{\mathsf{d}}})\approx\frac{1}{2}\frac{\mathrm{d}^{2}}{\mathrm{d}\epsilon^{2}}\bigg{|}_{\epsilon=0}\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star}+\epsilon\delta\underline{\boldsymbol{\mathsf{d}}})>0\,.$
(24)

Consequently, structural stability is associated with the second variation of
the incremental potential in the sense that

$\Delta\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}}^{\star},\delta\underline{\boldsymbol{\mathsf{d}}})=\frac{1}{2}\delta\underline{\boldsymbol{\mathsf{d}}}{}^{T}\big{[}\Pi{}^{\tau}(\underline{\boldsymbol{\mathsf{d}}}){}_{,\underline{\boldsymbol{\mathsf{d}}}\underline{\boldsymbol{\mathsf{d}}}}\big{]}_{\underline{\boldsymbol{\mathsf{d}}}^{\star}}\delta\underline{\boldsymbol{\mathsf{d}}}=\Lambda\delta\underline{\boldsymbol{\mathsf{d}}}^{T}\delta\underline{\boldsymbol{\mathsf{d}}}>0.$
(25)

The latter results in an eigenvalue problem in terms of the tangent of the
underlying system of equations according to

${\Lambda}=\min_{{\boldsymbol{\mathnormal{n}}}\in{\Ibb
N}^{3}}\min_{\underline{\boldsymbol{\mathsf{d}}}}\big{\\{}\Delta\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}},\delta\underline{\boldsymbol{\mathsf{d}}})\big{/}\|\delta\underline{\boldsymbol{\mathsf{d}}}\|^{2}\big{\\}}\begin{cases}{>0}\quad\text{for
stable state}\ \underline{\boldsymbol{\mathsf{d}}}^{\star}\\\ {\leq
0}\quad\text{for unstable state}\
\underline{\boldsymbol{\mathsf{d}}}^{\star}\end{cases}\,,$ (26)

where ${\boldsymbol{\mathnormal{n}}}\in{\Ibb N}^{3}$ denotes the number of
unit-cells RVEs with respect to three spatial directions [60]. The
minimization over ${\boldsymbol{\mathnormal{n}}}$ allows for the detection of
instabilities that would not be describable by means of a unit-cell RVE.
However, the above approach is computationally inefficient since it requires
the consideration of various realizations of RVEs that grow in size with
growing ${\boldsymbol{\mathnormal{n}}}$. Therefore, we use Bloch–Floquet
analysis to detect such critical structural instabilities on a unit-cell RVE
[24].

Based on Bloch–Floquet analysis, the variations
$\delta\underline{\boldsymbol{\mathsf{d}}}$ of the degrees of freedom can be
represented in terms of unit-cell variations
$\delta\underline{\boldsymbol{\mathsf{d}}}_{{\mathcal{B}}_{0}}$ as

$\delta\underline{\boldsymbol{\mathsf{d}}}({\boldsymbol{\mathnormal{X}}})=\delta\underline{\boldsymbol{\mathsf{d}}}_{{\mathcal{B}}_{0}}({\boldsymbol{\mathnormal{X}}})\exp[\mathrm{i}{\boldsymbol{\mathnormal{k}}}\cdot{\boldsymbol{\mathnormal{X}}}]\quad\text{with}\quad\delta\underline{\boldsymbol{\mathsf{d}}}^{+}_{{\mathcal{B}}_{0}}=\delta\underline{\boldsymbol{\mathsf{d}}}^{-}_{{\mathcal{B}}_{0}}\quad\text{on}\quad\partial{\mathcal{B}}_{0}=\partial{\mathcal{B}}_{0}^{+}\cup\partial{\mathcal{B}}_{0}^{-},$
(27)

where ${\boldsymbol{\mathnormal{k}}}$ is referred to as Bloch vector, which
characterizes the wavelength of the variations. Consequently, this expression
leads to the boundary conditions

$\delta\underline{\boldsymbol{\mathsf{d}}}({\boldsymbol{\mathnormal{X}}}^{+})=\delta\underline{\boldsymbol{\mathsf{d}}}({\boldsymbol{\mathnormal{X}}}^{-})\exp[\mathrm{i}{\boldsymbol{\mathnormal{k}}}\cdot({\boldsymbol{\mathnormal{X}}}^{+}-{\boldsymbol{\mathnormal{X}}}^{-})]$
(28)

for a unit-cell RVE. Here, ${\boldsymbol{\mathnormal{X}}}^{+}$ and
${\boldsymbol{\mathnormal{X}}}^{-}$ denote the reference placement of periodic
nodes, which are located at opposite boundaries of the unit-cell RVE, see Fig.
2. It follows a stability criterion related to the tangent of a unit-cell RVE
according to666For details related to the numerical implementation we refer to
[75, 3, 67, 66].

$\boxed{{\Lambda}=\min_{k_{i}\in[0,\pi]}\min_{\underline{\boldsymbol{\mathsf{d}}}}\big{\\{}\Delta\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}},\delta\underline{\boldsymbol{\mathsf{d}}})\big{/}\|\delta\underline{\boldsymbol{\mathsf{d}}}\|^{2}\big{\\}}\begin{cases}{>0}\quad\text{for
stable state}\ \underline{\boldsymbol{\mathsf{d}}}^{\star}\\\ {\leq
0}\quad\text{for unstable state}\
\underline{\boldsymbol{\mathsf{d}}}^{\star}\end{cases}}\,.$ (29)

Depending on the values of the Bloch vector components
$k_{i}\in{\mathcal{R}}^{3}$, we distinguish three types of instabilities: (i)
unit-cell periodic instabilities for $k_{i}=0$; (ii) short-wavelength
instabilities for finite values of $k_{i}$ usually spanning RVEs made of
multiple unit cells, where the explicit number of cells can be computed as
$n_{i}=2\pi/k_{i}$; (iii) long-wavelength instabilities for
${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$. We note that the
latter could also be determined from the effective mechanical moduli of
periodic hydrogels, see below.

###  3.2 Loss of strong ellipticity: Material stability analysis at macroscale

Having defined macroscopic mechanical boundary conditions in terms of
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}$,
we assume that at macroscale we need to fulfill the balance equation [65]

$\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Div}}[\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}{}(\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}},\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Div}}[\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\IKbb
H}}])]=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\rho}_{0}\ddot{\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\varphi}}}}\,.$
(30)

In order to check if the latter equation satisfies the strong-ellipticity
condition, we superimpose a plane wave of the form
$\delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\varphi}}}=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{n}}}}f(\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{X}}}}\cdot\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{N}}}}-\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{c}t)$,
such that for infinitesimal
$\delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\varphi}}}$
it holds that [31, 56]

$\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Div}}[\delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}{}]=\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Div}}[\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\Ibb
A}}{}:\delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}]=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\rho}_{0}\delta\ddot{\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\varphi}}}}\quad\text{with}\quad\delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}=\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Grad}}\delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\varphi}}}\,,$
(31)

where $\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\Ibb A}}$ are the
effective mechanical moduli of the microstructure. Note that since we do not
consider any macroscopic loading in terms of
$\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Div}}[\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\IKbb H}}]$, we do
not have to take into account a corresponding variation in the above equation.
After some modifications, (31) appears as

$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{n}}}}\cdot\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{Q}}}}\cdot\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{n}}}}=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\rho}_{0}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{c}{}^{2}\,,$
(32)

where
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{Q}_{ab}=(\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\Ibb
A}}{})\indices{{}_{a}^{A}{}_{b}^{B}}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{N}_{A}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{N}_{B}$
is the mechanical acoustic tensor defined in terms of the effective mechanical
moduli. It follows the criterion for the material stability at the macroscopic
scale

$\boxed{\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\Lambda}=\min_{\|\overset{\rule[-0.1507pt]{3.91806pt}{0.22603pt}}{{\boldsymbol{\mathnormal{N}}}}\|=1}\big{\\{}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{n}}}}\cdot\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{Q}}}}\cdot\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{n}}}}\big{\\}}\begin{cases}>0\quad\text{for
materially stable state}\\\ \leq 0\quad\text{for materially unstable state}\\\
\end{cases}}\,.$ (33)

###  3.3 Computation of the effective mechanical moduli of voided hydrogel
structures

We now describe the computation of the effective mechanical moduli of voided
hydrogel structures following the approach proposed by Miehe and Koch [54].
This approach allows to determine effective moduli based on surface integrals
along the external boundary of RVEs. Starting point is the definition of the
effective mechanical stresses (13)2 and the effective chemical potential [65]

$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}=\frac{1}{|{\mathcal{B}}_{0}|}\int_{\partial{\mathcal{B}}_{0}}({\boldsymbol{\mathnormal{P}}}{\boldsymbol{\mathnormal{N}}})\otimes{\boldsymbol{\mathnormal{X}}}\,\mathrm{d}A\quad\mbox{and}\quad\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\mu}=\frac{1}{3|{\mathcal{B}}_{0}|}\int_{\partial{\mathcal{B}}_{0}}(\mu{\boldsymbol{\mathnormal{N}}})\cdot{\boldsymbol{\mathnormal{X}}}\,\mathrm{d}A\,,$
(34)

where $\partial{\mathcal{B}}_{0}$ denotes the external boundary of an RVE. In
the discrete setting, the above equations can be given as follows (see [54]
for the purely mechanical case)

$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}=\frac{1}{|{\mathcal{B}}_{0}|}\sum_{I=1}^{P}{\boldsymbol{\lambda}}^{f}_{I}\otimes\left[\kern-1.66669pt\left[{\boldsymbol{\mathnormal{X}}}^{b}_{I}\right]\kern-1.66669pt\right]\quad\mbox{and}\quad\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\mu}=\frac{1}{3|{\mathcal{B}}_{0}|}\sum_{I=1}^{P}\lambda^{\mu}_{I}{\boldsymbol{\mathnormal{N}}}\cdot\left[\kern-1.66669pt\left[{\boldsymbol{\mathnormal{X}}}^{b}_{I}\right]\kern-1.66669pt\right]\,,$
(35)

where ${\boldsymbol{\lambda}}^{f}_{I}$, $\lambda^{\mu}_{I}$ and
${\boldsymbol{\mathnormal{X}}}^{b}_{I}$ denote the nodal forces, the chemical
potential and the coordinates associated with the node $I$ on the external
boundary of an RVE, respectively. The number $P=M/2$ is defined from the total
number of boundary nodes $M$. Introducing an array of nodal forces and
chemical potentials
${\boldsymbol{\lambda}}:=[{\boldsymbol{\lambda}}^{f}_{1},\lambda^{\mu}_{1},...,{\boldsymbol{\lambda}}^{f}_{P},\lambda^{\mu}_{P}]^{T}$
as well as an array containing the effective stresses (in vector notation) and
the effective chemical potential
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{P}}}:=[\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}},\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\mu}]^{T}$,
we can write (35) in a compact
form777$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{P}}}:=[\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11},\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{22},\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{12},\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{21},\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\mu}]^{T}$
is considered for the corresponding array in two dimensions.

$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{P}}}=\frac{1}{|{\mathcal{B}}_{0}|}\sum_{I=1}^{P}{\Ibb
Q}_{I}{\boldsymbol{\lambda}}_{I}=\frac{1}{|{\mathcal{B}}_{0}|}{\Ibb
Q}{\boldsymbol{\lambda}}\,,$ (36)

where ${\Ibb Q}:=[{\Ibb Q}_{1},...,{\Ibb Q}_{P}]$ is the array of projection
tensors with

${\Ibb
Q}_{I}^{T}:=\begin{bmatrix}\left[\kern-1.66669pt\left[({\boldsymbol{\mathnormal{X}}}^{b}_{I})_{1}\right]\kern-1.66669pt\right]&0&\left[\kern-1.66669pt\left[({\boldsymbol{\mathnormal{X}}}^{b}_{I})_{2}\right]\kern-1.66669pt\right]&0&0\\\\[5.69054pt]
0&\left[\kern-1.66669pt\left[({\boldsymbol{\mathnormal{X}}}^{b}_{I})_{2}\right]\kern-1.66669pt\right]&0&\left[\kern-1.66669pt\left[({\boldsymbol{\mathnormal{X}}}^{b}_{I})_{1}\right]\kern-1.66669pt\right]&0\\\\[5.69054pt]
0&0&0&0&\displaystyle\frac{1}{3}\left[\kern-1.66669pt\left[{\boldsymbol{\mathnormal{X}}}_{I}^{b}\right]\kern-1.66669pt\right]\cdot{\boldsymbol{\mathnormal{N}}}\end{bmatrix}$
(37)

in two dimensions.

In the following, we are only interested in the effective mechanical moduli.
Since these relate the incremental effective stresses to the incremental
effective deformation gradient according to
$\Delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\Ibb
A}}:\Delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}$,
we first determine the increment of the first Piola-Kirchhoff stress tensor in
vector notation from (36) as

$\Delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}={\Ibb
L}\Delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{P}}}=\frac{1}{|{\mathcal{B}}_{0}|}{\Ibb
L}{\Ibb Q}\Delta{\boldsymbol{\lambda}}\,,$ (38)

where ${\Ibb L}$ is a projection tensor with entries $0$ or $1$.

In order to obtain the final expression for the effective-moduli tensor, we
need to have a relation between the increment of the generalized nodal forces
$\Delta{\boldsymbol{\lambda}}$ and the effective deformation gradient
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}$.
We obtain such a relationship from the following variational principle in the
space-time discrete setting

$\widetilde{\Pi}{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}},{\boldsymbol{\lambda}})=\Pi{}^{\tau
h}(\underline{\boldsymbol{\mathsf{d}}})+{\boldsymbol{\lambda}}\cdot(\left[\kern-1.66669pt\left[\underline{\boldsymbol{\mathsf{d}}}^{b}\right]\kern-1.66669pt\right]-{\Ibb
Q}^{T}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{F}}})\,,$
(39)

where
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{F}}}:=[\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}},\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Div}}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\IKbb H}}]^{T}$ is an
array given by
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{F}}}=[\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{F}_{11},\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{F}_{22},\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{F}_{12},\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{F}_{21},\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Div}}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\IKbb H}}]^{T}$ in two
dimensions.

The second term on the right-hand side of (39) considers the Lagrange
multiplier method to enforce the jump boundary conditions for the deformation
map and the fluid-volume flux vector of the form

$\left[\kern-1.66669pt\left[{\boldsymbol{\varphi}}\right]\kern-1.66669pt\right]=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}\left[\kern-1.66669pt\left[{\boldsymbol{\mathnormal{X}}}\right]\kern-1.66669pt\right]\quad\mbox{and}\quad\left[\kern-1.66669pt\left[{\IKbb
H}\right]\kern-1.66669pt\right]\cdot{\boldsymbol{\mathnormal{N}}}=\frac{1}{3}\overset{\rule[-0.21529pt]{10.76385pt}{0.32289pt}}{\mathop{\operator@font
Div}}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\IKbb
H}}\left[\kern-1.66669pt\left[{\boldsymbol{\mathnormal{X}}}\right]\kern-1.66669pt\right]\cdot{\boldsymbol{\mathnormal{N}}}\,,$
(40)

refer to (12) as well as to [65] for more details. Furthermore,
${\boldsymbol{\mathnormal{X}}}^{b}:=[{\boldsymbol{\mathnormal{X}}}^{b}_{1},...,{\boldsymbol{\mathnormal{X}}}^{b}_{P}]^{T}$
in (39) contains the coordinates of the nodes at the external boundary of an
RVE. The Lagrange multiplier ${\boldsymbol{\lambda}}$ can be identified as the
array of nodal forces and chemical potentials. Its interpretation follows from
the Euler-Lagrange equations of the variational principle. Furthermore,
$\underline{\boldsymbol{\mathsf{d}}}^{b}:=[\underline{\boldsymbol{\mathsf{d}}}^{b}_{1},..,\underline{\boldsymbol{\mathsf{d}}}^{b}_{P}]^{T}$
denotes the degrees of freedom at the boundary of the unit-cell RVE. In (39),
we can express the jump of the degrees of freedom in terms of the latter as

$\left[\kern-1.66669pt\left[\underline{\boldsymbol{\mathsf{d}}}^{b}\right]\kern-1.66669pt\right]:=\mathop{\textbf{\Large{A}}}_{I=1}^{P}\begin{bmatrix}(\underline{\boldsymbol{\mathsf{d}}}^{b}_{\boldsymbol{\varphi}})_{I}^{+}-(\underline{\boldsymbol{\mathsf{d}}}^{b}_{\boldsymbol{\varphi}})_{I}^{-}\\\\[8.53581pt]
(\mathsf{d}^{b}_{h})_{I}^{+}+(\mathsf{d}^{b}_{h})_{I}^{-}\end{bmatrix}={\Ibb
P}\underline{\boldsymbol{\mathsf{d}}}^{b}\,,$ (41)

where ${\Ibb P}$ is a projection tensor with entries $\\{-1,0,1\\}$. Note
that, since we are considering Raviart–Thomas-type elements, the edge-based
solvent-volume flux degrees of freedom are anti-periodic at the external
boundary [65].

The microscopic problem is usually solved under constant macroscopic
deformation gradient for each time step. Therefore, the first variation of
(39) gives

$\displaystyle\delta\widetilde{\Pi}{}^{\tau
h}=\begin{bmatrix}\displaystyle\delta\underline{\boldsymbol{\mathsf{d}}}^{i}\\\
\delta\underline{\boldsymbol{\mathsf{d}}}^{b}\\\
\delta{\boldsymbol{\lambda}}\end{bmatrix}\cdot\begin{bmatrix}\displaystyle\partial_{\underline{\boldsymbol{\mathsf{d}}}^{i}}\Pi^{\tau
h}\\\ \partial_{\underline{\boldsymbol{\mathsf{d}}}^{b}}\Pi^{\tau h}+{\Ibb
P}{\boldsymbol{\lambda}}\\\ {\Ibb
P}\underline{\boldsymbol{\mathsf{d}}}^{b}-{\Ibb
Q}^{T}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{F}}}\end{bmatrix}=\begin{bmatrix}\displaystyle\delta\underline{\boldsymbol{\mathsf{d}}}^{i}\\\
\delta\underline{\boldsymbol{\mathsf{d}}}^{b}\\\
\delta{\boldsymbol{\lambda}}\end{bmatrix}\cdot\begin{bmatrix}\displaystyle\underline{\boldsymbol{\mathsf{R}}}^{i}\\\
\underline{\boldsymbol{\mathsf{R}}}^{b}+{\Ibb P}{\boldsymbol{\lambda}}\\\
{\Ibb P}\underline{\boldsymbol{\mathsf{d}}}^{b}-{\Ibb
Q}^{T}\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{F}}}\end{bmatrix}={\boldsymbol{\mathit{0}}}\,,$
(42)

where $\underline{\boldsymbol{\mathsf{d}}}^{i}$ denote the interior degrees of
freedom of a unit-cell RVE. Once the above system of equations is solved at
the microscale, the sensitivity of the system with respect to the macroscopic
deformation gradient can be determined from

$\displaystyle\Delta\delta\widetilde{\Pi}{}^{\tau
h}=\begin{bmatrix}\displaystyle\delta\underline{\boldsymbol{\mathsf{d}}}^{i}\\\
\delta\underline{\boldsymbol{\mathsf{d}}}^{b}\\\
\delta{\boldsymbol{\lambda}}\end{bmatrix}\cdot\begin{bmatrix}\displaystyle\underline{\boldsymbol{\mathsf{K}}}^{ii}&\underline{\boldsymbol{\mathsf{K}}}^{ib}&{\boldsymbol{\mathit{0}}}&{\boldsymbol{\mathit{0}}}\\\
\underline{\boldsymbol{\mathsf{K}}}^{bi}&\underline{\boldsymbol{\mathsf{K}}}^{bb}&{\Ibb
P}^{T}&{\boldsymbol{\mathit{0}}}\\\ {\boldsymbol{\mathit{0}}}&{\Ibb
P}&{\boldsymbol{\mathit{0}}}&-{\Ibb
Q}^{T}\end{bmatrix}\cdot\begin{bmatrix}\displaystyle\Delta\underline{\boldsymbol{\mathsf{d}}}^{i}\\\
\Delta\underline{\boldsymbol{\mathsf{d}}}^{b}\\\
\Delta{\boldsymbol{\lambda}}\\\
\Delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{F}}}\end{bmatrix}={\boldsymbol{\mathit{0}}}\,.$
(43)

First, from the variations with respect to
$\delta\underline{\boldsymbol{\mathsf{d}}}^{i}$ and
$\delta\underline{\boldsymbol{\mathsf{d}}}^{b}$ in (43)1 and (43)2, we obtain
elimination equations for $\Delta\underline{\boldsymbol{\mathsf{d}}}^{i}$ and
$\Delta\underline{\boldsymbol{\mathsf{d}}}^{b}$ as

$\Delta\underline{\boldsymbol{\mathsf{d}}}^{i}=-(\underline{\boldsymbol{\mathsf{K}}}^{ii})^{-1}\underline{\boldsymbol{\mathsf{K}}}^{ib}\Delta\underline{\boldsymbol{\mathsf{d}}}^{b}\quad\mbox{and}\quad\Delta\underline{\boldsymbol{\mathsf{d}}}^{b}=-(\underline{\boldsymbol{\mathsf{\widetilde{K}}}}{}^{bb})^{-1}{\Ibb
P}^{T}\Delta{\boldsymbol{\lambda}}\,,$ (44)

where
$\underline{\boldsymbol{\mathsf{\widetilde{K}}}}{}^{bb}=\underline{\boldsymbol{\mathsf{K}}}{}^{bb}-\underline{\boldsymbol{\mathsf{K}}}{}^{bi}(\underline{\boldsymbol{\mathsf{K}}}{}^{ii})^{-1}\underline{\boldsymbol{\mathsf{K}}}{}^{ib}$.
Next, we consider the variation with respect to $\delta{\boldsymbol{\lambda}}$
in (43) and obtain the searched expression of $\Delta{\boldsymbol{\lambda}}$
given by

$\Delta{\boldsymbol{\lambda}}=\big{[}{\Ibb
P}(\underline{\boldsymbol{\mathsf{\widetilde{K}}}}{}^{bb})^{-1}{\Ibb
P}^{T}\big{]}{}^{-1}{\Ibb
Q}^{T}\Delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{\boldsymbol{\mathsf{F}}}\,.$
(45)

Finally, we can compute the effective mechanical moduli as

$\boxed{\Delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{P}}}}=\frac{1}{|{\mathcal{B}}_{0}|}{\Ibb
L}{\Ibb Q}\big{[}{\Ibb
P}(\underline{\boldsymbol{\mathsf{\widetilde{K}}}}{}^{bb})^{-1}{\Ibb
P}^{T}\big{]}{}^{-1}{\Ibb Q}^{T}{\Ibb
L}^{T}\Delta\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}\quad\Rightarrow\quad\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\Ibb
A}}:=\frac{1}{|{\mathcal{B}}_{0}|}{\Ibb L}{\Ibb Q}\big{[}{\Ibb
P}(\underline{\boldsymbol{\mathsf{\widetilde{K}}}}{}^{bb})^{-1}{\Ibb
P}^{T}\big{]}{}^{-1}{\Ibb Q}^{T}{\Ibb L}^{T}}\,.$ (46)

##  4 Numerical investigation of instabilities in hydrogel structures

We now investigate computationally microscopic and macroscopic instabilities
in two- and three-dimensional periodic hydrogels, where all considered
boundary value problems are motivated from experimental studies (e.g., [84,
87, 78]). The two-dimensional studies of single- and two-phase voided hydrogel
microstructures under plane-strain conditions are summarized in Section 4.1.
Corresponding investigations of structrual instabilities of three-dimensional
hydrogel thin films are documented in Section 4.2. In all cases, we examine
the influence of physical properties and microscopic morphology on the onset
and the type of instabilities.

###  4.1 Two-dimensional hydrogel structures

In this section, we study instabilities in single-phase and two-phase voided
hydrogel microstructures. While the single-phase hydrogels consist of a
hydrogel matrix and voids, the two-phase hydrogels further contain stiff,
permeable coatings surrounding the voids as described in Fig. 3. The latter
structures are motivated from experimental studies on constrained plane
bilayer films, which reveal a rich set of wrinkling patterns [29, 27, 28]; see
also the analytical and numerical studies [47, 46, 9, 79, 16, 72]. Associated
wrinkling phenomena have not only been reported for plane bilayers, but also
for bilayer tubes [42, 72] and could be exploited in various engineering and
biomedical applications such as controlled formation of micro-gears [81],
generation of multi-cellular spheroids [86], controlled cell spreading [41],
see also [49, 26]. To the best of our knowledge, swelling-induced micro-
wrinkling at internal surfaces of periodic hydrogels has not been reported in
the literature yet. The sole experimental study coming close to the results
reported here is the one recently documented in [20].

In the following, we first discuss suitable boundary conditions and numerical
setups of boundary value problems to analyze the given hydrogel
microstructures. Next, we provide some numerical studies in order to give
further insight into the response of the voided hydrogel microstructures.
Eventually, we perform a set of computational investigations demonstrating
instabilities in two-dimensional hydrogel microstructures.

####  4.1.1 Boundary conditions in the two-dimensional setting

Figure 3: Unit-cell RVEs of two periodic hydrogel microstructures. We consider
(a) single-phase hydrogels (consisting of hydrogel matrix and voids) and (b)
two-phase hydrogels (consisting of hydrogel matrix, hydrogel coatings, and
voids). The coating is realized as a stiff film that is perfectly bonded to
the matrix (highlighted in green). It is further assumed that the solvent
enters the hydrogel through the internal boundary $\partial{\mathcal{H}}_{0}$,
where we prescribe suitable boundary conditions in terms of the chemical
potential. At the external boundary $\partial{\mathcal{B}}_{0}$, the jump
conditions for the deformation as well as for the solvent-volume flux are
applied as described in 2.3.

Since we consider periodic hydrogel microstructures, all computations are
reduced to a periodic unit-cell RVE as discussed in Section 2.3. Thus, the
jump conditions for the deformation map and the solvent-volume flux according
to (12) are considered at the external boundary of the RVE. In what follows,
we set
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}={\boldsymbol{\mathit{1}}}$
to mimic an effectively constrained hydrogel composite. We further assume that
a solvent diffuses into the hydrogel solely through the boundary of the voids
$\partial{\mathcal{H}}_{0}$ and thus consider the mechanical traction and
chemical potential on $\partial{\mathcal{H}}_{0}$ as

${\boldsymbol{\mathnormal{P}}}\cdot{\boldsymbol{\mathnormal{N}}}={\boldsymbol{\mathit{0}}}\quad\mbox{and}\quad\mu=\lambda(t)\mu_{0}\quad\text{on}\quad\partial{\mathcal{H}}_{0}\quad\text{with}\quad\lambda(t)=\begin{dcases}1-t\,&\text{for}\quad
t\leq 1\,\text{sec}\\\ 0\,&\text{for}\quad t>1\,\text{sec}\end{dcases}\,,$
(47)

where $\mu_{0}$ is the chemical potential of the reference state, see (11). As
can be seen in the latter equation, the chemical potential is linearly
increased from $\mu_{0}$ to zero within one second [7, 16].

In the present work, we consider a time step of $\tau=4\cdot
10^{-3}\,\text{sec}$ in order to determine instability points accurately. We
note that a suitable time-step size is related to the material parameters of
the boundary value problem as well as its numerical discretization [74].

In Fig. 3, we depict illustrations of microstructures considered in the first
part of our studies. The microstructures are discretized using Raviart–Thomas-
type finite elements of the type Q2RT0, where the displacement field is
interpolated using bi-quadratic Lagrangian shape functions associated with
node-based degrees of freedom and the solvent-volume flux field is
interpolated using linear Raviart–Thomas shape functions associated with edge-
based degrees of freedom [7, 65]. Depending on the geometry of the problem,
approximately $5,000$ to $9,000$ finite elements (corresponding to
approximately $20,000$ to $30,000$ nodes) are used for the discretization of a
unit-cell RVE. The used material parameters are summarized in Table 1.

####  4.1.2 Response of voided hydrogel microstructures

Table 1: Material parameters of the considered hydrogel microstructures. Parameter | Name, unit | Value | Eq.
---|---|---|---
$\gamma$ |  Shear modulus, N/mm2 | $0.01$ |  (8)
$\alpha$ |  Mixing modulus, N/mm2 | $20.0$ |  (8)
$\chi$ | Mixing parameter, – | $0.1$ |  (8)
$M$ |  Mobility parameter, mm4/(Ns) | $10^{-4}$ |  (9)
$\epsilon$ |  Penalty parameter, N/mm2 | $10\gamma$ |  (8)
$J_{0}$ | Pre-swollen Jacobian, – | $1.01$ |  (8) $\&$ (9)

$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$

$r_{A}/\text{mm}$

$t/\text{sec}$$10\%$$20\%$$30\%$$40\%$$50\%$macromicro

$0$$0.005$$0.01$$0.015$$0.02$$0.025$$0.03$$0.035$$0.04$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$

$\Delta r_{A}/\text{mm}$

$t/\text{sec}$$10\%$$20\%$$30\%$$40\%$$50\%$macromicro

(a) (b)

$0$$0.02$$0.04$$0.06$$0.08$$0.1$$0.12$$0.14$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$

$h_{void}/\text{(mm\textsuperscript{3}/s)}$

$t/\text{sec}$$10\%$$20\%$$30\%$$40\%$$50\%$macromicro

$-0.4$$-0.35$$-0.3$$-0.25$$-0.2$$-0.15$$-0.1$$-0.05$$0$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$

$\overset{\rule[-0.25835pt]{6.71666pt}{0.38748pt}}{P}{}_{11}/\gamma$

$t/\text{sec}$$10\%$$20\%$$30\%$$40\%$$50\%$macromicro

(c) (d)

Figure 4: _Response of voided hydrogel microstructures versus simulation
time._ a) Radius $r_{A}$ of voids and b) change in the radius $\Delta r_{A}$
of voids measured at the point $A$ as depicted in Fig. 3 with the initial
radius $r_{initial}=\sqrt{f_{0}/\pi}$; c) Solvent volume entering the hydrogel
matrix through the void boundary
$h_{void}=-\int_{\partial{\mathcal{H}}_{0}}{\IKbb
H}\cdot{\boldsymbol{\mathnormal{N}}}\,\mathrm{d}A$; d) Effective first Piola-
Kirchhoff stress $\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}$
normalized with respect to the shear modulus of the matrix $\gamma$. Figure 5:
_Snapshots of the deformed hydrogel microstructure with the void-volume
fraction $f_{0}=50\%$._ The diffusion of solvent volume takes place through
the boundary of the void $\partial{\mathcal{H}}_{0}$, where we prescribe the
chemical potential according to (47). The contour indicates the chemical
potential, where blue corresponds to the minimum and red to the maximum value,
respectively. We observe that at the time instant $t=2.12\,\text{sec}$, the
void boundary has developed creases. It is observed that for this
microstructure the strong ellipticity of the effective response is lost at
$t=1.312\,\text{sec}$.

$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$-3.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$

$r_{A}/\text{mm}$

$\mu^{N}_{void}/\alpha$$10\%$$20\%$$30\%$$40\%$$50\%$

$0$$0.005$$0.01$$0.015$$0.02$$0.025$$0.03$$0.035$$0.04$$-3.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$$0.004$$0.008$$0.012$$0.016$$-0.6$$-0.4$$-0.2$$0$

$\Delta r_{A}/\text{mm}$

$\mu^{N}_{void}/\alpha$

\makebox(0.0,0.0)[]{{}}

(a) (b)

$0$$0.05$$0.1$$0.15$$0.2$$0.25$$-3.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$

$h_{void}/\text{(mm\textsuperscript{3}/s)}$

$\mu^{N}_{void}/\alpha$$10\%$$20\%$$30\%$$40\%$$50\%$

$-0.1$$-0.08$$-0.06$$-0.04$$-0.02$$0$$-3.5$$-3$$-2.5$$-2$$-1.5$$-1$$-0.5$$0$

$\overset{\rule[-0.25835pt]{6.71666pt}{0.38748pt}}{P}{}_{11}/\gamma$

$\mu^{N}_{void}/\alpha$$10\%$$20\%$$30\%$$40\%$$50\%$

(c) (d)

Figure 6: _Response of hydrogel microstructures versus chemical potential
prescribed at the internal boundaries of the voids._ a) Radius $r_{A}$; b)
Change in the radius $\Delta r_{A}$ of the voids against the normalized
chemical potential $\mu^{N}/\alpha$ at point $A$; c) Solvent volume entering
the hydrogel matrix through the void boundary; d) Effective first Piola-
Kirchhoff stress $\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}$
normalized with respect to the shear modulus of the matrix $\gamma$.

To obtain a first insight into the behavior of voided hydrogel
microstructures, we analyze a sequence of RVEs with the void-volume fractions
$f_{0}\in[10\%-50\%]$ and illustrate their response over the simulation time
in Fig. 4. In Fig. 4a and b, we observe that the radius related to a given
void-volume fraction $f_{0}$ changes significantly due to the influx of
solvent, even after loading $t\geq 1\,\text{sec}$. In Fig. 4c, we illustrate
the solvent volume entering through the boundary of the voids over time
according to $h_{void}=-\int_{\partial{\mathcal{H}}_{0}}{\IKbb
H}\cdot{\boldsymbol{\mathnormal{N}}}\,\mathrm{d}A$ with
${\boldsymbol{\mathnormal{N}}}$ being a unit normal vector pointing outward
the hydrogel matrix. We observe that the influx of the solvent volume
increases as we decrease the chemical potential at the void boundary until it
finally reaches its maximum at $t=1\,\text{sec}$ and
$\mu^{N}_{void}=0.0\,\text{N/mm\textsuperscript{2}}$. In Fig. 4d, we have
plotted the effective first Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{22}$
normalized with respect to the matrix shear modulus $\gamma$, see (13)2.
Notably, the normalized effective stress response is almost identical for all
considered microstructures.

In Fig. 6, we depict associated response functions over the applied chemical
potential $\mu^{N}_{void}$ at the void boundary $\partial{\mathcal{H}}_{0}$.
While such an illustration gives insight into the chemically driven response,
is not capable of describing the behavior beyond $t>1\,\text{sec}$. Therefore,
the instability studies documented in the following sections will be discussed
with respect to time rather than applied chemical potential.

Creasing instabilities. Before we come to the analysis of structural
instabilities, we take a look at several snapshots of the unit-cell RVE with
50% void-volume fraction during its deformation, see Fig. 5. There we observe
the development of creasing during the swelling process. A subsequent analysis
of effective moduli indicates a loss of strong ellipticity at
$t=1.312\,\text{sec}$ (black bullet in Fig. 4), i.e., before creasing. In this
context we note that according to [58, 10] creasing is linked to quasi-
convexity at the boundary [2], which is a stronger condition than quasi-
convexity—and consequently rank-one convexity—in the bulk. Thus, there may be
a certain relation between the loss of strong ellipticity of effective moduli
and the development of creasing at internal boundaries of the microstructure;
see also [71, 80, 63]. However, an in-depth study of such a relation is beyond
the scope of the present work.

####  4.1.3 Instabilities in single-phase hydrogels with voids

$0.5$$1.0$$1.5$$2.0$$2.5$$3.0$$3.5$$4.0$$4.5$$5.0$$10$$15$$20$$25$$30$$35$$40$$45$$50$

$t_{crit}/\text{sec}$

$f_{0}/\%$$M=1.0\cdot 10^{-4}\,\text{mm}^{2}/\text{(Ns)}$$M=2.5\cdot
10^{-4}\,\text{mm}^{2}/\text{(Ns)}$$M=5.0\cdot
10^{-4}\,\text{mm}^{2}/\text{(Ns)}$$M=7.5\cdot
10^{-4}\,\text{mm}^{2}/\text{(Ns)}$macro
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$)

$-0.50$$-0.45$$-0.40$$-0.35$$-0.30$$-0.25$$-0.20$$-0.15$$-0.10$$-0.05$$0$$10$$15$$20$$25$$30$$35$$40$$45$$50$$-0.46\,\,$$-0.44\,\,$$-0.42\,\,$$-0.40\,\,$$-0.38\,\,$$10$$15$$20$

$\overset{\rule[-0.25835pt]{6.71666pt}{0.38748pt}}{P}_{11}/\gamma$

$f_{0}/\%$$M=1.0\cdot 10^{-4}\,\text{mm}^{2}/\text{(Ns)}$$M=2.5\cdot
10^{-4}\,\text{mm}^{2}/\text{(Ns)}$$M=5.0\cdot
10^{-4}\,\text{mm}^{2}/\text{(Ns)}$$M=7.5\cdot
10^{-4}\,\text{mm}^{2}/\text{(Ns)}$macro
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$)

\makebox(0.0,0.0)[]{{}}

(a) (b)

$0.8$$1.0$$1.2$$1.4$$1.6$$1.8$$2.0$$2.2$$2.4$$2.6$$10$$15$$20$$25$$30$$35$$40$$45$$50$

$t_{crit}/\text{sec}$

$f_{0}/\%$$\alpha=20\,\text{N/mm\textsuperscript{2}}$$\alpha=40\,\text{N/mm\textsuperscript{2}}$$\alpha=60\,\text{N/mm\textsuperscript{2}}$$\alpha=80\,\text{N/mm\textsuperscript{2}}$macro
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$)

$-0.50$$-0.45$$-0.40$$-0.35$$-0.30$$-0.25$$-0.20$$-0.15$$-0.10$$-0.05$$0$$10$$15$$20$$25$$30$$35$$40$$45$$50$$-0.46\,\,$$-0.44\,\,$$-0.42\,\,$$-0.40\,\,$$-0.38\,\,$$10$$15$$20$

$\overset{\rule[-0.25835pt]{6.71666pt}{0.38748pt}}{P}_{11}/\gamma$

$f_{0}/\%$$\alpha=20\,\text{N/mm\textsuperscript{2}}$$\alpha=40\,\text{N/mm\textsuperscript{2}}$$\alpha=60\,\text{N/mm\textsuperscript{2}}$$\alpha=80\,\text{N/mm\textsuperscript{2}}$macro
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$)

\makebox(0.0,0.0)[]{{}}

(c) (d)

Figure 7: _Onset of instabilities depending on the mobility $M$ and the
chemical parameter $\alpha$._ Results are plotted with respect to the critical
time $t$ (a and c) and the normalized effective first Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}/\gamma_{m}$ (b and
d) with $\gamma$ being the shear modulus of the matrix. We observe that while
the mobility and the chemical parameters have significant influence on the
critical time, the normalized critical stresses are all in the same range.
Note that next to short-wavelength (pattern-transforming) instabilities, we
observe long-wavelength instabilities for the void-volume fractions
$f_{0}=10\%$ and $f_{0}=15\%$.

We now study the onset of instabilities for the voided hydrogel
microstructures as described in Fig. 3a. In particular, we analyze the
influence of the void-volume fraction $f_{0}$, the mobility parameter $M$ and
the chemical parameter $\alpha$ of the matrix. As in the preceding studies,
the remaining material parameters are taken from Table 1.

The onset of instabilities for the considered initial boundary value problem
is depicted in Fig. 7a–d. In Fig. 7a and c, we plot the critical time versus
the void-volume fraction depending on the mobility parameter $M$ and the
chemical parameter $\alpha$, respectively. We observe that both parameters
have a significant impact on the time at which a critical instability occurs.
As expected, instabilities are initiated earlier for larger $M$ (that is, for
faster diffusion through the matrix) and larger $\alpha$ (that is, for an
enhanced chemo-mechanical interaction). In Fig. 7b and d, we plot the
normalized critical first Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{22}$
at the instability points for different values of the mobility and the
chemical parameter. We observe that the critical stresses associated with the
various hydrogel microstructures are all in the same range. This indicates
that the instability point is strongly correlated with the effective
normalized compressive stress state. We further observe that the short-
wavelength (microscopic) instabilities are primary for the microstructures
with the void-volume fractions $f_{0}\in\\{20\%-50\%\\}$ and that the long-
wavelength (macroscopic) instabilities are primary for the microstructures
with the void-volume fractions $f_{0}\in\\{10\%-15\%\\}$. As a result of the
microscopic instabilities, the hydrogel develops the so-called diamond-plate
pattern, resulting in a change of periodicity of the microstructure as
described in [84, 87, 78] and depicted in Fig. 9f.

We note that the macroscopic instabilities were detected by means of
Bloch–Floquet wave analysis with
${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$ and coincide with
the loss of strong ellipticity of the effective mechanical moduli as discussed
in Section 3.2.

####  4.1.4 Instabilities of two-phase hydrogels with voids

$0.7$$0.8$$0.9$$1.0$$1.1$$1.2$$0.02$$0.03$$0.04$$0.05$

$t_{crit}/\text{sec}$

$\omega/\text{mm}$$\gamma_{f}/\gamma_{m}=5.0$$\gamma_{f}/\gamma_{m}=20.0$$\gamma_{f}/\gamma_{m}=100.0$macro
$({\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}})$

$-1.0$$-0.9$$-0.8$$-0.7$$-0.6$$-0.5$$-0.4$$-0.3$$-0.2$$0.02$$0.03$$0.04$$0.05$

$\overset{\rule[-0.25835pt]{6.71666pt}{0.38748pt}}{P}_{11}/\gamma_{m}$

$\omega/\text{mm}$$\gamma_{f}/\gamma_{m}=1.0$$\gamma_{f}/\gamma_{m}=5.0$$\gamma_{f}/\gamma_{m}=20.0$$\gamma_{f}/\gamma_{m}=100.0$macro
$({\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}})$

(a) (b)

Figure 8: _Onset of instabilities for the composite hydrogels with void-volume
fraction $f_{0}=20\%$ for varying shear-modulus ratios $\gamma_{c}/\gamma_{m}$
of the coating ($c$) and the matrix ($m$) as well as coating thickness
$\omega$._ The results are plotted with respect to a) the critical time $t$
and b) the normalized effective first Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}/\gamma_{m}$. Again,
most instabilities induce pattern transformations from one unit-cell RVE to
$2\times 2$ unit-cell RVEs. Furthermore, we observe macroscopic instabilities
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$). Figure 9:
_Buckling modes of two-phase hydrogels depending on the shear modulus ratio of
the coating ( $c$) and the matrix ($m$)._ The periodic hydrogels are composed
of hydrogel matrix, hydrogel coatings with thickness $\omega=0.02\,\text{mm}$,
and voids with volume fraction $f_{0}=20\%$. Depending on the shear-modulus
ratio of coating and matrix $\gamma_{c}/\gamma_{m}$, pattern-transforming
instabilities with and without wrinkling of internal surfaces can be observed.

$0.80$$0.85$$0.90$$0.95$$1.00$$1.05$$1.10$$1.15$$1.20$$1.25$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

$t_{crit}/sec$

$\gamma_{f}/\gamma_{m}$$\omega=0.02\,\text{mm}$$\omega=0.03\,\text{mm}$$\omega=0.04\,\text{mm}$$\omega=0.05\,\text{mm}$macro
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$)

$-1.0$$-0.9$$-0.8$$-0.7$$-0.6$$-0.5$$-0.4$$-0.3$$-0.2$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

$\overset{\rule[-0.25835pt]{6.71666pt}{0.38748pt}}{P}_{11}/\gamma_{m}$

$\gamma_{f}/\gamma_{m}$$\omega=0.02\,\text{mm}$$\omega=0.03\,\text{mm}$$\omega=0.04\,\text{mm}$$\omega=0.05\,\text{mm}$macro
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$)

(a) (b)

Figure 10: _Onset of instabilities for the composite hydrogels with the void-
volume fraction $f_{0}=20\%$ depending on the shear modulus ratio
$\gamma_{c}/\gamma_{m}$ of the coating ($c$) and the matrix ($m$) as well as
the coating thickness $\omega$._ The results are plotted with respect to a)
the critical time $t$ and b) the normalized effective first Piola-Kirchhoff
stress $\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}/\gamma_{m}$.
The short-wavelength buckling patterns span $2\times 2$ unit-cell RVEs.

$1.00$$1.02$$1.04$$1.06$$1.08$$1.10$$1.12$$1.14$$1.16$$1.18$$0.02$$0.03$$0.04$$0.05$

$t_{crit}/\text{sec}$

$\omega/\text{mm}$$M_{f}/M_{m}=0.01$$M_{f}/M_{m}=0.1$$M_{f}/M_{m}=1.0$$M_{f}/M_{m}=10.0$$M_{f}/M_{m}=100.0$macro
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$)

$-0.44$$-0.42$$-0.40$$-0.38$$-0.36$$-0.34$$-0.32$$-0.30$$-0.28$$0.02$$0.03$$0.04$$0.05$

$\overset{\rule[-0.25835pt]{6.71666pt}{0.38748pt}}{P}_{11}/\gamma_{m}$

$\omega/\text{mm}$$M_{f}/M_{m}=0.01$$M_{f}/M_{m}=0.1$$M_{f}/M_{m}=1.0$$M_{f}/M_{m}=10.0$$M_{f}/M_{m}=100.0$macro
(${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$)

(a) (b)

Figure 11: _Onset of instabilities for the composite hydrogels with the shear-
modulus ratio $\gamma_{c}/\gamma_{m}=5$ and the void-volume fraction
$f_{0}=20\%$ depending on the mobility-parameter ratio $M_{m}/M_{c}$ of the
coating ($c$) and the matrix ($m$) as well as the coating thickness $\omega$._
The results are plotted with respect to a) the critical time $t$ and b) the
normalized effective first Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}/\gamma_{m}$. The
pattern-transforming instabilities span $2\times 2$ unit-cell RVEs.

We proceed with analyzing instabilities of two-phase voided hydrogel
microstructures as described in Fig. 3b. As can be seen in this figure, we
take into account microstructures composed of hydrogel matrix $(m)$, voids
($v$) and hydrogel coatings $(c)$ surrounding the voids.

In the following, we investigate the influence of the thickness of the coating
as well as of the ratios between the shear moduli and the mobility parameters
of the hydrogel coating and the hydrogel matrix on the formation of
instabilities. Here, we take into account a fixed void-volume fraction
$f_{0}=20\%$ and boundary conditions as described in Section 4.1. Material
parameters for the matrix and the coating are listed in Table 1.

In Fig. 8, we illustrate the critical time $t_{crit}$ and the effective first
Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{22}$
versus the coating thickness $\omega$ for various values of shear-modulus
ratio $\gamma_{c}/\gamma_{m}$ of the coating and the matrix. In all studies we
assume that the mechanical stiffness of the coating is higher than that of the
matrix $\gamma_{m}=0.01\,\text{N/mm\textsuperscript{2}}$. In the present
example the mobility parameter of both matrix and coating is set to
$M_{m}=M_{c}=10^{-4}\,\text{mm\textsuperscript{4}/(Ns)}$.

From Fig. 8a, we observe that the shear-modulus ratio $\gamma_{c}/\gamma_{m}$
has a significant influence on the critical instability point, such that
microscopic instabilities occur much earlier than in case of the corresponding
single-phase voided microstructures. Similarly, the critical effective
stresses strongly depend on the shear-modulus ratio, see Fig. 8b. In that
context we observe in particular that for $\gamma_{c}/\gamma_{m}=5.0$ and the
coating thickness $\omega<0.025$, instabilities occur much earlier at much
lower critical stresses.

We note that while most of the triggered instabilities show a critical
periodicity of $2\times 2$ unit cells, the observed microcopic patterns are
highly dependent on the thickness of the coating and the ratio of shear
moduli, see Fig. 9. In particular, two distinct kinds of pattern
transformations given by (i) micro-wrinkling and (ii) diamond-plate patterns
(mutual deformation of voids into ellipses) are recorded. In the first case,
pattern transformations due to wrinkling of the coating are detected for
$\gamma_{c}/\gamma_{m}<30$. Here, the number of micro-wrinkles can be tuned by
means of the shear moduli of the matrix and the coating. When the shear-
modulus ratio exceeds the value of $\gamma_{c}/\gamma_{m}=30$, diamond-plate
patterns with smooth internal boundaries as in case of the single-phase
hydrogel structures are observed.888We note that while the Bloch–Floquet
analysis indicated a long-wavelength instability for the shear-modulus ratio
$\gamma_{c}/\gamma_{m}=5.0$ (see black bullet in Fig. 8), the analysis of
effective moduli did not show a loss of strong ellipticity. A further analysis
of strong ellipticity beyond the critical point was not possible due to
missing convergence. Refining the time-step size did not help to overcome this
issue. We refer to [16], who observed a competition between short- and long-
wavelength instabilities for small shear-modulus ratios in plane bilayer
systems. Their results indicate that for a refined mesh, the long-wavelength
instabilities become the primary ones.

In Fig. 10, we further study instabilities for a wide range of shear-modulus
ratios and coating thicknesses
$\omega\in\\{0.02,\,0.03,\,0.04,\,0.05\\}\,\text{mm}$. The graphs describe the
critical time $t_{crit}$ and the effective first Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}$ depending on the
shear-modulus ratio $\gamma_{c}/\gamma_{m}$. Similar to the previous cases, we
observe that for small values $\gamma_{c}/\gamma_{m}<10$, the critical
instability time increases exponentially. As the stiffness of the coating
increases, the corresponding compressive stresses rise and induce earlier
buckling. A slightly different response is observed for the thicker and
stiffer coatings due to a stabilizing effect, see the intersection of
corresponding curves in Fig. 8a.

In Fig. 11, we study the influence of various mobility-parameter ratios
$M_{c}/M_{m}$ of the coating and the matrix on the onset and the type of
instabilities in consideration of a shear-modulus ratio of
$\gamma_{c}/\gamma_{m}=5.0$ with
$\gamma_{m}=0.01\,\text{N/mm\textsuperscript{2}}$ and a mobility parameter of
the coating $M_{c}=10^{-4}\,\text{mm\textsuperscript{4}/(Ns)}$. We illustrate
the critical time $t_{crit}$ and the effective first Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}=\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{22}$
depending on the thickness $\omega$ of the coating. Since the mobility
parameter is associated with the diffusivity of the hydrogel, the solvent
tends to accumulate in the coating when $M_{c}/M_{m}>1$. As a result, the
coating reaches the critical stress faster and the instabilities occur at
lower critical time, see Fig. 11a. In Fig. 11b, we depict the corresponding
effective critical stresses over the coating thickness $\omega$ depending on
the mobility-parameter ratio. In all considered problems, we observe the
buckling instabilities resulting in $2\times 2$ unit-cell RVEs. For a coating
thickness of $\omega=0.02\,\text{mm}$, the long-wavelength instabilities are
detected again.

###  4.2 Three-dimensional hydrogel films

The previous studies were limited to a two-dimensional plane-strain scenario.
We now extend the analysis to the three-dimensional case and investigate
structural instabilities in hydrogel thin films. The considered films comprise
cylindrical hydrogel inclusions that are embedded into a hydrogel matrix in a
periodic manner, see Fig. 12. We assume that the diffusivity of the inclusions
is low compared with the matrix, so that swelling of the inclusions is not
pronounced. The considered setup is motivated from the experimental works [77,
52].

Before we provide the results of a parametric study in consideration of
various geometrical properties of microstructures, we discuss suitable
boundary conditions for the numerical realization of three-dimensional
hydrogel films considering plane periodicity. In that context, we also comment
on relevant ingredients of the numerical discretization and finite-element
implementation.

####  4.2.1 Boundary conditions in the three-dimensional setting

Figure 12: Boundary value problem for the three-dimensional composite hydrogel
thin films. a) The periodic microstructure of a thin film is represented by a
unit-cell RVE, where the jump conditions for the deformation map and the
solvent-volume flux refer only to the side boundaries, see (48). In the three-
dimensional setting it is assumed that the solvent diffuses through the top
and the bottom surfaces, which is realized by applying associated boundary
conditions for the chemical potential. b) A top view of the thin film
indicated that the periodic microstructure can be described by means of unit-
cell RVEs given by right parallelogrammic prism with cylindrical inclusions.
In all considered examples, we assume that the inclusions have lower
diffusivity and higher mechanical stiffness than the matrix.

Because of the periodic nature of the problem, the following computations can
again be reduced to a selected unit-cell RVE. In this connection, we consider
the jump boundary conditions across the lateral boundaries of the RVE

$\left[\kern-1.66669pt\left[{\boldsymbol{\varphi}}\right]\kern-1.66669pt\right]=\widehat{\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}}\left[\kern-1.66669pt\left[{\boldsymbol{\mathnormal{X}}}\right]\kern-1.66669pt\right]\quad\mbox{and}\quad\left[\kern-1.66669pt\left[{\IKbb
H}\right]\kern-1.66669pt\right]\cdot{\boldsymbol{\mathnormal{N}}}=0\quad\text{on}\quad\partial{\mathcal{B}}_{0}^{\textrm{side}}\,,$
(48)

where
$\widehat{\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}}$
is the effective deformation gradient in the plane of periodicity, see also
Fig. 12a. We assume that the hydrogel is mechanically constrained in plane and
that no stresses arise in out-of-plane direction, that is

$\widehat{\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{{\boldsymbol{\mathnormal{F}}}}}=\left[\begin{array}[]{ccc}1&0&0\\\
0&1&0\\\ 0&0&\cdot\\\
\end{array}\right]\quad\mbox{and}\quad{\boldsymbol{\mathnormal{P}}}\cdot{\boldsymbol{\mathnormal{N}}}={\boldsymbol{\mathit{0}}}\
\mbox{on}\ \partial{\mathcal{B}}_{0}^{top,bottom}\,.$ (49)

On the chemical side, the solvent is assumed to enter the hydrogel only
through the top and the bottom surfaces of the film, which is reflected by the
chemical-potential boundary condition

$\mu=\lambda(t)\mu_{0}\quad\text{on}\quad\partial{\mathcal{B}}_{0}^{top,bottom}\quad\text{with}\quad\lambda(t)=\begin{dcases}1-t\,&\text{for}\quad
t\leq 1\,\text{sec}\\\ 0\,&\text{for}\quad t>1\,\text{sec}\end{dcases}\,,$
(50)

where $\mu_{0}$ is the chemical potential of the reference state, refer to
(11). In what follows, we consider a time-step size of $\tau=4\cdot
10^{-3}\,\text{sec}$.

The geometry of the boundary value problem and the corresponding unit-cell RVE
are described in Fig. 12b. Here, the film has the thickness
$\omega=0.02\,\text{mm}$ and plane dimensions of $a\times b$ with a fixed
length of $a=1.0\,\text{mm}$. The inclusions with radius $r=0.25\,\text{mm}$
have lower diffusivity and higher elastic stiffness than the matrix
($M_{i}/M_{m}=0.01$, $\gamma_{i}/\gamma_{m}=10$). The remaining material
parameters of both phases are taken from Table 1.

The microstructures described in Fig. 12b are discretized using hexahedral
finite elements of the type H1H1. Here, the displacement field and the
solvent-volume flux are assumed as nodal degrees of freedom resulting in six
degrees of freedom per node. Depending on the morphology, $5,500$ to $9,000$
finite elements ($8,500$ to $14,000$ nodes) are needed for the discretization
of a unit-cell RVE.

####  4.2.2 Instabilities of three-dimensional two-phase hydrogel thin films

$0.055$$0.060$$0.065$$0.070$$0.075$$0.080$$0.085$$0.090$$0.095$$0.100$$0.105$$6$$7$$8$$9$$10$

$t_{crit}/sec$

$b/\text{mm}$$\beta=45^{\circ}$$\beta=60^{\circ}$$\beta=75^{\circ}$${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$

$-0.0065$$-0.0060$$-0.0055$$-0.0050$$-0.0045$$-0.0040$$-0.0035$$-0.0030$$-0.0025$$6$$7$$8$$9$$10$

$\overset{\rule[-0.25835pt]{6.71666pt}{0.38748pt}}{P}_{11}/\gamma_{m}$

$b/\text{mm}$$\beta=45^{\circ}$$\beta=60^{\circ}$$\beta=75^{\circ}$${\boldsymbol{\mathnormal{k}}}\to{\boldsymbol{\mathit{0}}}$

(a) (b)

Figure 13: _Onset of instabilities for the three-dimensional composite
hydrogel films depending on their morphology._ A shear-modulus ratio of
$\gamma_{i}/\gamma_{m}=10$ between inclusions and matrix as well as a
morphology with $a=1.0\,\text{mm}$ and $r=0.25\,\text{mm}$ are considered.
Instabilities are depicted with respect to a) the critical time $t$ and b) the
normalized effective first Piola-Kirchhoff stress
$\overset{\rule[-0.21529pt]{5.59721pt}{0.32289pt}}{P}_{11}/\gamma_{m}$. Figure
14:  _Buckling modes of composite hydrogels with $\beta=45^{\circ}$ depending
on the width $b$ of the unit-cell RVE._ a)–c) Different buckling patterns
spanning $2\times 2$ unit-cell RVEs; d)–e) Contour plots of the critical
eigenvectors ($\delta v_{3}$; top view) reveal the change in the direction of
the buckling depending on the width of the unit-cell RVE. Figure 15:
_Buckling modes of composite hydrogels with $\beta=60^{\circ}$ depending on
the width $b$ of the unit-cell RVE._ a)–c) Different buckling patterns
spanning $2\times 2$ unit-cell RVEs; d)–e) Contour plots of the critical
eigenvectors ($\delta v_{3}$; top view) reveal the change in the direction of
the buckling depending on the width of the unit-cell RVE. Figure 16:
_Buckling modes of composite hydrogel microstructures composed of unit-cell
RVEs with width $b=0.8\,\text{mm}$._ a) We observe wrinkling of the
microstructure for $\beta=45^{\circ}$ and b) cooperative buckling for
$\beta=60^{\circ}$.

In Fig. 13, we illustrate the onset of instabilities depending on the
morphological parameters of the RVE. We observe that both the width of the
unit cell $b$ as well as the angle $\beta$ have a significant influence on the
critical time and the critical stress. Since the composite effectively becomes
softer with increasing $b$, instabilities occur earlier (cf. Euler buckling of
beams). A similar behavior can be observed for decreasing angle $\beta$, which
results in an effectively larger distance between inclusions in the
$X_{2}$-direction for a fixed value of $b$. While in most cases we observe
buckling into a configuration described by $2\times 2$ unit-cell RVEs, some
computations indicate a transition from short-wavelength to long-wavelength
instabilities. These points are indicated by black bullets in Fig. 13.

In Fig. 14, we illustrate representative buckling patterns for different
widths $b$ of the film for an angle $\beta=45^{\circ}$ (see Fig. 15 for
$\beta=60^{\circ}$). In case of $\beta=45^{\circ}$, we observe wrinkling of
the microstructure, see Fig. 14a–c, where the direction of the wrinkling
pattern can be influenced by adjusting the width $b$ of the unit-cell RVE (see
Fig. 14d–f). For $\beta=60^{\circ}$, we observe substantially different
buckling patterns in the form of cooperative saddle-point shapes, see Fig. 15.
For better insight into the overall buckling patterns, we refer to Fig. 16,
where we plot larger domains of the buckled films.

##  5 Conclusion

In the present work, we studied instabilities in two-dimensional and three-
dimensional composite hydrogels with periodic microstructures. In two
dimensions, we investigated the onset of instabilities for single-phase voided
microstructures as well as two-phase coated microstructures. We observed that
in case of single-phase hydrogels, the periodicity of the microstructure may
change from a unit-cell RVE to a larger RVE containing $2\times 2$ unit cells,
forming the well-known diamond-plate pattern frequently observed in
experiments. For two-phase composites with thin coating surrounding the voids,
we observed that the buckling modes can be substantially different from those
of the single-phase hydrogel microstructures. Here, the change in periodicity
and the development of patterns can be attributed either to wrinkling of the
coating or to cooperative deformation of voids into ellipses. A detailed
parametric study revealed that the observed patterns depend on the thickness
of the coating as well as on the material parameters of coating and matrix.
Our simulations of three-dimensional composite hydrogel films revealed various
buckling patterns depending on the arrangement of small reinforcing particles
within the hydrogel matrix of the microstructure. These findings could provide
access to the design of complex out-of-plane buckling patterns to be tuned by
adjustment of microstructures.

## Acknowledgements

Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
under Germany’s Excellence Strategy – EXC 2075 – 390740016. This funding is
gratefully acknowledged.

## References

  * Bai et al. [2019] Bai, R., Yang, J., Suo, Z., [2019]. _Fatigue of hydrogels_. European Journal of Mechanics - A/Solids 74, 337–370. (doi:https://doi.org/10.1016/j.euromechsol.2018.12.001)
  * Ball and Marsden [1986] Ball, J. M., Marsden, J. E., [1986]. _Quasiconvexity at the boundary, positivity of the second variation and elastic stability_. In: _The Breadth and Depth of Continuum Mechanics_. Springer, pp. 491–517.
  * Bertoldi et al. [2008] Bertoldi, K., Boyce, M., Deschanel, S., Prange, S., Mullin, T., [2008]. _Mechanics of deformation-triggered pattern transformations and superelastic behavior in periodic elastomeric structures_. Journal of the Mechanics and Physics of Solids 56, 2642–2668.
  * Bichara et al. [2012] Bichara, D. A., O’Sullivan, N.-A., Pomerantseva, I., Zhao, X., Sundback, C. A., Vacanti, J. P., Randolph, M. A., [2012]. _The tissue-engineered auricle: past, present, and future_. Tissue Engineering Part B: Reviews 18 (1), 51–61.
  * Bichara et al. [2010] Bichara, D. A., Zhao, X., Hwang, N. S., Bodugoz-Senturk, H., Yaremchuk, M. J., Randolph, M. A., Muratoglu, O. K., [2010]. _Porous poly (vinyl alcohol)-alginate gel hybrid construct for neocartilage formation using human nasoseptal cells_. Journal of Surgical Research 163 (2), 331–336.
  * Biot [1941] Biot, M. A., [1941]. _General theory of three-dimensional consolidation_. Journal of Applied Physics 12 (2), 155–164.
  * Böger et al. [2017] Böger, L., Nateghi, A., Miehe, C., [2017]. _A minimization principle for deformation-diffusion processes in polymeric hydrogels: constitutive modeling and fe implementation_. International Journal of Solids and Structures 121, 257–274.
  * Caló and Khutoryanskiy [2015] Caló, E., Khutoryanskiy, V. V., [2015]. _Biomedical applications of hydrogels: A review of patents and commercial products_. European polymer journal 65, 252–267.
  * Cao and Hutchinson [2012] Cao, Y., Hutchinson, J. W., [2012]. _Wrinkling phenomena in neo-hookean film/substrate bilayers_. Journal of applied mechanics 79 (3).
  * Cawte [2020] Cawte, J., [2020]. _Surface instabilities in nonlinear elastic materials_. Ph.D. thesis, University of Bath.
  * Chester and Anand [2010] Chester, S. A., Anand, L., [2010]. _A coupled theory of fluid permeation and large deformations for elastomeric materials_. Journal of the Mechanics and Physics of Solids 58 (11), 1879–1906.
  * Chester et al. [2015] Chester, S. A., Di Leo, C. V., Anand, L., [2015]. _A finite element implementation of a coupled diffusion-deformation theory for elastomeric gels_. International Journal of Solids and Structures 52, 1–18.
  * Coussy [2004] Coussy, O., [2004]. _Poromechanics_. John Wiley & Sons.
  * Coussy et al. [1998] Coussy, O., Dormieux, L., Detournay, E., [1998]. _From mixture theory to Biot’s approach for porous media_. International Journal of Solids and Structures 35 (34-35), 4619–4635.
  * Dervaux and Amar [2012] Dervaux, J., Amar, M. B., [2012]. _Mechanical instabilities of gels_. Annu. Rev. Condens. Matter Phys. 3 (1), 311–332.
  * Dortdivanlioglu and Linder [2019] Dortdivanlioglu, B., Linder, C., [2019]. _Diffusion-driven swelling-induced instabilities of hydrogels_. Journal of the Mechanics and Physics of Solids 125, 38–52.
  * Ehlers and Acartürk [2009] Ehlers, W., Acartürk, A., [2009]. _The role of weakly imposed dirichlet boundary conditions for numerically stable computations of swelling phenomena_. Computational Mechanics 43 (4), 545–557.
  * Ehlers et al. [2010] Ehlers, W., Acartürk, A., Karajan, N., [2010]. _Advances in modelling saturated soft biological tissues and chemically active gels_. Archive of Applied Mechanics 80 (5), 467–478.
  * Ehlers and Wagner [2019] Ehlers, W., Wagner, A., [2019]. _Modelling and simulation methods applied to coupled problems in porous-media mechanics_. Archive of Applied Mechanics 89 (4), 609–628.
  * Fan et al. [2022] Fan, X., Deng, C., Gao, H., Jiao, B., Liu, Y., Chen, F., Deng, L., Xiong, W., [2022]. _3d printing of nanowrinkled architectures via laser direct assembly_. Science Advances 8 (31), eabn9942.
  * Flory [1953] Flory, P. J., [1953]. _Principles of polymer chemistry_. George Banta Publishing Company.
  * Flory and Rehner [1943a] Flory, P. J., Rehner, J. J., [1943a]. _Statistical mechanics of cross-linked polymer networks i. rubberlike elasticity_. The Journal of Chemical Physics 11 (11), 512–520.
  * Flory and Rehner [1943b] Flory, P. J., Rehner, J. J., [1943b]. _Statistical mechanics of cross-linked polymer networks ii. swelling_. The Journal of Chemical Physics 11 (11), 521–526.
  * Geymonat et al. [1993] Geymonat, G., Müller, S., Triantafyllidis, N., [1993]. _Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity_. Archive of Rational Mechanics and Analysis 122, 231–290.
  * Gurtin et al. [2010] Gurtin, M. E., Fried, E., Anand, L., [2010]. _The mechanics and thermodynamics of continua_. Cambridge University Press.
  * Guvendiren and Burdick [2010] Guvendiren, M., Burdick, J. A., [2010]. _The control of stem cell morphology and differentiation by hydrogel surface wrinkles_. Biomaterials 31 (25), 6511–6518.
  * Guvendiren et al. [2010a] Guvendiren, M., Burdick, J. A., Yang, S., [2010a]. _Kinetic study of swelling-induced surface pattern formation and ordering in hydrogel films with depth-wise crosslinking gradient_. Soft Matter 6 (9), 2044–2049.
  * Guvendiren et al. [2010b] Guvendiren, M., Burdick, J. A., Yang, S., [2010b]. _Solvent induced transition from wrinkles to creases in thin film gels with depth-wise crosslinking gradients_. Soft Matter 6 (22), 5795–5801.
  * Guvendiren et al. [2009] Guvendiren, M., Yang, S., Burdick, J. A., [2009]. _Swelling-induced surface patterns in hydrogels with gradient crosslinking density_. Advanced Functional Materials 19 (19), 3038–3045.
  * Hill [1957] Hill, R., [1957]. _On uniqueness and stability in the theory of finite elastic strains_. Journal of the Mechanics and Physics of Solids 5, 229–241.
  * Hill [1962] Hill, R., [1962]. _Acceleration waves in solids_. Journal of the Mechanics and Physics of Solids 10 (1), 1–16.
  * Hoffman [1995] Hoffman, A. S., [1995]. _”intelligent” polymers in medicine and biotechnology_. Artificial organs 19 (5), 458–467.
  * Hoffman [2012] Hoffman, A. S., [2012]. _Hydrogels for biomedical applications_. Advanced drug delivery reviews 64, 18–23.
  * Hong et al. [2009] Hong, W., Liu, Z., Suo, Z., [2009]. _Inhomogeneous swelling of a gel in equilibrium with a solvent and mechanical load_. International Journal of Solids and Structures 46 (17), 3282–3289.
  * Hong et al. [2008a] Hong, W., Zhao, X., Suo, Z., [2008a]. _Drying-induced bifurcation in a hydrogel-actuated nanostructure_. Journal of Applied Physics 104 (8), 1–4.
  * Hong et al. [2008b] Hong, W., Zhao, X., Zhou, J., Suo, Z., [2008b]. _A theory of coupled diffusion and large deformation in polymeric gels_. Journal of the Mechanics and Physics of Solids 56 (5), 1779–1793.
  * Hong et al. [2008c] Hong, W., Zhao, X., Zhou, J., Suo, Z., [2008c]. _A theory of coupled diffusion and large deformation in polymeric gels_. Journal of the Mechanics and Physics of Solids 56 (5), 1779–1793.
  * Huggins [1941] Huggins, M. L., [1941]. _Solutions of long chain compounds_. The Journal of Chemical Physics 9 (5), 440.
  * Ilseng et al. [2019] Ilseng, A., Prot, V., Skallerud, B. H., Stokke, n. B. T., [2019]. _Buckling initiation in layered hydrogels during transient swelling_. Journal of the Mechanics and Physics of Solids 128, 219–238.
  * Ionov [2014] Ionov, L., [2014]. _Hydrogel-based actuators: possibilities and limitations_. Materials Today 17 (10), 494–503.
  * Jiang et al. [2002] Jiang, X., Takayama, S., Qian, X., Ostuni, E., Wu, H., Bowden, N., LeDuc, P., Ingber, D. E., Whitesides, G. M., [2002]. _Controlling mammalian cell spreading and cytoskeletal arrangement with conveniently fabricated continuous wavy features on poly (dimethylsiloxane)_. Langmuir 18 (8), 3273–3280.
  * Jin et al. [2018] Jin, L., Liu, Y., Cai, Z., [2018]. _Asymptotic solutions on the circumferential wrinkling of growing tubular tissues_. International Journal of Engineering Science 128, 31–43.
  * Kaessmair et al. [2021] Kaessmair, S., Runesson, K., Steinmann, P., Jänicke, R., Larsson, F., [2021]. _Variationally consistent computational homogenization of chemomechanical problems with stabilized weakly periodic boundary conditions_. International Journal for Numerical Methods in Engineering 122 (22), 6429–6454.
  * Kaessmair and Steinmann [2018] Kaessmair, S., Steinmann, P., [2018]. _Computational first-order homogenization in chemo-mechanics_. Archive of Applied Mechanics 88 (1-2), 271–286.
  * Kaneko et al. [2002] Kaneko, D., Gong, J. P., Osada, Y., [2002]. _Polymer gels as soft and wet chemomechanical systems–an approach to artificial muscles_. Journal of Materials Chemistry 12 (8), 2169–2177.
  * Kang and Huang [2010a] Kang, M. K., Huang, R., [2010a]. _Effect of surface tension on swell-induced surface instability of substrate-confined hydrogel layers_. Soft Matter 6 (22), 5736–5742.
  * Kang and Huang [2010b] Kang, M. K., Huang, R., [2010b]. _Swell-induced surface instability of confined hydrogel layers on substrates_. Journal of the Mechanics and Physics of Solids 58 (10), 1582–1598.
  * Kikuchi et al. [2022] Kikuchi, S., Matsubara, S., Nagashima, S., Okumura, D., [2022]. _Diversity of the bifurcations and deformations on films bonded to soft substrates: Robustness of the herringbone pattern and its cognate patterns_. Journal of the Mechanics and Physics of Solids 159, 104757.
  * Lam et al. [2006] Lam, M. T., Sim, S., Zhu, X., Takayama, S., [2006]. _The effect of continuous wavy micropatterns on silicone substrates on the alignment of skeletal muscle myoblasts and myotubes_. Biomaterials 27 (24), 4340–4347.
  * Larsson et al. [2010] Larsson, F., Runesson, K., Su, F., [2010]. _Variationally consistent computational homogenization of transient heat flow_. International Journal for Numerical Methods in Engineering 81 (13), 1659–1686.
  * Li et al. [2019] Li, C. Y., Hao, X. P., Wu, Z. L., Zheng, Q., [2019]. _Photolithographically patterned hydrogels with programmed deformations_. Chemistry–An Asian Journal 14 (1), 94–104.
  * Ma et al. [2019] Ma, P., Niu, B., Lin, J., Kang, T., Qian, J., Wu, Z. L., Zheng, Q., [2019]. _Sequentially controlled deformations of patterned hydrogels into 3d configurations with multilevel structures_. Macromolecular Rapid Communications 40 (3), 1800681.
  * Ma et al. [2016] Ma, S., Yu, B., Pei, X., Zhou, F., [2016]. _Structural hydrogels_. Polymer 98, 516–535.
  * Miehe and Koch [2002] Miehe, C., Koch, A., [2002]. _Computational micro-to-macro transitions of discretized microstructures undergoing small strains_. Archive of Applied Mechanics 72 (4), 300–317.
  * Miehe et al. [2015] Miehe, C., Mauthe, S., Teichtmeister, S., [2015]. _Minimization principles for the coupled problem of darcy–biot-type fluid transport in porous media linked to phase field modeling of fracture_. Journal of the Mechanics and Physics of Solids 82, 186–217.
  * Miehe et al. [1999] Miehe, C., Schotte, J., Schröder, J., [1999]. _Computational micro-macro transitions and overall moduli in the analysis of polycrystals at large strains_. Computational Materials Science 16, 372–382.
  * Miehe et al. [2016] Miehe, C., Vallicotti, D., Teichtmeister, S., [2016]. _Homogenization and multiscale stability analysis in finite magneto-electro-elasticity. application to soft matter EE, ME and MEE composites_. Computer Methods in Applied Mechanics and Engineering 300, 294–346. (doi:http://dx.doi.org/10.1016/j.cma.2015.10.013)
  * Mielke and Sprenger [1998] Mielke, A., Sprenger, P., [1998]. _Quasiconvexity at the boundary and a simple variational formulation of Agmon’s condition_. Journal of elasticity 51 (1), 23–41.
  * Miyoshi et al. [2021] Miyoshi, H., Matsubara, S., Okumura, D., [2021]. _Bifurcation and deformation during the evolution of periodic patterns on a gel film bonded to a soft substrate_. Journal of the Mechanics and Physics of Solids 148, 104272\.
  * Müller [1987] Müller, S., [1987]. _Homogenization of nonconvex integral functionals and cellular elastic materials_. Archive of Rational Mechanics and Analysis 99, 189–212.
  * Okumura et al. [2015] Okumura, D., Inagaki, T., Ohno, N., [2015]. _Effect of prestrains on swelling-induced buckling patterns in gel films with a square lattice of holes_. International Journal of Solids and Structures 58, 288–300. (doi:https://doi.org/10.1016/j.ijsolstr.2015.01.015)
  * Okumura et al. [2014] Okumura, D., Kuwayama, T., Ohno, N., [2014]. _Effect of geometrical imperfections on swelling-induced buckling patterns in gel films with a square lattice of holes_. International Journal of Solids and Structures 51 (1), 154–163. (doi:https://doi.org/10.1016/j.ijsolstr.2013.09.018)
  * Pandurangi et al. [2022] Pandurangi, S. S., Akerson, A., Elliott, R. S., Healey, T. J., Triantafyllidis, N., [2022]. _Nucleation of creases and folds in hyperelastic solids is not a local bifurcation_. Journal of the Mechanics and Physics of Solids 160, 104749\.
  * Pollmann et al. [2021] Pollmann, N., Larsson, F., Runesson, K., Lundgren, K., Zandi, K., Jänicke, R., [2021]. _Modeling and computational homogenization of chloride diffusion in three-phase meso-scale concrete_. Construction and Building Materials 271, 121558.
  * Polukhov and Keip [2020] Polukhov, E., Keip, M.-A., [2020]. _Computational homogenization of transient chemo-mechanical processes based on a variational minimization principle_. Advanced Modeling and Simulation in Engineering Sciences 7 (1), 1–26.
  * Polukhov and Keip [2021] Polukhov, E., Keip, M.-A., [2021]. _Multiscale stability analysis of periodic magnetorheological elastomers_. Mechanics of Materials 159, 103699.
  * Polukhov et al. [2018] Polukhov, E., Vallicotti, D., Keip, M.-A., [2018]. _Computational stability analysis of periodic electroactive polymer composites across scales_. Computer Methods in Applied Mechanics and Engineering 337, 165–197.
  * Rosiak and Yoshii [1999] Rosiak, J. M., Yoshii, F., [1999]. _Hydrogels and their medical applications_. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 151 (1-4), 56–64.
  * Rudzinski et al. [2002] Rudzinski, W. E., Dave, A. M., Vaishnav, U. H., Kumbar, S. G., Kulkarni, A. R., Aminabhavi, T. M., [2002]. _Hydrogels as controlled release devices in agriculture_. Designed Monomers and Polymers 5 (1), 39–65. (doi:10.1163/156855502760151580)
  * Savina et al. [2016] Savina, I. N., Ingavle, G. C., Cundy, A. B., Mikhalovsky, S. V., [2016]. _A simple method for the production of large volume 3d macroporous hydrogels for advanced biotechnological, medical and environmental applications_. Scientific reports 6 (1), 1–9.
  * Silling [1991] Silling, S. A., [1991]. _Creasing Singularities in Compressible Elastic Materials_. Journal of Applied Mechanics 58 (1), 70–74. (doi:10.1115/1.2897181)
  * Sriram et al. [2021] Sriram, S., Polukhov, E., Keip, M.-A., [2021]. _Transient stability analysis of composite hydrogel structures based on a minimization-type variational formulation_. International Journal of Solids and Structures 230-231, 111080. (doi:https://doi.org/10.1016/j.ijsolstr.2021.111080)
  * Tanaka et al. [1987] Tanaka, T., Sun, S.-T., Hirokawa, Y., Katayama, S., Kucera, J., Hirose, Y., Amiya, T., [1987]. _Mechanical instability of gels at the phase transition_. Nature 325 (6107), 796–798.
  * Teichtmeister et al. [2019] Teichtmeister, S., Mauthe, S., Miehe, C., [2019]. _Aspects of finite element formulations for the coupled problem of poroelasticity based on a canonical minimization principle_. Computational Mechanics, 1–32.
  * Triantafyllidis et al. [2006] Triantafyllidis, N., Nestorović, M. D., Schraad, M. W., [2006]. _Failure surfaces for finitely strained two-phase periodic solids under general in-plane loading_. ASME Journal of Applied Mechanics 73(3), 505–515. (doi:doi:10.1115/1.2126695)
  * Trujillo et al. [2008] Trujillo, V., Kim, J., Hayward, R. C., [2008]. _Creasing instability of surface-attached hydrogels_. Soft Matter 4 (3), 564–569.
  * Wang et al. [2017] Wang, Z. J., Zhu, C. N., Hong, W., Wu, Z. L., Zheng, Q., [2017]. _Cooperative deformations of periodically patterned hydrogels_. Science Advances 3 (9), e1700348.
  * Wu et al. [2014] Wu, G., Xia, Y., Yang, S., [2014]. _Buckling, symmetry breaking, and cavitation in periodically micro-structured hydrogel membranes_. Soft Matter 10 (9), 1392–1399.
  * Wu et al. [2013] Wu, Z., Bouklas, N., Huang, R., [2013]. _Swell-induced surface instability of hydrogel layers with material properties varying in thickness direction_. International Journal of Solids and Structures 50 (3-4), 578–587.
  * Yang et al. [2021] Yang, P., Fang, Y., Yuan, Y., Meng, S., Nan, Z., Xu, H., Imtiaz, H., Liu, B., Gao, H., [2021]. _A perturbation force based approach to creasing instability in soft materials under general loading conditions_. Journal of the Mechanics and Physics of Solids 151, 104401.
  * Yin et al. [2009] Yin, J., Bar-Kochba, E., Chen, X., [2009]. _Mechanical self-assembly fabrication of gears_. Soft Matter 5 (18), 3469–3474.
  * Zhang et al. [2009] Zhang, J., Zhao, X., Suo, Z., Jiang, H., [2009]. _A finite element method for transient analysis of concurrent large deformation and mass transport in gels_. Journal of Applied Physics 105 (9), 1–9.
  * Zhang et al. [2018] Zhang, X. N., Wang, Y. J., Sun, S., Hou, L., Wu, P., Wu, Z. L., Zheng, Q., [2018]. _A tough and stiff hydrogel with tunable water content and mechanical properties based on the synergistic effect of hydrogen bonding and hydrophobic interaction_. Macromolecules 51 (20), 8136–8146.
  * Zhang et al. [2008] Zhang, Y., Matsumoto, E. A., Peter, A., Lin, P.-C., Kamien, R. D., Yang, S., [2008]. _One-step nanoscale assembly of complex structures via harnessing of an elastic instability_. Nano letters 8 (4), 1192–1196.
  * Zhao et al. [2008] Zhao, X., Hong, W., Suo, Z., [2008]. _Stretching and polarizing a dielectric gel immersed in a solvent_. International Journal of Solids and Structures 45 (14–15), 4021–4031.
  * Zhao et al. [2014] Zhao, Z., Gu, J., Zhao, Y., Guan, Y., Zhu, X., Zhang, Y., [2014]. _Hydrogel thin film with swelling-induced wrinkling patterns for high-throughput generation of multicellular spheroids_. Biomacromolecules 15 (9), 3306–3312.
  * Zhu et al. [2012] Zhu, X., Wu, G., Dong, R., Chen, C.-M., Yang, S., [2012]. _Capillarity induced instability in responsive hydrogel membranes with periodic hole array_. Soft Matter 8 (31), 8088–8093.