11institutetext: Applied and Computational Mathematics Division, National
Institute of Standards and Technology (NIST), Gaithersburg, Maryland, U.S.A.,
11email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>22institutetext: Database Systems and Information
Management Group, Technische Universität, Berlin, Germany, 22email:
<EMAIL_ADDRESS>
33institutetext: Department of Mathematics and Computer Science, Virginia
State University, Petersburg, Virginia, U.S.A., 33email:
<EMAIL_ADDRESS>

# Digital Repository of Mathematical Formulae††thanks: The final publication
is available at http://link.springer.com.

Howard S. Cohl 11 Marjorie A. McClain 11 Bonita V. Saunders 11 Moritz Schubotz
22 Janelle C. Williams 33

###### Abstract

The purpose of the NIST Digital Repository of Mathematical Formulae (DRMF) is
to create a digital compendium of mathematical formulae for orthogonal
polynomials and special functions (OPSF) and of associated mathematical data.
The DRMF addresses needs of working mathematicians, physicists and engineers:
providing a platform for publication and interaction with OPSF formulae on the
web. Using MediaWiki extensions and other existing technology (such as
software and macro collections developed for the NIST Digital Library of
Mathematical Functions), the DRMF acts as an interactive web domain for OPSF
formulae. Whereas Wikipedia and other web authoring tools manifest notions or
descriptions as first class objects, the DRMF does that with mathematical
formulae. See http://gw32.iu.xsede.org/index.php/Main_Page.

##  1 Introduction

Compendia of mathematical formulae have a long and rich history. Many
scientists have developed such repositories as books and these have been
extremely useful to scientists, mathematicians and engineers over the last
several centuries (see [2, 3, 5, 6, 9, 13, 15] for instance). While there may
be some overlap of formulae in different compendia, one often needs to be
familiar with many different compendia to find a specific desired formula.
Online compendia of mathematical formulae exist, such as the NIST Digital
Library of Mathematical Functions (DLMF), subsets of Wikipedia, and the
Wolfram Functions Site. We hope to take advantage of the best aspects of these
online efforts while also incorporating powerful new features that a
community-arm of scientists should find beneficial. Our strategy is to start
with validated and trustworthy special function data from the NIST DLMF, while
adding Web 2.0 capabilities which will encourage community members to discuss
mathematical data associated with formulae. These discussions will include
internally hyperlinked proofs as well as mathematical connections between
formulae in the repository.

The online repository will be designed for a mathematically literate audience
and should (1) facilitate interaction among a community of mathematicians and
scientists interested in formulae data related to orthogonal polynomials and
special functions (OPSF); (2) be expandable, allowing the input of new
formulae; (3) be accessible as a standalone resource; (4) have a user
friendly, consistent, and hyperlinkable viewpoint and authoring perspective;
and (5) contain easily searchable mathematics and take advantage of modern
MathML tools for easy to read, scalably rendered mathematics. It is the desire
of our group to build a tool that brings the above features together in a
public website for mathematicians, scientists and engineers. We refer to this
web tool as the Digital Repository of Mathematical Formulae (DRMF).

Our project was motivated by the goal of creating an interactive online
compendium of mathematical formulae. This need was was addressed in SIAM
Activity Group OPSF-Net discussions, such as Dmitry Karp (OPSF-Net 18.4, Topic
#5). In that OPSF-Net edition, there were two related posts (OPSF-Net 18.4,
Topics #6,#7) with a follow-up post in OPSF-Net 18.6, Topic #3.

##  2 Implementation

In our project, we have taken advantage of the free and open source MediaWiki
wiki software as well as tools developed within the DLMF project [14], such as
LaTeXML and the DLMF LaTeX macros. DLMF macros (and extensions as necessary)
tie specific character sequences to unique mathematical objects such as
special functions, orthogonal polynomials, or to other mathematical symbols
associated with these. The DLMF macros are hence used to define OPSF within
DRMF and through LaTeXML, their corresponding rendered mathematical symbols.
Furthermore, the use of DLMF macros as linked to their definitions within the
DLMF, allows for easy access to precise OPSF definitions for the symbols used
within the LaTeX source for OPSF formulae. The committed use of DLMF macros
guarantees a mathematical and structural consistency throughout the DRMF. As a
web tool, the DRMF provides (1) formula interactivity, (2) formula home pages,
(3) centralized bibliography, and (4) mathematical search. The DRMF shares the
core DLMF component, LaTeXML, which (through the MediaWiki math extension)
processes Wikitext math markup written in LaTeX to produce XML and HTML. For
formula interactivity and menus linked to formulae, we have utilized the JOBAD
interactivity framework and are investigating the use of MathJax [4]. We have
also incorporated the MediaWiki: Math and MathSearch [16] extensions. Within
the DRMF, we will develop technology for users to interact with formulae using
a clipboard, which allows for easy copy/paste of formula source
representations (to include LaTeX with DLMF macros; presentation or content
MathML; as well as input formats for computer algebra systems such as
Mathematica, Maple, and Sage).

The DRMF treats formulae as first class objects, describing them in formula
home pages that currently contain: (1) a rendered description of the formula
itself (required); (2) bibliographic citation (required); (3) open section for
proofs (required); (4) list of symbols used and links to their definitions
corresponding to the DLMF macros (required); (5) open section for notes
relevant to the formula (e.g., formula name, if the formula is a
generalization or specialization of some other formula, growth or decay
conditions, links to errata pages, etc.); (6) open section for external links;
(7) substitutions with definitions required to understand the formula; and (8)
constraints the formula must obey. For each formula home page there is a
corresponding talk page, and we are incorporating a strategy for handling the
insertion of formula errata. A major resource in our ability to implement
effective and precise OPSF search will be the use of the DLMF macros in
building the LaTeX source for OPSF formulae and related mathematical data.
Next, we present an overview of the seed resources, which we plan to
incorporate within DRMF. We have been given permission and are seeding the
DRMF with data from the NIST DLMF [14]. We have also been given permission to
and are seeding LaTeX formulae data from [11] (KLS). We will also incorporate
Tom Koornwinder’s companion of recent arXiv published additions to KLS [12].
We have also been given permission to incorporate seed formula data from [5,
6] (BMP). Efforts to upload BMP data, as well as any book data without
existing LaTeX source, will prove extremely difficult, since this effort will
rely on the use of mathematical optical character recognition (OCR) software
such as InftyReader to produce LaTeX source for these formulae. Mathematical
OCR is still in its nascence and this effort is currently under consideration
for feasibility of use. We are in communication with other authors and
publishers to gain access and permission for other proven sources of
mathematical OPSF formulae such as [1, 7, 9, 10] and we are are excited about
the prospect of seeding proof data by Victor Moll and collaborators (see for
instance [8]). For LaTeX source where DLMF macros are not present (such as
KLS), we are developing tools which automate DLMF macro replacements. Seeding
and generating symbol lists are accomplished by converting LaTeX source into
Wikitext, in an automated fashion. We use Pywikibot to automate the uploading
of Wikitext pages to our demo site.

Acknowledgements 111The mention of specific products, trademarks, or brand
names is for purposes of identification only. Such mention is not to be
interpreted in any way as an endorsement or certification of such products or
brands by the National Institute of Standards and Technology, nor does it
imply that the products so identified are necessarily the best available for
the purpose. All trademarks mentioned herein belong to their respective
owners.
We are deeply indebted to Deyan Ginev for sharing with us his expansive vision
and especially for his support in the development of our proof of concept.
Without his guidance and coding, our present demonstration would not be
possible. We would also like to thank Bruce Miller at NIST for his invaluable
contributions regarding LaTeXML. We would also like to express our deep
gratitude to the KWARC group at Jacobs University, Bremen, Germany, and
especially to its group leader, Michael Kohlhase, for his advice and for
access to his group’s mathweb server for our initial DRMF development. We
would also like to thank Dan Lozier, Tom Koornwinder, Dmitry Karp, Dan
Zwillinger, Victor Moll, and Hans Volkmer for offering their advice and for
valuable discussions.

## References

  * Andrews et al., [1999] Andrews, G. E., Askey, R., and Roy, R. (1999).  Special Functions, volume 71 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge.
  * Brychkov, [2008] Brychkov, Y. A. (2008).  Handbook of Special Functions: Derivatives, Integrals, Series and Other Formulas.  Chapman & Hall/CRC Press, Boca Raton-London-New York.
  * Byrd and Friedman, [1954] Byrd, P. F. and Friedman, M. D. (1954).  Handbook of Elliptic Integrals for Engineers and Physicists.  Die Grundlehren der mathematischen Wissenschaften in Einzeldarstellungen mit besonderer Berücksichtigung der Anwendungsgebiete. Bd LXVII. Springer-Verlag, Berlin.
  * Cervone, [2012] Cervone, D. (2012).  Mathjax: A Platform for Mathematics on the Web.  Notices of the American Mathematical Society, 59(2):312–316.
  * [5] Erdélyi, A., Magnus, W., Oberhettinger, F., and Tricomi, F. G. (1981).  Higher Transcendental Functions. Vols. 1-3.
  * [6] Erdélyi, A., Magnus, W., Oberhettinger, F., and Tricomi, F. G. (1954).  Tables of Integral Transforms. Vols. 1-2.
  * Gasper and Rahman, [2004] Gasper, G. and Rahman, M. (2004).  Basic Hypergeometric Series, volume 96 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, second edition.  With a foreword by Richard Askey.
  * Glasser et al., [2012] Glasser, L., Kohl, K. T., Koutschan, C., Moll, V. H., and Straub, A. (2012).  The integrals in Gradshteyn and Ryzhik. Part 22: Bessel-$K$ functions.  Scientia. Series A. Mathematical Sciences. New Series, 22:129–151.
  * Gradshteyn and Ryzhik, [2007] Gradshteyn, I. S. and Ryzhik, I. M. (2007).  Table of Integrals, Series, and Products.  Elsevier/Academic Press, Amsterdam, seventh edition.
  * Ismail, [2005] Ismail, M. E. H. (2005).  Classical and Quantum Orthogonal Polynomials in one Variable, volume 98 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge.  With two chapters by Walter Van Assche, With a foreword by Richard A. Askey.
  * Koekoek et al., [2010] Koekoek, R., Lesky, P. A., and Swarttouw, R. F. (2010).  Hypergeometric Orthogonal Polynomials and their $q$-analogues.  Springer Monographs in Mathematics. Springer-Verlag, Berlin.  With a foreword by Tom H. Koornwinder.
  * Koornwinder, [2014] Koornwinder, T. H. (2014).  Additions to the formula lists in “Hypergeometric orthogonal polynomials and their $q$-analogues” by Koekoek, Lesky and Swarttouw.  arXiv:1401.0815.
  * Magnus et al., [1966] Magnus, W., Oberhettinger, F., and Soni, R. P. (1966).  Formulas and Theorems for the Special Functions of Mathematical Physics.  Third enlarged edition. Die Grundlehren der mathematischen Wissenschaften, Band 52. Springer-Verlag New York, Inc., New York.
  * Olver et al., [2010] Olver, F. W. J., Lozier, D. W., Boisvert, R. F., and Clark, C. W., editors (2010).  NIST Handbook of Mathematical Functions.  Cambridge University Press, Cambridge. Companion to http://dlmf.nist.gov.
  * Prudnikov et al., [1986] Prudnikov, A. P., Brychkov, Y. A., and Marichev, O. I. (1986).  Integrals and Series. Vols. 1-5.  Gordon & Breach Science Publishers, New York.
  * Schubotz, [2013] Schubotz, M. (2013).  Making Math Searchable in Wikipedia.  Conferences on Intelligent Computer Mathematics, abs/1304.5475.


##  1 Methods

###  1.1 Hurricane Sandy Test Case

The physical processes involved in the storm surge created by a hurricane are
essentially barotropic in nature and can therefore be modeled by a single-
layer ocean, the most important physical processes being forcing due to tidal
cycles, surface wind stress, and atmospheric pressure [pringle2021,
mandli2014]. Historically, MPAS-O has been used primarily for long-term (on
the order of hundreds of years) climate modeling; such models require a multi-
layered ocean that can resolve both barotropic and baroclinic motions. As
such, it is uncommon for MPAS-O to be run in a single layer configuration and
only recently was the ability to compute tidal forcing added to the model
[brus2022, barton2022].

Figure 1 shows observed sea-surface height (SSH) at a particular tidal gauge
versus time compared to the SSH predicted by MPAS-O in a multi-layer
configuration and a single-layer configuration. In Figure 1, the multi-layered
model uses a barotropic-baroclinic split-explicit time-stepping scheme
[higdon2005], and the single-layer model uses RK4. This figure shows that both
the single-layer and multi-layer models include the relevant physical
processes of tidal forcing, which produces semidiurnal oscillations, and the
surface winds and atmospheric pressure, which produce the storm surge.

Previous work on LTS in the MPAS framework has been devoted to the shallow
water equations (SWEs) [capodaglio2022], which are similar to the single-
layer, barotropic model used in this work. As a matter of fact, the present
model differs from the SWEs only in having realistic bathymetry, realistic
coastlines, and additional forcing terms needed to take tides, wind,
atmospheric pressure, and drag exerted on the fluid by the ocean floor into
account. The affinity between the SWEs and the present model, and the need for
high local resolution (generally with cells smaller than 1-5 kilometers in
width) for the accurate prediction of sea-surface height in coastal regions,
make storm surge modeling a good real-world test-case for LTS. Meshes for this
use case can be generated such that most of the global ocean is resolved by
relatively large cells and only areas of particular interest are highly
resolved by small cells (see Figure LABEL:fig:delaware_local).

In mandli2014 and pringle2021, where mesh design for storm surge modeling is
discussed, both works employ spatial grids that have a high resolution around
coastlines while the global ocean is generally kept at a low resolution. On
such meshes, the benefits of LTS are clear; one can use a small time-step on
the relatively small region(s) of high resolution, while being free to use a
much larger time-step elsewhere on the globe, with consequent savings in terms
of computational time.

Figure 1:  Sea-surface height (SSH) during Hurricane Sandy as predicted by a
multi-layer ocean model and single-layer model considered here compared to
observed data on the DelBay2km mesh (see Table 1). Split-explicit refers to
the first-order, explicit, split barotropic-baroclinic time-stepping scheme
that is the method of choice for a multi-layer model in MPAS-O [ringler2013].

####  1.1.1 Governing Equations

The momentum and thickness equations are given by eqn:model. Here the
thickness equation is conservation of volume for an incompressible fluid,
where the volume is normalized by the cell area, which is constant in time. In
these equations $\mathbf{u}$ is the horizontal velocity, $t$ is the time
coordinate, $f$ is the Coriolis parameter, $\mathbf{k}$ is the local vertical
unit vector, $K=\frac{\abs{\mathbf{u}}^{2}}{2}$ is the kinetic energy per unit
mass, $\rho_{0}$ is the (constant) fluid density, $p^{s}$ is the surface
pressure, $g$ is the gravitational constant, $\eta$ is the sea-surface height
perturbation, $\eta_{\text{EQ}}$ is the sea-surface height perturbation due to
equilibrium tidal forcing [barton2022], $\beta$ is the self-attraction and
loading coefficient [barton2022], $\chi$ is the spatially varying topographic
wave drag coefficient, $\frac{\mathcal{C}}{H}$ is an inverse time scale, $H$
is the resting depth of the ocean, $h$ is the total ocean thickness such that
$h=H+\eta$ $\mathcal{C}_{\text{D}}$ is the bottom drag coefficient,
$\mathcal{C}_{\text{W}}$ is the wind stress coefficient, and
$\mathbf{u}_{\text{W}}$ is the horizontal wind velocity.

$\cases{\dfrac}{\partial\mathbf{u}}{\partial
t}+\left(\nabla\times\mathbf{u}+f\mathbf{k}\right)\times\mathbf{u}=&-\nabla
K-\frac{1}{\rho_{0}}\nabla
p^{s}-g\nabla\left(\eta-\eta_{EQ}-\beta\eta\right)\\\
-\chi\dfrac{\mathcal{C}\mathbf{u}}{H}-\mathcal{C}_{\text{D}}\dfrac{\abs{\mathbf{u}}\mathbf{u}}{h}+\mathcal{C}_{\text{W}}\dfrac{\abs{\mathbf{u}_{\text{W}}-\mathbf{u}}\left(\mathbf{u}_{\text{W}}-\mathbf{u}\right)}{h}\,,\\\
\dfrac{\partial h}{\partial t}+\nabla\cdot(h\mathbf{u})=0\,.$ (1)

The wind velocity $\mathbf{u}_{\text{W}}$ is linearly interpolated from data
observed at one-hour increments between 10/10/2012 and 11/03/2012. Hurricane
Sandy hit the Eastern US coast on October 29, 2012. The wind and atmospheric
surface pressure data were obtained from the Climate Forecast System Version 2
(CFSv2) reanalysis product [saha2014].

The model equations are spatially discretized using a C-grid type finite
volume method called the TRiSK scheme wherein the fluid thickness $h$ is
stored at cell centers and the normal component of the fluid velocity
$\mathbf{u}$ is computed at cell edges. This method has been shown to conserve
the total energy of a system and robustly simulates potential vorticity
[thuburn2009, ringler2010]. The current code configuration was chosen to show
the performance improvements possible with LTS. Some processes, like wetting
and drying cells and spatially varying bottom friction, were not included here
because they are still under development, but would be expected to improve
model results relative to observations.

###  1.2 Time-stepping Schemes

The specific LTS scheme used in this paper was originally developed for use in
the MPAS framework by hoang2019. This scheme is based on so-called, strong
stability preserving Runge-Kutta (SSPRK) methods, which are explicit methods
for solving systems of ODEs resulting from the spatial discretization of
hyperbolic conservation laws. SSPRK methods satisfy the total variation
diminishing (TVD) property, which means that given a sufficiently small time-
step, the total variation of an SSPRK solution will not increase over time,
which implies that the solution will remain stable. For a more complete
discussion of SSPRK methods, see gottlieb1998 and gottlieb2001.

The LTS method used here is derived from a three-stage, third-order SSPRK
method (SSPRK3) and is itself third-order; we refer to this method as LTS3.
LTS3 was implemented in the MPAS framework by capodaglio2022. While a full
review of the LTS3 scheme is outside the scope of this paper, and can be found
in hoang2019 and capodaglio2022, we give a short description of how the scheme
works for completeness. Consider Figure 2, which serves as a hypothetical mesh
wherein there are four classes of cells. The red cells are referred to as the
coarse cells, the yellow as interface 2 cells, the pink as interface 1 cells,
and the blue as fine cells. The coarse, interface 1, and interface 2 cells all
advance with the coarse time-step $\Delta t_{\text{coarse}}$ and the fine
cells advance with the fine time-step $\Delta t_{\text{fine}}$ such that
$\Delta t_{\text{fine}}=\frac{\Delta t_{\text{coarse}}}{M}$ for some positive
integer $M$. For practical applications, the labels fine and coarse denote the
regions by their cell size, but the algorithm simply advances two different
time steps in the two regions, so the regions could have other cell sizes—for
example, identically-sized cells for testing purposes.

Figure 2: A mesh with cells labeled for LTS. Red cells are coarse, yellow
cells are interface 2, pink cells are interface 1, and blue are fine.

The LTS3 scheme proceeds as follows:

  1. 1.

Starting from time $t^{n}$, advance the solution on the coarse, interface 1,
and interface 2 cells with the first two stages of SSPRK3 using $\Delta
t_{\text{coarse}}$.

  2. 2.

Sub-cycle on the fine cells; advance to $t^{n+1}$ by repeating all three
stages of SSPRK3 $M$ times using $\Delta t_{\text{fine}}$.

     * •

This requires that we know values for $\mathbf{u}$ and $h$ at intermediate
time-levels on interface layer 1 cells. Obtain predicted values for these with
an appropriate prediction step using data from time $t^{n}$ and the first two
stages of SSPRK3 that were obtained in step 1.

  3. 3.

Advance to time $t^{n+1}$ on the coarse cells with the final stage of SSPRK3.

  4. 4.

Correct the values for $\mathbf{u}$ and $h$ at time $t^{n+1}$ on interface 1
and interface 2 cells by accounting for fluxes coming from the fine cells
during the sub-cycling in step 2.

######  Remark 1.1.

(LTS3 Step 1)  In order to advance with the second stage of SSPRK3 on the
interface cells that border the fine region, we need the first stage SSPRK
data on the three layers of fine cells adjacent to the interface region. Three
layers are required here because the spatial discretization of the model
equations mandates that a given cell’s advancement depends on the three layers
of adjacent cells. This data is not used to advance the fine cells themselves,
but this fact does require that these specific fine cells, referred to as the
interface-adjacent fine cells, be large enough that they admit the coarse
time-step.

######  Remark 1.2.

(LTS3 Step 2)  The prediction of values for $\mathbf{u}$ and $h$ at
intermediate time-levels on interface 1 cells can be thought of as an
interpolation of the existing data from $t^{n}$ and the first two stages of
SSPRK3. This is a substantial simplification and we refer the reader to
hoang2019 for a full description of this step.

This scheme is $\mathcal{O}\left((\Delta t)^{3}\right)$ everywhere, including
on the interface layers. Figures 5 and 5 show the convergence in time for both
$\mathbf{u}$ and $h$ using RK4 and LTS3 respectively, on the spatially
discretized versions of the model equations eqn:model. The root-mean-square
(RMS) error is defined as

$E_{\text{RMS}}=\sqrt{\frac{\sum_{i=1}^{N}(s_{i}-m_{i})^{2}}{N}},$ (2)

where $\\{s_{i}\\}_{i=1}^{N}$ is the discrete reference solution,
$\\{m_{i}\\}_{i=1}^{N}$ is the discrete model solution, and $N$ is the number
of discretization points.

0.65
0.65

Figure 3:  Figure 4:  Figure 5:  Temporal convergence for layer thickness $h$
and normal velocity $\mathbf{u}$ and using RK4 (a) and LTS3 (b) on the
discretized version of the model equations eqn:model. Convergence tests were
conducted on the DelBay2km mesh (see Table 1 and Figure 6). For LTS3, we used
$M=4$. The errors were calculated against a reference solution generated by
RK4 with $\Delta t=0.1$. The run duration is 1 hour. The RMS error is as
defined in Eq. 2.

######  Remark 1.3.

One should note that RK4 uses one more Runge-Kutta stage per time-step than
LTS3 and therefore has a lessened restriction on the size of the largest time-
step it can use on high-resolution cells, e.g. see Table 1 and note that on
the mesh labeled DelBay2km, RK4 uses a global time-step of 30 seconds, while
LTS3 uses a fine time-step of 18 seconds. Consider that, on the smallest
cells, RK4 advances on average $\frac{30}{4}=7.5$ seconds per stage while LTS3
advances on average $\frac{18}{3}=6$ seconds per stage, which means that RK4
would be 25% faster were LTS3 using a global time-step. This puts LTS3 at an
inherent disadvantage when comparing computational performance on these
meshes. Nevertheless, we will show that LTS3 is still capable of considerably
outperforming RK4 in almost all of our case studies. We have chosen to use RK4
as our point of comparison rather than a third order method because RK4 is the
state-of-the-art scheme in MPAS-O for a single-layer model.

###  1.3 Meshes and LTS Configuration

Here we consider five different meshes; the relevant parameters for each are
given in Table 1. Of particular interest is the physical distribution of
regions of low and high resolutions. Each mesh is divided into five regions,
each of which is populated by cells of a different size. These regions are
Delaware Bay, the area around the Delaware coast, the area around the eastern
coast of the US, the western Atlantic ocean, and the rest of the globe. In the
Resolutions section of Table 1, we give the cell width in kilometers for each
of these five regions, for all meshes considered. In particular, note that the
highest resolutions in these meshes range from 2 kilometers down to 125 meters
and the ratio of the lowest to highest resolution ranges from 60 to 240. It
should be noted that even our lowest resolution mesh using 2 kilometer cells
in Delaware Bay is of extremely high resolution for MPAS-O. Simulations with
these ultra-high-resolution meshes are possible because the single-layer
configuration is much faster than standard MPAS-O climate domains, which use
60 to 100 layers, vertical advection and diffusion, and additional
parameterizations.

All meshes were created with JIGSAW [engwirda2017], a library that can quickly
generate high-quality variable-resolution meshes. Global meshes are designed
by specifying the cell width distribution across the sphere (Figure 6), which
is then passed into JIGSAW. The five regions of the mesh were used as a
straightforward means of providing regional refinement in the Delaware Bay
estuary within a global mesh. The Delaware Bay region is the highest level of
refinement, while the Delaware coast, and eastern US coast regions provide a
smooth transition to the coarser global resolution. The western Atlantic Ocean
region uses a higher level of refinement than the global ocean to better
resolve the hurricane wind and pressure fields, while the rest of the Earth
has the coarsest resolution to reduce the total number of cells in the mesh.

| DelBay2km | DelBay1km | DelBay500m | DelBay250m | DelBay125m
---|---|---|---|---|---
Resolutions: Grid cell width (km) |  |  |  |  |
Global background | 120 | 60 | 30 | 30 | 30
Western Atlantic | 30 | 30 | 15 | 15 | 15
Eastern US coast | 10 | 5 | 2.5 | 1.25 | 0.625
Delaware coast | 5 | 2.5 | 1.25 | 0.625 | 0.3125
Delaware Bay | 2 | 1 | 0.5 | 0.25 | 0.125
Mesh Parameters |  |  |  |  |
Number of cells | 58,240 | 198,776 | 794,172 | 1,591,416 | 4,617,565
Number of MPI ranks | 2 | 8 | 32 | 64 | 178
Thousands of cells per MPI rank | 29 | 25 | 25 | 25 | 26
LTS Parameters |  |  |  |  |
Eastern US Coast (EC) Fine Region |  |  |  |  |
Number of interface layers | 2 | 5 | 10 | 20 | 50
Count ratio | 1.88 | 1.71 | 1.71 | 0.46 | 0.12
Resolution ratio | 15 | 30 | 30 | 60 | 120
$\Delta t_{\text{RK4}}$ (s) | 30 | 15 | 7.5 | 3.75 | 1.875
$\Delta t_{\text{fine}}$ (s) | 18 | 8 | 4 | 2 | 1
$\Delta t_{\text{coarse}}$ (s) | 72 | 72 | 24 | 24 | 24
$M$ | 4 | 9 | 6 | 12 | 24
Western Atlantic (WA) Fine Region |  |  |  |  |
Number of interface layers | 2 | 5 | 10 | 20 | 50
Count ratio | 0.92 | 1.28 | 1.28 | 0.39 | 0.11
Resolution ratio | 60 | 60 | 60 | 120 | 240
$\Delta t_{\text{RK4}}$ (s) | 30 | 15 | 7.5 | 3.75 | 1.875
$\Delta t_{\text{fine}}$ (s) | 18 | 8 | 4 | 2 | 1
$\Delta t_{\text{coarse}}$ (s) | 306 | 152 | 80 | 80 | 80
$M$ | 17 | 19 | 20 | 40 | 80
|  |  |  |  |

Table 1:  Relevant parameters for each mesh and LTS configuration used in performance experiments.  |
---|---
| Atlantic Ocean | Eastern US Coast
DelBay2km |  |
DelBay1km |  |
DelBay500m |  |
DelBay250m |  |
DelBay125m |  |
Figure 6:  Global cell width in kilometers for each mesh described in Table 1.
Each subplot shares the same color scale, which is in log-space, so
differences in high-resolution areas can be easily distinguished.

####  1.3.1 Choice of Time-steps

The time-steps used on each mesh and each LTS configuration for performance
tests are given in Table 1. These time-steps were obtained experimentally and
are the largest values that are admissible to guarantee stability of the
model.

Because there are two different time-steps for the LTS scheme and these must
differ by a factor of $M$ such that $\Delta t_{\text{fine}}=\frac{\Delta
t_{\text{coarse}}}{M}$, they can be obtained in different orders; one could
find the largest time-step admittable on the coarse region of the mesh then
find the smallest value of $M$ that gave an admittable fine time-step, or
first find the largest admittable fine time-step then the largest value of $M$
that gives an admittable coarse time-step. This has the result of the time-
step in one region or the other not technically being maximal.

There are advantages to both methods that depend on the mesh itself. For
instance, on a mesh where there are sufficiently more coarse cells than fine
cells, one might wish to first maximize the time-step in the coarse region and
take a small penalty to the size of the fine time-step. In this paper however,
we have opted to find admittable time-steps by first maximizing the fine time-
step, then by finding the largest value of $M$ that gives a valid coarse time-
step with the intention on minimizing the work done in the fine region.

####  1.3.2 LTS Mesh Parameters

When designing a mesh for use with LTS there are two parameters of particular
interest, which we refer to as the count ratio and the resolution ratio. Both
of these parameters depend on the number of cells where the fine time-step is
used, referred to as the fine cells, and the number of cells where the coarse
time-step is used, referred to as the coarse cells. Note that, in this
discussion, the label of fine or coarse does not necessarily reference the
size of a cell, but rather which time-step is used to advance it. The region
made up of fine cells is called the fine region while the rest of the globe is
called the coarse region.

The count ratio is the ratio of the number of cells in the coarse region to
the number of cells in the fine region, i.e.

$\text{countratio}=\frac{\text{numberofcoarsecells}}{\text{numberoffinecells}}.$

The resolution ratio is the ratio of the cell width of the coarse cells to the
cell width of the fine cells, i.e.

$\text{resolutionratio}=\frac{\text{cellwidthofcoarsecells}}{\text{cellwidthoffinecells}}.$

In the case where there are cells of multiple resolutions in either region, as
is the case in our meshes, we consider the smallest value of cell width in a
given region as it is the smallest cell that restricts the size of the time-
step admittable in that region. For example, in DelBay2km with the EC
configuration the smallest fine cells have a cell width of 2 kilometers and
the smallest coarse cells have a cell width of 30 kilometers so we would say
that the resolution ratio was 15.

These parameters are the primary drivers behind the performance of LTS. The
higher the resolution ratio is, the larger the coarse time-step can be,
allowing LTS to take significantly less time-steps on coarse cells compared to
a global method, which must use a small time-step everywhere. Similarly, the
higher the count ratio is, the more coarse cells there are in a mesh where LTS
has to do less work than a global method.

In a preliminary analysis, we investigated how varying the count and
resolution ratios affected the performance of LTS in an idealized test case.
Namely, in the shallow-water core of MPAS-O, we created a test case following
test case five from williamson1992 and generated a series of meshes that had
varying values for the count ratio and the resolution ratio. We produced the
plots in Figure 9, which are relevant to the work done here for two reasons.
First, the shallow-water equations used to generate this data are identical to
those used in our Hurricane Sandy test case, except for additional forcing
terms used here. That is, we can expect for the results shown in Figure 9 to
be similar to the results we would get with the hurricane Sandy test case.
Second, this data clearly demonstrates the importance of both the count and
resolution ratios when considering the performance of LTS. We can see that as
either parameter increases, the speedup increases, and that this increase is
drastic close to zero. The curves in both plots eventually level off as the
speedup becomes limited by the other parameter, which is fixed.

0.8
0.8

Figure 7:  Figure 8:  Figure 9:  Speedup versus mesh design parameters: (a)
varying count ratio while resolution ratio is held constant at 32; (b) varying
resolution ratio as count ratio is held near one. In both cases, speedup
increases and then levels off at 50–70%. Speedup is calculated as in Eq.
LABEL:eqn:speedup.

Figure 12 shows two configurations for the placement of the fine and coarse
regions. The configuration in Figure 12 will be referred to as the EC
configuration and the configuration in Figure 12 will be referred to as the WA
configuration. We consider these different configurations because moving the
fine region to include the western Atlantic cells greatly increases the size
of our coarse time-step, with the penalty of adding a non-trivial number of
cells to our fine region where we must use the fine time-step restricted by
the smallest cell in the region. This trade-off between the count ratio and
the resolution ratio is a significant part of our investigation and will be
discussed in Section LABEL:subsec:performance.

0.45  0.45

Figure 10: East US coast (EC) fine region Figure 11: Western Atlantic (WA)
fine region Figure 12: The two configurations of the fine region used for
DelBay meshes: the smaller East US Coast fine region (a) and the larger
Western Atlantic fine region (b).

####  1.3.3 Load Balancing, Communication, and Additional Interface Layers

Performing the LTS3 algorithm in parallel requires careful consideration of
load balancing. LTS requires that cells in a given mesh be sorted into
different groups, some are coarse cells, some are fine cells, and some are
interface cells. Each type of cell requires a different amount of work, so
load balancing among MPI processes is non-trivial; an efficient way to address
this has been a primary concern of capodaglio2022. For effective load-
balancing, each MPI rank is given three sets of cells, one set of fine cells,
one set of coarse cells, and one set of interface cells. The different groups
of cells are spread across each rank so that each process has an approximately
equal number of each type. This helps to ensure that no process is idle while
it waits for the others to finish work. An issue with this scheme is that in
most mesh configurations, there are very few interface cells relative to the
number of fine or coarse cells. As a result, when running with more than a few
processes, there may not be enough interface cells in the mesh for an equal
number to be spread across processes in a way that also respects the expense
of MPI communication calls.

Figure 13:  An example of a mesh needing additional layers of interface cells
for load balancing in DelBay1km. Here, there are 10 interface layers.

To avoid this issue, additional interface layers are added to the mesh (see
Figure 13). Because LTS3 is $\mathcal{O}\left((\Delta t)^{3}\right)$
everywhere on the mesh, these additional interface layers have no negative
impact on the model’s accuracy. Adding these additional interface layers
allows us to insure that there are enough interface cells to spread across a
given number of processes. It is shown in capodaglio2022 that, as a rule of
thumb, each process should own at least 100 interface cells to achieve
sufficient load-balancing. On each mesh and for each LTS fine region
configuration, we have set a number of interface layers that ensures this
requirement is met. The number of interface layers used in each case is
reported in Table 1.


# Direct 3D Printing of Soft Fluidic Actuators with Graded Porosity

Nick Willemstein
Department of Biomechanical Engineering
University of Twente
Enschede, The Netherlands
<EMAIL_ADDRESS>
Herman van der Kooij
Department of Biomechanical Engineering
University of Twente
Enschede, The Netherlands
<EMAIL_ADDRESS>
Ali Sadeghi
Department of Biomechanical Engineering
University of Twente
Enschede, The Netherlands
<EMAIL_ADDRESS>

###### Abstract

New additive manufacturing methods are needed to realize more complex soft
robots. One example is soft fluidic robotics, which exploits fluidic power and
stiffness gradients. Porous structures are an interesting type for this
approach, as they are flexible and allow for fluid transport. Within this
work, the Infill-Foam (InFoam) is proposed to print structures with graded
porosity by liquid rope coiling (LRC). By exploiting LRC, the InFoam method
could exploit the repeatable coiling patterns to print structures. To this
end, only the characterization of the relation between nozzle height and coil
radius and the extruded length were necessary (at a fixed temperature). Then
by adjusting the nozzle height and/or extrusion speed the porosity of the
printed structure could be set. The InFoam method was demonstrated by printing
porous structures using styrene-ethylene-butylene-styrene (SEBS) with
porosities ranging from 46% to 89%. In compression tests, the cubes showed
large changes in modulus (more than 200 times), density (-89% compared to
bulk), and energy dissipation. The InFoam method combined coiling and normal
plotting to realize a large range of porosity gradients. This grading was
exploited to realize rectangular structures with varying deformation patterns,
which included twisting, contraction, and bending. Furthermore, the InFoam
method was shown to be capable of programming the behavior of bending
actuators by varying the porosity. Both the output force and stroke showed
correlations similar to those of the cubes. Thus, the InFoam method can
fabricate and program the mechanical behavior of a soft fluidic (porous)
actuator by grading porosity.

##  1 Introduction

Soft robotics utilize softness and smart design to realize robots that are
inherently safe for human-robot interaction and compliant to unforeseen
disturbances. Soft fluidic robotics are among the most popular type due to
their flexibility, wide range of applications, and (relative) simplicity for
actuation. They have been used in a broad range of applications from grippers1
to autonomous systems2.

Soft fluidic robotics work by transferring fluidic power (flow and pressure)
from one place to another. This power can be used for actuation1, sensing3,
and control2, 4. Soft fluidic actuators (SFAs) consist of mechanically
compliant chambers that can be deformed using fluidic power. The deformation
of these chambers is defined by the stiffness gradient of the structure. This
gradient allows for motions such as translation, bending, and twisting 5.
Researchers employ a broad number of methods to realize this stiffness
gradient, which includes exploiting material properties and/or geometry. An
approach is to add a second material such as fibers5 textile1, and/or paper6
in a dedicated geometry. Other approaches use a single material and exploit
geometrical features, such as origami7 and kirigami8. A prominent example of
exploiting internal geometry is porous structures, such as foam.

An advantage of porous structures, is that they can, inherently, allow for
fluidic transport and are flexible. In addition, porous structures are
compressible. This feature has been exploited to realize vacuum-driven
actuators. Examples thereof include: a continuum soft robot9, high-force
sensorized foam actuators10, and bending actuators11. Lastly, positive
pressure actuation is also feasible with porous structures for both
actuation12 and control13.

Existing methods for fabricating porous structures span a broad spectrum of
approaches. A popular approach is to use commercial foams combined with semi-
automatic/manual processing, such as laser-cutting, gluing, and coating with
silicone or laminating (for airtightness)9, 10, 14. Another method is to use
casting with a sacrificial scaffold/material such as salt 12, 13, sugar15 or
Polyvinyl Alcohol (PVA)16 that can be removed afterwards.

To realize the multifunctional nature of soft robotic systems (including SFAs
based on porous structures), additive manufacturing (AM) is considered a
promising fabrication method. Literature has already demonstrated that AM can
realize (sensorized) SFAs17, 18, 19, 20, 21. However, the fabrication of
porous structures for soft robots by AM methods is less developed.

Recently, a custom ink22 for realizing local porosity was developed by adding
a porogen. By regulating the ratio of ammonium bicarbonate (the porogen) and
silicone rubber the porosity can be changed locally. This approach requires
modification of the material itself.

Besides using a chemical reaction, a porous structure can also be realized by
the fabrication of miniature patterns during the AM process. Liquid rope
coiling (LRC) is an interesting candidate for the fabrication of these
patterns. LRC is the coiling of a viscous fluid, such as honey, due to
buckling when falling from an elevated height23.

LRC has already been exploited to fabricate porous structures24, 25, 26, 27,
28. Results in the literature indicate that LRC can realize stiffness changes
of more than one order of magnitude24, 28. In addition, it has been shown to
be feasible to print different patterns by proper selection of process
parameters. 26, 25, 28 Such an approach has been applied to AM processes
extrude molten materials that solidify in the process, such as glass25 and
thermoplastics27, 28. Although these initial results indicate that porous
structures can be realized with LRC and AM, its application to soft robotics
has, to the knowledge of the authors, not yet been explored. In addition, we
also explore how to use LRC to create multiple levels of porosity in a single
structure to realize mechanical programming with porosity. Within this work, a
new printing approach that exploits LRC for the fabrication of porous
structures using hyperelastic materials is proposed. This new approach is
referred to as the Infill Foam (InFoam) approach. We demonstrate how LRC’s
coiling behavior can be incorporated into the fused deposition modeling (FDM)
process to fabricate porous hyperelastic structures. To use hyperelastic
materials within the FDM process a screw extruder (see Fig. 1(a)) was used29,
20. Subsequently, the InFoam method is used for mechanical programming of
stiffness, density, and energy dissipation. Then the extension to graded
porosity is demonstrated as a tool for programming the deformation and to
program the behavior of soft bending vacuum actuators.

##  2 Materials and Methods

Figure 1: (a) Screw extruder, (b) liquid rope coiling setup with important
parameters, (c) coil geometry, (d) grading approach, (e) cubes with different
porosity under the same load, and (f) porosity gradients for soft vacuum
actuators.

###  2.1 From Coiling to Printing

Liquid rope coiling is a phenomenon that can be observed when extruding a
viscous liquid above the printing (see Fig. 1(b)). This situation is also
observed in extrusion-based AM methods such as FDM. When there is a gap
between the nozzle and printing surface, the deposited material acts as a
flexible (liquid) rope that buckles, which leads to a pattern of coils on the
surface (see Fig. 1(b)). A wide variety of shapes can be achieved through LRC,
such as a meander and figure of eight30, 31. Within this work, the focus is
purely on the circular coil pattern. The repetitive nature of the circular
coiling pattern makes it an interesting tool for 3D printing.

Literature has shown that the shape of the deposited coils (called coiling
pattern) depends on three parameters, namely: (I) the viscosity of the
extruded material, (II) the height ($H$) of the nozzle above the platform,
(III) the ratio of extrusion ($V_{e}$) and printhead speed ($V_{F}$)30, 32.
All three parameters can be set by the user in the FDM process. The viscosity
can (primarily) be set using the printing temperature and the height by moving
the nozzle. Whereas the ratio of extrusion and printhead speed is set by the
extrusion multiplier ($\alpha$). The variable $\alpha$ is the screw rotation
(in radians) per millimeter moved of the nozzle. In essence, $\alpha$ scales
the distance traveled by the printhead to compute the required screw rotation,
which defines the ratio of speeds.

These parameters can be set to control the coiling pattern and to plan the
desired shape. However, modeling the coiling pattern is not straightforward
due to the changes in material properties (such as viscosity) due to the
solidification of the thermoplastic during deposition. Within this work, an
experimental method, similar to25, is used to characterize the coiling
pattern.

Similar to conventional printing, the InFoam method needs the dimensions
(width and height) of the coiling pattern for path planning. The coiling
pattern consists of coils that stack on top of each other in line with the
movement of the printhead (see Fig. 1(b)). The geometry of this coiling
pattern can be used to compute the width and height. As observed in Fig. 1(c),
the width $W$ of the coiling pattern is equal to the diameter of the coil,
which is equal to

$W=2R_{c}+d$ (1)

Within this equation, $R_{c}$ is the coil radius and $d$ the nozzle diameter.
To compute the coil height $h_{c}$, a second parameter of the coiling pattern
needs to be quantified, namely the $N$-value. The $N$-value represents the
coil density within a row of coils, which is defined by the number of coils
within one (outer) coil diameter ($W$). This value can be computed after
deposition by measuring the distance between two coils $\Delta x$ and using
the relation:

$N=\sqrt{\frac{W^{2}}{\Delta x^{2}-d^{2}}}$ (2)

The $N$-value can also be computed from the coiling behavior by examining the
ratio of the extruded length over the linear distance of one outer coil
diameter ($W$). This modelling approach allows for a relation between the
(user set) parameters $\alpha$ and $H$ (controls the value of $R_{c}$) and the
$N$-value. After rewriting the coil density ($N$-value) can be expressed as
(see derivation in Supplementary Information)

$N=\frac{G\alpha W}{2\pi R_{c}}$ (3)

The variable $G$ (mm/rad) represents the length of material extruded per
radian turned by the screw. The exact value of $G$ is defined by the material
and geometry of the screw and barrel33. Within this work, this value is
acquired by fitting experimental data. Thus, in Equation 3 moving the linear
distance $W$ the screw rotates $\alpha W$ radians leading to an extruded
length of $G\alpha W$. This length is used to compute the coil density by
dividing it by the circumference of the coil ($2\pi R_{c}$). Thus, the coil
density can be adjusted by changing the extrusion multiplier ($\alpha$) under
the assumption of constant $R_{c}$. Adjustment of the extrusion multiplier for
the desired density at a different $R_{c}$ requires scaling by the factor
$W/(2\pi R_{c})$. The scaled extrusion multiplier $\alpha W/(2\pi R_{c})$ is
referred to as the coil density multiplier.

The coil density and radius can be used to compute the coil height $h_{c}$.
The coil can be modelled as a rectangle with rounded corners (see Fig. 1(c))
under an angle of $\theta=\arctan\left(\frac{Nd}{W}\right)$. The height is
then equal to

$h_{c}=W\sin\left(\theta\right)+\left(1-\sin\left(\theta\right)\right)d$ (4)

The preceding equations can be used for path planning and characterization of
the coiling pattern. Specifically, Equations 1 and 4 can be used to compute
the width and height of the coiling pattern. Whereas Equation 2 and 3 can be
used to measure and set the desired coil density ($N$-value), respectively. To
use these equations the value of $R_{c}$ needs to be characterized to relate
it to $H$. The geometry of the deposited coiling patterns will introduce
porosity. This porosity can be used for mechanical programming of a
structure34. The possibility of grading the porosity allows for localized
mechanical programming. Combining LRC and the normal plotting of the 3D
printer enables such functionality by using the deposited coils as a scaffold
for the normal plotting. A solution for printing non-porous zones is to move
the nozzle down after extruding a layer of coils (i.e. the scaffold) and then
locally plot additional material(see Fig. 1(d)). This simple approach enables
porosity grading and ensures proper force transfer by welding layers of coils
together.

###  2.2 Mechanical Programming with Porosity

The InFoam method can produce (graded) coiling patterns with different coiling
radii and $N$-values by adjusting the extrusion multiplier ($\alpha$) and/or
the height ($H$). The realized coiling patterns can be used to print porous
structures. The porosity will have a major influence on the resulting
mechanical properties34. The repetitive nature of the coiling pattern was
exploited to estimate the porosity $\phi$ (in percent), which is equal to (see
Supplementary Information for the derivation)

$\phi=100\left(1-\frac{G\alpha\pi d^{2}}{4Wh_{c}}\right)$ (5)

This equation estimates the porosity by dividing the extruded volume
(numerator) by the spanned volume (denominator). The spanned volume is based
on a rectangular bar of coil height $h_{c}$ and outer coil diameter $W$
whereas the length is canceled out (see Supplementary Information). This
equation shows that the user can decrease the porosity nonlinearly by both
$\alpha$ and $H$. The latter is implicit as changing the $H$ will change the
coil radius $R_{c}$, which is related to $W$ and $h_{c}$.

The porosity can be used to program a wide variety of mechanical properties,
which include elastic modulus and density. The elastic modulus decreases
quadratically with the porosity34. The large change in stiffness is visualized
by the samples shown in Fig. 1(e) wherein the same load lead to deformations
ranging from minor ($\phi=46$%) to large ($89$%).

In general, the effect of porosity on mechanical properties can be captured in
a power-law relation through the phenomenological model34:

$\frac{p_{p}}{p_{s}}=C\left(1-\phi/100\right)^{n}$ (6)

Wherein the $p_{p}$ and $p_{s}$ represent the mechanical property of interest
of the porous structure and solid, respectively. Whereas $C$ and $n$ are two
fitting values. The relation between density and porosity for air (and other
low-density fluids), is an example of such a relation (with $n=1$ and $C=1$)

$\frac{\rho_{p}}{\rho_{b}}=1-\phi/100$ (7)

The variables $\rho_{p}$ and $\rho_{b}$ represent the densities of the porous
structure and bulk material, respectively. This equation approximates the
porosity well as low-density fluids (such as air) have densities that are
multiple orders of magnitude below that of the solid material.

###  2.3 Experimental Setup & Methodology

####  2.3.1 3D Printing Setup

For printing with the InFoam method, a modified Creality Ender 5 Plus
(Shenzhen Creality 3D Technology Co., Ltd., China) was used (see Figure S1 in
the Supplementary Information). The printer was modified to use a screw
extruder (see schematic in Fig. 1(a)) to allow for the use of thermoplastic
elastomers. To magnify the torque a 22:63 gear ratio between the stepper motor
and screw was added. However, the extrusion multiplier ($\alpha$) will always
refer to the screw’s rotation. As printing material G1657 styrene-ethylene-
butylene-styrene (SEBS) pellets (Kraton Corporation, USA) with a Shore
hardness of 47A and bulk density of 900 kg/m3 were used. In addition, all
experiments were conducted with a printing temperature of 230∘C and with a
nozzle diameter of $d=0.4$ mm. Lastly, the G-code for printing with the InFoam
method was generated using a custom script in MATLAB (The Mathworks, Inc.,
USA). Printing with the InFoam method is shown in Supplementary Movie 1 in the
Supplementary Information.

####  2.3.2 Characterization of Coiling Pattern

Using the InFoam method for printing complex structures requires
characterization of the coiling pattern (i.e. the coil radius ($R_{c}$) and
density ($N$-value)). The LRC effect indicates that the coil radius and the
$N$-value are controlled by the viscosity, extrusion multiplier ($\alpha$),
and height ($H$). Within this experiment, the $H$ and $\alpha$ were varied at
a fixed temperature but due to the shear-thinning behavior of the
thermoplastic, the extrusion speed will also affect the viscosity (and
therefore the coil radius).

For this characterization, several lines were printed. Afterwards, the printed
lines were scanned with a USB camera, and the acquired images were analyzed
using ImageJ35. To relate the pixel distance to real-world distances dots on
the printbed were used as a reference distance. Specifically, the width ($W$)
and linear distance ($\Delta x$) were measured. These measurements were used
to determine $R_{c}$ and the $N$-value using Equations 1 and 2, respectively.

####  2.3.3 Mechanical Programming of Soft Cubes

To investigate the effect of the coil radius ($R_{c}$) and density ($N$-value)
on the porosity, stiffness, and damping a set of 25x25x25 mm3 cubes were
printed. The cubes were printed (in triplicate) for different heights ($H$)
and $N$-values. The corresponding extrusion multiplier ($\alpha$) and $R_{c}$
were based on the results of the coiling characterization to use for process
planning (i.e Equations 1, 3, and 4). The cubes were divided into two groups
of six (with one the same in both $[H,N]=[6,3]$) with either a change in
$N$-value$=[2.2,3,6.3,8.5,10.75,12.75]$ or $H=[2,4,6,8,10,15]$. All cubes were
printed with the coils in the same direction for each layer.

After printing the cubes, the mass and height were measured using a scale and
caliper, respectively. The resulting porosity was then computed using the
measured density by dividing the mass by the cube’s volume. The resulting
porosity (in percent) is then computed using Equation 7.

Subsequently, the stiffness and damping were investigated using a universal
tensile tester (Instron 3343, Instron, USA). The testing was done in four
phases. Firstly, the cubes were pre-loaded by compressing to 5 mm at a speed
of 5 mm/min. Secondly, it returned to its original position at 60 mm/min.
These two load steps were performed to have a similar pre-load for all cubes.
In the third phase, the cubes were compressed for 21 mm at 10 mm/min. This
loading step was slow to have a primarily elastic response. Lastly, the
tensile tester returned to its original position with a speed of 60 mm/min.
The speed was higher in this step to increase the viscous damping component.
These loading steps were performed for all three axes (each in triplicate).

The stiffness was approximated by computing the segment modulus at different
levels of strain. A quadratic curve was fitted at specific levels of strain
using the twelve preceding points. The result of this fitting was used to
compute the segment modulus. In addition, the dissipated energy by hysteresis
was investigated to evaluate the damping. The dissipated energy was computed
by taking the ratio of the stress-strain curve’s inner shape (i.e. between the
forward and backward motion) and the overall curve.

####  2.3.4 Graded Porosity for Soft Fluidic Actuators

By moving from a single level of porosity to a graded porosity, the
deformation can be programmed. To demonstrate this grading capability, a set
of soft vacuum actuators with different porosity gradients were printed. These
included the contraction pattern of Fig. 1(f), which used high (80+%) and low
porosity discs stacked on top of each other to realize contraction. The low
porosity discs had a low porosity ring (<20%) on the outside and high porosity
inside to improve contraction but preserve airflow (the exceptions are the top
and bottom parts, which are completely low porosity). In addition,bending
actuators were printed (see Fig. 1(f)) with and without spacer (with a
porosity of 83.8% and spacers with 15%). A zero porosity layer was added as a
stiff substrate, such that when a vacuum is applied the structure would
collapse in a bending motion. Similarly, a twisting actuator was printed with
the spacers at a 45-degree angle to make a twisting motion the preferred
deformation. Lastly, a screw pattern (see Fig. 1(f) with high porosity of 81%
and spacers of 15%) had two zero-porosity sections perpendicular to each
other. By spacing these two sections in succession, actuation would lead to a
screw-like motion. The outer dimensions were set the same as the bending
actuators. After printing, all printed actuators were put in a 0.4 mm thick
SEBS (heat-sealed) sleeve and were actuated using a vacuum.

The bending actuators with spacers were further investigated to characterize
the effect of porosity on the output force and stroke. To this end, four soft
bending actuators were printed with different levels of porosity
($\phi=[77.6,79.1,83.8,85.2]$). The nominal dimensions of these actuators were
15 mm wide, 75 mm long and 8 mm high. These actuators used a zero porosity
strip as a (relatively) stiff substrate (1 mm thick). Subsequently, the
actuators were put in a heat-sealed SEBS sleeve (0.4 mm thick). Both the
sleeves and actuators were fabricated in triplicate. The setup used for
characterization of the output force and stroke evaluation is shown in Fig. 2.
The swing acts as a grasping object with a bearing to reduce friction. An
LSB200 loadcell (Futek, USA) was integrated into the swing to measure the
output force, which was connected to an Arduino Uno (Arduino AG, Italy)
through an HX711 loadcell amplifier (Avia Semiconductor, China). The loadcell
was tared right before starting a measurement to offset the mass of the
actuator. For actuation, a vacuum source was operated using a manually
triggered valve, and the vacuum pressure was set using a flow valve. The
pressure was read by the Arduino Uno using an MXP5011DP pressure sensor (NXP
Semiconductors, The Netherlands). The Arduino and a USB webcam (to record the
actuator’s deformation) were connected through USB to MATLAB. This experiment
was performed in triplicate for each actuator for 30 seconds.

Figure 2: Actuator characterization setup.

For analysis of the data, a Maxwell viscoelastic model36 with one term was
fitted to the force data, which is defined as:

$F(t)=(K_{0}+K_{1}\exp(-t/\tau_{1}))F_{max}$ (8)

Wherein $t$ is the time in seconds after the maximum force $F_{max}$ with
$K_{0}$ and $K_{1}$ fitting parameters and $\tau_{1}$ the relaxation time.
This equation, was used to predict the steady-state force
($F_{ss}=K_{0}F_{max}$) and settling time. The settling time was set to the
time needed to be within 5% of final value (i.e.
$0.05K_{0}=K_{1}\exp(-t_{s}/\tau_{1})$).

In addition, the bending actuators were characterized in terms of stroke using
the same setup. For these experiments (repeated thrice), the swing was removed
and the actuator pressurized. Afterwards, the captured images were used to
compute the stroke. Specifically, MATLAB was used to fit a circle along the
deformed structure using three points to compute the maximum curvature.

##  3 Results and Discussion

###  3.1 Characterization of Coiling Pattern

The effect of the height ($H$) on the coil radius ($R_{c}$) is shown in Fig.
3(a). The plotted values of $R_{c}$ are plotted for multiple values of screw
rotational speed ($\alpha V_{f}$ with $V_{f}$ the printhead speed) at the same
height. It can be observed that $R_{c}$ is linearly related to the height
($H$) in the range of 2.5 mm to 15 mm, with almost one order of magnitude of
increase. In contrast, the increase of $R_{c}$ seems to stall between 15 and
20 mm. This stalling implies that increasing $H$ will no longer increase
$R_{c}$ (and could even decrease)32. This behavior makes further increasing
the height not interesting for the InFoam method as it does not allow for
larger coils. In addition, it further cools down the thermoplastic, which can
lead to bonding issues making it even less appealing.

Figure 3: (a,b) Coil radius versus the height $H$ (a) and (b) rotational speed
of screw ($\alpha V_{F}$), and (c) The coil density ($N$-value) versus coil
density multiplier $(\alpha W(2\pi R_{c}))$ used to estimate $G$ (see Equation
3).

Besides $H$, the literature indicates that $R_{c}$ decreases when the
viscosity decreases32. Due to the shear-thinning behavior of the
thermoplastic, a decrease in viscosity (for constant temperature) should
therefore coincide with an increased extrusion speed. The extrusion speed is
equivalent (under the assumption of linearity) to an increased rotational
speed of the screw $\alpha V_{F}$, which is plotted against $R_{c}$ in Fig.
3(b). It can be observed that there is a nonlinear and downwards trend with
increasing extrusion speed. This behavior is in line with a power-law shear-
thinning relationship. To compare the shear-thinning behavior, the viscosity
of the SEBS material was first measured using a Rheocompass MC501 (Anton-Paar
GmbH, Germany) with the data shown in Fig. S2 of the Supplementary
Information. The data were then fitted to a power-law normalized by $R_{c}$
(using the exponent of the viscosity ($n-1=-0.09$), average coil radius
$\hat{R}_{c}$ per $H$, and rotational speed of the screw $\alpha V_{F}$):

$R_{c}=a(\alpha V_{F})^{n-1}\hat{R}_{c}$ (9)

Fitting the value of $a$ for all heights and averaging leads to the curve seen
in Fig. 3(b). The fits had errors of less than 10%. Thereby implying that
there is a correlation between the change in coiling radius $R_{c}$ and the
extrusion speed, which can explain the spread of coiling radii. This extrusion
speed-induced change provides an additional opportunity to fine-tune $R_{c}$
at a fixed $H$. However, increasing the extrusion speed will also affect other
properties such as the diameter (by increasing die swell).

The $N$-value was plotted against the coil density multiplier $\alpha W/(2\pi
R_{c})$ (see Equation 3) in Fig. 3(c) (using the mean value of $R_{c}$ per
$H$). It can be observed that there is indeed a linear relation for the
$N$-value. The linearity is also observed by fitting a line with $G=0.17$
mm/mrad (Equation 3). The accuracy thereof implies that the proposed
relationship is valid. In addition, when extruding in the air the value of $G$
was found to be in the range of 0.154-0.165 mm/mrad, which is in the same
range as the fitted value (see Supplementary Information). Lastly, it can be
observed that lower $H$ do not achieve high $N$-values. This discrepancy is
due to the coils stacking up too much. When this occurs the coiling pattern is
replaced with an accumulation of material31. This accumulation phase is more
chaotic and not used by the InFoam method.

Lastly, the coil height ($h_{c}$) estimation (Equation 4) was validated using
the coiling patterns to be on average within 4% (see Supplementary
Information).

###  3.2 Mechanical Programming of Soft Cubes

The measured and computed porosity of the soft cubes versus the height ($H$)
and the coil density ($N$-value) are shown in Fig. 4(a). It can be observed
that the porosity ranged from 45 to 89% with similar levels of porosities were
attained by changing $H$ and the $N$-value. Increases in the $N$-value/$H$
lead to a linear decreasing and nonlinear increase of the porosity,
respectively. The predicted porosity of Equation 5 (with the estimate of $G$
from the previous section) correlated quite well. Even though the accuracy for
the $N$-value is lower than the coil radius ($R_{c}$) estimate, the mean
absolute error was less than 4% for both ($H$ (2.6%) and the $N$-value 3.9%).
This graph demonstrates that the $N$-value can be used to tune the level of
porosity independent of $H$. This change in porosity is advantageous as it
means the InFoam method can adjust the $N$-value by setting the extrusion
multiplier ($\alpha$) accordingly (see Equation 3). By doing so the InFoam
method can use both large and small coiling radii to increase accuracy
(smaller) or printing speed (larger), while still achieving the desired
porosity.

Figure 4: (a) Porosity versus height ($H$) and coil density ($N$-value), (b,c)
stress-strain curve at 46% and 89% porosity, (d,e) segment modulus versus
porosity for different axes at maximum strain (d) and different levels of
strain (e), and (f) the energy dissipation versus porosity.

The results of the compression tests on the cubes were analyzed with the axes
defined as: in line with the coils $x$, perpendicular to the coils but in-
plane $y$, and layer-wise $z$. The stress-strain behavior of the soft cubes
are shown in Fig. 4(c) and (d) for the $x$\- and $z$-axis ($y$-axis is similar
but not shown for clarity) for both low (45.8%,(c)) and high (89%,(d))
porosity, respectively. Their general shape is similar with a clear hysteretic
behavior but with much lower magnitudes for the higher porosity cube. It can
also be noted that both have a linear region and a nonlinear region. The
nonlinearity occurs due to the densification of the porous structure (i.e.
when all air has been pushed out). This densification limits the effective
deformation one can practically achieve37, which is reached at an earlier
stage in the higher porosity than the lower porosity. The stress in $x,z$
seems linear up to 0.2 and 0.5 strain, for the high and low porosity cube,
respectively. The compression of three foam cubes with different porosity is
shown in Supplementary Movie 2 in the Supplementary Information.

The segment modulus at maximum strain (0.84) is plotted in Fig. 4(d) for all
three axes. The raw data is shown for the $x$ and $z$ axis (see Fig. S3 of the
Supplementary Information for the $y$). Error bars were not added for clarity
but the standard errors were on average within 13%. It can be observed that
the trend in both $N$-value and $H$ are similar. This observation implies that
the magnitude of the porosity is important and not the geometry of the
underlying coiling pattern (i.e. the specific $R_{c}$ and $N$-value
combination). In addition, it can be observed that all three axes have a
power-law relation with the porosity, as expected. The fitted line shows that
a (nearly) quadratic function describes all three directions reasonably well,
which is in line with the decay of the elastic modulus seen in normal foams34.
In addition, different moduli can be observed for the three axes, which
indicates that the layout of the coiling pattern will induce some anisotropy.
A possible method to reduce the anisotropy between $x$ and $y$ is to change
the coil trajectory by 90 degrees between layers. This approach leads to
similar stress in the $x$ and $y$-direction (see Fig. S4 in the Supplementary
Information).

The segment modulus was also evaluated at multiple levels of strain for the
$x$-axis (see Fig. 4(e)). It can be observed that before densification (low
strain) the moduli are more similar. After densification the slope at which
moduli changes increases significantly (from 11.25 to 19.09). Interestingly,
the modulus decreases faster but still in a very similar magnitude to normal
foams (i.e. around quadratically)34. This change allows for mechanical
programming to have parts that behave very differently under the same load. To
compare the magnitudes to the bulk material, the modulus was approximated
using the Shore hardness (47A) and the equation38 $E=0.486\exp{(0.0345\cdot
47)}=2.46$ MPa. It can be observed that even at a strain of 0.4 that the 45.8%
and 89% cubes were 3.26 and 246 times as small. Thereby indicating that the
InFoam method can reduce the modulus by over two orders of magnitude. This
value exceeds the 66 times already demonstrated in literature when printing
with LRC24.

Lastly, the dissipated energy (in percentage) versus porosity is plotted in
Fig. 4(f) with the standard deviation. It can be observed that the
$x$-direction (with 31-40%) dissipates the most energy. A power-law
relationship with the porosity can be observed for the $x$ and $y$ directions
but is nearly constant in $z$. We expect that is due to the buckling of coils
in both $x$ and $y$, whereas the coils are compressed/bent in the
$z$-direction. In addition, the effect of $H$ and the $N$-value are similar
(just as for the stiffness) providing more evidence that not the coil geometry
but porosity is the important value.

###  3.3 Graded Porosity for Soft Fluidic Actuators

The printed structures with porosity gradients are shown in Fig. 5 and
Supplementary Movie 3 (see Supplementary Information). It can be seen that the
porosity gradient can deform to the screw motion, which is shown without a
sleeve (Fig. 5(a)) and actuated (5(b)). In addition, the twisting (5(c)) and
contraction (5(d)) are shown in both actuated and unactuated form. Lastly,
Fig. 5(e) and (f) show that the absence (5(e)) and presence of a porosity
gradient with spacers (5(f)) can significantly increase the bending. In
general, screw, twisting, and two bending actuators all have similar external
geometry but using different porosity gradients allows for a wide variety of
deformation patterns. These examples show that the InFoam method can change
the deformation behavior significantly. This grading capability during AM
enables designers to design the deformation behavior purely by the porosity
gradient.

Figure 5: Deformation pattern for different porosity patterns. Figure 6: (a)
Maximum and steady state force, (b) settling time, and (c) curvature.

The maximum and steady-state force of the bending actuators (of Figure 5(f))
are plotted in Fig. 6(a) (and in Supplementary Movie 4 of the Supplementary
Information). For clarity of the image, the error bars are omitted but these
were lower than 10% for most (the data is provided with error bars in Fig. S5
in the Supplementary Information). It can be observed that an increase in
pressure leads to an increase in both steady-state and maximum force. The
maximum forces were in the range of 67.2 to 21.1 grams for actuators, which is
10.2 to 3.8 times their weight (6.6 and 5.5 grams, respectively). In addition,
a decrease in force for higher porosity actuators is observed. This behavior
is expected to be due to the lower stiffness. The lower stiffness would allow
the actuator to deform such that the output force is reduced. The steady-state
value of the force increases less with increasing pressure for lower porosity
samples. This trend is expected to be linked to densification, which decreases
the achievable deformation of lower porosity samples. Fitting a curve for the
maximum and steady-state force curve gives exponents that are similar to each
other. However, there is a discrepancy between the maximum and steady-state
values in terms of their exponent. This discrepancy is expected to be due to
the viscous component of the thermoplastic, which impacts the maximum force
only. Interestingly, the power-law fits for both maximum and steady-state
force correlate well with data, which indicates that porosity could be an
interesting metric to predict the mechanical performance of an SFA.

The settling time, i.e. when the force is within 5% of the steady-state value,
is plotted in Fig. 6(b). Although there is a significant spread, the trend
does indicate that the settling time decreases with increasing porosity. This
decrease is expected to be due to the combination of decreased density and
increased energy dissipation associated with higher porosity, as seen in the
soft cubes. Fits based on the average settling time for 20 and 75 kPa, follow
a power-law behavior with the porosity. Thereby supporting that the porosity
could be an interesting metric to predict the mechanical performance.

Lastly, the maximum curvature versus porosity for different pressures is shown
in Fig. 6(c). A nonlinear relation can be observed between the curvature and
porosity. This behavior is expected as the higher porosity actuators will
deform more under the same load. Similar to the force metrics discussed
previously, these curves can be described with a power-law relation. Fitting
the data shows that the curvature increases with the porosity with an exponent
around the cubic root of the porosity. Thereby indicating again that porosity
could be a useful tool for programming deformation behavior. In addition, a
nonlinear relation between the curvature and pressure can be observed.
Independent of the porosity the actuators seem to plateau at 60 and 75 kPa.
This plateau is expected to be due to the densification of the porous
structure. The densification means that exponentially more pressure would be
required to increase curvature. The higher densification strain of the higher
porosity structures allows them to deform more before reaching this point.

##  4 Conclusion

Soft fluidic actuators are a popular choice in soft robotics due to their
versatility, accessibility, and (relative) simplicity. To further expand their
capabilities, fabrication methods for structures that simultaneously allow for
the transfer of fluidic power and have a stiffness gradient are required. Our
InFoam method, as proposed in this work, can provide this capability by
allowing for the direct printing of porous structures. The printed porous
structures allow for fluid transport and mechanical flexibility. In addition,
the structures used a hyperelastic thermoplastic to further reduce the
stiffness due to the usage of a hyperelastic thermoplastic.

The InFoam method is capable of mechanical programming by a single level or
gradient of porosity (by combining it with normal plotting). To achieve this
capability it only requires the characterization of the coil radius (with
respect to the height) and the length of extruded material per unit of
movement ($G$) at the printing temperature. Subsequently, the InFoam method
can use the developed equations to adjust the porosity of a structure by
setting the height and/or extrusion speed while taking the coiling pattern’s
width and height into account.

With a single level of porosity, the InFoam method can program the density,
stiffness, and energy dissipation. Large changes could be realized such as an
89% lower density and scaling the modulus by a factor of 246. In addition, the
energy dissipation in $x$ and $y$ could be programmed without affecting the
$z$-direction. All of these properties were shown to correlate well through a
power-law relationship with the porosity.

By printing a porosity gradient, the InFoam method allowed for a wide range of
motion (including twisting, contraction, and bending). In addition, the
grading was shown to be capable of programming the behavior of soft bending
actuators. Specifically, it was shown that the porosity could be used to
program both the output force (magnitude and settling time) and stroke.

Thus, the InFoam method can print graded porous structures that are
mechanically programmed by porosity. The InFoam method allows this capability
by solely adjusting the height of the nozzle and the ratio of extrusion and
movement speed. An interesting observation of the printed porous structures is
that the mechanical behavior of both the cubes and bending actuators could be
correlated through a power-law relationship with the porosity. In addition,
they showed similar trends in their changes in behavior. Further investigation
of this correlation could enable tuning of the mechanical performance of soft
(porous) actuators using a simple compression experiment. In addition, the
InFoam method’s stiffness grading can be further investigated. A structure
with high porosity has low stiffness and allows fluid transport, which can be
interesting for applications such as scaffolds for more complex SFAs and
support. Whereas graded porosity can be a tool for soft fluidic
actuators/sensors with complex deformation patterns.

## References

  * 1 B. Mosadegh, P. Polygerinos, C. Keplinger, S. Wennstedt, R. F. Shepherd, U. Gupta, J. Shim, K. Bertoldi, C. J. Walsh, and G. M. Whitesides, “Pneumatic networks for soft robotics that actuate rapidly,” _Advanced functional materials_ , vol. 24, no. 15, pp. 2163–2170, 2014.
  * 2 M. Wehner, R. L. Truby, D. J. Fitzgerald, B. Mosadegh, G. M. Whitesides, J. A. Lewis, and R. J. Wood, “An integrated design and fabrication strategy for entirely soft, autonomous robots,” _nature_ , vol. 536, no. 7617, pp. 451–455, 2016.
  * 3 N. Cheng, J. H. Low, B. W. Ang, A. J. Goh, and C.-H. Yeow, “Soft fabric-based pneumatic sensor for bending angles and contact force detection,” _IEEE Sensors Journal_ , vol. 19, no. 4, pp. 1269–1279, 2018.
  * 4 S. T. Mahon, A. Buchoux, M. E. Sayed, L. Teng, and A. A. Stokes, “Soft robots for extreme environments: Removing electronic control,” in _2019 2nd IEEE international conference on soft robotics (RoboSoft)_. IEEE, 2019, pp. 782–787.
  * 5 F. Connolly, P. Polygerinos, C. J. Walsh, and K. Bertoldi, “Mechanical programming of soft actuators by varying fiber angle,” _Soft Robotics_ , vol. 2, no. 1, pp. 26–32, 2015.
  * 6 R. V. Martinez, C. R. Fish, X. Chen, and G. M. Whitesides, “Elastomeric origami: programmable paper-elastomer composites as pneumatic actuators,” _Advanced functional materials_ , vol. 22, no. 7, pp. 1376–1384, 2012.
  * 7 L. Paez, G. Agarwal, and J. Paik, “Design and analysis of a soft pneumatic actuator with origami shell reinforcement,” _Soft Robotics_ , vol. 3, no. 3, pp. 109–119, 2016.
  * 8 M. A. Dias, M. P. McCarron, D. Rayneau-Kirkhope, P. Z. Hanakata, D. K. Campbell, H. S. Park, and D. P. Holmes, “Kirigami actuators,” _Soft matter_ , vol. 13, no. 48, pp. 9087–9092, 2017.
  * 9 M. A. Robertson and J. Paik, “New soft robots really suck: Vacuum-powered systems empower diverse capabilities,” _Science Robotics_ , vol. 2, no. 9, p. eaan6357, 2017.
  * 10 S. P. Murali Babu, F. Visentin, A. Sadeghi, A. Mondini, F. Meder, and B. Mazzolai, “Sensorized foam actuator with intrinsic proprioception and tunable stiffness behavior for soft robots,” _Advanced Intelligent Systems_ , p. 2100022, 2021.
  * 11 Y. Yamada and T. Nakamura, “Laminated foam-based soft actuator for actuatable flexible structure,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 4359–4364.
  * 12 B. C. Mac Murray, X. An, S. S. Robinson, I. M. van Meerbeek, K. W. O’Brien, H. Zhao, and R. F. Shepherd, “Poroelastic foams for simple fabrication of complex soft robots,” _Advanced Materials_ , vol. 27, no. 41, pp. 6334–6340, 2015.
  * 13 C. C. Futran, S. Ceron, B. C. Mac Murray, R. F. Shepherd, and K. H. Petersen, “Leveraging fluid resistance in soft robots,” in _2018 IEEE International Conference on Soft Robotics (RoboSoft)_. IEEE, 2018, pp. 473–478.
  * 14 Y. Yamada and T. Nakamura, “Actuatable flexible large structure using a laminated foam-based soft actuator,” in _2020 IEEE/SICE International Symposium on System Integration (SII)_. IEEE, 2020, pp. 74–79.
  * 15 S. Bilent, T. H. N. Dinh, E. Martincic, and P.-Y. Joubert, “Influence of the porosity of polymer foams on the performances of capacitive flexible pressure sensors,” _Sensors_ , vol. 19, no. 9, p. 1968, 2019.
  * 16 S. Bilent, E. Martincic, and P.-Y. Joubert, “Fabrication of porous polymer based pressure sensors using sacrificial pva 3d scaffolds for the refined tuning of sensor performances,” in _2020 Symposium on Design, Test, Integration & Packaging of MEMS and MOEMS (DTIP)_. IEEE, 2020, pp. 1–4.
  * 17 J. Morrow, S. Hemleben, and Y. Menguc, “Directly fabricating soft robotic actuators with an open-source 3-d printer,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 1, pp. 277–281, 2016.
  * 18 H. K. Yap, H. Y. Ng, and C.-H. Yeow, “High-force soft printable pneumatics for soft robotic applications,” _Soft Robotics_ , vol. 3, no. 3, pp. 144–158, 2016.
  * 19 S. Walker, U. Daalkhaijav, D. Thrush, C. Branyan, O. D. Yirmibesoglu, G. Olson, and Y. Menguc, “Zero-support 3d printing of thermoset silicone via simultaneous control of both reaction kinetics and transient rheology,” _3D Printing and Additive Manufacturing_ , vol. 6, no. 3, pp. 139–147, 2019\.
  * 20 M. A. H. Khondoker and D. Sameoto, “Direct coupling of fixed screw extruders using flexible heated hoses for fdm printing of extremely soft thermoplastic elastomers,” _Progress in Additive Manufacturing_ , vol. 4, no. 3, pp. 197–209, 2019.
  * 21 A. Georgopoulou, L. Egloff, B. Vanderborght, and F. Clemens, “A sensorized soft pneumatic actuator fabricated with extrusion-based additive manufacturing,” in _Actuators_ , vol. 10, no. 5. Multidisciplinary Digital Publishing Institute, 2021, p. 102\.
  * 22 O. D. Yirmibesoglu, L. E. Simonsen, R. Manson, J. Davidson, K. Healy, Y. Menguc, and T. Wallin, “Multi-material direct ink writing of photocurable elastomeric foams,” _Communications Materials_ , vol. 2, no. 1, pp. 1–14, 2021.
  * 23 N. M. Ribe, “Coiling of viscous jets,” _Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences_ , vol. 460, no. 2051, pp. 3223–3239, 2004.
  * 24 J. I. Lipton and H. Lipson, “3d printing variable stiffness foams using viscous thread instability,” _Scientific reports_ , vol. 6, no. 1, pp. 1–6, 2016.
  * 25 P.-T. Brun, C. Inamura, D. Lizardo, G. Franchin, M. Stern, P. Houk, and N. Oxman, “The molten glass sewing machine,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 375, no. 2093, p. 20160156, 2017.
  * 26 X. Tian, J. Plott, H. Wang, B. Zhu, and A. J. Shih, “Silicone foam additive manufacturing by liquid rope coiling,” _Procedia Cirp_ , vol. 65, pp. 196–201, 2017.
  * 27 S. Zou, D. Therriault, and F. P. Gosselin, “Spiderweb-inspired, transparent, impact-absorbing composite,” _Cell Reports Physical Science_ , vol. 1, no. 11, p. 100240, 2020.
  * 28 B. Emery, D. Revier, and J. Lipton, “Applied viscous thread instability for manufacturing 3d printed foams,” in _2021 International Solid Freeform Fabrication Symposium_. University of Texas at Austin, 2021.
  * 29 M. Saari, M. Galla, B. Cox, P. Krueger, A. Cohen, and E. Richer, “Additive manufacturing of soft and composite parts from thermoplastic elastomers,” in _2015 International Solid Freeform Fabrication Symposium_. University of Texas at Austin, 2015.
  * 30 S. Chiu-Webster and J. Lister, “The fall of a viscous thread onto a moving surface: a’fluid-mechanical sewing machine’,” _Journal of Fluid Mechanics_ , vol. 569, p. 89, 2006.
  * 31 H. Yuk and X. Zhao, “A new 3d printing strategy by harnessing deformation, instability, and fracture of viscoelastic inks,” _Advanced Materials_ , vol. 30, no. 6, p. 1704028, 2018.
  * 32 N. M. Ribe, M. Habibi, and D. Bonn, “Liquid rope coiling,” _Annual review of fluid mechanics_ , vol. 44, pp. 249–266, 2012.
  * 33 R. J. Crawford and P. Martin, _Plastics engineering_. Butterworth-Heinemann, 2020.
  * 34 I. Gibson and M. F. Ashby, “The mechanics of three-dimensional cellular materials,” _Proceedings of the royal society of London. A. Mathematical and physical sciences_ , vol. 382, no. 1782, pp. 43–59, 1982.
  * 35 J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid _et al._ , “Fiji: an open-source platform for biological-image analysis,” _Nature methods_ , vol. 9, no. 7, pp. 676–682, 2012.
  * 36 C. W. Macosko, _Rheology Principles, Measurements, and Applications_. Wiley-VCH, 1994.
  * 37 B. C. Mac Murray, C. C. Futran, J. Lee, K. W. O’Brien, A. A. Amiri Moghadam, B. Mosadegh, M. N. Silberstein, J. K. Min, and R. F. Shepherd, “Compliant buckled foam actuators and application in patient-specific direct cardiac compression,” _Soft robotics_ , vol. 5, no. 1, pp. 99–108, 2018.
  * 38 K. Larson, “Can You Estimate Modulus From Durometer Hardness for Silicones?” https://www.dow.com/content/dam/dcc/documents/en-us/tech-art/11/11-37/11-3716-01-durometer-hardness-for-silicones.pdf, 2016, [Accessed 3rd of March 2022].