# Third medium finite element contact formulation for pneumatically actuated
systems

Ondřej Faltus<EMAIL_ADDRESS>Martin Horák<EMAIL_ADDRESS>Martin
Doškář<EMAIL_ADDRESS>Ondřej Rokoš<EMAIL_ADDRESS>
###### Abstract

Mechanical metamaterials are artificially engineered microstructures that
exhibit novel mechanical behavior on the macroscopic scale. Active
metamaterials can be externally controlled. For instance, pneumatically
actuated metamaterials can change their mechanical, acoustic, or other types
of effective behavior in response to applied pressure with possible
applications ranging from soft robotic actuators to phononic crystals. To
facilitate the design of such pneumatically actuated metamaterials and
structures by topology optimization, a robust way of their computational
modeling, capturing both pneumatic actuation of internal voids and internal
contact, is needed. Since voids in topology optimization are often modeled
using a soft material model, the third medium contact formulation lends itself
as a suitable stepping stone. In this manuscript we propose a single
hyperelastic material model capable of (i) maintaining exactly a prescribed
hydrostatic Cauchy stress within a void in the pre-contact phase while (ii)
simultaneously acting as a third medium to enforce frictionless contact, which
contrasts existing third medium approaches focused typically solely on
contact. To achieve this goal, we split the overall third-medium energy
density into contact, regularization, and pneumatic pressure contributions,
all of which can be individually controlled and tuned. To prevent distortions
of the compliant third medium, we include curvature penalization in our model.
This improves on existing formulations in terms of compliant third medium
behavior, leading ultimately to better numerical stability of the solution.
Since our formulation is energetically consistent, we are able to employ more
advanced finite element solvers, such as the modified Cholesky algorithm to
detect instabilities. We demonstrate the behavior of the proposed formulation
on several examples of traditional contact benchmarks, including a standard
patch test, and validate it with experimental measurement.

###### keywords:

Contact , third medium , pneumatic actuation , second-order continuum
formulation

††journal: Elsevier

[ctumech] organization=Department of Mechanics, Faculty of Civil Engineering,
Czech Technical University in Prague,addressline=Thákurova 7, postcode=166 29,
city=Prague 6, country=Czech Republic

[tuemom] organization=Mechanics of Materials, Department of Mechanical
Engineering, Eindhoven University of Technology,addressline=P.O. Box 513,
postcode=5600MB, city=Eindhoven, country=The Netherlands

[casita] organization=Institute of Information Theory and Automation, Czech
Academy of Sciences,addressline=Pod Vodárenskou věží 4, postcode=182 00,
city=Prague 8, country=Czech Republic

##  1 Introduction

The concept of a metamaterial has its origins in electromagnetic and optical
applications [1, 2], being defined as a man-engineered material with
properties pushed beyond the limits of materials that typically appear in
nature. Mechanical metamaterials specifically focus on artificially designed
microstructures leading to interesting mechanical characteristics on the macro
scale [3], such as auxeticity (negative Poisson ratio) [1, 4], variable
stiffness [5] and programmable motion [6, 7].

Pattern-forming 2D metamaterials are a specific subclass of mechanical
metamaterials originated from the behavior of a two-dimensional polymer sheet
perforated in regular intervals by voids. These materials are based on a
process of internal in-plane buckling of ligaments between the voids, which
leads to the development of an internal pattern modifying the macroscopic
behavior of the microstructure. The most prominent among these materials are
honeycomb microstructures with circular voids in a hexagonal arrangement,
which buckle into three different patterns depending on the biaxiality ratio
of the applied macroscopic load [8, 9, 10], or square-stacked void
microstructures [11], which exhibit auxeticity in compression [12], possibly
even in a programmable manner [13].

The field of mechanical metamaterials is increasingly focused on design of
actuated microstructures. Active control of metamaterials can be achieved by a
variety of methods, such as electromagnetic or pneumatic actuation [14].
Possible applications of pneumatic actuation are wide-ranging from changes in
electromagnetic properties of sandwiched material designs [15, 16] to soft
robots capable of walking motion [17]. For pattern-forming metamaterials,
pneumatic actuation has been studied by Chen et al. [18] to describe the
relationship between the actuation and the patterning behavior. The use of a
simple actuated microstructure as a gripper with programmable motion has been
demonstrated by Yang et al. [7].

The natural design progression leads from intuition-based designs towards the
use of generative methods such as topology optimization, which traditionally
parameterizes material distribution via a density variable [19]. This shift in
the design paradigm thus necessitates a robust simulation method capable of
capturing pneumatic actuation as well as internal contact in the voids upon
closing, enabling multi-switching designs exploiting contact in the
microstructure [20]. Even though these requirements can be met with
traditional contact methods and simulation of pressure by follower loads [21],
such a formulation leads to a significant increase in computational complexity
due to the need to identify interfaces and their normal vectors (not known a
priori) from the density variable. To address this issue, topology
optimization of pneumatically actuated structures performed recently by
Caasenbrood et al. [22] incorporated pneumatic actuation in the form of a
third medium, in which an arbitrary eigenstrain is prescribed to regions with
a low value of the density variable, leading to an approximate simulation of a
pneumatically actuated void.

The concept of a third-medium method more traditionally appears as an
alternative to classical methods for contact, such as penalty or Lagrange
multiplier methods [23]. The third medium formulation relies on the
introduction of a virtual material, the so-called third-medium material, into
the space between the contacting solid bodies. Properties of this material
need to be chosen such that penetration is prevented without significantly
influencing the solution. Among the advantages of this approach are the
simplicity of meshing without the need to define contact interfaces and
universal enforcement of contact regardless of which parts of the surface come
into contact with each other, which is especially convenient for topology
optimization. Inclusion of a third medium contact method is generally
straightforward for tasks where a void space is already meshed. For instance,
Dev et al. [24] demonstrated that meshing free space is necessary for certain
tasks of topology optimization of electroactive polymers. Disadvantages of the
third medium method may include additional computational expenses due to
meshed voids, overt reliance on material parameters of both bulk and third
medium, and the need for additional regularization to prevent excessive
distortions within void regions.

The third medium method for handling computational contact mechanics was first
proposed by Wriggers et al. [25], even though the fictitious domain of Pagano
and Alart [26] might be considered its precursor. Wriggers et al.’s approach
is largely influenced by traditional contact methods, as the third medium here
is modeled as a material with variable anisotropy, the direction of which
changes with the direction of the contact normal identified from the direction
of the principal stretch. Despite the promising behavior demonstrated in [25],
further research into the area of third medium contact remains, to the
authors’ knowledge, sparse. Among the rare follow-ups, we mention the work by
Bog et al. [27], and an extension to an isogeometric formulation by Kruse et
al. [28]. Recently, a related concept has emerged that revolves around the use
of an Eulerian framework to handle solid-to-solid contact [29], which is
tailored towards easy inclusion of fluid-solid interactions and growth
phenomena.

The appealing applicability of the third medium contact for topology
optimization was recently showcased by Bluhm et al. [30]. Relying on the same
hyperelastic material model for bulk material as well as voids, the authors
proposed a simplified version of a third medium material with a regularization
based on a locally-computed second gradients of displacement at the element
level to prevent oscillations and excessive distortions of the third medium. A
follow-up work with further improvement was recently presented by Frederiksen
et al. [31]. The design principles behind both models inspired the model of
the presented in this paper.

The main contribution of this manuscript lies in integrating third medium
contact with pneumatic actuation and improving on the aforementioned
approaches in several ways: First, we pair an improved third medium
regularization method with precise pneumatic actuation. The second-gradient
regularization employed in our approach avoids unnecessarily penalizing
stretch components of the deformation gradient, leading to an overall more
compliant behavior of the third medium. The pneumatic actuation is represented
exactly as a prescribed Cauchy stress in the material voids. Moreover, our
formulation is energetically consistent. Within the present paper, we limit
the scope to two-dimensional plane strain analysis, as is traditional for
pattern-forming metamaterials. We comment however, on possible challenges
related to extension to 3D.

The remainder of this manuscript is structured as follows. Section 2 provides
an overview of the model with its component terms illustrated by examples
throughout. Section 3 then presents additional numerical examples to showcase
the workings of the method, while Section 4 offers concluding remarks together
with an outline of possible future work.

##  2 Proposed model

The desired third medium model is based on hyperelasticity and tackles both
pneumatic actuation and contact. Its strain energy density function $W$ can be
expected to contain three distinct terms:

$W=\psi_{\mathrm{p}}W_{\mathrm{p}}+\psi_{\mathrm{c}}W_{\mathrm{c}}+\psi_{\mathrm{r}}W_{\mathrm{r}}$
(1)

where the first energy term $W_{\mathrm{p}}$ ensures that the hydrostatic
Cauchy stress across the third medium material is exactly equal to a
prescribed pressure value. The neo-Hookean hyperelastic term $W_{\mathrm{c}}$
is then intended to enforce contact by its rapid response to material volume
approaching zero, ensuring sudden stiffening of the initially compliant third
medium. As pointed out in literature [30], this term is not sufficient on its
own; an additional regularization term $W_{\mathrm{r}}$ based on second
gradients of displacement is necessary to stabilize the response of the
compliant third medium, especially in cases in which a free surface of the
third medium is present. The magnitude of the different strain energy density
terms is governed by the corresponding scalar multipliers $\psi_{\bullet}$,
$\bullet\in\\{p$, $c$, $r\\}$. In the sections below, we discuss the rationale
and derivation of individual regularization terms; the pressure actuation is
discussed in Section 2.1, contact in Section 2.2, whereas the second-gradient
regularization is detailed in Section 2.3.

###  2.1 Pneumatic term

Pneumatic actuation in a microstructural void means an introduction of a
uniform pressure acting on the boundary of the void. The magnitude of this
pressure is externally prescribed as a loading parameter. While naturally,
this would be simulated by deformation-dependent loading, it can also be
performed using a third medium.

This method has recently been discussed by Caasenbrood et al. [22] within a
context of topology optimization in soft robotics. The geometric shape of a
structure in topology optimization is dictated by the density variable $\rho$,
which interpolates between a bulk material ($\rho=1$) and a void
($\rho=\varepsilon$, where $\varepsilon$ is a small number, e.g.
$\varepsilon=10^{-3}$). Finite elements in regions with vanishing density
variable $\rho=\varepsilon$ are selected by Caasenbrood et al. to represent
pneumatically actuated voids. To introduce internal pressure, the global force
vector is enhanced for each of those elements with

$\displaystyle\bm{t}_{e}=\int_{\Omega_{e}}\mathsf{B}_{NL}^{T}\mathsf{D}_{e}\mathsf{e}_{\mathrm{v}}\;\mathrm{d}\Omega$
(2)

where $\Omega_{e}$ is the element domain in 2D, $\mathsf{B}_{NL}$ is the
nonlinear strain-displacement matrix that relates displacements to the Green-
Lagrange strain, and $\mathsf{D}_{e}$ is a constitutive elasticity tensor in
the matrix format calculated using Yeoh strain energy density in [22].
Finally, $\mathsf{e}_{\mathrm{v}}$ is a prescribed constant volumetric strain
that is to be interpreted as the Green-Lagrange strain. The product
$\mathsf{D}_{e}\mathsf{e}_{\mathrm{v}}$ can therefore be interpreted as the
volumetric part of the second Piola-Kirchhoff stress. Even though this
approach may seem appropriate, it is not exactly consistent with the pressure
loading, as shown below.

Herein we propose a new formulation that captures the pressure loading
exactly. To this end, we introduce the following strain energy density term of
the pneumatic part $W_{\mathrm{p}}$:

$\displaystyle W_{\mathrm{p}}=pJ(\bm{F})$ (3)

where the prescribed hydrostatic pressure $p$ works on the volumetric change
$J=\det{\bm{F}}$, with $\bm{F}=\bm{I}+\nabla\bm{u}$ being the deformation
gradient tensor, $\bm{I}$ the second order identity tensor, $\nabla$ the
material gradient, and $\bm{u}$ the displacement vector. Integrating the
strain energy density over the third medium domain $\Omega$ gives the strain
energy

$\Theta_{\mathrm{p}}(\bm{u})=\int_{\Omega}\underbrace{p(\bm{x}){\rm
det}(\nabla\bm{x})}_{W_{p}}{\rm d\Omega}$ (4)

For the sake of generality and to demonstrate the full consistency of our
approach with the follower force approach, we assume temporarily that the
prescribed pressure $p$ can vary in the spatial coordinates
$\bm{x}=\bm{X}+\bm{u}$. Taking the first variation of the strain energy and
exploiting integration by parts leads to

$\displaystyle\delta\Theta_{\mathrm{p}}\left(\bm{u},\delta\bm{u}\right)$
$\displaystyle=$
$\displaystyle\int_{\Omega}J\nabla_{x}p\cdot\delta\bm{u}+p\;{\rm
cof}\bm{F}:\nabla\delta\bm{u}{\rm d\Omega}$ $\displaystyle=$
$\displaystyle\left[\int_{A}\left(p\;{\rm
cof}\bm{F}\cdot\bm{N}\right)\cdot\delta\bm{u}\;{\rm
d}A+\int_{\Omega}\left(J\nabla_{x}p-{\rm Div}\left(p{\rm
cof}\bm{F}\right)\right)\cdot\delta\bm{u}\;{\rm d\Omega}\right]$

where $\delta\bm{F}=\delta(\bm{I}+\nabla\bm{u})=\nabla\delta\bm{u}$ and we
have used the fact that the derivative of the determinant of $\bm{F}$ is its
cofactor. In addition, we have introduced the unit normal in the reference
configuration $\bm{N}$, the reference configuration boundary area $A$, spatial
gradient $\nabla_{x}$, material divergence ${\rm Div}$, and the variation of
the displacement field $\delta\bm{u}$. Moreover, using Nanson’s formula
establishing a relationship between the normal vector in the reference
configuration $\bm{N}$ and the normal vector in the deformed configuration
$\bm{n}$, and Piola’s identity, stating that the divergence of the cofactor is
zero, we arrive at

$\displaystyle\delta\Theta_{\mathrm{p}}$ $\displaystyle=$
$\displaystyle\Bigg{[}\int_{A}p\;\underbrace{\left({\rm
cof}\bm{F}\bm{N}\right)\cdot\delta\bm{u}\;{\rm
d}A}_{\bm{n}\cdot\delta\bm{u}\;{\rm d}a}$ (6) $\displaystyle+$
$\displaystyle\int_{\Omega}\left(\underbrace{J\nabla p\cdot\bm{F}^{-1}-\;{\rm
cof}\bm{F}\cdot\nabla p}_{\bm{0}}-p\underbrace{{\rm Div}({\rm
cof}\bm{F})}_{\bm{0}}\right)\cdot\delta\bm{u}\;{\rm d\Omega}\Bigg{]}$
$\displaystyle=$ $\displaystyle\int_{a}p\;\bm{n}\cdot\delta\bm{u}\;{\rm d}a$

which is the well-known formula for the pressure follower load boundary term
[32, 33], proving the exact representation of pneumatic pressure by our
improved approach. Note that in Equation (6) $a$ denotes the deformed boundary
area.

To further compare the presented approach with that of Caasenbrood et al., the
first Piola-Kirchhoff stress tensor $\bm{P}$ can be derived as

$\bm{P}=\frac{\partial W_{p}}{\partial\bm{F}}=pJ\bm{F}^{-T}$ (7)

which can be transformed to Cauchy stress $\bm{\sigma}$, yielding

$\bm{\sigma}=\frac{1}{J}\bm{P}\cdot\bm{F}^{T}=p\bm{I}$ (8)

resulting in a hydrostatic Cauchy stress with magnitude $p$. From this we can
also compute the second Piola-Kirchhoff stress $\bm{S}$ as

$\displaystyle\bm{S}=\bm{F}^{-1}\cdot\bm{P}=pJ\bm{F}^{-1}\cdot\bm{F}^{-T}$ (9)

It may be now apparent from Equation (2) that through prescribing a constant
Green-Lagrange strain $\bm{\varepsilon}_{v}=\bm{I}$, Caasenbrood et al. [22]
essentially propose to prescribe the hydrostatic part of the second Piola-
Kirchhoff stress tensor, rather than the Cauchy stress tensor, as an arbitrary
constant

$p_{S}=\frac{1}{3}{\rm
tr}(\bm{D}_{e}:\underbrace{\bm{\varepsilon}_{v}}_{\bm{I}})=\frac{1}{3}{\rm
tr}(\bm{S})=pJ{\rm tr}(\bm{F}^{-1}\cdot\bm{F}^{-T})$ (10)

For a deformation-independent constitutive tensor $\bm{D}_{e}$ this however
does not lead to constant hydrostatic pressure $p=p_{S}$, but rather to a
deformation-dependent relation

$p=\frac{p_{S}}{J{\rm tr}(\bm{C}^{-1})}$ (11)

where $\bm{C}=\bm{F}^{T}\cdot\bm{F}$ is the right Cauchy-Green deformation
tensor. Hence, using Equation (2) necessitates control of $\bm{D}_{e}$ such
that $p$ in Equation (11) remains constant.

Note that in the following, we expect pressurization by air pumps. Thus, we
use $\Delta p$ instead of $p$ to represent a pressure difference of the air in
the void as compared to the atmospheric pressure. Since we assume undeformed
voids under atmospheric pressure, $\Delta p=0$ corresponds to the initial
state. Consequently suction corresponds to negative pressure difference while
inflation is a result of positive pressure difference.

The contribution of the pneumatic term to the overall energy density is
weighed by the scalar parameter $\psi_{\mathrm{p}}$. Outside of topology
optimization, this material parameter is unnecessary and should always be set
to $1$. In topology optimization, where the transition from void to bulk
material is smooth and only the void portion of the domain is actuated,
$\psi_{\mathrm{p}}$ depends on the density variable $\rho$.

It is worth noting that while the presented pneumatic energy density term of
Equation (3) indeed results in hydrostatic Cauchy stress within the third
medium, it cannot be used alone for numerical reasons, as it leads to a
material with no shear stiffness. An additional energy term is needed for
regularization. Its presence then distorts the perfect representation of the
Cauchy stress, hence its magnitude needs to be kept comparatively low to
ensure that the activated pneumatic term dominates the material energy and the
contribution of other terms is negligible.

###  2.2 Contact terms

The contact-enforcing part of the energy density $W_{\mathrm{c}}$ is based on
a compressible neo-Hookean material model [34, 35] that is split into its
volumetric $W_{\mathrm{c,vol}}$ and isochoric $W_{\mathrm{c,iso}}$ parts as

$W_{\mathrm{c}}=\underbrace{\ln^{2}{J}}_{W_{\mathrm{c,vol}}}+\underbrace{\left(J^{-2/3}I_{1}-3\right)}_{\mathrm{W_{c,iso}}}$
(12)

with $I_{1}=\mathrm{tr}\left(\bm{C}\right)$ the first invariant of the right
Cauchy-Green deformation tensor.

The primary criterion for selecting a third medium contact material lies in
its overall high compliance, characterized by stiffness coefficients several
orders of magnitude lower than that of the bulk material constituting the
contacting bodies. It is crucial, however, that initially high compliance is
complemented by a swift transition to high stiffness upon total compression of
the third medium, effectively preventing any form of penetration. This is
satisfied, in a general 3D case, e.g., by the logarithmic nature of the
volumetric term in Equation (12). Under plane strain conditions, however,
where the out-of-plane stretch is fixed to $1$, no purely volumetric
deformation is ever possible (apart from the trivial case $\bm{F}$ =
$\bm{I}$). Thus, the isochoric part by itself is active in response to in-
plane compression and approaches infinity; moreover, it does so at a larger
rate and more abruptly than the logarithmic term under the same conditions.
This behavior is illustrated in Figure 1 in comparison of plane strain and
full 3D cases. As a consequence, the inclusion of the volumetric term
$\ln^{2}J$ of Equation (12) is not strictly necessary in plane strain analysis
as presented below in this paper, and hence will be left out. The volumetric
part becomes necessary for 3D applications, as well as in cases where it is
desirable for the third medium to correspond in its material formulation to a
neo-Hookean bulk material, such as for topology optimization.

Figure 1: Behavior of the neo-Hookean strain energy density terms in response
to in-plane biaxial and full triaxial compression. Under in-plane compression,
due to the $J^{-2/3}$ term, the isochoric energy term tends to infinity more
rapidly than the volumetric term. The isochoric term alone can therefore act
as the contact term for plane strain geometries.

The neo-Hookean terms of the third medium model are closely related to the
material model of the bulk, i.e., the contacting bodies. Their treatment
differs for applications in and outside of topology optimization.

Outside of topology optimization, the material law of the bulk is arbitrary,
and it is only necessary to ensure that the third medium is comparatively
compliant. The material parameter $\psi_{\mathrm{c}}$ governing the magnitude
of the neo-Hookean terms’ contribution is thus set inside the third medium to
a small constant positive value $\gamma$, which can be understood as the
contact stiffness of the material. The choice of $\gamma$ in relation to the
material parameters of the bulk contacting bodies influences the accuracy of
the solution, with larger values leading to a larger gap being left in the
geometry of the converged solution. Smaller values, conversely, can lead to
numerical difficulties with the solution because of badly conditioned global
stiffness matrix. Natural comparison invites itself to the penalty contact
method, where a similar behavior is observed: the larger the penalty
parameter, the more accurate the solution and the larger the issues with
stiffness conditioning. On the other hand, while the gap arising from the
inaccuracy of the penalty method is negative (i.e., overlap), the gap caused
by the inaccuracy of the third medium method is positive.

In topology optimization, it is advantageous to smoothly transition between
bulk and void (third medium). The assignment of a region to either the void or
bulk material is typically governed by the density $\rho$, i.e., a smooth
field representing a gradual transition between the two materials. In a
context in which the bulk material is neo-Hookean, $\psi_{\mathrm{c}}$ can
thus be a function of the density variable or directly the density variable
itself, reducing the stiffness of the base hyperelastic model in the void
region and increasing it in the bulk region. Thus, a smooth transition between
bulk and void is achieved. To satisfy this relationship between bulk and void,
other hyperelastic material models can also supply the "neo-Hookean term" in
the third medium model. The only condition is that the response to compression
to zero volume needs to rapidly approach infinity, compare Figure 1.

The aforementioned terms represent a way to simulate pneumatic actuation with
the pneumatic term, a way to stabilize the pneumatic term with the isochoric
neo-Hookean term, and finally even a way to enforce contact, since the
isochoric neo-Hookean term ensures that as well under plane strain conditions.
With these abilities, it is now possible for the third medium model to be
tested. In the following example, loosely motivated by soft robotics similarly
to Caasenbrood et al. [22], we are going to see the application of pneumatic
suction to a rectangular box with a closed void, examining the influence of
the material parameters of the third medium on the behavior of the solution. A
rectangular void is positioned inside of a rectangular box, as pictured in
Figure 2. The dimensions of the finite sample are $B\times
H=2\;$\mathrm{m}$\times 0.5\;$\mathrm{m}$$ with the walls
$t=$0.1\text{\,}\mathrm{m}$$ thick. Both bulk material and void are
discretized with a regular mesh of 8-node quadratic quadrilateral elements
with 9 Gauss integration points. For this and all following examples, unless
specified otherwise, the bulk material is represented by a neo-Hookean solid
with the bulk and shear moduli of $K=$2000\text{\,}\mathrm{MPa}$$ and
$G=$10\text{\,}\mathrm{MPa}$$, respectively.

Figure 2: Pneumatic box subjected to suction within the internal void:
geometry and boundary conditions.

In Figure 3, results of three simulations of the aforementioned example are
illustrated. They differ in the parameters used for the third medium. In the
first case, pictured in Figure 3a (left), only the pneumatic and isochoric
terms were used, with $\psi_{\mathrm{c}}=\gamma=10$ and $\gamma_{p}=1$. The
calculation converged in the whole loading range, however this relatively
large contact stiffness resulted in imprecisely identified contact onset, as
can be seen from the smooth shape of the associated pressure-gap diagram in
Figure 3b (red dashed line). Reducing the contact stiffness $\gamma$ further
to $1$ leads on the other hand to a loss of convergence of the numerical
solver due to the oscillating third medium. The last converged state can be
seen in Figure 3a (center) and the pressure-gap diagram cut short in Figure 3b
(black line with a black cross for the last converged state).

The problem of the oscillating third medium can be averted by introducing a
regularization term, explained in the next section, where a new regularization
stiffness parameter $k_{r}$ will be introduced to stabilize the third medium,
see Equation (19). From the results we can already see that such a
regularization leads to good accuracy while stabilizing the third medium. The
satisfactory deformed shape and sharp pressure-gap diagram are visible in
Figure 3a (right) and from its associated blue dotted line in Figure 3b,
respectively. Ignoring this regularization leads either to a too stiff
response or a failure of the numerical solver.

---
(a)
(b)
Figure 3: Simulation of monotone suction in a single pneumatic cell. (a) Mesh
and deformed configuration for $\gamma=10$ (left), for $\gamma=1$ (center),
where only the last converged configuration is shown, with contact penetration
and distorted third medium, and for $\gamma=1$ with the use of additional
regularizing terms (right). (b) Comparison of pressure-gap diagrams for
different values of $\gamma$ and regularization (only relevant area shown).

Note, in addition, that the demonstrated behavior with the third medium’s
computational stability depending either on sufficiently high contact
stiffness at a cost of accuracy, or on additional regularization, does not
only apply to pneumatic actuation; it is very similar for pure enforcement of
contact conditions, too.

To show this, the sample from Figure 3 is now loaded by a force acting in the
center of the top boundary with no internal pressure applied (compare Figure
4), leading to a large enough deformation to enforce contact between the top
and bottom walls. Results of such a simulation are pictured in Figure 5. For
$\gamma=10$ we can again see a smooth force-displacement diagram (red dashed
line in Figure 5b), capturing the contact imprecisely. Reducing $\gamma$ to
$1$ leads to oscillations, see Figure 5a center and the inset, with the force
displacement diagram cut short by a loss of convergence (black line). Addition
of regularization, see Section 2.3 below, solves this issue (Figure 5a right
and blue dotted line in the force-displacement diagram).

Figure 4: Geometry and boundary conditions of the second loading case, i.e., a
box with an internal void subjected to displacement of its top boundary.
---
(a)
(b)
Figure 5: A rectangular cell loaded by a prescribed displacement, without
pneumatic actuation. (a) Mesh and deformed configuration for $\gamma=10$
(left), for $\gamma=1$ (center), where only the last converged configuration
is shown, with contact penetration and distorted third medium (see the inset),
and for $\gamma=1$ with the use of additional regularizing terms (right). (b)
Comparison of force-displacement diagrams for different values of $\gamma$ and
regularization (only vicinity of contact initiation at $u_{y}=0.3$\mathrm{m}$$
shown).

###  2.3 Regularization term

The necessity of additional regularizing energy has been identified previously
by Bluhm et al. [30] for cases where the third medium has a free surface.
However, as the previous two examples have demonstrated, regularization is
needed even for enclosed void regions. Due to the generally very small
stiffness of the third medium in comparison to materials of the surrounding
solid bodies, it tends to deform rather excessively in the pre-contact phase.
The third medium is fictional and does not describe any real material or
structure, which means that this is not a problem in the physical sense; it
can, however, become a problem in the computational sense. Excessive
distortions limit the ability of the third medium to be volumetrically
compressed while properly capturing contact. Additionally, oscillating
behavior can also take place in the third medium displacements, especially for
structured finite element discretizations, where multiple energetically
equivalent deformations exist for the third medium.

One way to ensure reasonable behavior of the third medium in the pre-contact
phase is to penalize the shape distortion of its finite elements. Bluhm et al.
[30] achieved this by introducing a regularization term that penalizes second
gradients of the displacement, formulating a regularizing energy term in the
form of

$W_{\mathrm{r}}=\frac{1}{2}c\nabla\bm{F}\mathbin{\vdots}\nabla\bm{F}$ (13)

where $c$ is a material parameter. This leads to a much better-behaved third
medium, which deforms uniformly and does not tend to bulge out of free
surfaces nor oscillates, see Figure 6 for a comparison on a simple benchmark
(the detail setting of which is explained later in Section 3.2). A
disadvantage, however, is that gradients of stretches are also penalized by
this term, and thus with deformation the penalization term stiffens the third
medium as a whole. While this can be addressed by recursive adjustments of
material parameters, it is, in the authors’ experience, rather cumbersome.
Additionally, the formulation by Bluhm et al. introduces a switching term
serving to restrict the regularization to compression only. Moreover, this
switch, in the form of the $e^{-5J}$ multiplier, only appears in the first
displacement variation of strain energy density, causing the model to become
energetically inconsistent, i.e., without an energy potential. This may lead,
in turn, to undesirable energy dissipation under cyclic loading and to a non-
symmetric stiffness tangent. For this reason, this particular term has already
been removed in a follow-up work [31].

Figure 6: Effect of second gradient regularization introduced by Bluhm et al.
[30] on pre-contact deformation of third mediums with free surface.

In this contribution, we opt for a different approach based on a penalization
of curvatures only, expressed as material gradients of rotations. While
rotations $\bm{R}$ can be obtained from the deformation gradient $\bm{F}$ by
polar decomposition $\bm{F}=\bm{R}\cdot\bm{U}$, this expression is not easily
differentiable either by material coordinates $\bm{X}$ to obtain the gradients
$\nabla\bm{R}$, or by deformation gradient $\bm{F}$ to calculate stresses and
material stiffness necessary for any finite element implementation. A possible
approach to this problem has been recently identified in the works of Smith et
al. [36] and Poya et al. [37]. For the purposes of penalizing a fictitious
third medium, we find it substantially easier to just penalize the gradient of
an alternative rotational measure associated with the material spin, which we
have denoted $\bm{Q}$ to avoid confusion with the rotation tensor $\bm{R}$.
The tensor $\bm{Q}$ is defined by the incremental relation in pseudo-time $t$

$\displaystyle\bm{Q}(t=0)=\bm{I}$ (14)
$\displaystyle\dot{\bm{Q}}=\bm{w}_{X}\cdot\bm{Q}$ (15)

where $\bm{I}$ is the second order identity tensor, $\dot{\bm{Q}}$ is the
pseudo-time derivative of $\bm{Q}$, and $\bm{w}_{X}$ is the material spin,
defined as the skew-symmetric part of a material gradient of pseudo-velocity
$\bm{v}$:

$\bm{w}_{X}=\frac{1}{2}\left(\nabla\bm{v}-\nabla\bm{v}^{T}\right)=\frac{1}{2}(\dot{\bm{F}}-\dot{\bm{F}}^{T})$
(16)

with $\dot{\bm{F}}$ being the rate of the deformation gradient. To simplify
the implementation even further, it is possible to substitute for $\bm{Q}$ its
logarithm $\ln\bm{Q}$ without a significant change in effect, according to our
numerical experiments. With this modification, it is now possible to construct
the term $\nabla\ln\bm{Q}$ by a simple additive incremental relation. This
term then serves as an easier-to-compute proxy of $\nabla\bm{R}$. Clearly, the
simplifications included in the derivation of this term represent a compromise
between computational convenience and physical accuracy in capturing
curvature, but we believe it is acceptable due to the fictitious nature of the
third medium. A detailed discussion of this reasoning and the description of
the regularization term’s implementation can be found in A. The penalization
thus follows from a quadratic strain energy density term

$W_{r}=\frac{1}{2}c\left(\nabla\ln\bm{Q}\mathbin{\vdots}\nabla\ln\bm{Q}\right)$
(17)

The regularization term based on the $\bm{Q}$-rotation indeed leads to element
shape regularization without affecting the volumetric stiffness of the third
medium. Consider Figure 7, where two simulations are compared with the same
material parameters on the C-shape specimen, showing the additional stiffness
exerted by the original $\nabla\bm{F}$ regularization compared to the
$\nabla\ln\bm{Q}$ regularization, visible in the larger contact gap.

Figure 7: Comparison of the effect of regularization terms $\nabla\bm{F}$ and
$\nabla\ln\bm{Q}$. The regularization based on rotation ($\nabla\ln\bm{Q}$)
leads, for the same material parameters, to a lower stiffness of the third
medium and hence smaller contact gap ($\varepsilon<<\delta$).

There is one drawback, however: due to the larger freedom in element
deformation caused by the lack of stretch and shear penalization, elements may
approach a deformed state with locally zero volume in certain geometries of
the third medium. The resistance of the contact part of the third medium model
to this condition then leads either to numerical instabilities (consider
Figure 8) or, for stiffer material parameters, to the creation of an
artificially stiff band across the geometry of the third medium and
consequently premature enforcement of contact constraints. Fortunately, these
issues are easily remedied by an additional penalty on the gradient of $J$,
forcing uniformity of volume change across each element. Thus, the final
regularization term we propose consists of a sum of two quadratic forms:

$W_{\mathrm{r}}=\frac{1}{2}c\left(\nabla\ln\bm{Q}\mathbin{\vdots}\nabla\ln\bm{Q}+\nabla
J\cdot\nabla J\right)$ (18)

Figure 8: Comparison of the effect of regularization terms. The regularization
based on rotation ($\nabla\ln\bm{Q}$, left), can lead to numerical instability
due to flipped elements and hence negative finite element Jacobians. This is
avoided by adding an enforcement of uniform volume change across element
(penalizing $\nabla J$, right).

The proposed regularization term meets expectations and allows for
regularization without influencing the contact stiffness and with stable
numerical behavior. Uniform rotation, i.e., zero curvature and uniform volume
change, are enforced. Contrary to Bluhm et al.’s approach of penalizing second
displacement gradients as a whole, however, stretch and certain modes of
symmetric extension or contraction of the elements are left unpenalized, which
improves the compliant behavior of the third medium. Detailed comparison of
the different approaches to regularization with numerical examples can be
found in Section 3.2. Apart from being much less dependent on the contact
stiffness, the present approach also tends to perform better with more
compliant third mediums, which greatly increases the accuracy of contact
constraint enforcement.

Because of this regularization term, any finite element use of the third
medium model necessitates discretization into at least quadratic elements. It
is also worth noting that such an implementation is local only in the
rotational and Jacobian gradients, with no continuity required of the second
gradients across element boundaries. For practical reasons, higher continuity
is not enforced, since this would require using splines as interpolation
functions [38] or introduction of an additional gradient field coupled weakly
with displacements [39]. While this is a physical violation in a model of a
realistic material, it is of no consequence here due to the fictitious nature
of the third medium.

The material parameter $\psi_{\mathrm{r}}$ should reflect the contact
stiffness already included in the neo-Hookean terms’ multiplier. In all
examples here, it is set to

$\psi_{\mathrm{r}}=\gamma{k}_{r}$ (19)

where $\gamma$ is the previously discussed contact stiffness appearing in
$\psi_{c}$ and $k_{r}$ is a regularization stiffness coefficient set usually
to a value between $10^{2}$ and $10^{6}$. To maintain unit consistency of
material parameters, it is possible to include some characteristic stiffness
and length parameters of the bulk material within $c$. Within topology
optimization, a switch based on the density variable is necessary to disable
the regularization term for bulk regions.

##  3 Numerical experiments

We demonstrate efficacy of the proposed third medium model on three benchmark
and test problems: a patch test (Section 3.1), a self-contact example inspired
by Bluhm et al. [30] (Section 3.2), and a pneumatic actuation of a pattern-
forming metamaterial sample (Section 3.3). In the spirit of topology
optimization and pattern-forming metamaterials, we use a neo-Hookean material
model for the solid material, which could represent a silicone rubber sample
produce by casting in a 3D-printable mold. The bulk and shear moduli
$K=$2000\text{\,}\mathrm{MPa}$$ and $G=$10\text{\,}\mathrm{MPa}$$ are used.
The model is defined with the strain energy density in the form

$W_{\mathrm{NH,bulk}}=\frac{K}{2}\ln^{2}{J}+\frac{G}{2}\left(J^{-2/3}I_{1}-3\right)$
(20)

The discretization of each problem consist of quadratic 8-node quadrangular
finite elements with 9 Gauss integration points. We found that a Lobatto
integration rule (i.e., integration points at element nodes) improves the
precision of contact enforcement. In the authors’ experience, this improvement
is particularly significant for coarser meshes. With standard integration
rules, interpenetration of contact surfaces might occur, its magnitude being
approximately the distance between the element boundary and the integration
point of the third medium element. For meshing, the open-source GMSH software
has been used [40]. All presented examples have been calculted under plane
strain conditions.

For the patch test benchmark, the standard Newton solver has been used to
solve the system. For the subsequent examples, a modified Newton algorithm
[41, 42] has been employed, using Cholesky decomposition to modify the Hessian
to repel the solution from local maxima and saddle points towards a local
minimum. Additionally, a smart step length reduction scheme was in effect to
allow the simulation to proceed in cases where large deformations,
instabilities, and oscillations within the third medium deformed geometry may
occur. When the solver had not converged within a specified maximum of Newton
iterations, step length was halved, this process repeating for a maximum of 14
times before concluding the solver had indeed failed. Successful convergence
of a Newton step lead to an increase of step length by the factor of $1.5$,
conversely.

###  3.1 Patch test

Performance of a computational contact method is traditionally verified using
a standard contact patch test. Here, we adopted a patch test proposed
originally by Crisfield [43], cf. Figure 9, in which a block of an elastic
material of dimensions $1$\mathrm{m}$\times 1$\mathrm{m}$$, cut in two by a
horizontal contact interface, is loaded by a uniform downward displacement on
the top boundary of the upper layer, thus pressing the top part against the
bottom one. This should lead to a uniform stress distribution in the material,
identical on both sides of the contact interface. For small strain elasticity
in a two-dimensional domain under plane strain conditions, an analytical
solution of $\sigma_{22}=E\Delta/(1-\nu^{2})$ can be derived, with
$\sigma_{22}$ being the vertical normal component of the stress tensor, $E$
the elasticity modulus, $\Delta$ the imposed displacement (corresponding to
vertical strain due to the unit dimensions of the block), and $\nu$ the
Poisson’s ratio. If the interface nodes of the finite element mesh are aligned
across the contact interface, the test is typically fulfilled for all contact
methods. If they are misaligned, however, some commonly used methods may fail
unless sufficiently modified, resulting in a horizontally non-uniform pressure
distribution, see, e.g., [44].

For the purposes of third medium contact, a slight modification is necessary:
an initial gap has to be introduced in the contact interface to accommodate a
third medium mesh inbetween the two elastic blocks, see Figure 9. We set the
gap to be $0.1\text{\,}\mathrm{m}$ high, with the blocks placed above and
below the gap without being reduced in size. Total prescribed displacement on
the top boundary is $\Delta=-$0.2\text{\,}\mathrm{m}$$; that is,
$-$0.1\text{\,}\mathrm{m}$$ to close the gap, and $-$0.1\text{\,}\mathrm{m}$$
to induce a vertical strain of $\varepsilon_{22}=0.1$ in the elastic block.
Parameters of the third medium model used for this benchmark are $\gamma=1$
and $k_{r}=2\times 10^{3}$, recall Equation (19). Since a large strain
hyperelastic material is used for the bulk, instead of relying on a small-
strain analytic solution, the comparison will be made with a numerical result
from an undivided elastic block without any contact interface influence (i.e.,
a contact-free solution).

Figure 9: Geometry and boundary conditions of the standard patch test [43] as
adapted for third medium contact.

Figure 10 shows results of two patch test variants, one with aligned interface
nodes and the other with misaligned interface nodes. Their meshes and
corresponding deformed states can be seen in Figures 10a-b. Figure 10c,
combined for both cases, reveals a complete overlap among the uniform stress
distributions of the contact-free solution and of the benchmark results. Those
were sampled along the top and bottom boundaries of the block, as well as
along the gap. The third medium method thus fulfills the patch test for both
the aligned and misaligned node cases. Admittedly, a detailed look at the
resulting stress distribution reveals slight imperfections of the solutions,
see the inset of Figure 10c. Oscillations in the stress distributions in the
vicinity of the gap can be observed, while all the stress curves tend to rise
slightly close to the vertical edge of the domain. This effect is caused by
the small free surface of the third medium on both sides of the gap, where the
third medium elements have the freedom to bulge out. It can be reduced by a
more aggressive second gradient regularization, i.e. by setting a larger
$k_{r}$, or completely avoided by a boundary condition preventing horizontal
movement of the free surface. Both this free surface effect and the
oscillations, however, remain in magnitudes negligible compared to the
solution.

|
---|---
(a) | (b)

---
(c)

Figure 10: Patch tests of the presented third medium contact formulation
conducted on an elastic block with a gap with uniform displacement applied to
the top boundary. Illustration of the undeformed and the deformed
configurations (bulk material in blue, third medium in gray) for two test
variants: (a) aligned interface nodes, (b) misaligned interface nodes. (c)
Uniform stress propagation through the third medium, resulting in a uniform
stress distribution in the elastic block, with identical magnitudes on the top
and bottom; comparison with contact-free benchmark solution. Note that a
slight stress increase and oscillations around the gap near the interface are
present, negligible compared to its average.

###  3.2 Self-contacting C-shape example

Bluhm et al. [30] have recently developed a particular benchmark for third
medium contact, which has proven useful to demonstrate the characteristic
behavior of third medium contact models. It consists of a rectangular
structure in the shape of the letter C elongated along the horizontal axis, as
pictured in Figure 11. The dimensions in the examples presented here forth are
$B=$1\text{\,}\mathrm{m}$$, $H=$0.5\text{\,}\mathrm{m}$$ and
$t=$0.2\text{\,}\mathrm{m}$$. The entire vertical left edge of the C-shape
domain is fixed. The rightmost corner of the upper horizontal cantilever is
loaded with a prescribed displacement $u_{y}$, forcing it to bend up to the
point of initial contact with the lower cantilever and beyond. To capture this
contact, the inner area of the C-shape specimen has to be filled with the
third medium material. The third medium, therefore, has a relatively long free
surface along the right-hand side vertical boundary of the domain. The
particular usefulness of this simple benchmark is twofold: first, the
rectangular geometry of the domain allows for perfectly regular meshes, on
which all distortions of the third medium and its boundary are clearly
visible; second, the contact point itself lies on the free surface of the
third medium.

Figure 11: Geometry and boundary conditions of the C-shape benchmark example
originally proposed by Bluhm et al. [30].

Among the peculiarities of this benchmark problem is the necessity to include
at least one extra column of third medium elements along the free surface.
This has already been demonstrated by Bluhm et al. [30], and remains true for
the present material model as well. The reason behind this requirement is
purely numerical, caused by the behavior of the deformed third medium free
surface. This is illustrated in Figure 12, where a very coarse third medium
mesh, consisting only of three elements per gap height without this extra
column is shown pressed up to the point of contact and beyond. If the vertical
free surface bulges ever so slightly inwards, penetration occurs, with the
contact point passing through the boundary of a non-neighboring third medium
element. Thus, the third medium is prevented from arranging itself in the
contact gap and activating the resistance to compression of the hyperelastic
law.

Figure 12: Results of a tentative C-shape benchmark without an extra column of
third medium elements. With loading, the free surface (the black dashed line
in the insets) bulges inward (left column), leading to interpenetration of
non-neighboring finite elements (center column) and finally failure to enforce
contact constraints (right column). Note the overlapping colored elements of
the third medium in the detailed insets.

A correctly behaving C-shape benchmark simulation, i.e., with the additional
layer of elements, is depicted in Figure 13, where the newly proposed approach
to regularization of the third medium (denoted herein as $\nabla\ln\bm{Q}$,
see Section 2.3) is compared to the original approach of Bluhm et al. [30]
(denoted here as $\nabla\bm{F}$).

The parameters of the third medium for the $\nabla\ln\bm{Q}$ regularization
are $\gamma=1,k_{r}=2\times 10^{3}$, whereas for the $\nabla\bm{F}$
regularization the regularization stiffness $k_{r}=2\times 10^{3}$ and
$\gamma\in\\{0.1,1\\}$ is used. Since the contact stiffness $\gamma=1$ leads
to premature contact, we include also a lower setting of $\gamma=0.1$.

The newly proposed $\nabla\ln\bm{Q}$ regularization proves to be more robust,
allowing for much larger deformation (see Figure 13c) before the simulation
stops due to numerical oscillations and instabilities caused by large sliding
along the contact surface (as the surface has rotated to be loaded almost
entirely in a tangential direction). This is mainly due to the ability of the
new formulation to utilize a larger value of contact stiffness without
experiencing premature contact, thus achieving better performance, with the
computational cost of both formulations being comparable.

Figure 13: Self contact example and a comparison of the newly proposed
regularization $\nabla\ln\bm{Q}$ to a $\nabla\bm{F}$ formulation of Bluhm et
al. [30]. In the pre-contact phase, $\nabla\bm{F}$ with $\gamma=1$ exhibits
premature contact (in blue), unlike $\nabla\ln\bm{Q}$ (in red) or
$\nabla\bm{F}$ with $\gamma=0.1$ (in cyan). The final configuration plot and
its associated force-displacement curves with highlighted points of failure
demonstrate the higher robustness of the $\nabla\ln\bm{Q}$ regularization. The
$\nabla\ln\bm{Q}$ formulation only becomes numerically unstable when the
contact surface starts being loaded in a tangential direction.

###  3.3 Buckling of an internally pressurized metamaterial

The last example showcases the third medium model in simulations of
pneumatically actuated pattern-forming metamaterials. Specifically, we
consider a silicone rubber sample with four circular voids arranged in a
square lattice. This arrangement is known upon compression to change its
geometry from initially circular holes to elliptical holes with alternating
horizontal and vertical major axes as a result of an internal instability,
resulting in auxetic behavior [11, 12]. Alternatively to mechanical
compression, this behavior can also be introduced by pneumatic pressure in the
voids [18].

Here, the numerical model is compared against experimental data. To this end,
a square specimen has been manufactured and tested with the dimensions of
$B\times H=$40\text{\,}\mathrm{mm}$\times$40\text{\,}\mathrm{mm}$$ and four
voids of $d=$15\text{\,}\mathrm{mm}$$ diameter located
$t=$3.25\text{\,}\mathrm{mm}$$ from the edge, see Figure 14.

Figure 14: Geometry and boundary conditions of the symmetrical metamaterial
rubber sample.

In the experiment, the sample was placed between a pair of lubricated
plexiglass sheets, see Figure 15a. Through four drilled holes therein, air
pumps were attached to introduce a constant suction pressure into all four
voids. The evolution of pressure difference with respect to the ambient
atmospheric pressure was measured during the experiment to determine a
critical value of the pressure difference leading to internal buckling. In a
repeated pressure loading and unloading procedure, the sample was observed to
buckle upon loading and then return to an almost undeformed configuration upon
unloading. Numerical sensitivity analysis with regards to geometric
imperfections shows that this effect of material memory should have a
negligible influence on the results. The buckling point was consistently
achieved for the critical pressure of $(\Delta
p)_{\mathrm{crit}}=$-8.45\text{\,}\mathrm{kPa}$\pm$0.3\text{\,}\mathrm{kPa}$$.

To reflect the conditions of the experiment, a square-shaped finite structure
with free edges except for preventing rigid body motion is considered in the
numerical simulation. The considered plane strain model corresponds well to
reality, since the out-of-plane deformation is prevented by the fixed pair of
plexiglass sheets. For the silicone rubber, yet again the neo-Hookean model is
used with its strain energy density as defined in Equation (20). The material
parameters $K$ and $G$ have been determined from an assumption of near
incompressibility and uniaxial compression tests on samples of the silicone
rubber, resulting in the values of $K=$91.111\text{\,}\mathrm{MPa}$$ and
$G=$0.182\text{\,}\mathrm{MPa}$$. Further details of the material parameter
identification are included in B.

The third medium parameters in this case are chosen as $\gamma=1\times
10^{-3}$ and $k_{r}=2\times 10^{3}$. The $\nabla\ln\bm{Q}$ model with
pneumatic actuation is used. The contact stiffness is somewhat lower than the
previous examples to maintain comparative compliance of the third medium in
relation to the parameters of the solid material.

In both the simulation and experiment, the negative pressure difference leads
to an internal instability causing a pattern to emerge in the shape of the
voids as a result of bifurcation. Photos of the experiment are compared to
deformed configurations resulting from the simulation in Figure 15a for three
values of pressure difference: one pre-bifurcation, one post-bifurcation, and
one at the end of the experiment where voids are in contact. The critical
pressure value of the instability found experimentally is highlighted in a
pressure-volume diagram in Figure 15b. In the simulation, the critical
pressure value obtained was slightly lower, namely $(\Delta
p)_{\mathrm{crit}}=$-8.96\text{\,}\mathrm{kPa}$$. This slight error is
explained by measurement imprecision, viscous effects in the silicone rubber
material not captured by the hyperelastic model, general geometrical
imperfection of the experiment, and friction between the plexiglass sheets and
the specimen. Nevertheless, it can be said that the numerical results agree
satisfactorily with the experimental observations, with the model accurately
capturing the experimentally determined behavior.

---
(a)
(b)
Figure 15: Simulation and experimental data of a pneumatically actuated sample
of a pattern-forming metamaterial. (a) Comparison of the deformed states from
simulation and experiment for various levels of internal pressure loading. (b)
A corresponding pressure-volume diagram, with an apparent point of internal
instability. The horizontal line highlights the experimentally determined
buckling pressure of $(\Delta
p)_{\mathrm{crit}}=$-8.45\text{\,}\mathrm{kPa}$$.

##  4 Conclusions

In this manuscript, we have introduced a novel material model for a
computational fictitious third medium enabling concurrent modeling of
pneumatic actuation and contact. Combined pneumatic actuation and contact can
be a strong design tool in research, topology optimization, and design within
the soft robotics field and mechanics of secondary stiffening or densification
of pattern-forming metamaterials. The proposed model is formulated as
hyperelastic, albeit with inclusion of second gradient terms. Its strain
energy density function comprises three distinct parts representing the
different functionality of the third medium (contact, regularization of pre-
contact behavior, pneumatic pressurization). Individual parts can be switched
on and off independently depending on the use case. Noteworthy is the utility
of this formulation for topology optimization of both metamaterials and
macrostructures. In that case, the presence of a third medium mesh is
implicit, which eliminates the main disadvantage of all third medium methods:
the need for extra computational effort due to additional finite elements
within void regions.

The main contributions of this manuscript can be summarized as follows.
Several limitations of existing third medium formulations for contact and
pneumatic actuation have been eliminated. First, our pneumatic actuation model
describes pressure loads precisely, maintaining consistency with a follower-
load approach. Next, a second gradient regularization term, which governs the
behavior of the compliant third medium, significantly improves computational
stability. Our modification is based on an energy potential, meaning that the
deformation of the third medium is thus fully reversible without an energy
dissipation and the entire model has a symmetric tangent stiffness. This
further allows for the employment of advanced mathematical programming methods
and finite element solvers, such as buckling analysis. Finally, only material
gradients of selected deformation mode fields are penalized, greatly improving
the compliant behavior of the third medium.

The proposed methodology has been demonstrated on a number of numerical
examples, demonstrating the usefulness of the model for both free surface and
closed void configurations including simulations of pneumatically actuated
structures. In all cases, our formulation exhibits in improved behavior
compared to existing third medium contact formulations. It is furthermore
shown that presented contact method fulfills the standardized benchmark patch
test.

Although the model formulation has been presented for plane strain geometries,
its extension to 3D is straightforward, albeit technical. The incremental
additive method for the calculation of rotational gradient terms relies on
commutation between rotation tensors, generally not present in 3D geometries;
the gradient terms can, however, be reformulated in a more general exponential
form to avoid this limitation. Desired future extensions concern
implementation within a topology optimization framework and inclusion of
frictional contact and, which is an ongoing effort [25].

#### Acknowledgement

This work has been supported by the Czech Science Foundation [grant no.
GA19-26143X]. The authors would like to thank Jan Havelka and Pavel Rychnovský
at the Czech Technical University in Prague for their help with physical
experiments.

##  Appendix A Implementation of second gradient regularization

Third medium is prone to excessive deformation of its highly compliant finite
elements, which can lead to numerical instabilities. To prevent this issue, it
is necessary to penalize sharp changes in rotations, i.e., to penalize
curvatures. To this effect, a regularization term is introduced, penalizing
the gradient of large rotations. The rotation tensor $\bm{R}$ can be obtained
by polar decomposition of the deformation gradient $\bm{F}$ into stretch and
rotations as

$\bm{F}=\bm{R}\cdot\bm{U}$ (21)

However, calculation of $\nabla\bm{R}$ is computationally challenging,
requiring spectral decomposition of $\bm{F}$, which is process that is
difficult to differentiate. Since the desired penalization is not physical but
purely numerical, we can opt for penalizing related quantities instead. Thus
we propose to penalize the gradient of the logarithm of an alternative
rotational measure $\bm{Q}$, i.e. $\nabla\ln\bm{Q}$, which can be in the given
plane strain finite element context, calculated more easily with the use of
incremental linearization.

The $\bm{Q}$-rotations can be incrementally calculated in pseudo-time
derivative from the material spin tensor $\bm{w}_{X}$ as

$\dot{\bm{Q}}=\bm{w}_{X}\cdot\bm{Q}$ (22)

The material spin tensor $\bm{w}_{X}$ is defined as the skew-symmetric part of
the _material_ gradient of the spatial velocity $\nabla\bm{v}$:

$\bm{w}_{X}=\frac{1}{2}\left(\nabla\bm{v}-\left(\nabla\bm{v}\right)^{T}\right)$
(23)

The spatial velocity $\bm{v}$ is defined as the change of displacement
$\bm{u}$ within a computational timestep $\Delta t$:

$\bm{v}=\frac{\Delta\bm{u}}{\Delta t}$ (24)

The use of the material gradient in Equation (23) distinquishes the material
spin from the commonly known spin $\bm{w}$, which is defined using the
_spatial_ gradient instead [45], and is linked to the rotational tensor
$\bm{R}$ through the relation [46]

$\dot{\bm{R}}=\bm{w}\cdot\bm{R}-\bm{R}\cdot{\rm
Skw}\left(\dot{\bm{U}}\cdot\bm{U}^{-1}\right)$ (25)

where ${\rm Skw}\left(\bm{A}\right)$ denotes the skew symmetric part of a
tensor $\bm{A}$. A comparison between Equations (25) and (22) reveals the
simplifications made in the use of $\bm{Q}$ instead of $\bm{R}$, that is
namely the exchange of spatial spin for material spin and neglect of the
correction term removing the rotated skew symmetric parts of stretch from the
result. While these simplifications introduce errors in physical description
of curvature, we do not find them to hinder the purpose of the regularization
term, i.e., stabilization of behavior of the fictional third medium.

A linearization of Equation (22) with the use of exponential mapping [45]
leads to

$\bm{Q}(t_{n+1})={\rm exp}(\Delta t\bm{w}_{X})\cdot\bm{Q}(t_{n})$ (26)

where $t_{n}$ and $t_{n+1}$ are arbitrarily chosen pseudo-times of the finite
element solver and $\Delta t=t_{n+1}-t_{n}$ is their difference.

A logarithm can now be applied to equation (26) and since in the two-
dimensional case, the tensors $\exp(\bm{w}_{X})$ and $\bm{Q}$ commute (for 3D
see A.1), the rule on sum of logarithms can be applied, yielding

$\ln\bm{Q}(t_{n+1})=\Delta t\bm{w}_{X}+\ln\bm{Q}(t_{n})$ (27)

Finally, applying the gradient operator to both sides of the equation, a
simple additive incremental formulation for $\nabla\ln\bm{Q}$ is obtained:

$\nabla\ln\bm{Q}(t_{n+1})=\Delta t\nabla\bm{w}_{X}+\nabla\ln\bm{Q}(t_{n})$
(28)

The material spin tensor gradient $\nabla\bm{w}_{X}$ can be calculated at a
Gauss point of a standard second-order finite element from nodal displacements
and shape functions derivatives, as detailed below. The incremental nature of
the formula requires a basic implementation of material status memory, which
is otherwise uncommon for hyper-elasticity, but common, e.g., for materials
with plastic behavior. As an initial condition, the value
$\nabla\ln\bm{Q}(t_{0}=0)=\bm{0}$ is assumed for the undeformed state of the
third medium.

The definition of spatial velocity (recall Equation (24)) can be used,
introducing the incremental displacement gradient
$\bm{f}=\nabla(\Delta\bm{u})$, to move the time step length to the left side
of Equation (23):

$\Delta t\bm{w}_{X}=\frac{1}{2}\left(\bm{f}-\bm{f}^{T}\right)$ (29)

To obtain the time step increment to the gradient $\bm{Q}$-rotation logarithm
term (recall Equation (28)), only an application of the gradient operator
remains:

$\Delta
t\nabla\bm{w}_{X}=\frac{1}{2}\left(\nabla\bm{f}-\left(\nabla\bm{f}^{T}\right)\right)$
(30)

Upon discretization of the domain into finite elements, it would be
advantageous to express this relation in dependence on the nodal displacements
vector $\mathsf{d}$, its time increment $\Delta\mathsf{d}$ and their variation
$\delta\mathsf{d}$ (note that $\delta\Delta\mathsf{d}=\delta\mathsf{d}$
because $\mathsf{d}(t_{n+1})=\mathsf{d}(t_{n})+\Delta\mathsf{d}$). So it
follows for the incremental term and its variation that

$\displaystyle\Delta t\nabla\bm{w}_{X}$ $\displaystyle=$
$\displaystyle\mathsf{H}_{w}\cdot\Delta\mathsf{d}$ (31)
$\displaystyle\delta\left(\Delta t\nabla\bm{w}_{X}\right)$ $\displaystyle=$
$\displaystyle\mathsf{H}_{w}\cdot\delta\mathsf{d}$ (32)

where $\mathsf{H}_{w}$ is an appropriate linear mapping of second derivatives
of element shape functions to nodal displacements. In an engineering matrix
notation, Equation (31) takes on the form

$\footnotesize\Delta t\begin{bmatrix}w^{X}_{11,1}\\\ w^{X}_{11,2}\\\
w^{X}_{12,1}\\\ w^{X}_{12,2}\\\ w^{X}_{21,1}\\\ w^{X}_{21,2}\\\
w^{X}_{22,1}\\\
w^{X}_{22,2}\end{bmatrix}=\frac{1}{2}\begin{bmatrix}0&0&0&0&&...&0&0\\\
0&0&0&0&&...&0&0\\\
N_{1,21}&-N_{1,11}&N_{2,21}&-N_{2,11}&&...&N_{i,21}&-N_{i,11}\\\
N_{1,22}&-N_{1,12}&N_{2,22}&-N_{2,12}&&...&N_{i,22}&-N_{i,12}\\\
-N_{1,21}&N_{1,11}&-N_{2,21}&N_{2,11}&&...&-N_{i,21}&N_{i,11}\\\
-N_{1,22}&N_{1,12}&-N_{2,22}&N_{2,12}&&...&-N_{i,22}&N_{i,12}\\\
0&0&0&0&&...&0&0\\\ 0&0&0&0&&...&0&0\end{bmatrix}\begin{bmatrix}\Delta
d_{1}\\\ \Delta d_{2}\\\ \Delta d_{3}\\\ \Delta d_{4}\\\ ...\\\ \Delta
d_{2i-1}\\\ \Delta d_{2i}\end{bmatrix}\normalsize$ (33)

where $X_{1}$ and $X_{2}$ are components of the material coordinate tensor
$\bm{X}$, the subscript $\bullet_{,k}$ denotes a derivative with regards to
$X_{k}$, $N_{1}$ to $N_{i}$ are the shape functions of an $i$-node element
dependent on the material position $\bm{X}$, and $\Delta d_{1}$ to $\Delta
d_{2i}$ are the incremental values of displacements in the degrees of freedom
corresponding to those nodes.

The relevant part of the strain energy density function, as discussed in
Section 2.3, has a quadratic form:

$W_{\mathrm{r,\ln(Q)}}=\frac{1}{2}c\;\nabla\ln\bm{Q}\mathbin{\vdots}\nabla\ln\bm{Q}$
(34)

with $c$ being a material parameter. Its first and second variations read

$\displaystyle\delta W_{\mathrm{r,Q}}$ $\displaystyle=$ $\displaystyle
c\left(\nabla\ln\bm{Q}\right)\mathbin{\vdots}\delta\left(\nabla\ln\bm{Q}\right)$
(35) $\displaystyle\delta^{2}W_{\mathrm{r,Q}}$ $\displaystyle=$ $\displaystyle
c\;\delta\left(\nabla\ln\bm{Q}\right)\mathbin{\vdots}\delta\left(\nabla\ln\bm{Q}\right)$
(36)

Due to the additive incremental nature of Equation (28) the variation
$\delta\left(\nabla\ln\bm{Q}\right)$ only depends on the variation of the
material spin gradient $\nabla\bm{w}_{X}$ and does not depend on any material
status history, i.e.,

$\displaystyle\nabla\ln\bm{Q}(t_{n+1})$ $\displaystyle=$
$\displaystyle\nabla\ln\bm{Q}(t_{n})+\mathsf{H}_{w}\cdot\Delta\mathsf{d}$ (37)
$\displaystyle\delta\left(\nabla\ln\bm{Q}\right)$ $\displaystyle=$
$\displaystyle\mathsf{H}_{w}\cdot\delta\mathsf{d}$ (38)

###  A.1 Alternative energy formulation with direct exponential mapping

The introduction of logarithm to the $\bm{Q}$-rotation term is not strictly
necessary; an alternative approach is possible using directly the exponential
mapping defined in Equation (26). This would be especially important in a 3D
implementation, where the rule on sum of logarithms (consider Equation 27)
cannot be used.

Upon application of the gradient operator, the expression

$\nabla\bm{Q}(t_{n+1})=\nabla\left(\exp(\Delta
t\bm{w}_{X})\right)\cdot\bm{Q}(t_{n})+\exp(\Delta
t\bm{w}_{X})\cdot\nabla\bm{Q}(t_{n})$ (39)

is obtained, with $t_{n}$ denoting the $n$-th time step.

Values of $\bm{Q}(t_{n})$ and $\nabla\bm{Q}(t_{n})$ are always known from the
previous time step, seeded as $\bm{I}$ and $\bm{0}$, respectively, at
$t_{0}=0$. For the calculation of the exponential term and its gradient, the
method proposed by Ortiz et al. [47] can be applied. The exponential is
expressed as a sum of an infinite series

$\displaystyle\exp(\Delta t\bm{w}_{X})$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\frac{1}{k!}(\Delta t\bm{w}_{X})^{k}$ (40)

The $(k+1)$-th member of this series and its gradient can be calculated from
the $k$-th member and its gradient:

$\displaystyle\exp(\Delta t\bm{w}_{X})^{(k+1)}$ $\displaystyle=$
$\displaystyle\frac{1}{k+1}\exp(\Delta t\bm{w}_{X})^{(k)}\cdot\left(\Delta
t\bm{w}_{X}\right)$ (41) $\displaystyle\nabla\exp(\Delta t\bm{w}_{X})^{(k+1)}$
$\displaystyle=$ $\displaystyle\frac{1}{k+1}\Big{(}\nabla\exp(\Delta
t\bm{w}_{X})^{(k)}\cdot\left(\Delta t\bm{w}_{X}\right)$ $\displaystyle+$
$\displaystyle\exp(\Delta t\bm{w}_{X})^{(k)}\cdot\left(\Delta
t\nabla\bm{w}_{X}\right)\Big{)}$

considering the zeroth member as

$\exp(\Delta t\bm{w}_{X})^{(0)}=\bm{I}\quad\quad\quad\nabla\exp(\Delta
t\bm{w}_{X})^{(0)}=\bm{0}$ (43)

This allows construction of the strain energy density function without the
logarithm in the form

$W_{\mathrm{r,Q}}=\frac{1}{2}c\;\nabla\bm{Q}\mathbin{\vdots}\nabla\bm{Q}$ (44)

The recursive formulae used in this approach present an increase both in
computational time and in implementation difficulty. Differentiation of the
recursive formulas for the calculation of energy variations, i.e., forces and
stiffness tangent, is certainly possible, but cumbersome. Alternatively, a
finite difference scheme may be applied on the strain energy function at a
further cost to computational time, as was done in our testing. In our
experience, the effect of this more precise approach to penalization of
curvature is negligible and does not compensate for this increase in
complexity of the material model. It is thus only recommended to use this
scheme in 3D context, where it is necessary.

###  A.2 Penalization of Jacobian gradients

Due to instabilities and/or stiffening arising in certain geometries of third
medium elements penalized by curvature only, as discussed in Section 2.3,
enforcement of uniform volume change across an element is required. For this
reason, another term is introduced into the energy density function,
penalizing gradients of deformation gradient Jacobian $J$:

$W_{\mathrm{r,J}}=\frac{1}{2}c\;\left(\nabla J\cdot\nabla J\right)$ (45)

The Jacobian is defined as the determinant of the deformation gradient:

$\displaystyle J=\det{\bm{F}}$ (46)

and thus its gradient can be expressed as

$\nabla J=\frac{\partial\det{\bm{F}}}{\partial\bm{F}}:\nabla\bm{F}={\rm
cof}\bm{F}:\nabla\bm{F}$ (47)

with ${\rm cof}\bm{F}=\det(\bm{F})\bm{F}^{-T}$ signifying the tensor cofactor
operation. An efficient use can be made of an alternative definition of the
tensor cofactor [48], leading to

$\displaystyle J$ $\displaystyle=$
$\displaystyle\frac{1}{6}\bm{F}\times\bm{F}:\bm{F}$ (48) $\displaystyle\nabla
J$ $\displaystyle=$ $\displaystyle\frac{1}{2}\bm{F}\times\bm{F}:\nabla\bm{F}$
(49)

where $\times$ denotes the tensor cross product operator. In the indicial
notation the expressions amount to

$\displaystyle J$ $\displaystyle=$
$\displaystyle\frac{1}{6}\epsilon_{ikp}\epsilon_{jlq}F_{kl}F_{pq}F_{ij}$ (50)
$\displaystyle J_{,m}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{ikp}\epsilon_{jlq}F_{kl}F_{pq}F_{ij,m}$
(51)

where $\epsilon_{ijk}$ is the Levi-Civita symbol and the subscript
$\bullet_{,m}$ denotes the $m$-th coordinate of the gradient, i.e., a
derivative with regards to the reference coordinate $X_{m}$.

Computing the necessary variations of the strain energy density function
$W_{\mathrm{r,J}}$, it is now possible to derive the contributions to the
first Piola-Kirchhoff stress $\bm{P}_{\mathrm{r,J}}$, second order stress
$\bm{T}_{\mathrm{r,J}}$, elastic stiffness tangent
$\mathbf{D}_{\mathrm{r,J}}$, mixed order stiffnesses
$\mathbf{C}_{\mathrm{PS,r,J}}$ and $\mathbf{C}_{\mathrm{SP,r,J}}$, and second
order stiffness $\mathbf{C}_{\mathrm{SS,r,J}}$. In indicial notation, these
can be written as

$\displaystyle W_{\mathrm{r,J}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}cJ_{,m}J_{,m}$ (52) $\displaystyle
P^{\mathrm{r,J}}_{rs}$ $\displaystyle=$ $\displaystyle cJ_{,m}\frac{\partial
J_{,m}}{\partial F_{rs}}$ (53) $\displaystyle T^{\mathrm{r,J}}_{rst}$
$\displaystyle=$ $\displaystyle cJ_{,m}\frac{\partial J_{,m}}{\partial
F_{rs,t}}$ (54) $\displaystyle D^{\mathrm{r,J}}_{rsuv}$ $\displaystyle=$
$\displaystyle c\frac{\partial J_{,m}}{\partial F_{uv}}\frac{\partial
J_{,m}}{\partial F_{rs}}+cJ_{,m}\frac{\partial^{2}J_{,m}}{\partial
F_{rs}\partial F_{uv}}$ (55) $\displaystyle C^{\mathrm{PS,r,J}}_{rsuvw}$
$\displaystyle=$ $\displaystyle c\frac{\partial J_{,m}}{\partial
F_{uv,w}}\frac{\partial J_{,m}}{\partial
F_{rs}}+cJ_{,m}\frac{\partial^{2}J_{,m}}{\partial F_{rs}\partial F_{uv,w}}$
(56) $\displaystyle C^{\mathrm{SP,r,J}}_{rstuv}$ $\displaystyle=$
$\displaystyle c\frac{\partial J_{,m}}{\partial F_{uv}}\frac{\partial
J_{,m}}{\partial F_{rs,t}}+cJ_{,m}\frac{\partial^{2}J_{,m}}{\partial
F_{rs,t}\partial F_{uv}}$ (57) $\displaystyle C^{\mathrm{SS,r,J}}_{rstuvw}$
$\displaystyle=$ $\displaystyle c\frac{\partial J_{,m}}{\partial
F_{uv,w}}\frac{\partial J_{,m}}{\partial
F_{rs,t}}+cJ_{,m}\frac{\partial^{2}J_{,m}}{\partial F_{rs,t}\partial
F_{uv,w}}$ (58)

where

$\displaystyle\frac{\partial J_{,m}}{\partial F_{rs}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{irp}\epsilon_{jsq}F_{pq}F_{ij,m}+\frac{1}{2}\epsilon_{ikr}\epsilon_{jls}F_{kl}F_{ij,m}=\epsilon_{irp}\epsilon_{jsq}F_{pq}F_{ij,m}$
(59) $\displaystyle\frac{\partial J_{,m}}{\partial F_{rs,t}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{rkp}\epsilon_{slq}F_{kl}F_{pq}\delta_{mt}$
(60) $\displaystyle\frac{\partial^{2}J_{,m}}{\partial F_{rs}\partial F_{uv}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{iru}\epsilon_{jsv}F_{ij,m}+\frac{1}{2}\epsilon_{iur}\epsilon_{jvs}F_{ij,m}=\epsilon_{iru}\epsilon_{jsv}F_{ij,m}$
(61) $\displaystyle\frac{\partial^{2}J_{,m}}{\partial F_{rs}\partial
F_{uv,w}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{urp}\epsilon_{vsq}F_{pq}\delta_{mw}+\frac{1}{2}\epsilon_{ukr}\epsilon_{vls}F_{kl}\delta_{mw}=\epsilon_{urp}\epsilon_{vsq}F_{pq}\delta_{mw}$
(62) $\displaystyle\frac{\partial^{2}J_{,m}}{\partial F_{rs,t}\partial
F_{uv}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{rup}\epsilon_{svq}F_{pq}\delta_{mt}+\frac{1}{2}\epsilon_{rku}\epsilon_{slv}F_{kl}\delta_{mt}=\epsilon_{rup}\epsilon_{svq}F_{pq}\delta_{mt}$
(63) $\displaystyle\frac{\partial^{2}J_{,m}}{\partial F_{rs,t}\partial
F_{uv,w}}$ $\displaystyle=$ $\displaystyle 0$ (64)

with $\delta_{ij}$ as the Kronecker delta symbol.

With these first and second order stresses and stiffnesses, the first and
second variation of $W_{\mathrm{reg,J}}$ can be constructed as

$\displaystyle\delta W_{\mathrm{r,J}}$ $\displaystyle=$
$\displaystyle\bm{P}_{\mathrm{r,J}}:\delta\bm{F}+\bm{T}_{\mathrm{r,J}}:\delta\nabla\bm{F}$
(65) $\displaystyle\delta^{2}W_{\mathrm{r,J}}$ $\displaystyle=$
$\displaystyle\delta\bm{F}:\mathbf{D}_{\mathrm{r,J}}:\delta\bm{F}+\delta\nabla\bm{F}:\mathbf{C}_{\mathrm{PT,r,J}}:\delta\bm{F}$
$\displaystyle+$
$\displaystyle\delta\bm{F}:\mathbf{C}_{\mathrm{TP,r,J}}:\delta\nabla\bm{F}+\delta\nabla\bm{F}:\mathbf{C}_{\mathrm{TT,r,J}}:\delta\nabla\bm{F}$

leading after finite element discretization and transfer to engineering
notation to

$\displaystyle\delta W_{\mathrm{r,J}}$ $\displaystyle=$
$\displaystyle\mathsf{P}_{\mathrm{r,J}}\mathsf{B}\delta\mathsf{d}+\mathsf{T}_{\mathrm{r,J}}\;\mathsf{H}\delta\mathsf{d}$
(67) $\displaystyle\delta^{2}W_{\mathrm{r,J}}$ $\displaystyle=$
$\displaystyle\delta\mathsf{d}^{\mathrm{T}}\mathsf{B}^{\mathrm{T}}\;\mathsf{D}_{\mathrm{r,J}}\;\mathsf{B}\delta\mathsf{d}+\delta\mathsf{d}^{\mathrm{T}}\mathsf{H}^{\mathrm{T}}\;\mathsf{C}_{\mathrm{PT,r,J}}\;\mathsf{B}\delta\mathsf{d}$
$\displaystyle+$
$\displaystyle\delta\mathsf{d}^{\mathrm{T}}\mathsf{B}^{\mathrm{T}}\;\mathsf{C}_{\mathrm{TP,r,J}}\;\mathsf{H}\delta\mathsf{d}+\delta\mathsf{d}^{\mathrm{T}}\mathsf{H}^{\mathrm{T}}\;\mathsf{C}_{\mathrm{TT,r,J}}\;\mathsf{H}\delta\mathsf{d}$

with $\mathsf{B}$ a linear mapping of element shape function derivatives to
nodal displacements and $\mathsf{H}$ likewise a linear mapping of second
derivatives of element shape functions to nodal displacements.

This approach provides a universal derivation valid in all cases. It is,
however, notable that due to the tensor cross product being only defined in
three dimensions, a two-dimensional deformation gradient tensors have to be
brought to the three-dimensional space (and later back) to allow for use of
the formulas. Particularly in a plane strain case, a large simplification is
possible, if a direct definition for a determinant of a second order tensor is
used:

$\displaystyle J=F_{11}F_{22}-F_{12}F_{21}$ (69)

Now it is possible to simply differentiate for the gradient:

$\displaystyle
J_{,m}=F_{11,m}F_{22}+F_{11}F_{22,m}-F_{12,m}F_{21}-F_{12}F_{21,m}$ (70)

Equations (59)-(64) thus greatly simplify to:

$\displaystyle\frac{\partial J_{,m}}{\partial F_{rs}}$ $\displaystyle=$
$\displaystyle
F_{11,m}\delta_{2r}\delta_{2s}+\delta_{1r}\delta_{1s}F_{22,m}-F_{12,m}\delta_{2r}\delta_{1s}-\delta_{1r}\delta_{2s}F_{21,m}$
(71) $\displaystyle\frac{\partial J_{,m}}{\partial F_{rs,t}}$ $\displaystyle=$
$\displaystyle\delta_{1r}\delta_{1s}\delta_{mt}F_{22}+F_{11}\delta_{2r}\delta_{2s}\delta_{mt}$
$\displaystyle-$
$\displaystyle\delta_{1r}\delta_{2s}\delta_{mt}F_{21}-F_{12}\delta_{2r}\delta_{1s}\delta_{mt}$
$\displaystyle\frac{\partial^{2}J_{,m}}{\partial F_{rs}\partial F_{uv}}$
$\displaystyle=$ $\displaystyle 0$ (73)
$\displaystyle\frac{\partial^{2}J_{,m}}{\partial F_{rs}\partial F_{uv,w}}$
$\displaystyle=$
$\displaystyle\delta_{1u}\delta_{1v}\delta_{mw}\delta_{2r}\delta_{2s}+\delta_{1r}\delta_{1s}\delta_{2u}\delta_{2v}\delta_{mw}$
$\displaystyle-$
$\displaystyle\delta_{1u}\delta_{2v}\delta_{mw}\delta_{2r}\delta_{1s}-\delta_{1r}\delta_{2s}\delta_{2u}\delta_{1v}\delta_{mw}$
$\displaystyle\frac{\partial^{2}J_{,m}}{\partial F_{rs,t}\partial F_{uv}}$
$\displaystyle=$
$\displaystyle\delta_{1r}\delta_{1s}\delta_{mt}\delta_{2u}\delta_{2v}+\delta_{1u}\delta_{1v}\delta_{2r}\delta_{2s}\delta_{mt}$
$\displaystyle-$
$\displaystyle\delta_{1r}\delta_{2s}\delta_{mt}\delta_{2u}\delta_{1v}-\delta_{1u}\delta_{2v}\delta_{2r}\delta_{1s}\delta_{mt}$
$\displaystyle\frac{\partial^{2}J_{,m}}{\partial F_{rs,t}\partial F_{uv,w}}$
$\displaystyle=$ $\displaystyle 0$ (76)

Note that in this plane strain case, all the tangent contributions are
generally independent of the in-plane members of the deformation gradient,
making them constant.

##  Appendix B Identification of material parameters for the silicone rubber
used in physical experiment

In Section 3.3, we present a verification of the proposed third medium
material model against a physical experiment of pneumatic actuation of a
pattern-transforming metamaterial. The experimental sample was cast from
silicone rubber [49] in a 3D-printed mold. To match simulation and
experimental data, material properties need to be found for this silicone
rubber material.

The chosen material model is the neo-Hookean hyperelastic model with the
following strain energy density

$W_{\mathrm{NH,bulk}}=\frac{K}{2}\ln^{2}{J}+\frac{G}{2}\left(J^{-2/3}I_{1}-3\right)$
(77)

with $K$ and $G$ being the material parameters that need to be determined.
They are a representation of (and when linearized around the reference
configuration they exactly correspond to) the bulk and shear moduli,
respectively, of linear small-strain elasticity. Consequently, they can be
equivalently expressed in terms of Young’s modulus $E$ and the initial
Poisson’s ratio $\nu$:

$K=\frac{E}{3(1-2\nu)}\quad\quad G=\frac{E}{2(1+\nu)}$ (78)

Polymer-based elastic materials such as silicone rubbers are generally assumed
to be incompressible with a Poisson’s ratio approaching $\nu=0.5$ [50].
Introducing an assumption of near incompressibility into our model and setting
$\nu=0.499$, we are left with two parameters which both are just a linear
function of an unknown Young’s modulus $E$. No specific material treatment is
required for this value of $\nu$.

To determine this unknown value, uniaxial compression tests were conducted on
samples of the silicone rubber material. Three separate cyllindrical samples
of length $L=$25\text{\,}\mathrm{mm}$$ and diameter
$d=$13\text{\,}\mathrm{mm}$$ were compressed in the axial direction in an
open-hardware Thymos111http://thymos.cz/ loading frame. The experimental setup
with different stages of loading can be seen in Figure 1. The loading
continued up to the point at which each individual sample slid out of its
fixings. The force-displacement curve was measured at displacement increments
of $\Delta u=$0.05\text{\,}\mathrm{mm}$$. Converting the resulting data to
axial stretch $\lambda_{1}=(L-u)/L$ and the normal component of the first
Piola-Kirchhoff stress $P_{11}=N/(\pi d^{2}/4)$ where $N$ is the measured
normal force response, gives three stress-strain curves shown in Figure 2
(dotted lines).

|  |  |
---|---|---|---
Figure 1: Uniaxial compression of silicone rubber cylinders, various stages of
loading (compressive strain increases left to right). The samples are fixed in
the loading machine between lubricated plexiglass plates.

The conditions of the uniaxial compression test can be described by the neo-
Hookean model. Starting from the strain energy density in Equation (77) and
differentiating with respect to the deformation gradient, the first Piola-
Kirchhoff stress is obtained in the indicial notation as

$P_{ij}=\frac{\partial W_{\mathrm{NH,bulk}}}{\partial
F_{ij}}=K\ln{J}F^{-1}_{ji}-\frac{G}{3}J^{-2/3}F^{-1}_{ji}\;I_{1}+GJ^{-2/3}F_{ij}$
(79)

Assuming laterally unconstrained uniaxial compression, the matrix
representation of the deformation gradient takes the form of

$\mathsf{F}=\begin{bmatrix}\lambda_{1}&0&0\\\ 0&\lambda_{2}&0\\\
0&0&\lambda_{3}\end{bmatrix}$ (80)

where $\lambda_{i}$ are the principal stretches with $\lambda_{1}$ oriented
along the axis of the cylinder. The invariants $J$ and $I_{1}$ then read

$J=\lambda_{1}\lambda_{2}\lambda_{3}\quad\quad
I_{1}=\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}$ (81)

Substituting into Equation (79), normal components of the first Piola-
Kirchhoff stress in the direction of $\lambda_{i}$ can be expressed as

$\displaystyle P_{ii}$ $\displaystyle=$
$\displaystyle\frac{K}{\lambda_{i}}\left(\ln{\lambda_{1}}+\ln{\lambda_{2}}+\ln{\lambda_{3}}\right)$
$\displaystyle-$
$\displaystyle\frac{G}{3}\left(\lambda_{1}\lambda_{2}\lambda_{3}\right)^{-2/3}\left(\frac{1}{\lambda_{i}}\left(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}\right)-3\lambda_{i}\right)$

Thus it is possible to compute for a given axial stretch $\lambda_{1}$ the
corresponding lateral stretches $\hat{\lambda}=\lambda_{2}=\lambda_{3}$, using
the condition of zero normal stress in the lateral direction

$P_{22}=P_{33}=0=\frac{K}{\hat{\lambda}}\left(\ln{\lambda_{1}}+2\ln{\hat{\lambda}}\right)-\frac{G}{3}\left(\lambda_{1}\hat{\lambda}^{2}\right)^{-2/3}\left(\frac{1}{\hat{\lambda}}\left(\lambda_{1}^{2}+2\hat{\lambda}^{2}\right)-3\hat{\lambda}\right)$
(83)

and then substitute back into Equation (B) to obtain the normal component of
the first Piola-Kirchhoff stress $P_{11}$.

The curve computed as described above has been fitted to the experimental data
using a least square fit procedure for the value of $E$. Errors with respect
to all three experimental datasets have been weighted equally. The resulting
fit is pictured in Figure 2. The optimal Young’s modulus value found is
$E=$0.547\text{\,}\mathrm{MPa}$$, which substituting into Equations (78) gives
the material parameters of

$K=$91.111\text{\,}\mathrm{MPa}$\quad\quad G=$0.182\text{\,}\mathrm{MPa}$$
(84)

Figure 2: Stress-strain diagrams for uniaxial compression tests on cylindrical
silicone rubber samples shown in Figure 1. Experimental data from three
independent tests and a response curve of a near incompressible neo-Hookean
material model with parameters determined from a least-square fit procedure.

## References

  * [1] X. Ren, R. Das, P. Tran, T. D. Ngo, Y. M. Xie, Auxetic metamaterials and structures: A review, Smart Materials and Structures 27 (2) (2018) 023001.  doi:10.1088/1361-665X/aaa61c.
  * [2] S. A. Ramakrishna, Physics of negative refractive index materials, Reports on Progress in Physics 68 (2005) 449–521.  doi:10.1088/0034-4885/68/2/R06.
  * [3] J.-H. Lee, J. P. Singer, E. L. Thomas, Micro-/nanostructured mechanical metamaterials, Advanced Materials 24 (36) (2012) 4782–4810.  doi:10.1002/adma.201201644.
  * [4] R. Lakes, Foam structures with a negative Poisson’s ratio, Science 235 (1987) 1038–1040.  doi:10.1126/science.235.4792.1038.
  * [5] X. Yu, J. Zhou, H. Liang, Z. Jiang, L. Wu, Mechanical metamaterials associated with stiffness, rigidity and compressibility: A brief review, Progress in Materials Science 94 (2018) 114–173.  doi:10.1016/j.pmatsci.2017.12.003.
  * [6] D. Goswami, S. Liu, A. Pal, L. G. Silva, R. V. Martinez, 3D-Architected soft machines with topologically encoded motion, Advanced Functional Materials 29 (24) (2019) 1808713.  doi:10.1002/adfm.201808713.
  * [7] D. Yang, B. Mosadegh, A. Ainla, B. Lee, F. Khashai, Z. Suo, K. Bertoldi, G. M. Whitesides, Buckling of elastomeric beams enables actuation of soft machines, Advanced Materials 27 (2015) 6323–6327.  doi:10.1002/adma.201503188.
  * [8] S. D. Papka, S. Kyriakides, Biaxial crushing of honeycombs: Part I: Experiments, International Journal of Solids and Structures 36 (29) (1999) 4367–4396.  doi:10.1016/S0020-7683(98)00224-8.
  * [9] S. D. Papka, S. Kyriakides, In-plane biaxial crushing of honeycombs: Part II: Analysis, International Journal of Solids and Structures 36 (29) (1999) 4397–4423.  doi:10.1016/S0020-7683(98)00225-X.
  * [10] N. Ohno, D. Okumura, H. Noguchi, Microscopic symmetric bifurcation condition of cellular solids based on a homogenization theory of finite deformation, Journal of the Mechanics and Physics of Solids 50 (2002) 1125–1153.  doi:10.1016/S0022-5096(01)00106-5.
  * [11] T. Mullin, S. Deschanel, K. Bertoldi, M. C. Boyce, Pattern transformation triggered by deformation, Physical Review Letters 99 (8) (2007) 084301.  doi:10.1103/PhysRevLett.99.084301.
  * [12] K. Bertoldi, P. M. Reis, S. Willshaw, T. Mullin, Negative Poisson’s ratio behavior induced by an elastic instability, Advanced Materials 22 (2010) 361–366.  doi:10.1002/adma.200901956.
  * [13] B. Florijn, C. Coulais, M. V. Hecke, Programmable mechanical metamaterials, Physical Review Letters 113 (17) (2014) 175503.  doi:10.1103/PhysRevLett.113.175503.
  * [14] S. Xiao, T. Wang, T. Liu, C. Zhou, X. Jiang, J. Zhang, Active metamaterials and metadevices: A review, Journal of Physics D: Applied Physics 53 (50) (2020) 503002.  doi:10.1088/1361-6463/abaced.
  * [15] X. Su, C. Feng, Y. Zeng, H. Yu, A dual-band tunable metamaterial absorber based on pneumatic actuation mechanism, Optics Communications 459 (2020) 124885.  doi:10.1016/j.optcom.2019.124885.
  * [16] I. Khodasevych, W. Rowe, A. Mitchell, Reconfigurable fishnet metamaterial using pneumatic actuation, Progress in Electromagnetics Research B 38 (2012) 57–70.  doi:doi.org/10.2528/PIERB11102505.
  * [17] Y. Matia, G. H. Kaiser, R. F. Shepherd, A. D. Gat, N. Lazarus, K. H. Petersen, Harnessing nonuniform pressure distributions in soft robotic actuators, Advanced Intelligent Systems 2023 (2023) 2200330.  doi:10.1002/aisy.202200330.
  * [18] Y. Chen, L. Jin, Geometric role in designing pneumatically actuated pattern-transforming metamaterials, Extreme Mechanics Letters 23 (2018) 55–66.  doi:10.1016/j.eml.2018.08.001.
  * [19] M. P. Bendsøe, O. Sigmund, Topology Optimization, Springer Berlin Heidelberg, 2004.  doi:10.1007/978-3-662-05086-6.
  * [20] C. Coulais, A. Sabbadini, F. Vink, M. van Hecke, Multi-step self-guided pathways for shape-changing metamaterials, Nature 561 (2018) 512–515.  doi:10.1038/s41586-018-0541-0.
  * [21] V. B. Hammer, N. Olhoff, Topology optimization of continuum structures subjected to pressure loading, Structural and Multidisciplinary Optimization 19 (2000) 85–92.  doi:10.1007/s001580050088.
  * [22] B. Caasenbrood, A. Pogromsky, H. Nijmeijer, A computational design framework for pressure-driven soft robots through nonlinear topology optimization, in: 2020 3rd IEEE International Conference on Soft Robotics (RoboSoft), IEEE, 2020, pp. 633–638.  doi:10.1109/RoboSoft48309.2020.9116010.
  * [23] P. Wriggers, Computational Contact Mechanics, Springer Berlin Heidelberg, 2006.  doi:10.1007/978-3-540-32609-0.
  * [24] C. Dev, G. Stankiewicz, P. Steinmann, On the influence of free space in topology optimization of electro-active polymers, Structural and Multidisciplinary Optimization 66 (8 2023).  doi:10.1007/s00158-023-03634-5.
  * [25] P. Wriggers, J. Schröder, A. Schwarz, A finite element method for contact using a third medium, Computational Mechanics 52 (2013) 837–847.  doi:10.1007/s00466-013-0848-5.
  * [26] S. Pagano, P. Alart, Self-contact and fictitious domain using a difference convex approach, International Journal for Numerical Methods in Engineering 75 (2008) 29–42.  doi:10.1002/nme.2241.
  * [27] T. Bog, N. Zander, S. Kollmannsberger, E. Rank, Normal contact with high order finite elements and a fictitious contact material, Computers & Mathematics with Applications 70 (7) (2015) 1370–1390.  doi:10.1016/j.camwa.2015.04.020.
  * [28] R. Kruse, N. Nguyen-Thanh, P. Wriggers, L. D. Lorenzis, Isogeometric frictionless contact analysis with the third medium method, Computational Mechanics 62 (2018) 1009–1021.  doi:10.1007/s00466-018-1547-z.
  * [29] F. Lorez, M. Pundir, D. S. Kammer, Eulerian framework for contact between solids represented as phase fields, Computer Methods in Applied Mechanics and Engineering 418 (2024) 116497.  doi:10.1016/j.cma.2023.116497.
  * [30] G. L. Bluhm, O. Sigmund, K. Poulios, Internal contact modeling for finite strain topology optimization, Computational Mechanics 67 (4) (2021) 1099–1114.  doi:10.1007/s00466-021-01974-x.
  * [31] A. H. Frederiksen, O. Sigmund, K. Poulios, Topology optimization of self-contacting structures, Computational Mechanics (2023).  doi:10.1007/s00466-023-02396-7.
  * [32] J. Bonet, A. J. Gil, R. D. Wood, Nonlinear Solid Mechanics for Finite Element Analysis: Statics, Cambridge University Press, 2016.  doi:10.1017/CBO9781316336144.
  * [33] M. Kružík, T. Roubíček, Mathematical Methods in Continuum Mechanics of Solids, Springer International Publishing, 2019.  doi:10.1007/978-3-030-02065-1.
  * [34] R. S. Rivlin, Large elastic deformations of isotropic materials. I. Fundamental concepts, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 240 (822) (1948) 459–490.  doi:10.1098/rsta.1948.0002.
  * [35] T. J. Pence, K. Gou, On compressible versions of the incompressible neo-Hookean material, Mathematics and Mechanics of Solids 20 (2015) 157–182.  doi:10.1177/1081286514544258.
  * [36] B. Smith, F. D. Goes, T. Kim, Analytic eigensystems for isotropic distortion energies, ACM Transactions on Graphics 38 (2 2019).  doi:10.1145/3241041.
  * [37] R. Poya, R. Ortigosa, A. J. Gil, Variational schemes and mixed finite elements for large strain isotropic elasticity in principal stretches: Closed-form tangent eigensystems, convexity conditions, and stabilised elasticity, International Journal for Numerical Methods in Engineering 124 (2023) 3436–3493.  doi:10.1002/nme.7254.
  * [38] P. Fischer, M. Klassen, J. Mergheim, P. Steinmann, R. Müller, Isogeometric analysis of 2d gradient elasticity, Computational Mechanics 47 (2011) 325–334.  doi:10.1007/s00466-010-0543-8.
  * [39] M. Horák, M. Kružík, Gradient polyconvex material models and their numerical treatment, International Journal of Solids and Structures 195 (2020) 57–65.  doi:10.1016/j.ijsolstr.2020.03.006.
  * [40] C. Geuzaine, J. F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (2009) 1309–1331.  doi:10.1002/nme.2579.
  * [41] P. E. Gill, W. Murray, Newton-type methods for unconstrained and linearly constrained optimization, Mathematical Programming 7 (1974) 311–350.  doi:10.1007/BF01585529.
  * [42] R. Fletcher, T. L. Freeman, A modified Newton method for minimization, Journal of Optimization Theory and Applications 23 (1977) 357–372.  doi:10.1007/BF00933446.
  * [43] M. A. Crisfield, Re-visiting the contact patch test, International Journal for Numerical Methods in Engineering 48 (2000) 435–449.  doi:10.1002/(SICI)1097-0207(20000530)48:3<435::AID-NME891>3.0.CO;2-V.
  * [44] G. Zavarise, L. de Lorenzis, A modified node-to-segment algorithm passing the contact patch test, International Journal for Numerical Methods in Engineering 79 (2009) 379–416.  doi:10.1002/nme.2559.
  * [45] E. A. de Souza Neto, D. Perić, D. R. J. Owen, Computational Methods for Plasticity: Theory and Applications, Wiley, 2008.  doi:10.1002/9780470694626.
  * [46] M. E. Gurtin, E. Fried, L. Anand, The Mechanics and Thermodynamics of Continua, Cambridge University Press, 2010.  doi:10.1017/CBO9780511762956.
  * [47] M. Ortiz, R. A. Radovitzky, E. A. Repetto, The computation of the exponential and logarithmic mappings and their first and second linearizations, International Journal for Numerical Methods in Engineering 52 (2001) 1431–1441.  doi:10.1002/nme.263.
  * [48] J. Bonet, A. J. Gil, R. Ortigosa, On a tensor cross product based formulation of large strain solid mechanics, International Journal of Solids and Structures 84 (2016) 49–63.  doi:10.1016/j.ijsolstr.2015.12.030.
  * [49] Havel Composites, P20 Liquid Silicone Data Sheet, revised 10 2023 (11 2019).
  * [50] P. H. Mott, C. M. Roland, Limits to Poisson’s ratio in isotropic materials, Physical Review B - Condensed Matter and Materials Physics 80 (13) (2009) 132104.  doi:10.1103/PhysRevB.80.132104.