SIMULATION AND DESIGN OF SOFT ROBOTIC SWIMMERS WITH ARTIFICIAL

MUSCLE

By

Andrew Michael Hess

A DISSERTATION

Michigan State University

in partial fulﬁllment of the requirements

Submitted to

for the degree of

Mechanical Engineering – Doctor of Philosophy

Computational Mathematics, Science, and Engineering – Dual Major

2019

ABSTRACT

SIMULATION AND DESIGN OF SOFT ROBOTIC SWIMMERS WITH ARTIFICIAL

MUSCLE

By

Andrew Michael Hess

Soft robots take advantage of rich nonlinear dynamics and large degrees of freedom to perform
actions often by novel means beyond the capability of conventional rigid robots. Nevertheless, there
are considerable challenges in analysis, design, and optimization of soft robots due to their complex
behaviors. This is especially true for soft robotic swimmers whose dynamics are determined by
highly nonlinear ﬂuid-structure interactions. We present a holistic computational framework that
employs a multi-objective evolutionary method to optimize feedback controllers for maneuvers of
a soft robotic ﬁsh under artiﬁcial muscle actuation.

The resultant ﬂuid-structure interactions are fully solved by using a novel ﬁctitious domain/active
strain method, developed to numerically study the swimming motion of thin, light-weight soft
robots composed of smart materials that can actively undergo reversible large deformations. We
assume the elastic material to be neo-Hookean, and behave like an artiﬁcial “muscle” which, when
stimulated, generates a principal stretch of contraction. Instead of imposing active stresses, here we
adopt an active strain approach to impose contracting strains to drive elastic deformation following
a multiplicative decomposition of the deformation gradient tensor. The hydrodynamic coupling
between the ﬂuid and the solid is then resolved by using the ﬁctitious domain method where the
induced ﬂow ﬁeld is virtually extended into the solid domain. Pseudo body forces are employed
to enforce the interior ﬁctitious ﬂuid motion to be the same as the structural movement. To
demonstrate this method we carry out a series of numerical explorations for soft robotic swimmers
of both 2D and 3D geometries that prove these robot prototypes can eﬀectively perform undulatory
and jet-propulsion locomotion when active contracting strains are appropriately distributed on
elastica.

To complete our framework we then demonstrate the design and optimization of a maneuverable,
undulatory, soft robotic ﬁsh.
In particular, we consider a two-dimensional elastic plate with
ﬁnite thickness, subjected to active contractile strains on both sides of the body. Compared
to the conventional approaches that require specifying the entire-body curvature variation, we
demonstrate that imposing contractile active strains locally can produce various swimming gaits,
such as forward swimming and turning, using far fewer control parameters. The parameters of
a pair of proportional-integral-derivative (PID) controllers, used to control the amplitude and the
bias of the active strains, respectively, are optimized for tracking a moving target involving diﬀerent
trajectories and Reynolds numbers, with three objectives, tracking error, cost of transport, and
elastic strain energy. The resulting Pareto fronts of the multi-objective optimization problem reveal
the correlation and trade-oﬀ among the objectives and oﬀer key insight into the design and control
of soft swimmers.

Furthermore, we expand upon these core studies by investigating a variety processes and
phenomena that are of potential importance to the design of soft robotic ﬁsh. Rich and complex
phenomena occur at the interface between ﬂuids and elastic solids that can inﬂuence the behavior
of the overall system. This is especially true in the case of soft robots where wrinkling, folding,
and similar wavelike behavior can occur at the surface and aﬀect the momentum transfer which is
necessary to produce locomotion. To begin to understand this behavior we analytically investigate
the linear stability of a viscoelastic gel-like ﬁlm subjected to Newtonian Couette ﬂow in the limit
of vanishing Reynolds number.

The ﬁnal avenue of study attempts to answer the question of how to design and optimize the
body shape in conjunction with the actuation scheme of a soft robotic ﬁsh. The ﬂexibility of soft
materials and actuators lead to large degrees of freedom creating a wide array of possibilities but
at the same time overwhelming the traditional design process. We demonstrate a viable solution
to this issue by using evolutionary optimization to simultaneously optimize both the shape and
actuation of a simple 2D soft robotic ﬁsh in a simple case study.

Copyright by
ANDREW MICHAEL HESS
2019

Dedicated to my parents, whose support and patience made all of this possible.

v

ACKNOWLEDGEMENTS

First and foremost I would like to thank my advisor, Dr. Tong Gao, for his expertise and guidance
throughout the projects that make up this dissertation. As his ﬁrst doctoral student, this process
was a learning experience for us both and I believe one that left us both the better for it. I hope that
the feeling is mutual and that our continued collaborations will bear fruit for years to come.

I would like to thank my committee members, Dr. Farhad Jaberi, Dr. Ahmed Naguib, and Dr.
Mohsen Zayernouri for their support, suggestions, and imparted wisdom from the beginning as I
sat in their classrooms to the very end as I stood before them to defend my work.

I gratefully acknowledge my collaborators such as Dr. Xiaobo Tan and Dr. Zhaowu Lin among
innumerable others who have and continue to share their knowledge and expertise without which
none of this would have been possible.

Last but not least, I thank those who supported me through thick and thin. Mom & Dad and
my grandparents who were always there for me regardless of how rocky the path, at times, became.
Brittany, my longtime partner in crime, who stuck with me from undergrad through all the trials
and tribulations that lead to us each becoming some variety of “Dr”. And friends new and old
whose companionship kept me sane along the way, Krissy, Coco, Mattei, Seward, Kevin, Chris,
Rachel, Joe, Chad, Taylor, Connor, Sanders, ...

vi

TABLE OF CONTENTS

.
.

.
.

.
.

.
.

.
.

.
.

.
.

ix

x

1

2.1 Mathematical Model

LIST OF TABLES .
LIST OF FIGURES .

.
.
INTRODUCTION .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
CHAPTER 1
CHAPTER 2 STABILITY OF COUETTE FLOW PAST A GEL FILM . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Constitutive equations
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Universal Eulerian description . . . . . . . . . . . . . . . . . . . . . . . .
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
5
5
7
9
2.2.1 Base state .
9
2.2.2 Normal modes and interfacial conditions . . . . . . . . . . . . . . . . . . . 10
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 12
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.1 Growth Rate .
.
2.3.2 Mechanism . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.3 Marginal Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
. 16

2.2 Linear Stability Analysis
. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.4 Concluding Remarks

2.3 Results .

.

.
.
.

.

.

.

. .

. .

. .

.

.

CHAPTER 3 A FLUID–STRUCTURE INTERACTION STUDY OF SOFT ROBOTIC

SWIMMER USING A FICTITIOUS DOMAIN/ACTIVE-STRAIN METHOD 19
3.1 Mathematical Model and Numerical Method . . . . . . . . . . . . . . . . . . . . . 20
3.1.1 Hill’s Muscle Model
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.2 Active Strain Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.3 Distributed Lagrange Multiplier/Fictitious Domain Method (DLM/FD) . . . 23
3.1.3.1
Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.3.2 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.3.3 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2.1 Two-Dimensional Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.2 Three-Dimensional Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.2.1 Deformation Control . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.2.2 Undulatory Swimming . . . . . . . . . . . . . . . . . . . . . . . 38
Jet-propulsion Swimming . . . . . . . . . . . . . . . . . . . . . 40
3.2.2.3
. 44

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2 Results and Discussion .

3.3 Concluding Remarks

.

CHAPTER 4 CFD-BASED MULTI-OBJECTIVE CONTROLLER OPTIMIZATION

FOR SOFT ROBOTIC FISH WITH MUSCLE-LIKE ACTUATION . . . . . 47
4.1 Mathematical Model
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2 Controller Design and Optimization . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.2.1 Control Scheme .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.2.2 Multi-Objective Optimization . . . . . . . . . . . . . . . . . . . . . . . . 54

.

.

vii

4.3 Results and Discussion .

.

4.3.1 Undulatory Swimming . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Optimization and Analysis . . . . . . . . . . . . . . . . . . . . . . . . .
Trajectory Comparison . . . . . . . . . . . . . . . . . . . . . .
Pareto Fronts of Optimized Solutions . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
. 56
. 58
. 58
4.3.2.1
4.3.2.2
. 60
4.3.2.3 Objective Correlations and Parameter Trends . . . . . . . . . . . 62
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.1 Body Shape Optimization for Muscle Driven Undulatory Swimming . . . . . . . . 67
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
. 72
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 6 CONCLUSIONS .
.
BIBLIOGRAPHY .

4.4 Concluding Remarks

. .
CHAPTER 5 FUTURE WORK .

.

.

.

.

.

.

.

viii

LIST OF TABLES

Table 3.1: The amplitude of the central point on the free end and the Strouhal number

and the average drag coeﬃcient for a ﬂexible beam . . . . . . . . . . . . . . . . 29

Table 3.2: The amplitude of the central point on the plate trailing edge and the Strouhal

number and the average drag coeﬃcient for a ﬂag ﬂapping in uniform ﬂow.

. . . 30

Table 4.1: Summary of the variables (PID tuning parameters) and objectives considered

in the optimization. . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

ix

LIST OF FIGURES

Figure 2.1: Schematic for a Couette ﬂow past a gel ﬁlm. Nonlinear solid models: (a)

neo-Hookean, (b) Kelvin-Voigt, (c) Zener.

. . . . . . . . . . . . . . . . . . . .

5

Figure 2.2: The growth rate Re(α) as a function of k for the neo-Hookean solid: (a)
H = 0.4,T = 0; (b) H = 4.0,T = 0; (c) H = 0.4,T = 10; (d) H = 4.0,T = 10.
The open circles are the results obtained by Gkanis & Kumar [33].

. . . . . . . 12

Figure 2.3: Comparison of the growth rate Re(α) as a function of k for the neo-Hookean,
Kelvin-Voigt, and Zener solids: (a) H = 0.4,T = 0, G = 5; (b) H = 4.0,T =
0, G = 5; (c) H = 0.4,T = 10, G = 11; (d) H = 4.0,T = 10, G = 5. . . . . . . . . 13

Figure 2.4:

(a) Schematic for a sinusoidal surface perturbation. The contributions of
the ﬁrst normal-stress diﬀerence and the solid viscous eﬀect are marked by
the black and red solid arrows, respectively. The open arrows represent the
stabilizing eﬀect due to surface tension. (b,c): Critical values of the imposed
strain Gc and the critical wavenumber kc as a function of H. The black dotted
lines are drawn at the locations where either Gc or kc tends to inﬁnity. Inset
in (c): Growth-rate curves near H ≈ 1 where short-wave instabilities recur.

. . 14

Figure 3.1: Schematic of the 1D three element Hill’s muscle model. It is composed of
a contractile element (CE) in series with a nonlinear spring (SE) which are
both in parallel with another nonlinear spring (PE).

. . . . . . . . . . . . . . . 21

Figure 3.2: Time evolution of the displacements in both x− and y−directions of the
central point on the plate trailing edge (point A, see inset of (a)) of a thin
ﬂexible plate in uniform ﬂow.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Figure 3.3: Time evolution of the transverse displacement of the central point on the plate

trailing edge of a thin ﬂexible plate in uniform ﬂow.

. . . . . . . . . . . . . . . 30

Figure 3.4:

(a) Schematic of a 2D beam under contracting actuation. The active segment
is highlighted in orange. The microscopic mechanical model of the active
contraction is illustrated by the zoom-in schematic on bottom. (b) Either a
constant (i.e., a step-like) or a sinusoidal contraction is applied during the ﬁrst
half period 0 ≤ t < T/2 on the two surfaces alternatively. Then the actuation
is turned oﬀ for the second half T/2 ≤ t < T. (c) Convergence study for
the 2D beam under a sinusoidal actuation by varying domain size, mesh
resolution, and time step. The other parameters are Re = 500, G = 1000,
δ = 0.4, α = 0.15, and T = 2.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 31

x

Figure 3.5: Swimming motions of a 2D elastic beam under the constant and sinusoidal
actuation when ﬁxing Re = 500 and T = 2.0. (a) Instantaneous velocity as a
function time. (b) Strain energy E as a function of time. Inset in (b): Envelope
trajectories of the plate midline during one actuation period. Snapshots of
the corresponding vorticity ﬁelds are shown in panels (c) (constant) and (d)
(sinusoidal).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

. .

. .

.

Figure 3.6: Comparisons of the free swimming motion of a 2D elastic beam at Re =
100 and Re = 500 when ﬁxing T = 1.0 for sinusoidal actuation following
Eq. (3.20). (a) Instantaneous velocities as a function time. (b) Mean elastic
energy densities e as a function of time. Inset in (b): Envelope trajectories of
the plate midline during one actuation period. Snapshots of the corresponding
vorticity ﬁelds are shown in panels (c) at Re = 100 and (d) at Re = 500,
respectively.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

. .

. .

Figure 3.7: Measurements of the averaged swimming speed (cid:104)U(cid:105) in (a-c), and power P
(solid blue lines) and eﬃciency η (dashed blue lines) in (d-f) of an swimming
beam under a sinusoidal actuation following Eq. (3.20) when ﬁxing Re = 500,
α = 0.15, and varying T, G and δ, respectively. . . . . . . . . . . . . . . . . . . 34

Figure 3.8:

(a) Examples of diﬀerent actuation patterns (horizontal vs. diagonal) on the
top surface of a square plate, and the corresponding bending deformations in
(b) and (c). .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

.

.

.

.

.

Figure 3.9: Snapshots of circular plate under (a-c) radial and (d-e) azimuthal actuation. . . . 36

Figure 3.10: Snapshots of an oblate bell under (a-c) radial and (d-e) azimuthal actuation.

. . 37

Figure 3.11: Snapshots of 3D vortex-core structures (a) and the projection of 2D vortices
on the x − y plane (b) of an elastic plate that performs undulation under
a sinusoidal actuation following Eq. (3.20). The computation parameters
are chosen as α = 0.2, δ = 0.5, T = 1.0, G = 1000, and Re = 1000.
(c) Instantaneous swimming velocity as a function of time. The color bar
represents the magnitude of the 3D vorticity, i.e., |ω f | (b) Drag coeﬃcient
(red line) and wake-induced force (i.e., Lamb vector)
2 (cid:104)U(cid:105)2 A
(blue line) as functions of time. (e) Strain
energy E as a function of time.

(cid:17)
Cd = Fy/(cid:16) 1
coeﬃcient Cw = −FL/(cid:16) 1

2 (cid:104)U(cid:105)2 A

(cid:17)

. . . . . . . . . . . . . . . . . . . . . . . . . . 38

Figure 3.12: Simulation results of an undulating swimming motion by varying δ while
(a) Instantaneous
ﬁxing α = 0.2, T = 1.0, G = 1000 and Re = 1000.
(b) Time-averaged
swimming speed U at δ = 0.3,0.5,0.7, respectively.
swimming speed (cid:104)U(cid:105) and strain energy (cid:104)E(cid:105) as a function of δ for periodic
locomotion. (c) Power P and eﬃciency η as a function of δ. . . . . . . . . . . . 39

xi

Figure 3.13: Typical vortex structures of an elastic bell undergoing jellyﬁsh-like swimming
in the moving coordinates, when choosing α = 0.8, T = 2.0, G = 100,
Re = 500, ε = 0.1, δ = 0.21, and h = 0.03. (a) Vortex core overlaid on the
solid mesh. The color bar represents the magnitude of the 3D vorticity, i.e.,
|ω f |. Inset: Schematic of the swimming elastic bell with a thickened top and
azimuthal actuation. (b-i) In-plane ﬂuid velocity vector ﬁeld (u, v) overlaid
on the colormap of vorticity on a 2D cross section plane during one period. . . . 41

Figure 3.14: Highlights of bell’s rowing motion accompanying whole-body vibrations
through the kinematics of the bell arm during the contracting and reﬁlling
process (a). The instantaneous strain energy (b) and COM velocity(c) as
function of time.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

. .

Figure 3.15: (a) Phase diagram of (α,Re) for a swimming bell undergoing jet-propulsion
locomotion. The while and black symbols represent the forward and backward
locomotion, respectively. The background color represents the time-averaged
steady-state swimming velocity (cid:104)U(cid:105) at late times. (b) Typical instantaneous
COM velocities when choosing α(cid:48)s in diﬀerent regimes at Re = 250.

. . . . . . 42

Figure 3.16: Comparison of the dynamics of soft robotic baby-jellyﬁsh by distributing
active contractions (marked by the blue lines) on (a) arms (i.e., radially)
and partial of the circular disk (i.e., azimuthally), and (b) arms only. The
four sequential snapshots in (a) highlight a stable directional locomotion in
one actuation period. The snapshot in (b) shows an example of unstable
asymmetric deformations without applying the azimuthal contraction in the
disk to stabilize the structure.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 45

Figure 4.1:

(a) Schematic of a 2D beam under contractile actuation (view from above).
The active segment is of length δ. The microscopic mechanical model of
the active contraction is illustrated by the zoom-in schematic on bottom. (b)
Contractile actuation applied alternatively on the surfaces of two sides, which
exponentially decays in the thickness (y) direction (not shown in the ﬁgure).

. . 49

Figure 4.2: Schematic illustrating the control scheme based on the translational and rota-

tional motion of a soft swimmer in relation to a moving target xt.

. . . . . . . . 52

Figure 4.3: Free-swimming of a 2D elastic beam at Re = 500 and T = 2.0. (a) Instan-
taneous snapshot of vorticity ﬁeld.
(b) Envelope trajectories of the beam
midline during one actuation period: (top) numerical results; (bottom) ﬁt-
ted results by Eq. (4.14). (c) Mean values of U, E, and Cot as function of
contraction strength α. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Figure 4.4: Simulation results for the soft swimmer at α0 = 0.2, Re = 500, and G = 1000:
(a) Center of mass (CoM) paths at diﬀerent values of β; (b) Turning curvature
R−1 as a function of β.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

xii

Figure 4.5: CoM trajectories of a soft swimmer tracking a moving target that is moving
along both the (a) straight and (b) sinusoidal paths (black lines), when the
swimmer is subjected to the PID control. The (red) and (blue) lines correspond
to the controllers that have the lowest and a moderate error ¯e (but lowest elastic
energy ¯E), respectively. The black arrows mark the target’s moving directions.

. 60

Figure 4.6: Pareto fronts in the space of tracking error ( ¯e), cost of transport (Cot), and
elastic strain energy ( ¯E) for various target speeds v with Re = 500. Blue-
square: v = 0.1. Red-triangle: v = 0.2. Green-circle: v = 0.3. Straight
path: (a-d). Sinusoidal path: (e-h). The leftmost panels (a) and (e) are 3D
views of the objective spaces. The panels to their right are 2D projections
of the respective distribution that visualize the trends between each pair of
objectives. (b) and (f): ¯e vs Cot. (c) and (g): ¯e vs ¯E. (d) and (h): Cot vs ¯E. . . . 62

Figure 4.7: Pareto fronts in the space of tracking error ( ¯e), cost of transport (Cot), and
elastic strain energy ( ¯E) for various Re with v = 0.2. Blue-square: Re = 500.
Red-triangle: Re = 1000. Green-circle: Re = 1500. Yellow-diamond: Re =
2000. Straight path: (a-d). Sinusoidal path: (e-h). The leftmost panels (a)
and (e) are 3D views of the objective spaces. The panels to their right are
2D projections of the respective distribution that visualize the trends between
each pair of objectives. (b) and (f): ¯e vs Cot. (c) and (g): ¯e vs ¯E. (d) and (h):
Cot vs ¯E. .

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Figure 4.8: Comparison of the Pearson correlation coeﬃcients (ρ) of the objectives for
varying v (left) and varying Re (right). Blue-square = ¯e to Cot. Red-triangle
= ¯e to ¯E. Yellow-circle = Cot to ¯E.

. . . . . . . . . . . . . . . . . . . . . . . . 64

Figure 4.9: The trends of kp,position vs.

(c-d) Cases of
varying Re.
(Left) panels are from the straight path and (right) from the
sinusoidal path. The black dashed lines delineate the borders between the
meaningful and trivial solutions.

. . . . . . . . . . . . . . . . . . . . . . . .

(a-b) Cases of varying v.

¯e.

. 65

Figure 5.1: An example of a randomly generated swimmer body. The blue line is the
body surface and the dots are interpolation points for the Catmull-Rom spline.
Muscle actuation is limited to within the active region. . . . . . . . . . . . . . . 68

Figure 5.2: Sample objective space population distribution for generation 100 of a vari-

able body shape optimization. . . . . . . . . . . . . . . . . . . . . . . . . . .

. 68

Figure 5.3: Sample convergence history of the extreme (most optimal) values for each

objective in a variable body shape optimization.

. . . . . . . . . . . . . . . .

. 69

Figure 5.4: Sample conﬁgurations over progressive generations: (a) fastest swimming

speed, (b) highest swimming eﬃciency, and (c) lowest average elastic energy. . . 70

xiii

CHAPTER 1

INTRODUCTION

There has been tremendous interest in the development and design of soft robots in recent years,
owing to rapid advancement of technology in soft actuators, sensing, and additive manufacturing.
Compared to the conventional “rigid” robots, soft robots can take advantage of a large number of
degrees of freedom to achieve versatile locomotion and dexterous manipulation as many biological
organisms do. These abilities open up the possibility of soft robots that are smaller, lighter, and more
eﬃcient than their traditional counterparts [4, 60, 84, 23, 65]. While many promising examples
of soft robots have appeared, including soft swimmers emulating rays [74], juvenile jellyﬁsh [70],
and eel larvae [15], most soft robots use pneumatic actuation or cable-driven actuation [84], which
presents challenges in terms of system complexity, noise, and footprint, among other concerns.
Liquid crystal elastomers (LCE), dielectric elastomer (DEA) [15], and shape memory alloy [32],
oﬀer the promise of compact, muscle-like actuation for soft robots. These so-called smart materials
have attracted growing attention because of the intriguing shape or volume recovery properties under
diﬀerent external stimuli, such as electric, magnetic, thermal, chemical, light and pneumatic [83, 7,
50, 72, 89, 96]. However, understanding of the interplay between soft body and such muscle-like
actuation remains limited. In particular, a systematic design tool for soft robotic systems is lacking
that accounts for model dynamics and multiple objectives of interest.

Designing soft swimming robots that undergo active deformations in a ﬂuid is considerably
challenging. First of all, fast swimming motions are typically resultant from a signiﬁcant amount of
momentum exchange between the ﬂuid and solid structures to overcome viscous drag force in the
ﬂuid, which require robots to generate rapid and stable structural deformations reversibly. Mean-
while, eﬃcient locomotion of a deformable object requires the employment of speciﬁc swimming
patterns (or swimming gaits) to take advantage of thrust forces from the resultant ﬂuid drag and
wake structures behind [57], which is critical especially in the small or ﬁnite Reynolds number
regime where the viscous eﬀect is important [79]. To take all the factors into account, the dynamical

1

performances of soft robots with various geometries, material properties, as well as the imposed
active control schemes, need to be determined jointly with the induced ﬂuid motions.
In addi-
tion, dynamic instabilities may occur when light-weight structures move in ﬂuid [81, 24], which
adds complexities in understanding the control mechanisms of how soft robots perform a stable
swimming motion. In general, while various diﬀerent types of soft robots have been manufactured
and tested, it is desired to understand their precise swimming mechanisms, which requires the
combination of experimental studies with accurate modeling and simulations in design, analysis,
and optimization.

In Chapter 4 we present a novel study on multi-objective optimization of feedback controllers
for soft robots based on high-ﬁdelity CFD simulation that explicitly accounts for muscle-like
biomimetic actuation. Oftentimes, instead of fully resolving such nonlinear ﬂuid/elastic-structure
interactions (FESI), people oversimplify the problem by solving a reduced ﬂuid-solid coupling sys-
tem in mathematical modeling and numerical simulations. For example, the majority of numerical
studies of ﬁsh-like swimmers assume that the undulatory motion (or swimming gait) approximately
follows a traveling-wave form which can be experimentally measured for a free-swimming biolog-
ical ﬁsh [32, 48, 99]. Typically, one treats the ﬁsh body as deformable but prescribes the curvature
of its backbone deformation to a certain traveling-wave form [25, 93]. While convenient, these
studies may oversimplify the nonlinear FESI to be a (quasi) one-way coupling problem, and, in
the meantime, introduce too many free parameters to deﬁne the body actuation, which adds signif-
icant complexities in designing control schemes later. More importantly, soft swimmers will not
always follow a speciﬁc swimming gait, which only approximate the solitary swimming motion at
a quasi-steady state. In reality, they can constantly switch between straight-swimming and turning
motions with acceleration/deceleration, especially when coordinating with many other swimmers
during schooling [76, 31, 47].

To model soft actuation materials, an alternative to prescribing the entire body deformation
or curvature is to build micro-mechanical models that take into account the intrinsic muscle-like
behaviors locally. Generally speaking, there are typically two diﬀerent ways of modeling the

2

so-called artiﬁcial muscles. One way is to decompose the total deformation of the material into
two parts: elastic deformation caused by mechanical stress and active deformation/strain caused
by other stimuli [14, 95, 58]. For example, deformation of hydrogel can be decomposed into
the elastic part and the swelling part [14]; deformation of liquid crystal elastomer can also be
decomposed into the elastic strain and the stimuli-induced strain. Another approach is the active
stress method, where the total stress is decomposed into a mechanical part and an active part, both
of which can induce deformation [91].
In many cases, the above two methods can be ﬂexibly
used in achieving arbitrary deformations. In fact they are mathematically identical, due to the fact
that the active deformation introduced here is a ﬁrst-order approximation of nonlinear elasticity,
and is independent of mechanical stress. These two methods have been widely used in modeling
active soft elastic structures in biological systems (e.g., cardiovascular mechanics) and synthetic
active soft materials [2]. While these elasticity models provide much more ﬂexible yet simple local
actuation schemes, solving for the resultant nonlinear deformation is computationally expensive
when coupled with another (nonlinear) ﬂow solver. To address these challenges, we have recently
implemented the active strain approach in a ﬁctitious domain (FD) method, and demonstrated that
the FD/active-strain method is capable of solving the fully-coupled FESI for an arbitrary shape
subjected to distributed contractive active strains [58], which is presented in Chapter 3.

When a soft material interacts with ﬂuid ﬂows, the coupling between the ﬂuid force and material
elasticity can generate waves propagating at the ﬂuid/solid interface. Understanding this behavior
is critical to the study of biological swimmers and their artiﬁcial analogs as it can greatly eﬀect
the viability and eﬃciency of a potential swimming mode. The soft materials of which soft robots
are made are often viscoelastic with nonlinear behavior. Previous studies suggest that nonlinearity
has a destablizing eﬀect on the ﬂuid-elastic solid system thus further amplifying the importance
of understanding this behavior [33]. Chapter 2 presents an analytical study of the stability of the
interface between a nonlinear, viscoelastic gel ﬁlm and a Newtonian ﬂuid undergoing Couette ﬂow.

3

CHAPTER 2

STABILITY OF COUETTE FLOW PAST A GEL FILM

The interaction of viscous ﬂuid and soft objects is of considerable importance in a wide range of
problems, such as rheology of complex ﬂuids, coating, biological locomotion, and soft lubrication
[108, 90, 29, 87]. Kumaran et al. [51] ﬁrst studied the stability of an incompressible viscoelastic gel
ﬁlm in a Newtonian Couette ﬂow by ignoring inertia, where a linear model is adopted to describe
deformation. They found the ﬂuid/solid interface becomes unstable when the imposed shear goes
beyond a certain critical number, and the critical value of the imposed ﬂuid shear strength varies
inversely with the ﬁlm thickness for suﬃciently thick solids, which is veriﬁed by the following
experiments by Kumaran & Muralikrishnan [52, 69]. For suﬃciently thin solids, however, the
linear elastic model overpredicts the critical values of the ﬂuid shear that drives the interfacial
instability. Gkanis & Kumar [33] studied the similar problem by employing a neo-Hookean solid
model which admits ﬁnite/large deformation. While observing similar behaviors for thick gels,
their model predicts a much smaller critical shear for thin gels, which suggests that incorporating
solid nonlinearity can eﬀectively destabilize the system. It has been identiﬁed that such interfacial
instability under shear is mainly due to the ﬁrst normal-stress diﬀerence appearing at the base state
solutions, which is known as the Poynting eﬀect in nonlinear solids [78, 64], and is also similar to
the situation of two coupled viscoelastic liquids [80, 11]. When surface tension is incorporated into
the model, it changes the short-wave instability in thin ﬁlms to be ﬁnite-wavelength. Later Gkanis
& Kumar [34] investigated how the ﬂow ﬁeld and combined pressure gradient impact stability of a
neo-Hookean solid.

Although such elastohydrodynamic instability has been studied for simple neo-Hookean solids,
in practice soft materials often exhibit more complicated constitutive behaviors than hyperelasticity.
Especially for gel-like (e.g., hydrogel) materials that are typically composed of a large amount of
solvent such as water and long chain polymers which can form a complex network by chemical
crosslinkers such as covalent bond, or physical crosslinkers such as ionic bond or van der waals

4

Figure 2.1: Schematic for a Couette ﬂow past a gel ﬁlm. Nonlinear solid models: (a)
neo-Hookean, (b) Kelvin-Voigt, (c) Zener.

interaction. The combined eﬀects of elasticity, viscosity, and surface tension may generate new
dynamic behaviors of the material [53], as well as some intriguing interfacial instability phenomena
[68, 10, 102].

Here we study stability of an incompressible, impermeable gel ﬁlm when subjected to a Newto-
nian Couette ﬂow; see schematic in ﬁgure 2.1. We employ an Eulerian representation in both ﬂuid
and solid phases, and solve for the velocity, pressure and stress ﬁelds simultaneously. Compared
to the previous works [51, 33, 34], we examine the contributions from both the viscous eﬀect and
the nonlinear elasticity to the ﬂuid/solid interface instability, especially for the thin ﬁlms where
material nonlinearity is manifested.

2.1 Mathematical Model

2.1.1 Constitutive equations

We ﬁrst derive equivalent “diﬀerential” models for the evolution equations of solid stress in an
Eulerian frame. In this way, the governing equations in two phases are consistently deﬁned in terms
of velocity, pressure and stress [27, 28, 29]. In characterizing material constitutive relations, virtual
models of springs and dashpots are often employed to describe the contributions from the elastic
and viscous eﬀects, respectively. Diﬀerent combinations of springs and dashpots can eﬀectively
represent complex viscoelastic behaviors as highlighted by the schematics for the three models

5

x1x2UwFluidx2 = Rx2 = 0x2 = -HRvSSvη12totv′v′′v1S2Sηtot12Gel(a)           (b)                 (c)n(a-c) as shown in ﬁgure 2.1. The constitutive relation for the simplest nonlinear hyperelastic solid,
i.e., neo-Hookean solid, can be characterized by a linear stress-strain relation

τnh = S(B − I)

(2.1)

in a spring, where τnh is the elastic stress tensor, S is the shear modulus, and

B = F · FT

(2.2)
is the Finger tensor where F (or index form Fi j = ∂xi/∂Xj) is the deformation gradient tensor
serving as a mapping from the initial (X) to the current (x) reference frame. It is straightforward to
show that the rate of change of the Finger tensor satisﬁes

(cid:219)B − L · B − B · LT = 0,

(2.3)

where

L = (cid:219)F · F−1 = (∇v)T

(2.4)
(here Li,j = ∂vj/∂xi ) and dot represents the material time derivative (cid:219)Bi j = ∂Bi j/∂t + vk Bi j,k [45].
Then the diﬀerential form of the neo-Hookean model can be written as [28]:

neo-Hookean :

∇
τnh = (cid:219)τnh − L · τnh − τnh · LT = S(∇v + (∇v)T),

where

∇
τ = (cid:219)τ − L · τ − τ · LT

(2.5)

(2.6)

is the so called upper-convected time derivative.

To proceed, the simplest viscoelastic solid, a Kelvin-Voigt solid, can be characterized by
connecting a spring with a dashpot in parallel. When deformed, each element exerts stress (τ) and
undergoes deformation rate (v) accordingly. The total stress of the element is hence deﬁned by

τtot = τ1 + τ2,

(2.7)

where τ1 and τ2 are the stresses in the respective branches. By deﬁning τ1 = τnh from (2.5) and
the viscous stress as

τ2 = η(∇v + (∇v)T)

6

(2.8)

(η is a solid viscosity), we derive a two-equation model for a Kelvin-Voigt material:

Kelvin-Voigt : τtot = τ1 + η(∇v + (∇v)T),

∇
τ1 = S(∇v + (∇v)T).

(2.9)

However, this model is known to produce reasonable values of creep but incorrectly predicts the
behavior of stress relaxation [63], which can be improved by adding a second neo-Hookean spring
element (with modulus S2) in series with the dashpot to provide a Zener model. Note that on the
right branch, we need to consider individual velocities associated with the spring v(cid:48) and the dashpot
v(cid:48)(cid:48), which are related to the total velocity v through
(cid:48) + v

(2.10)

v = v

(cid:48)(cid:48)

.

Therefore, the neo-Hookean stress in the right branch (spring 2) now satisﬁes

the viscous stress becomes

∇
τ2 = S2(∇v

(cid:48) + ∇v

(cid:48)T);

τ2 = η(∇v

(cid:48)(cid:48) + ∇v

(cid:48)(cid:48)T).

(2.11)

(2.12)

Again, by making use of the fact that τtot = τ1 + τ2, after some algebra we derived the governing
equations for the total stress (τtot) and the neo-Hookean (τ1) stress on the left branch (spring 1) in
the Zener model:

(cid:16)∇v + ∇vT(cid:17)

,

(2.13)

Zener : τtot + λ1

∇
τtot = τ1 + λ2(∇v + ∇vT),

∇
τ1 = S1

where λ1 = η
S2

and λ2 = η
S2

(S1 + S2).

2.1.2 Universal Eulerian description

Consider an incompressible Newtonian ﬂuid past a gel ﬁlm as shown in ﬁgure 2.1 with material
constitutive relations constructed in Section 2.1.1. The ﬂuid layer with thickness R is bounded by
the solid at the bottom and a rigid wall on the top moving with a constant velocity Uw. The solid
gel with thickness HR is ﬁxed on a rigid plate. Following Kumaran et al. [51] and Gkanis &

7

Kumar [33], we use R as the length scale, neo-Hookean modulus S1 as the pressure/stress scale,
µ f /S1 as the time scale, and RS1/µ f as the velocity scale. The Reynolds number is then deﬁned
by Re = ρ f S1R2/µ2
f . We ignore the inertia eﬀect (Re (cid:28) 1), and assume the ﬂuid ﬂow is governed
by the Stokes equation. The dimensionless governing equations can be written as:

∇ · v f = 0, −∇p f + ∇2v f = 0,

(2.14)

where p f and v f are the pressure and the velocity ﬁeld in the ﬂuid, respectively. The solid is
assumed to be incompressible. The conservation of mass and momentum equations are written in
the Eulerian frame as:

∇ · vs = 0, −∇ps + ∇ · τs = 0,

(2.15)

where vs and ps are the dimensionless solid velocity and pressure, respectively. And τs is a total
stress tensor. The constitutive equations for the Zener model become

τs + ˆλ1

∇
τs = τnh + ˆλ2(∇vs + ∇vT
s ),

∇
τnh = ∇vs + ∇vT
s ,

(2.16)

(cid:16) S1+S2

(cid:17) ηs

µ f

S2

ηs
µ f

and ˆλ2 =

where ˆλ1 = S1
respectively represent the relaxation and viscous eﬀects,
S2
and τnh represents the dimensionless neo-Hookean stress. Note that ˆλ1 and ˆλ2 are related to each
other by ˆλ2 = ˆλ1 + ηs/µ f , and hence ˆλ1 < ˆλ2. Furthermore, in the limit ˆλ1 → ˆλ2 = λ, it is easy
to show that the above equations reduce to

λ(τs − τnh)∇ + (τs − τnh) = 0,

(2.17)

or in the Lagrangian frame,

(cid:104)

F−1 · (τs − τnh) · F−T(cid:105)

(cid:104)

F−1 · (τs − τnh) · F−T(cid:105)

+

= 0.

(2.18)

λ

d
dt

When there is no pre-stress applied on the gel, it is trivial to show that τs = τnh, i.e., reduction from
the full Zener model to the neo-Hookean model. Therefore, the neo-Hookean model ( ˆλ1 = ˆλ2 = 0)
and the Kelvin-Voigt model ( ˆλ1 = 0, ˆλ2 (cid:44) 0) appear to be the two asymptotic limits of the Zener
model.

8

To supplement the governing equations with boundary conditions, we assume periodicity in
x1−direction, and apply no-slip conditions on the top (moving) and bottom (stationary) rigid walls.
The dimensionless wall velocity is v f = Gˆe1 at x2 = 1, where G = µ f Uw/RS1 characterises the
strength of the applied shear ﬂow ﬁeld on the solid. At the ﬂuid/solid interface, the velocity and
traction are continuous:

v f = vs, σf · n + T κn = σs · n,

(2.19)

where n is the normal vector to the interface, T is the dimensionless surface tension scaled by
S1R and κ is the curvature of the interface scaled by R−1. The Eulerian models derived above
are equivalent to the conventional Lagrangian models which typically solve for the displacement
ﬁeld, and simplify the mathematical treatment of the ﬂuid-solid coupling at the interface to permit
a more straightforward implementation of the velocity and traction continuity conditions [27].

2.2 Linear Stability Analysis

2.2.1 Base state

To perform a linear stability analysis, we need to derive the base-state solutions when the interface
is ﬂat. The steady-state velocity solution in the ﬂuid is simply a planar Couette ﬂow v f = (v1, v2) =
(Gx2,0), and the ﬂuid stress is given by

σf = −p f I + τf =(cid:169)(cid:173)(cid:173)(cid:171) −p f G

G − p f

(cid:170)(cid:174)(cid:174)(cid:172) .

In the solid, the base-state solutions can be derived by solving the initial transient dynamics during
which a fully developed shear ﬂow start exerting on the gel. For the neo-Hookean model (τs = τnh),
the (symmetric) total stress tensor in the solid is deﬁned as

σs = −ps I + τnh =(cid:169)(cid:173)(cid:173)(cid:171) −ps + τnh

τnh
12

11

(cid:170)(cid:174)(cid:174)(cid:172) .

τnh
12

− ps + τnh
22

(2.20)

(2.21)

It is also reasonable to assume that during the initial transient, the induced velocity ﬁeld is unidi-
rectional, i.e., ∂/∂x1 = 0 and vs = (v1(x2),0). Then the constitutive equations in (2.16) reduce

9

to

∂τnh
11
∂t

∂τnh
12
∂t
When there is no pre-stress in the solid, i.e. τnh
i j
integrating both sides of the remaining two stress equations leads to

∂τnh
22
∂t

= 2τnh
12

∂v1
∂x2

∂v1
∂x2

=

,

,

= 0.

(2.22)

= 0 at t = 0, it is apparent that τnh

22 = 0. Next,

∂v1
∂x2

dt =

∂u1
∂x2

0
0 v1dt is the displacement. Thus the total stress tensor at equilibrium reads

0

0

(cid:16)

(cid:17)2

∂

τnh
12
∂t

dt =G2.

Þ t
where u1 =Þ t

τnh
12 =

(2.23)

(2.24)

dt =

= G,

τnh
12

∂v1
∂x2

τnh
11 = 2

Þ t
(cid:170)(cid:174)(cid:174)(cid:172) =(cid:169)(cid:173)(cid:173)(cid:171) G2 G
τnh =(cid:169)(cid:173)(cid:173)(cid:171) τnh
σs = −ps I + τs =(cid:169)(cid:173)(cid:173)(cid:171) −ps + G2 G

τnh
12
τnh
22

11
τnh
12

G

G

0

Þ t
(cid:170)(cid:174)(cid:174)(cid:172) .
(cid:170)(cid:174)(cid:174)(cid:172) .

In the Kelvin-Voigt model and the Zener model, since vs = 0 at equilibrium, equation (2.24) also
serves as the base-state solution for the three solid models. Therefore, the base-state total stress
can be written as

Note that a ﬁrst normal stress diﬀerence appears, i.e., τs
n = (0,1)), equations (2.19) reduce to σf · n = σs · n and p f = ps.

− ps
11 − τs
22 = G2. At the ﬂat interface (i.e.,

(2.25)

2.2.2 Normal modes and interfacial conditions

A standard linear stability analysis is employed to study the growth of disturbances [9]. All physical
disturbance variables f (cid:48) can be expanded with respect to normal modes

f (cid:48)(x1, x2, t) = ˜f(x2) exp(ik x1 + αt),

(2.26)

where ˜f is the complex-valued amplitude functions, k is the wavenumber, and α is the complex-
valued growth rate. The governing equations in (2.14)-(2.16) are linearized, and then transformed
into the Fourier space for the disturbances of (˜v f , ˜p f ) in ﬂuid and (˜vs, ˜ps, ˜τs, ˜τnh) in solid. After

10

some algebra, we are able to derive two fourth-order ordinary diﬀerential equations for the vertical
velocity components ˜v

2 which admit the following general solutions:

f
2 and ˜vs

where

f

2 = A1 exp(k x2) + A2 exp(−k x2) + A3x2 exp(k x2) + A4x2 exp(−k x2),
˜v
2 = B1 exp(k x2) + B2 exp(−k x2) + B3 exp(ξ1x2) + B4 exp(ξ2x2),
˜vs
(cid:17)2

(cid:114)

(cid:16)

(cid:17)

(cid:16)

(cid:16)

(cid:17)(cid:16) ˆλ1 − ˆλ2

−G2α

1 + α ˆλ1

(cid:17) ± k

−iGk

1 + α ˆλ1

ξ1,2 =

1 + α ˆλ2

(2.27)

(2.28)

+

1 + α ˆλ2

,

(2.29)

Ai and Bi (i=1,2,3,4) are eight unknown coeﬃcients, four of which can be eliminated by boundary
conditions at the two walls. For the neo-Hookean solid, we obtain ξ1,2 = (−iG ± 1) k independent
of α, and again, recover the solutions obtained by Gkanis & Kumar [33].
The two solutions above are coupled through the kinematic condition

∂δ
∂t

+ v1

∂δ
∂x1

= v2,

(2.30)

imposed at the perturbed interface x2 = δ(x1) where |δ| (cid:28) 1.

= v(cid:48)

2 around the ﬂat plane at x2 = 0, and v(cid:48)
∂δ
∂t
velocity and traction continuity conditions lead to

It can be further linearized as
2 is the ﬂuid (or solid) velocity disturbance. The

,

s

,

(cid:48)
v
1

(cid:48)
f = v
2
f

(cid:48)
f + Gδ = v
1
(cid:48)
s
v
2
∂v(cid:48)
2
∂x1
∂v(cid:48)
f
+ T ∂2δ
2
∂x2
∂x2
1

+ G2 ∂δ
∂x1
= −p(cid:48)

+

(cid:48)
12,

= τ

(cid:48)
22.

s + τ

f

∂v(cid:48)
1
∂x2

−p(cid:48)
f

+ 2

(2.31)

(2.32)

(2.33)

(2.34)

The ﬁrst normal stress diﬀerence from the base state (G2) now appears at the left-hand-side of
(2.33). The above equations are transformed using the normal modes to the Fourier space where
we solve the growth rate α numerically for given values of H, k,T, G, ˆλ1, and ˆλ2.

11

Figure 2.2: The growth rate Re(α) as a function of k for the neo-Hookean solid: (a)
H = 0.4,T = 0; (b) H = 4.0,T = 0; (c) H = 0.4,T = 10; (d) H = 4.0,T = 10. The open circles are
the results obtained by Gkanis & Kumar [33].

2.3 Results

2.3.1 Growth Rate
In ﬁgure 2.2, we show the real part of the growth rate of disturbance, Re(α), by varying thickness
ratio (H), surface tension (T), and applied shear strength (G). When there is no surface tension at the
interface (T = 0), for a thin solid (H = 0.4, panel(a)), Re(α) reaches a plateau at high wavenumber
asymptotically, indicating a short wave instability. For a relatively thick ﬁlm (H = 4.0, panel(b)),
the maximum growth rate occurs at a ﬁnite k, hence indicating a ﬁnite-wavelength instability.
Panels(c) and (d) show Re(α) as a function of k when surface tension is included, and chosen
as T = 10. In both cases, the surface tension can eﬀectively eliminate instability growth at high
wavenumber. The obtained growth rates of the neo-Hookean model are identical to those obtained
by Gkanis & Kumar [33].

Although the instabilities in panels(b)-(d) all appear to be ﬁnite-wavelength, their proﬁles are
very diﬀerent, leading to distinctive features shown later in the corresponding marginal-stability
curves. In panel(b) where T = 0, while instability is induced at ﬁnite k and then more and more
enhanced as the applied ﬂuid shear G increases, the short-wave feature still remains as Re(α)
saturates at large k. In panels(c) and (d) when T (cid:44) 0, there are typically two separated regimes
of peaks visible: One occurs at small k, and the other occurs at relatively large k. It is seen that

12

kRe(α)0246810-0.4-0.200.20.4(a)G=7.0G=5.0G=4.0G=3.0k0246810-1-0.8-0.6-0.4-0.200.20.4(b)G=5.0G=4.0G=2.0G=1.0k0246810-0.6-0.4-0.200.2(c)G=13.0G=11.0G=9.0G=7.0k0246810-1.2-1-0.8-0.6-0.4-0.20(d)G=5.0G=4.0G=2.0G=1.0Figure 2.3: Comparison of the growth rate Re(α) as a function of k for the neo-Hookean,
Kelvin-Voigt, and Zener solids: (a) H = 0.4,T = 0, G = 5; (b) H = 4.0,T = 0, G = 5; (c)
H = 0.4,T = 10, G = 11; (d) H = 4.0,T = 10, G = 5.

with surface tension, the disturbance grows much faster at ﬁnite/large k regimes as G increases.
Conversely, for thick ﬁlms, instabilities appear to be more strengthened at small k.

Next, we consider the two viscoelastic models and make comparisons with the neo-Hookean
model. The ﬁve cases presented in ﬁgure 2.3 follow a progression from the basic neo-Hookean
model (black) to the Kelvin-Voigt model (red), and then to a full Zener model (blue). In the Kelvin-
Voigt model where ˆλ1 = 0 and ˆλ2 (cid:44) 0, again, we ﬁnd that increasing the eﬀective solid viscosity
(i.e., increasing ˆλ2) suppresses instability at large k in both thin and thick gels, which makes the
ﬁnite-wavelength instabilities dominate. On the other hand, in panels(a,c) we also observe that the
instabilities are enhanced at ﬁnite wavenumbers to form local maximum for thin gels.

To examine the relaxation eﬀect in the Zener model, we ﬁx ˆλ2 = 0.1 and vary the dimensionless
relaxation time ˆλ1. We observe that increasing ˆλ1 shifts Re(α) upwards. As ˆλ1 approaches ˆλ2, this
eﬀect becomes increasingly pronounced and the curve gradually returns to that of the neo-Hookean
model, consistent with the analytical predictions as ˆλ1 → ˆλ2 in Section 2.1.2.

2.3.2 Mechanism

As noted by Gkanis & Kumar [33, 34], the interfacial instability in the neo-Hookean solid is
due to the coupling of the ﬁrst normal stress diﬀerence and the surface ﬂuctuation δ appearing
in the interfacial condition. Here we adopt their explanation by considering a sinusoidal surface

13

kRe(α)0246810-0.200.2Neo-HookeanKelvin-Voigt, λ2=0.05Kelvin-Voight, λ2=0.1Zener, λ1=0.05, λ2=0.1Zener, λ1=0.09, λ2=0.1(a)k0246810-0.200.2(b)k0246810-1.2-0.8-0.400.4(c)k0246810-0.8-0.6-0.4-0.20(d)(cid:16) 2π

(cid:17)

Figure 2.4: (a) Schematic for a sinusoidal surface perturbation. The contributions of the ﬁrst
normal-stress diﬀerence and the solid viscous eﬀect are marked by the black and red solid arrows,
respectively. The open arrows represent the stabilizing eﬀect due to surface tension. (b,c):
Critical values of the imposed strain Gc and the critical wavenumber kc as a function of H. The
black dotted lines are drawn at the locations where either Gc or kc tends to inﬁnity. Inset in (c):
Growth-rate curves near H ≈ 1 where short-wave instabilities recur.

k x1

as shown in the schematic in ﬁgure 2.4(a). The traction continuity
perturbation δ = ε sin
condition in (2.33) essentially suggests the coupling term G2δx1 has a 90o phase lag compared to
velocity (and shape δ) ﬂuctuations in (2.31), which tends to compress the wave crests and extend
the valleys (marked by the black solid arrows), and hence destabilize the interface.

To understand the solid viscous eﬀects in the Kelvin-Voigt model, we rewrite the shear stress

component in equation (2.33) in the Fourier space as:

(cid:34)

(cid:33)(cid:35)

(cid:32)

D2
k2

+ 1

˜τelast
12

− ˜τ

f

12 = G2 ˜δx1 − ˜τvis
12 =

G2 − α ˆλ2

˜δx1 .

(2.35)

In the above we split the solid stress τs into the contributions from the hyperelastic (due to the
neo-Hookean elasticity, and denoted by “elast”) and viscous (due to the viscous term in the Kelvin-
Voigt model, and denoted by “vis”) eﬀects. It is clearly seen that on the right-hand-side, as k → ∞
the viscous term provides an asymptotic limit −α ˆλ2 acting against the ﬁrst coupling term due to the
normal stress diﬀerence. As illustrated by the red arrows in ﬁgure 2.4(a), this relaxes the pinching
eﬀect at large k to stabilize the surface though its behavior at small k is somewhat complex. A
similar argument can be made based on the traction balance in the normal direction when examining

14

112sinxxk(a)1x2xHGc10-110010110-1100101Neo-Hookean, T=0Neo-Hookean, T=10Kelvin Voigt, λ2=0.1, T=0Kelvin Voigt, λ2=0.1, T=10Zener, λ1=0.05, λ2=0.1, T=0Zener, λ1=0.05, λ2=0.1, T=10(b)-1Hkc10-110010110-1100101102(c)viscoelastic, T=0Small HH ≈ 1Large HRe(α)kcRe(α)kckRe(α)kc(cid:18)

(cid:19)

˜σelast
22

− ˜σ

f

22 = T ˜δx1x1 − ˜σvis
22 =

T + 2α ˆλ2

D
k2

˜δx1x1 .

the surface tension eﬀect. Following the same procedure we rewrite equation (2.34) as

The second derivative on δ results in a 180o phase shift so there is a normal force acting to restore the
perturbed interface. Note that the right-hand-side of (2.36) can be written as
suggesting that the contribution from the surface tension is scaled by k2, and hence the inclusion
of surface tension causes a signiﬁcant damping eﬀect to suppress disturbance growth, especially at
large k [51, 33].

(2.36)

(cid:16)−k2T − 2α ˆλ2D

(cid:17) ˜δ,

2.3.3 Marginal Instability

We show the marginal stability curves for the critical shear stress Gc (ﬁgure 2.4(b)) and the
corresponding critical wavenumber kc (ﬁgure 2.4(c)) as functions of the ﬁlm thickness ratio H.
One obvious observation is that as H becomes large, all the models produce almost the same results
as the linear elastic model [51] due to smaller deformation. We also observe that the marginal
stability curves of the Zener model again vary between those of the neo-Hookean model and the
Kelvin-Voigt model.

In panel(b) when surface tension is lacking, the three models in fact predict very similar values
of Gc, suggesting that adding viscous eﬀects to the solid does not necessarily stabilize the interface.
In contrast, the ﬁnite-wavelength instabilities observed in ﬁgure 2.3(a) and (b) for the Kelvin-Voigt
model may even cause a slight reduction in Gc, making the interface even more unstable under
shear. The critical stresses for thin ﬁlms are found to be around 2.93 while vary inversely with H
for thick ﬁlms, with transitions occurring for ﬁlms with ﬁnite thickness (H ≈ 1).

Adding surface tension eﬀectively stabilizes the interface for all values of H, and yields ﬁnite-
wavelength instabilities. Compared to the neo-Hookean model, the solid viscous eﬀects signiﬁ-
cantly change the critical values for thin ﬁlms, leading to yet more complex behaviors. When H is
suﬃciently small, we ﬁrst identify a narrow regime where the system is stable to the ﬂuid shear,
i.e., Gc → ∞ (marked by black dotted lines). At relatively large H, ﬁnite-wavelength instabili-

15

ties appear, and predict lower critical values than the neo-Hookean case. As H further increases
close to the transition regime, Gc rises again in both the Kelvin-Voigt and Zener models, up to
maximum values around H = 0.7. Such non-monotonic behaviors are due to the ampliﬁcation
(suppression) of instability in the viscoelastic ﬁlms at small (large) k as discussed in ﬁgure 2.3(c),
which correspondingly picks relatively low (high) critical values of Gc for thin (thick) ﬁlms.

While the Gc−H curves are very similar, the marginal curves of kc in panel(c) exhibit distinctive
features between the three models at T = 0. At small thickness, the ﬁnite-wavelength instabilities in
the two viscoelastic models are visible by the critical values of kc beyond the stable regime; while
kc → ∞ for the neo-Hookean solid due to its short-wave nature. Interestingly, we observe that
in both the Kelvin-Voigt and Zener models, discontinuities of kc occur near the transition regimes
when H ≈ 0.7, in accordance with the peaks in panel(b). As illustrated by the growth-rate curves
in the inset of panel(c), in the transition regimes, increasing ﬁlm thickness H drives a secondary
ﬁnite-wavelength instability arising at small wavenumber. As the growth-rate maximum shifts
from large towards small k, short-wave instabilities may recur, and hence kc saturates at inﬁnity,
resembling those observed in the thin neo-Hookean ﬁlms in ﬁgure 2.2. Moreover, inclusion of
surface tension eﬀectively reduces kc as shown in the lower branch of marginal curves in panel(c).
Also a stable regime at small H is identiﬁed for viscoelastic gels in accordance with panel(b) where
Gc tends to inﬁnity, and then ﬁnite-wavelength instabilities dominate as H increases.

2.4 Concluding Remarks

In this work, we investigate the interfacial stability of a soft gel ﬁlm when subjected to a
Newtonian Couette ﬂow at zero Reynolds number. We have constructed three diﬀerential models
for the solid stress including the neo-Hookean, Kelvin-Voigt, and Zener models, and solve for the
velocity, pressure and stress ﬁelds in both ﬂuid and solid phases simultaneously in the Eulerian
frame. We have performed linear stability analysis to study the interfacial instability by exploring
the parameter spaces, and compared the critical behaviors of the three solid models. We focus
our attention in the thin-ﬁlm regime where the material nonlinearity is pronounced; while for

16

suﬃciently thick ﬁlms, all the models produce very similar results as those predicted by linear
models due to small deformation. We ﬁnd the interfacial instability is driven by the ﬁrst-normal
stress diﬀerence in the base-state solution in the solid phase. Compared with the neo-Hookean
model, inclusion of solid viscous eﬀects leads to ﬁnite-wavelength instabilities in both thin and
thick ﬁlms, although we still observe short-wave instabilities recur in a narrow transition regime
at ﬁnite ﬁlm thickness. When surface tension is included, the ﬂuid/solid interface becomes more
stable, and ﬁnite-wavelength instabilities are found to dominate in all models.
In addition, for
viscoelastic gels we have observed stable regimes for suﬃciently thin ﬁlms, as well as intriguing
non-monotonic features on the marginal stability curves near the transition regimes where the ﬁlm
thickness is ﬁnite.

It is interesting to observe that when surface tension is negligible, the critical values of the
applied shear strength (Gc ≈ 2.93) seem to be well-characterized by the neo-Hookean model, i.e.,
material’s hyperelasticity. Further analysis reveals that the results obtained above can be directly
extended to more general hyperelastic Mooney-Rivlin model which is used for constitutive relations
in a wide range of materials. Recall the constitutive equation for a Mooney-Rivlin solid is given as

τs = g1B + g2B−1,

(2.37)

where g1 and g2 are scalar functions of the invariants of B [59]. Thus the neo-Hookean model is
a special case of the Mooney-Rivlin model by choosing g1 = 1 and g2 = 0. It is easy to show that
B−1 satisﬁes a lower-convected derivative:

(cid:219)B

−1 + LT B−1 + B−1L = 0

(2.38)

[45]. Following the same procedure in Section 2.2, we ﬁnd that a lower-convected model for the
elastic stress τs ∝ I − B−1 yields the base-state solutions

σs =(cid:169)(cid:173)(cid:173)(cid:171) −ps

(cid:170)(cid:174)(cid:174)(cid:172)

G

G − ps − G2

p f = ps + G2,

17

and

(2.39)

(2.40)

which lead to identical stability analysis results as the neo-Hookean model where τs ∝ B − I.
Hence the Mooney-Rivlin model will generate the same results as the neo-Hookean model as well.
With this study, it is straightforward to construct more elaborated solid models (e.g., poroelasticity),
and apply similar methodologies to investigate their stabilities when coupled with ﬂuid ﬂows. It is
also desired to perform direct simulations to resolve new physics and phenomena associated with
nonlinear deformations in ﬂuid/elastic-structure interactions.

18

CHAPTER 3

A FLUID–STRUCTURE INTERACTION STUDY OF SOFT ROBOTIC SWIMMER

USING A FICTITIOUS DOMAIN/ACTIVE-STRAIN METHOD

In the past decades, numerous computational methods have been developed to simulate interactions
between ﬂuid and moving elastic objects. One example is the so-called boundary conforming
method, such as the Arbitrary Lagrangian-Eulerian (ALE) ﬁnite element method [40, 27] and
the space-time ﬁnite element method [62].
In these methods, the governing equations of both
the ﬂuid and solid phases are solved on one set of mesh that is adapted to the deforming solid
objects or geometric boundaries. Although the sharp-interface methods can accurately resolve the
complex geometries of moving boundaries, the computation cost is high, especially in 3D. This
is because for the strongly coupled ﬂuid-elastic system, additional subiterations between the ﬂuid
and solid solvers are necessary in order to achieve better convergence [46, 92]. Moreover, frequent
re-meshing might be necessary in order to follow boundary movement. When the moving-mesh
technique is employed (e.g., ALE method), the associated mesh velocities have to be solved as
additional unknowns together with the variables in the ﬂuid and solid phases [40, 27]. In addition,
unstructured meshes (e.g., triangular, quadrilateral, etc.) are often employed in these methods,
which adds complexities in parallelization.

Instead of employing body-ﬁtted meshes, the so-called Cartesian grid methods typically use
either separate or non-boundary ﬁtted meshes. This approach simpliﬁes the initial mesh generation
and avoids the need for moving and/or adaptive meshing to resolve the changing interface which
reduces the overall cost and complexity. For example, the immersed boundary method employs
overlaid Eulerian and Lagrangian meshes to solve ﬂuid/elastic structure interactions: The Navier-
Stokes (N-S) equations are solved on a ﬁxed Eulerian mesh; while the embedded boundaries are
tracked by a set of freely moving Lagrangian points. To account for the no-slip conditions at
ﬂuid-solid interfaces, appropriate forcing density terms are added to the N-S equations by using
proper interpolations between the solution variables at the ﬁxed grid points in the vicinity of the

19

moving boundary and the nearest Lagrangian points [49, 77, 67, 66, 6, 94]. Instead of solving
the N-S equations in the ﬂuid phase, both the particle-based methods (e.g., Smoothed Particle
Hydrodynamics (SPH) method [3]) and the kinetic methods (e.g., Lattice Boltzmann method
(LBM) [37, 13, 30, 41, 42]) can be ﬂexibly combined with the immersed boundary method to
handle complex geometries, and has become more and more popular in dealing with ﬂuid-structure
interactions (FSI).

Inspired by the recent experimental designs [70, 74], in this chapter we numerically study
the swimming motion of soft robotic swimmers that perform large reversible elastic deformations
driven by the imposed actuation. We assume that the elastic material is hyper-elastic, and can be
described by the neo-Hookean constitutive relation that permits ﬁnite/large nonlinear deformations.
To model the muscle-like behaviors, we adopt an active strain model to distribute contracting strains
on soft elastica, which is then coupled with the ﬁctitious domain (FD) method [104, 105] to resolve
the FSIs. The central idea of the active strain model is based on the multiplicative decomposition
where the (total) deformation gradient tensor is split into two parts that represent the sequential
contributions from the imposed contractions and the elastic response of the material [54, 55]. It
is diﬀerent from the so-called active stress approaches where an additional stress term (diﬀerent
from applied mechanical stress) is introduced into the material model [2]. By appropriately
applying contracting strains on elastic plates of diﬀerent (e.g., square, circular, or bell-like) shapes
and rigidities, we show that soft robot prototypes of simple geometries can eﬀectively perform
undulatory or jet-propulsion swimming motions.

3.1 Mathematical Model and Numerical Method

3.1.1 Hill’s Muscle Model

The three element Hill’s muscle model is the standard basis for modeling the mechanical behavior
of skeletal muscle [39, 61]. In its original form it is a 1D model composed of two nonlinear springs
and a contractile element. The contractile element (CE) is arranged in series with one spring and
both of which are in parallel with another spring. Generally speaking, the series spring (SE) models

20

the elastic behavior of the active components of the muscle while the parallel spring models the
passive elastic behavior, even while the contractile element is inactive. The contractile element is
the active part of the muscle that contracts while activated and typically is freely extensible while
not.

Figure 3.1: Schematic of the 1D three element Hill’s muscle model. It is composed of a
contractile element (CE) in series with a nonlinear spring (SE) which are both in parallel with
another nonlinear spring (PE).

Many variations exist and for our study we ﬁnd that a two element model, excluding the parallel
spring element, can model the types of active materials that we are investigating. These materials
tend to be uniform in nature so we expect the elastic behavior to be reasonably deﬁned by a single
element. As opposed to biological muscle which may be composed of many smaller structures and
materials.

3.1.2 Active Strain Method

To begin with, we seek mathematical models that describe the nonlinear mechanics of soft actuating
material. Generally speaking, there are typically two diﬀerent ways of modelling soft actuating
material or so-called artiﬁcial muscle. One way is decomposing the total deformation of the material
into two parts: elastic deformation caused by mechanical stress and active deformation/strain caused
by other stimuli such as light, electric ﬁeld, thermal ﬁeld etc. For example, deformation of hydrogel
can be decomposed to the elastic part and swelling part [14]; deformation of liquid crystal elastomer
can also be decomposed to the elastic strain and stimuli-induced strain. In experiments, active strain

21

CESEPEcan be easily measured. Without any mechanical load, deformation of a free-standing soft material
can be measured as a function of the intensity of external stimuli [95]. Here we employ an active
strain approach which applies principal contractions in a manner similar to artiﬁcial muscle. To
apply actuation the deformation gradient tensor F is decomposed into an active deformation tensor
Fa and an elastic deformation tensor Fe following multiplicative decomposition [54, 55] such that

F = Fe · Fa.

(3.1)

(cid:113)

Here Fa is eﬀectively an arbitrary function applied to the reference conﬁguration that can be
designed in terms of the desired location, direction, sign and strength. For incompressible solids
that are considered here, we apply an incompressible restriction on actuation such that det(F) = 1.
For simplicity Fa can be deﬁned in the principal coordinates and transformed to the desired
orientation by a standard solid body rotation coordinate transformation. For an artiﬁcial muscle
that contracts uniaxially when activated, we deﬁne Fa = diag[λ1, λ2, λ3] where λ1 < 1 represents
λ−1
1 > 1 are the correspondingly strains in the other two
a principal stretch, and λ2 = λ3 =
directions. With total F being mapped appropriately, the resultant elastic stress can be calculated
through the constitutive relation τ = τ (Fe) = τ
for various diﬀerent kinds of materials
(hyperelastic, viscoelastic, composite, etc.). In this study, we assume the active strain tensor to
be either a constant or time-dependent. Nevertheless, it is important to mention that the active
strain of real muscle can also depend on its loading condition. Therefore, to better characterize
the behaviors of artiﬁcial muscles (e.g., electroactive polymer or liquid crystal elastomer) under
actuation, systematic experimental measurements are required to obtain the relationship between
the active strain tensor and the imposed loading conditions.

(cid:16)F · F−1

(cid:17)

a

Alternatively, another widely used approach is the active stress method where the total stress is
decomposed into a mechanical part and an active part, and both of which can induce deformation.
For instance, Maxwell stress is usually introduced in the constitutive models of dielectric elastomer
[91].
In many cases, the above two methods are mathematically identical, due to the fact that
the active deformation introduced here is a ﬁrst order approximation of nonlinear elasticity, and is
independent of mechanical stress. There are also other approaches to implement (soft) actuation

22

without using stress or strain decompositions. For example, to mimic swimming motions of slender
ﬂexible swimmers such as sperms or nematodes using the immersed boundary method, active forces
are directly added to the right-hand-side of the N-S equation by taking the variational derivative of
the discrete elastic energy that is constructed based on the targeted local stretches and curvatures
[25, 93].

3.1.3 Distributed Lagrange Multiplier/Fictitious Domain Method (DLM/FD)

3.1.3.1 Formulation

To resolve the FSI of soft robotic swimmers, we implement the above active strain model in
a ﬁctitious domain (FD) method proposed by Yu [104], who extended the FD formulation of
Glowinski et al. [35] for the rigid particles to the case of ﬂexible bodies. The key idea of the FD
method is that the interior of the solid is assumed to be ﬁlled with a ﬁctitious ﬂuid that is constrained
to move at the same velocity with the solid by a pseudo body force (i.e. Lagrange multiplier). The
advantage of the FD method is that the ﬂow ﬁelds can be solved eﬃciently with a ﬁxed Cartesian
grid. The reader is referred to [104] for the details on the derivation of the FD formulation and the
numerical schemes.

In the following, we brieﬂy introduce the formulation and the numerical schemes of the FD
method. Suppose that a deformable body of density ρs is immersed in the incompressible Newtonian
ﬂuid of viscosity µ and density ρ f . Let Ω denote the entire computational domain containing both
solid and ﬂuid domains, and P(t) represent the solid domain. It is noted that tracking swimming
objects in the ﬁxed coordinates requires using a very large computational domain, which makes
computation very expensive especially in 3D. Alternatively, here we employ an instantaneous
inertial frame Ω that co-moves with the swimmer at a certain reference speed U [56]. Then the
dimensionless FD governing equations in weak form for the incompressible Newtonian ﬂuid in the
co-moving relative references as

Þ

(cid:18) ∂u f

∂t

Ω

(cid:19)

Þ

(cid:18)

Ω

+ ˆu f · ∇u f

· v f dx +

: ∇v f dx =

λ · v f dx,

(3.2)

(cid:19)

−pI +

1
Re

(∇u f )T

23

Þ

P

(3.3)
where ˆu f = u f − U. The dimensionless governing equations for neo-Hookean solid material
(τ = G(B − I)) are solved in the absolute references as

Ω

q∇ · u f dx = 0,

(cid:18) dus

Þ

(cid:20)
(cid:18)

P
(∇vs)T :

Þ

P

−

− Fr

g
g

dt
[∇u f + (∇u f )T]

(cid:19)

· vsdx +

P
dx = −

(ρr − 1)
1
−pI +
Re

(∇vs)T : [λ0 ln JI + G(B − I)]dx
λ · vsdx,

Þ

Þ

P

Þ

(cid:19)(cid:21)

Þ

P

(u f − us) · ζ dx = 0

(3.4)

(3.5)

Equations (4.3)-(4.6) represent the ﬂuid momentum equation, the ﬂuid continuity equation, the
solid momentum equation, and the velocity constraint in the solid domain, respectively. In these
equations, u f is the ﬂuid velocity, us the solid velocity, p the ﬂuid pressure, and λ the pseudo body-
force (i.e., the Lagrange multiplier). v f , vs, q and ζ are the corresponding variations, respectively.
For the original FD method for the passive deformation model, B = F· FT is the left Cauchy-Green
deformation tensor, here F being the deformation gradient tensor deﬁned as: F = ∂x/∂X, in which
x and X are the solid current and reference conﬁgurations, respectively. J is the determinant of F,
and J =1 for the incompressible solid. In contrast, for the active strain model studied,

B = Fe · FT

e = (F · F−1

a ) · (F · F−1

a )T,

(3.6)

where Fe is the elastic deformation tensor which causes the elastic stress, and Fa is an input
deformation tensor without generating elastic stress. Also g denotes the gravitational acceleration.
The following characteristic scales are used for the non-dimensionlization scheme: Lc for length,
c/Lc for Lagrange multiplier λ. The
Uc for velocity, Lc/Uc for time, ρ f U2
c for pressure p, and ρ f U2
following dimensionless control parameters are also introduced:

24

Density ratio : ρr =

ρs
ρ f

Material parameters

: λ0 =

Reynolds number

: Re =

Froude number

: Fr =

, G =

λ0
ρ f U2
c
ρ f UcLc

G
ρ f U2
c

µ
gLc
U2
c

(3.7)

where g is the gravity constant, λ0 is related to the compressibility property of the material and G
represents the shear modulus. It should be mentioned that for robotic swimmers, the characteristic
velocity is deﬁned as Uc = fcLc where a reference actuation frequency is chosen as fc = 1Hz. In
addition, in order to enforce the incompressibility constraint of the solid material, i.e., J = 1, we
impose a penalty function-like approximation by setting a large enough value for λ0 (104 ∼ 105).

3.1.3.2 Numerical scheme

The problem (4.3)-(4.6) is decomposed into ﬂuid, solid and Lagrange-multiplier sub-problems in
a partitioned coupling manner as follows [104]:

1. Fluid sub-problem (for u∗

Þ
(cid:32)

(cid:32)u∗

f and p)
f − un
∆t
(∇u∗

f

Þ

+

Ω

Ω

−pI +

(cid:33)

)

· v f dx

f − ˆun−1

f

· ∇un−1

f

: ∇v f dx =

Þ

Pn

(cid:33)

1
2

+

(3ˆun
f · ∇un
f )T + (∇un
f )T
2Re
Þ

q∇ · u∗

Ω

f dx = 0.

λn · v f dxn,

(3.8)

(3.9)

where the vector Lagrange multiplier λ serves as a pseudo body force, and needs to be
interpolated from the Lagrangian mesh via a Dirac delta function δ via:

25

Þ

Pn

λn(x) =

λn(X)δ(x − X)dX.

(3.10)

In the above the Eulerian and the Lagrangian mesh are denoted by x and X, respectively. We
approximate the delta function as a tri-linear function to transfer the quantities between the
unstructured Lagrangian mesh and the uniform Eulerian mesh (with grid size ∆x = ∆y = ∆z)
[77]

δ(x − X) =

1

∆x

∆x3 φ( x − X
φ(r) =

1 − |r|

0

)φ( y − Y

)φ( z − Z

∆x

∆x
|r| ≤1
|r| >1

.

).

(3.11)

(3.12)

This sub-problem is essentially the solution of the N-S equation. An eﬃcient ﬁnite-diﬀerence-
based projection method on a homogeneous half-staggered grid is employed [107]. All
spatial derivatives are discretized with the second-order central diﬀerence scheme. Here we
choose the reference speed U as the swimmer’s COM speed at previous time step, and the
computational domain position x is updated at the end of the ﬂuid sub-problem step as

xn+1 = xn + U∆t.

(3.13)

a )TFT − F−1)(cid:105)n+1

(cid:33)(cid:35)

+

Po

P0

=

+

a (F−1

(∇0vs)T :

P0
+ (ρr − 1)

λ0 ln JF−1 + G(F−1

· vsdX +
u∗
f (xn)
∆t

ρrxn+1
∆t2
ρrxn
∆t2
(∇0vs)T :
(3.14)
P0
where σ∗
f (xn) denotes the viscous stress of the ﬁctitious ﬂuid inside the solid domain. The
velocity ﬁeld u∗

f (xn)(cid:105)

λn · vsdxn,

· vsdX

+ Fr

dX

g
g

∗

f is interpolated from the Eulerian mesh via
u∗
f (x)δ(x − X)dx.

u∗
f (X) =

(3.15)

2. Solid sub-problem (for xn+1)

Þ
Þ
Þ

(cid:34)

Þ
(cid:104)(Fn)−1 · σ

(cid:104)
(cid:32)xn − xn−1
Þ

∆t2
dX −

Pn

Þ

Ω

26

This sub-problem has been re-formulated to enhance the computational robustness at ρr = 1
without sacriﬁcing accuracy [104]. A Lagrangian ﬁnite element method is used for the
solution of this problem. We use the eight-node brick element (i.e., tri-linear interpolant) for
the spatial discretization of the solid conﬁguration and the Lagrange multiplier. In equation
3.14, the penalty function term (involving the material expansion modulus) is integrated
using Gaussian rule with only one Gaussian point, and all other terms except the Lagrange
multiplier term are integrated with 8 Gaussian points (2 × 2 × 2). The Lagrange multiplier
term is integrated using the trapezoidal rule. The resulting non-linear algebraic equations
are solved with Newton iterations and the linearized equations in each iteration are solved
with the conjugate gradient iterative method. Note that Fa is a function of the reference
conﬁguration and time, and is independent of the current conﬁguration, thus it is not more
diﬃcult to derive the Jacobian matrix for the active strain model compared to the original
passive deformation model.

3. Lagrange multiplier sub-problem (for un+1

and λn+1)

Þ

Ω

(cid:169)(cid:173)(cid:171)un+1
f − u∗
(cid:32)
Þ
∆t

f

Þ

f

(cid:170)(cid:174)(cid:172) · v f dx =

f − xn+1 − xn
un+1

∆t

Pn

· ζ dxn = 0.

(λn+1 − λn) · v f dxn,

Pn

(cid:33)

(3.16)

(3.17)

This sub-problem is solved with the direct-forcing scheme [107, 106], instead of the original
Uzawa iterations [104]. It has been shown that two schemes produced the same results [106].

3.1.3.3 Validation

Considering that no benchmark problem for the active strain model is available and the original FD
method for the passive deformation model has not been validated with quantitative comparisons,
we verify the FD method with two benchmark problems: i) ﬂow-induced vibration of a ﬂexible 2D
beam that is attached to the downstream side of a stationary cylinder, and ii) ﬂow-induced ﬂapping
motion of a 3D ﬂag.

27

Figure 3.2: Time evolution of the displacements in both x− and y−directions of the central point
on the plate trailing edge (point A, see inset of (a)) of a thin ﬂexible plate in uniform ﬂow.

We ﬁrst consider the case of a 2D ﬂow past a ﬂexible beam that is attached to a rigid cylinder; see
inset in Fig. 3.2(a). The cylinder has a diameter D, and the attached beam has a length L = 3.5D and
a thickness h = 0.2D. The ﬂuid domain is Lx × Ly = 11D×4.1D We impose are no-slip conditions
on the top and bottom, and a far-ﬁeld condition (i.e., ∂u/∂x = ∂p/∂x = 0) on the right boundaries
of the computational domain. On the left boundary, we specify a parabolic velocity proﬁle with an
average velocity U0 as: U(y) = 6U0y(Ly − y)/L2
y. In this case, we choose Re = 100 by choosing
the velocity scale Uc = U0. For elastic material, its rigidity is chosen as E∗ = E/(ρ f U2
0) = 1400.
The density ratio between the solid and the ﬂuid phase is chosen as ρs/ρ f = 10. The ﬂuid mesh is
discretized into uniform grids in both directions with a resolution of 1024 × 256, while the solid
domain is discretized into brick elements with a resolution of 210 × 12. The time step is chosen
as ∆t = 0.002. The time-varying displacements in both directions of the central point on the free
end of the beam (point A) is shown in Fig 3.2, and is in an excellent agreement with the previous
studies. The amplitudes of the y-coordinate of the free end, the Strouhal number and the average

drag coeﬃcient (Cd = (cid:104)Fx(cid:105) /(cid:16) 1

(cid:17)

2 ρ f U2
0 D

) are also summarized in Table 1.

We then consider the second case of a 3D ﬂexible plate of size (L, h,W) = (1,0.01,1) in a
uniform ﬂow whose leading edge is clamped (h is the plate thickness). The computational domain
is rectangular, and its streamwise, transverse and spanwise directions are denoted by X, Y, and

28

timexA02468105.75.85.96PresentTurek & Hron(2006)(a)timeyA024681011.522.53PresentTurek & Hron(2006)(b)Table 3.1: The amplitude of the central point on the free end and the Strouhal number and the
average drag coeﬃcient for a ﬂexible beam

Amplitude Strouhal number Drag coeﬃcient

Works
Turek and Hron (2006) [97]
Tian et al. (2014) [94]
Present

0.83
0.78
0.81

0.19
0.19
0.19

4.13
4.11
4.10

(cid:17)

ρ f L

mass ratio ρsh/(cid:16)

to be 1.0, the dimensionless bending rigidity Eh3/[12(1 − ν2

Z, respectively. Currently the domain size is chosen as (Lx, Ly, Lz) = (8L,8L,2L) based on the
plate length L. The uniform velocity U0 is imposed on the inlet, top and bottom boundaries,
respectively. The far-ﬁeld and periodic conditions are imposed on the outlet and the spanwise
direction, respectively. The ﬂuid mesh resolution is 512 × 512 × 128, the number of the solid
elements is 80 × 4 × 80, and the time step is ∆t = 0.001. Here we choose Re to be 200, the
0 L3]
to be 0.0001. The Young’s modulus is computed correspondingly for the given plate size and
the Poisson’s ratio (νs = 0.4), and then the shear and volume moduli can be obtained via the
relationships G = E/[2(1 + νs)] and λ0 = E/[3(1 − 2νs)]. Again, in Fig. 3.3 and Table 2, our
results of both tip displacement and time-averaged drag coeﬃcients show good agreements with the
previous studies by Huang and Sung [43], and by Tian et al. [94]. In addition, we have examined the
impacts of diﬀerent interpolating functions on numerical accuracy by choosing either the tri-linear
or the smoothed 2-point function [103],

s)ρ f U2

φ(r) =

3/4 − r2
9/8 − 3|r|/2 + r2/2

0

|r| ≤0.5
0.5< |r| ≤1.5
1.5< |r|

,

(3.18)



and found that they indeed yield very similar results.

3.2 Results and Discussion

Note that slender biological swimmers (e.g., ﬁsh and jellyﬁsh) not only adopt some unique
shapes to enhance thrust generation but also distribute muscle inhomogeneously across their thin

29

Figure 3.3: Time evolution of the transverse displacement of the central point on the plate trailing
edge of a thin ﬂexible plate in uniform ﬂow.

Table 3.2: The amplitude of the central point on the plate trailing edge and the Strouhal number
and the average drag coeﬃcient for a ﬂag ﬂapping in uniform ﬂow.

Amplitude Strouhal number Drag coeﬃcient

Works
Huang and Sung (2010) [43]
Tian et al. (2014) [94]
Present(tri-linear function)
Present(smoothed 2-point function)

0.780
0.806
0.762
0.760

0.260
0.266
0.266
0.260

0.543
0.579
0.586

body to facilitate generation of desired deformation. Therefore, in order for soft robots to achieve
eﬀective swimming motions, it is critical to design appropriate swimming strategies (i.e. swimming
gaits) by taking into account their geometrical characteristics, material properties, and deformation
generation due to the applied active strains. In the following, we show a few successful designs of
soft robotic swimmer prototypes that undergo undulatory swimming and jet propulsion by applying
active (contracting) strains on diﬀerent locations of thin elastic plates with diﬀerent shapes. To
seek certain optimized swimming strategies, we systematically explore the parameter spaces by
performing direct simulations to resolve the nonlinear ﬂuid/elastic-structure interactions at the
moderate Reynolds numbers regime where both the viscous and the inertia forces are important.
In all simulations, we ﬁx the density ratio ρr = 1 for light-weight structures in both 2D and 3D.

30

timeyA02468-0.500.5PresentHuang & Sung (2010)Tian et al. (2014)Figure 3.4: (a) Schematic of a 2D beam under contracting actuation. The active segment is
highlighted in orange. The microscopic mechanical model of the active contraction is illustrated
by the zoom-in schematic on bottom. (b) Either a constant (i.e., a step-like) or a sinusoidal
contraction is applied during the ﬁrst half period 0 ≤ t < T/2 on the two surfaces alternatively.
Then the actuation is turned oﬀ for the second half T/2 ≤ t < T. (c) Convergence study for the 2D
beam under a sinusoidal actuation by varying domain size, mesh resolution, and time step. The
other parameters are Re = 500, G = 1000, δ = 0.4, α = 0.15, and T = 2.

3.2.1 Two-Dimensional Case

We begin by manipulating a 2D beam of length L = 1 and a uniform thickness h = 0.03 as shown
in Fig. 4.1(a). Here we consider undulation swimming motion for slender biological swimmers that
perform forward locomotion by generating wavy deformations to produce thrust force. As suggested
by the simulation results that use the immersed boundary method [36, 98], following the classical
Hill muscle model by distributing three-element spring units on the top and bottom surfaces
can eﬀectively generate oscillating-switching bending deformations which resemble undulatory
swimming motions. Without using complicated models to mimic the biological tissues’ mechanical
properties (e.g., viscoelasticity, electrophysiology eﬀects [100]), we assume the soft active material
is hyperelastic which, at the microscopic level, is driven by a contracting element with a length l and
= 1− λ0 locally. It connects with
an initial length l0, yielding an eﬀective contractile strain λa = l
l0
a neo-Hookean spring that generate elastic stresses in response to the active input of the contraction
ﬁeld which, in 2D, is simply chosen as a homogeneous ﬁeld Fa = diag
applied on the

λa, λ−1

(cid:16)

(cid:17)

a

31

time|u|0246800.10.20.30.40.50.6∆x=1/256, dt=2.5×10-4,  Lx×Ly=16×2∆x=1/512, dt=1.25×10-4, Lx×Ly=16×2∆x=1/256, dt=2.5×10-4, Lx×Ly=16×4(c)active segment δ.

Figure 3.5: Swimming motions of a 2D elastic beam under the constant and sinusoidal actuation
when ﬁxing Re = 500 and T = 2.0. (a) Instantaneous velocity as a function time. (b) Strain
energy E as a function of time. Inset in (b): Envelope trajectories of the plate midline during one
actuation period. Snapshots of the corresponding vorticity ﬁelds are shown in panels (c)
(constant) and (d) (sinusoidal).

Figure 3.6: Comparisons of the free swimming motion of a 2D elastic beam at Re = 100 and
Re = 500 when ﬁxing T = 1.0 for sinusoidal actuation following Eq. (3.20). (a) Instantaneous
velocities as a function time. (b) Mean elastic energy densities e as a function of time. Inset in
(b): Envelope trajectories of the plate midline during one actuation period. Snapshots of the
corresponding vorticity ﬁelds are shown in panels (c) at Re = 100 and (d) at Re = 500,
respectively.

(cid:16) 2π

We follow the geometry adopted by Hamlet et al. [36] by choosing the ﬁrst 10% of body length
from the left to be passive, and then connecting an active section of length δ. To activate the
beam, as shown in panel(b), we apply the constant (λa ∼ 1− α) or time-dependent (e.g., sinusoidal
λa ∼ 1 − α sin
) contractile strain ﬁeld alternatively on the both sides within a time period
T. Here α is an amplitude which characterizes the contraction strength.
In the meantime, the
maximum strain is imposed on the surface with an exponential decay in the thickness direction.

(cid:17)

T t

32

Mathematically, it is convenient to deﬁne λa for a constant actuation as

λa =

1 − α exp

1 − α exp
1 − α sin(2πt

1 − α sin(2πt

T ) exp
T ) exp

d0

(cid:17)
(cid:16)− h−y
(cid:17)
(cid:16)− y
(cid:16)− h−y
(cid:16)− y
(cid:17)

d0

d0

d0

(cid:17)

Or,

λa =

0 ≤ t ≤ T/2
,
, T/2 < t ≤ T.

0 ≤ t ≤ T/2,
,
, T/2 < t ≤ T

(3.19)

(3.20)

for a sinusoidal actuation. In the above, d0 controls the steepness of the decay, and is chosen as
h/3 in all simulations in both 2D and 3D. Typical resultant bending deformations are illustrated in
Fig. 4.1(a) by a few snapshots of the midline position (grey dashed lines). As shown in Fig. 4.1(c),
we have performed a convergence study for the swimming beam under a sinusoidal actuation by
varying domain size, mesh resolution, and time step. In the following simulations, we have chosen
the Eulerian domain to be Lx × Ly = 16 × 2, which is uniformly partitioned in the x − y plane with
the grid size ∆x = 1/256. The solid domain is partitioned by using rectangular brick elements with
the resolution 160 × 8 so that the area ratio of the Eulerian element to the Lagrangian element is
about 1.5 [104].

In Fig. 4.3, we compare the swimming motions due to diﬀerent actuation methods (also see
movie S1) when choosing Re = 500, G = 1000, δ = 0.4, α = 0.15, and T = 2. We choose
the computation mesh Panel(a) shows the magnitude of the instantaneous velocity, i.e., |u|, of
the center-of-mass (COM) positions. The spikes suggest sudden changes in the swimming speed
due to the way that the actuation abruptly switches between the constant actuation applied two
surfaces, while the velocity varies smoothly when the two sides are under a periodic sinusoidal
actuation. When projecting the velocity on its swimming direction, the two actuation methods result
in approximately the same time-averaged swimming speed (cid:104)U(cid:105) ≈ 0.13. In panel(b), we calculate
is the energy
density for incompressible materials in terms of the ﬁrst invariant of the elastic stress tensor [71].
Apparently the sinusoidal actuation requires much gentler elastic deformations compared to the
constant one. The snapshots in panels(c) and (d) highlight the drastic diﬀerences of the wake

the instantaneous elastic strain energy E(t) =Þ

(cid:16)Fe · FT
e − I(cid:17)

P wdV where w = 1

2Gtr

33

structures generated the diﬀering actuation schemes as characterized by the envelop trajectories of
the midline (see Inset of panel(b)). As highlighted by the black solid lines in the inset of (b), the
beam exhibits larger deformations under a constant actuation. A pair of vortices of opposite sign is
generated during each tail upward/downward beating, leading to two vortex streets that are widely
separated; see panel(c). In comparison, a sinusoidal actuation yields a smooth swimming motion,
shedding vorticies to form a reverse von Kármán vortex street.

Figure 3.7: Measurements of the averaged swimming speed (cid:104)U(cid:105) in (a-c), and power P (solid blue
lines) and eﬃciency η (dashed blue lines) in (d-f) of an swimming beam under a sinusoidal
actuation following Eq. (3.20) when ﬁxing Re = 500, α = 0.15, and varying T, G and δ,
respectively.

Next, we examine the Reynolds number eﬀect by simulating the swimming motions under a
sinusoidal actuation when choosing computational parameters δ = 0.4, α = 0.15, T = 1, G = 1000
at Re = 100,500 respectively. As shown in Fig. 3.6(a), the swimming plate moves faster as inertia
increases even if the elastic deformations at the two Reynolds numbers are similar as shown in
panel(b). Panel(c) suggests that at Re = 100, the viscous eﬀect still dominates without signiﬁcant
ﬂow separation and vorticity generation. The undulating deformation is useful to overcome the
viscous eﬀect by taking advantage of the anisotropic drag forces exerted upon the beam. As shown

34

T<U>0.511.522.500.050.10.150.20.25(a)Re = 500α = 0.15G = 1000δ = 0.4G<U>05001000150000.050.10.15(b)Re = 500α = 0.15T = 2δ = 0.4δ<U>0.10.20.30.40.50.600.060.120.18(c)Re = 500α = 0.15T = 2G = 1000TPη0.511.522.500.10.20.300.01(d)Re = 500α = 0.15G = 1000δ = 0.4GPη05001000150000.020.040.0600.005(e)Re = 500α = 0.15T = 2δ = 0.4δPη0.10.20.30.40.50.600.100.004(f)Re = 500α = 0.15T = 2G = 1000in panel(d), apparently the inertial eﬀect at Re = 500 enhances the momentum exchange to generate
a reverse von Kármán vortex street, which eﬀectively generates thrust forces to propel the swimmer.
After choosing a particular actuation method (e.g., sinusoidal), we may study the 2D undulating
swimming motion by systematically exploring the parameter space of (Re, α,T, G, δ). Since simul-
taneously optimizing all these parameters can be very challenging, here we examine the diﬀerences
of the 2D swimming motions by varying individual parameters. As shown in Fig. 3.7(a)-(c),
we vary parameter T, G, or δ individually to seek the fastest moving speed (cid:104)U(cid:105) and maximum
swimming eﬃciency. In panel(a), such motion occurs when T ≈ 1.3 and the other parameters are
ﬁxed. Asymmetric and non-periodic undulations have been observed when actuation is either very
fast (T < 0.5) or slow (T > 2.5), which barely generates a directional motion. In panel(b), when
varying G only, the corresponding fastest swimming motion occurs at G ≈ 900. Interestingly, we
have observed that at small G (i.e., G < 200), while swimming slowly, the beam in fact can move
back and forth as actuation switches sides during one period. When G becomes larger, the beam
exhibits a directional motion with its COM velocity sharply increasing before saturating around
U = 0.13 at G > 500. In panel(c), the maximum speed occurs when approximately half of the
body is actuated. As δ goes beyond 0.6, we have again observed that non-periodic and asymmetric
swimming motions.

In panels(d-f) we estimate the power done by the beam upon the ﬂuid as [99]

Þ

Þ

Ω f

(3.21)

(3.22)

P =

1
2

d
dt

|u f |2dV +

2
Re

Ω f

|∇u f |2dV,

based on which we deﬁne the swimming eﬃciency as

Ek

,

E f + Ek

η =

2 M f (cid:104)U(cid:105)2 is the kinetic energy of the elastic beam, and E f =Þ T

where Ek = 1
0 Pdt represents the
total energy delivered to the ﬂuid due to the mechanical work performed by the swimming motion.
We have found that P varies approximately monotonically with the varying parameters, leading to
the maximum eﬃciencies around T = 1.1, G = 400 and δ = 0.3, respectively.

35

3.2.2 Three-Dimensional Cases

3.2.2.1 Deformation Control

Figure 3.8: (a) Examples of diﬀerent actuation patterns (horizontal vs. diagonal) on the top
surface of a square plate, and the corresponding bending deformations in (b) and (c).

Figure 3.9: Snapshots of circular plate under (a-c) radial and (d-e) azimuthal actuation.

In 3D, the active strain ﬁeld can be deﬁned as Fa = diag(λ1, λ2, λ3), deﬁned in the principal
coordinates in the un-deformed conﬁguration. We align λ1 = λa < 1 with the local contracting
direction while deﬁning the remaining principal stretches as λ2 = λ3 = λ
. This assumes
that the deformation is homogeneous in the directions of λ2 and λ3. However, design and control
of elastic deformations in 3D can prove to be much more complicated than in 2D, especially
when dealing with complex geometries. As demonstrated in Fig. 3.8(a)-(c) for a rectangular plate,

−1/2
a

36

Figure 3.10: Snapshots of an oblate bell under (a-c) radial and (d-e) azimuthal actuation.

the resultant bending deformation can be adjusted by distributing the constant contracting ﬁeld
of λ1 along arbitrary directions (e.g., horizontal in (b) and diagonal in (c)) in the un-deformed
conﬁguration.

Another example is shown in Fig. 3.9, where we change the plate shape from rectangular to
circular, and apply isometric actuation in either radial or azimuthal direction. Apparently, when
radially activated (i.e., λaer, see panel(a)), the plate apparently experiences a two-step deformation
during which an intermediate equilibrium shape occurs in panel(b) due to a wrinkling instability,
and then relaxes and switches to become a steady-state saddle shape in panel(c). Active deformation
involving complex transient dynamics like this is diﬃcult to control, and hence the radial actuation
is not preferred. On the other hand, when a constant contraction is applied in the azimuthal direction
(i.e., λaeθ, see panel(d)), the circular plate can generate stable reversible bending deformation as
shown in panels(e,f). Similar actuation strategies may be used for other thin soft structures with
azimuthal symmetry. As shown in Fig. 3.10, when a radial actuation is applied on the inner surface
of an oblate bell (see panel(a)), an opposing expansion occurs in the azimuthal direction due to
the material incompressibility, driving a transient wrinkling deformation (see panel(b)). Then the
whole body relaxes to become a steady ﬂattened disk as shown in panel(c). Alternatively, as shown

37

in Fig. 3.9(d-f), applying an azimuthal actuation results in eﬃcient contraction of the entire-body,
which, as shown below, can eﬀectively be used to generate reversible contraction motions to power
a jellyﬁsh robotic swimmer.

3.2.2.2 Undulatory Swimming

Figure 3.11: Snapshots of 3D vortex-core structures (a) and the projection of 2D vortices on the
x − y plane (b) of an elastic plate that performs undulation under a sinusoidal actuation following
Eq. (3.20). The computation parameters are chosen as α = 0.2, δ = 0.5, T = 1.0, G = 1000, and
Re = 1000. (c) Instantaneous swimming velocity as a function of time. The color bar represents
(red line)

the magnitude of the 3D vorticity, i.e., |ω f | (b) Drag coeﬃcient Cd = Fy/(cid:16) 1
and wake-induced force (i.e., Lamb vector) coeﬃcient Cw = −FL/(cid:16) 1

(cid:17)
2 (cid:104)U(cid:105)2 A

functions of time. (e) Strain energy E as a function of time.

(cid:17)

2 (cid:104)U(cid:105)2 A

(blue line) as

Similar to the 2D cases discussed in Section 3.2.1, we ﬁrst examine the undulatory swimming
motion of a 3D rectangular plate. Here we adopt the 2D shape in Fig. 4.1(a), and extrude the
beam of length L and thickness h along its third dimension for a certain width W. In Fig. 3.11, we
show a case of a swimming plate of size (L, h,W) = (1,0.03,0.25) that undergoes a quasi-steady

38

Figure 3.12: Simulation results of an undulating swimming motion by varying δ while ﬁxing
α = 0.2, T = 1.0, G = 1000 and Re = 1000. (a) Instantaneous swimming speed U at
δ = 0.3,0.5,0.7, respectively. (b) Time-averaged swimming speed (cid:104)U(cid:105) and strain energy (cid:104)E(cid:105) as a
function of δ for periodic locomotion. (c) Power P and eﬃciency η as a function of δ.

undulatory swimming motion, due to sinusoidal contractions (Eq. (3.20)) applied alternatingly on
both sides. Again, a uniform mesh is used for the Eulerian domain Lx × Ly × Lz = 2.8 × 2.8 × 0.7
with the resolution 256 × 256 × 64. The Lagrangian domain L × h × W = 1 × 0.03 × 0.25 is
partitioned into brick elements with the solution 84 × 3 × 21, leading to the volume ratio (i.e.,
Eularian vs. Lagrangian) about 1.0. The other computation parameters are chosen as T = 1.0,
δ = 0.5, G = 1000 and Re = 1000.

In Fig. 3.11(a), we visualize the 3D vortex structures by using the so-called λci criterion [109,
− ∂v
12], together with the distribution of the vorticity component ωz = ∂u
∂x on a 2D cross-section on
∂ y
the (x-y) plane in panel (b). In contrast to the 2D case where the shedding vorticies form a reverse
von Kármán vortex street, here ﬂow separation generates two toroidal structures during each beating
motion, leading to two arrays of vortex rings. Apparently, the dynamics quickly become periodic as
(cid:17)
shown by the COM velocity in panel(c), as well as the strain energy in panel(e), suggesting a quasi-
2 (cid:104)U(cid:105)2 A
in
2 (cid:104)U(cid:105)2 A
,
where A = LW measures the dimensionless area of plate. Here the drag force is deﬁned by net
Ωc ω f × u f dV from

steady free swimming motion. In panel(d), we show the drag coeﬃcient Cd = Fy/(cid:16) 1
the swimming (y) direction, as well as the wake-induced force coeﬃcient Cw = −FL/(cid:16) 1
force exerted on the deformable body. The wake-induced thrust force FL =Þ

the vortex structures behind the swimmer, where ω f × u f =
vector [101]. We have found that the time-averaged net force is approximately zero, i.e., (cid:104)Cd(cid:105) = 0,

(cid:17)
(cid:17) × u f is the so-called Lamb

(cid:16)∇ × u f

39

timeU024600.20.40.6δ = 0.3δ = 0.5δ = 0.7(a)δηP00.20.40.600.0060.0120.01800.010.020.030.040.05ηP(c)δ<E><U>00.20.40.600.050.10.1500.10.20.30.40.5(b)<E><U>since the active strain inputs do not generate extra forces. The oscillations of the Cd − t curve
are largely due to the abrupt switches of the actuation on the two sides. In contract, the vortex
rings induce negative hydrodynamic forces that vary periodically in the y direction to push ﬂuid
backward, which hence contributes to the thrust force generation for the swimming plate.

In Fig. 3.12 we show an example of varying the active segment length δ using a sinusoidal
actuation when choosing α = 0.2, T = 1.0, G = 1000, and Re = 1000. According to the
instantaneous swimming velocity U in panel(a), the plate typically starts a steady-steady navigation
after performing 5-6 strokes, with a signiﬁcant increase in the speed when δ goes beyond 0.3. In
panel(b), it is observed that the time-averaged swimming speed U increases with δ, and reaches
a maximum value when δ ≈ 0.5. The time-averaged strain energy (cid:104)E(cid:105) grows monotonically as δ
increases, driven by more inputs from the contracting strain ﬁeld through Fa. In panel(c), compared
to the 2D case where the swimming power P grows monotonically with δ, it approximately saturates
around δ = 0.6, leading to the maximum eﬃciency at δ ≈ 0.5.

3.2.2.3

Jet-propulsion Swimming

Next we examine jet-propulsion based swimming (e.g., jellyﬁsh) that propels the swimmer by
squirting out water from the muscular cavity through periodic contractions. We follow the design
reported in previous studies [85, 38, 1, 73] to construct an oblate 3D bell with its initial shape a
partial prolate ellipse which can be described as

(cid:19)2

(cid:18) x

R1

(cid:19)2

(cid:18) y

R2

(cid:19)2

(cid:18) z

R3

+

+

= 1

(3.23)

where R1 = R3 = 0.5 and R2 = 0.25, with a thickness h = 0.03. As discussed in Section 4.2, in
order to generate robust whole-body contractions for a thin circular or bell-like plate, it is preferred
to apply an active strain ﬁeld with the principal contraction occurring in the azimuthal direction to
avoid complex transient dynamics and instabilities of thin structures. We employ the same strategy
here to apply the active contractions on the outer region of the inner surface marked by the blue
circles in the inset of Fig. 3.13(a). The active segment is taken approximately 1/3 of the total

40

Figure 3.13: Typical vortex structures of an elastic bell undergoing jellyﬁsh-like swimming in the
moving coordinates, when choosing α = 0.8, T = 2.0, G = 100, Re = 500, ε = 0.1, δ = 0.21, and
h = 0.03. (a) Vortex core overlaid on the solid mesh. The color bar represents the magnitude of
the 3D vorticity, i.e., |ω f |. Inset: Schematic of the swimming elastic bell with a thickened top and
azimuthal actuation. (b-i) In-plane ﬂuid velocity vector ﬁeld (u, v) overlaid on the colormap of
vorticity on a 2D cross section plane during one period.

41

Figure 3.14: Highlights of bell’s rowing motion accompanying whole-body vibrations through the
kinematics of the bell arm during the contracting and reﬁlling process (a). The instantaneous
strain energy (b) and COM velocity(c) as function of time.

Figure 3.15: (a) Phase diagram of (α,Re) for a swimming bell undergoing jet-propulsion
locomotion. The while and black symbols represent the forward and backward locomotion,
respectively. The background color represents the time-averaged steady-state swimming velocity
(cid:104)U(cid:105) at late times. (b) Typical instantaneous COM velocities when choosing α(cid:48)s in diﬀerent
regimes at Re = 250.

42

timeU051015-0.200.20.40.60.81(c)timeE051015051015(b)(a)t = 2.0t = 0refilling(a)∆t = 0.1contractingarm-length of the bell. In addition, we thicken the bell top (denoted by ε, and we ﬁx ε = 0.1) in
order to stabilize the elastic deformation [86, 75], which is important when both α and Re are large.
Here we employ a constant contraction following Eq. (3.19), instead of a gentle sinusoidal one
which is not eﬀective in generating rapid deformations as required by jet propulsion. Here we choose
parameters based on the in vivo measurements of jellyﬁsh tissue (e.g., mesoglea) [26]. By using the
Young’s modulus G ≈ 100Pa (the experimental data are approximately between 10Pa and 100Pa)
as well as choosing the characteristic length Lc ≈ 1cm [70] and velocity Uc ≈ 1cm/s, we choose
the estimated the range of the dimensionless parameters as G = 10 ∼ 100 and Re = 10 ∼ 1000. An
example of bell swimming through jet-propulsion is shown in Fig. 3.13. The periodic axisymmetric
contractions of the oblate elastic bell is reminiscent of the “rowing” motion of a jellyﬁsh [16, 17,
18]. As shown in panel(a) for the 3D vortex cores overlaid on the solid mesh, after the initial
transient, the elastic bell quickly starts to perform a stable forward swimming motion as the bell
arm paddles, leaving behind toroidal vortices [19, 20, 73, 75]. The typical vortex dynamics during
one period are illustrated in sequential snapshots of a 2D cross section shown in panels(b-i). In
the power stroke corresponding to the ﬁrst half-period (0 < t < T/2, see panels(b-e)), the applied
active strains drive an inward contracting deformation. A shear-induced starting vortex is generated
at the bell tip, interacting with the stopping vortex from the previous recovery stroke. After the
starting vortex pair is separated from the bell arm, it is quickly advected with the jet ﬂow due to
the entrained ﬂuid mass that is ejected downward (panels(c-e)). In the 2D plane, it is visualized as
a vortex dipole. In the next half-period (T/2 < t < T, see panels(f-i)), turning oﬀ the active strains
relaxes the elastic bell, driving a recovery stroke during which the ﬂuid mass reﬁlls the bell. In the
meantime, a stopping vortex with an opposite sign is generated and then adhered on the bell tip in
the expansion process.

Details of the bell kinematics during one period is shown in Fig. 3.14. In panel(a), the sequential
snapshots of the swimmer’s contour are plotted. Note that while actively interacting with the ﬂuid
ﬂows, imposing a large step-like contracting deformation of thin soft structures may easily trigger
resonances during which the entire body oscillates quickly. For example, the bell top bounces back

43

and force during contraction (also see Fig. 3.13(b-e) as well as movie S2), and then overshoots
the equilibrium position and swings back during relaxation. Such oscillatory behaviors can be
quantitatively characterized by examining the variation of the total strain energy as a function of
time in panel(b) where small amplitudes ﬂuctuations are observed, immediately after the active
strains being turned on and oﬀ. This is also reﬂected on the instantaneous average velocity in
panel(c). As the bell starts to contract, its COM velocity quickly reaches a maximum, and then
relaxes following the structural oscillations.

Systematically varying α and Re reveals various diﬀerent regimes of jellyﬁsh-like swimming
motion as a result of the complex ﬂuid-thin structure interactions. The dynamic behaviors are
summarized in the phase diagram in Fig. 3.15(a) where the background colors represents the
time-averaged COM velocities at steady state. As explained above, in the regime of large α and
Re, the bell exhibits a fast forward swimming motion (see COM velocity marked by the red line
in panel(b)), due to the strong positive feedbacks obtained from the fast ﬂuid jet that carries ﬂuid
momentum backward.
In the regime of smaller α and lower Re, however, while the bell arms
are rowing, the locomotion is signiﬁcantly weakened, largely due to the reduced amount of inputs
from active contractions which hence leads to much weaker elastic deformations. Therefore, we
ﬁnd that the contracting deformation and the induced whole-body vibration are comparable. Their
interplay eventually leads to complex arm kinematics as well as the momentum exchange between
the ﬂuid and the solid phase. Interestingly, we have observed that while swimming slowly (see the
COM velocities marked by the blue and black colors in panel(b)), the bell can exhibit both forward
(marked by solid squares) and backward (marked by black triangles, see movie S3) movement.

3.3 Concluding Remarks

In this work we have developed a FD-based computational framework to simulate swimming
motion of soft robots by distributing active contracting stains on the elastic solid phase. The active
strain model is based on the multiplicative decomposition of the deformation gradient tensor. As
a ﬁrst order approximation, the induced active deformation is independent of mechanical stress,

44

and can be easily measured experimentally by varying the external stimuli without imposing any
mechanical load. From the biomimetic perspective, the active strain approach makes it straightfor-
ward to distribute active contractions in the muscle-ﬁber directions by following biological muscle
architectures. Moreover, it is useful to take advantage of ﬁnite thickness of thin elastic structures
by imposing an inhomogeneous contraction ﬁeld to achieve desired elastic deformation.

Figure 3.16: Comparison of the dynamics of soft robotic baby-jellyﬁsh by distributing active
contractions (marked by the blue lines) on (a) arms (i.e., radially) and partial of the circular disk
(i.e., azimuthally), and (b) arms only. The four sequential snapshots in (a) highlight a stable
directional locomotion in one actuation period. The snapshot in (b) shows an example of unstable
asymmetric deformations without applying the azimuthal contraction in the disk to stabilize the
structure.

Using the FD/active-strain method, we have been able to virtually design 2D/3D thin, light-
weight soft robots with simple geometries (e.g., square, circle, bell) that perform undulation and
jet-propulsion that are commonly employed by biological swimmers.
It is also straightforward
to activate soft swimmers of more complex geometries by distributing active strains with certain
strength and direction on the desired locations. As shown in Fig. 3.16 and movie S4, we have
followed the design by Nawroth et al. [70] to simulate a baby jellyﬁsh (e.g. early-stage Ephyrae)
that has a circular-disk body and multiple attached “arms”. In particular, we have found that it
is important to apply the active contracting stains to follow the muscle architecture observed in
jellyﬁsh. As shown in Fig. 3.16(a), we make the contractions occur in the radial direction to be
aligned with the arms, together with a layer of contraction simultaneously imposed in the azimuthal

45

direction near the rim of the circular disk, which turns out to be critical in stabilizing the thin-
structure deformation in forward swimming. A typical failure case is shown in Fig. 3.16(b) (also
see movie S4) where only the radial contractions are applied.
It is seen that the baby jellyﬁsh
cannot generate a clear directional locomotion due to very unstable deformations as shown by the
snapshot. We have also tried varying the area covered by radial actuation but without success (not
reported here).

While we only consider some relatively thin elastic structures in this paper, the FD/active strain
model can be generally applied to manipulate elastic structures of arbitrary shapes and sizes. In
particular, it is important to incorporate more elaborate constitutive models for viscoelastic or
composite materials in order to mimic the behaviors of real biological tissues. Moreover, the
active strain model can be ﬂexibly implemented in other computational frameworks, either using
boundary conforming methods or Cartesian grid methods, to handle coupling of soft actuation and
the induced ﬂuid motion. Last but not least, in order to optimize soft robotic systems with a large
number of parameters (e.g., shape, stiﬀness, and actuation distribution) and possibly multiple design
objectives (e.g., swimming speed, power, and eﬃciency), it is desirable to couple the existing CFD
tools with multi-objective evolutionary algorithms and packages such as the Uniﬁed Non-dominated
Sorting Genetic Algorithm (U-NSGA) II/III [22, 88, 21, 44]. The integrated framework like this
will facilitate design of high-performance soft robotic machines for a wide range of applications in
engineering and applied sciences.

46

CHAPTER 4

CFD-BASED MULTI-OBJECTIVE CONTROLLER OPTIMIZATION FOR SOFT

ROBOTIC FISH WITH MUSCLE-LIKE ACTUATION

In this chapter, we develop a holistic computational framework for the design, control, and opti-
mization of a soft robotic swimmer. With the aid of the high-ﬁdelity FD/active-strain simulator, we
implement feedback control strategies for the fully-coupled ﬂuid/elasticity systems under actuation.
Here we model the robotic ﬁsh as a 2D swimming elastic beam of a ﬁnite thickness, with contractive
strains being imposed on two sides (with sharp decays in the thickness direction) periodically. When
coupled with ﬂuid ﬂows, we demonstrate that the resultant undulatory free-swimming motions can
be tuned by changing the active strain magnitude and frequency. Moreover, turning motions of the
swimmer can be eﬀectively accomplished by imposing asymmetric active strains on the two sides
of the body. To achieve feedback control of the swimming motion, we couple the FD simulator
with model-free proportional-integral-derivative (PID) control, which, arguably, is the most widely
used feedback control scheme [82]. We seek optimal parameters of the PID controllers by using
a genetic-algorithm-based multi-objective evolution method, U-NSGA-III [88], which produces a
Pareto set of optimal solutions that clearly illustrate the tradeoﬀs among the objectives, with which
the designer is able to make informed decisions.

To put the study in context, we consider a series of tasks where the soft swimmer tracks a
moving target involving diﬀerent trajectories and velocities, where a pair of PID controllers are
used to adjust the amplitude and the bias of the active contractile strains imposed on the swimmer.
The parameters of the PID controllers are optimized with three objectives: (1) tracking error, (2)
cost of transport, and (3) the average elastic strain energy of the morphing body. In particular,
the tracking error refers to the error between the swimmer position and the moving target position
and is thus a metric on tracking performance, while the cost of transport represents a measure of
locomotion eﬃciency. The elastic strain energy is also of importance because it is not only closely
related to the required performance of the actuator, but also to the rate of fatigue failure to which

47

soft materials can be prone. The optimization is repeated at several Reynolds numbers for cases
with diﬀerent velocities and trajectories of a moving target. The resulting Pareto fronts of the
multi-objective optimization problem reveal the correlations and trade-oﬀs among the objectives
and oﬀer key insight into the design and control of soft swimmers. For example, the tracking
error shows strong inverse correlation with the elastic energy. As another example, the optimized
proportional gain parameters are also inversely related to the tracking error. Finally, the controller
parameters produced via the proposed multi-objective optimization method demonstrate robustness
across diﬀerent conditions for the moving target and for the ﬂuid environment.

4.1 Mathematical Model

In this section, we brieﬂy review the major components of the FD/active strain method, which
is the core of our CFD solver. The reader is referred to [58] for more details, especially on
the derivation of the FD formulation and the numerical schemes. Mathematically, the solid
deformation can be described by the deformation gradient tensor F = ∂x/∂X, which projects the
current deformed state to the initial reference state. To capture the local actuation, we decompose
the deformation gradient tensor F into an active deformation tensor Fa and an elastic deformation
tensor Fe following multiplicative decomposition [54, 55] such that

F = Fe · Fa.

(4.1)

Here Fa is eﬀectively an arbitrary function applied to the reference conﬁguration that can be de-
signed in terms of the desired location, direction, sign and strength. For incompressible solids that
are considered here, we apply an incompressible restriction on actuation such that det(F) = 1, where
“det” represents the determinant. For simplicity, Fa can be deﬁned in the principal coordinates and
transformed to the desired orientation by a standard rigid body rotation coordinate transformation.
For an artiﬁcial muscle that contracts uniaxially when activated, we deﬁne Fa = diag[λ1, λ2, λ3],
λ−1
1 > 1 are correspond-
where λ1 < 1 represents a principal compression ratio, and λ2 = λ3 =
ingly the homogeneous deformation in the other two directions. With the total F being mapped
appropriately, the resultant elastic stress can be calculated through certain constitutive relation

(cid:113)

48

(cid:16)F · F−1

a

(cid:17)

τ = τ (Fe) = τ
conditions.

for various kinds of material (hyperelastic, viscoelastic, composite, etc.)

It is well understood that during an undulatory motion, ﬁsh use the anisotropic viscous drag
force exerted in the longitudinal and transverse directions along the wavy body to create a net
propulsive force [57]. As shown in Fig. 4.1(a), we adopt a simple 2D rectangular beam of length
L and uniform thickness h to characterize the slender ﬁsh body shape. Without using complicated
interconnected spring models to mimic the biological tissues’ mechanical properties (e.g., the
artiﬁcially coupled viscoelasticity and electrophysiology eﬀects as what have been done by using
the immersed boundary method [36, 98]), we assume that the soft active material is a continuous
hyperelastic elastica which, at the microscopic level, is driven by a contracting element with an
= 1 − λ0. At the microscopic
initial length l0, yielding an eﬀective local contraction ratio λa = l
l0
scale, it connects with a neo-Hookean spring that generate elastic stresses in response to the
active input of the contraction ﬁeld. The latter, in 2D, is simply chosen as a homogeneous ﬁeld
Fa = diag

applied on the active segment.

(cid:16)

(cid:17)

λa, λ−1

a

Figure 4.1: (a) Schematic of a 2D beam under contractile actuation (view from above). The active
segment is of length δ. The microscopic mechanical model of the active contraction is illustrated
by the zoom-in schematic on bottom. (b) Contractile actuation applied alternatively on the
surfaces of two sides, which exponentially decays in the thickness (y) direction (not shown in the
ﬁgure).

Here we follow the geometry adopted by Hamlet et al. [36] by choosing the ﬁrst 10% of body
length from the left to be passive, and then connecting an active section of length δ. To activate the
beam, as shown in Fig. 4.1(b), we apply a constant (λa ∼ 1−α, for a ﬁxed α) or time-dependent (e.g.,

49

htime0(a)(b)Lyxl0l/2T1a(cid:16) 2π

(cid:17)

T t

sinusoidal, λa ∼ 1 − α sin
) contractile strain ﬁeld alternatively on both sides within a time
period T. Here magnitude α characterizes the contraction strength. In the meantime, the maximum
strain is imposed on the left (denoted by “L”) or right (denoted by “R”) side with an exponential
decay in the thickness direction. Mathematically, it is convenient to deﬁne the sinusoidal actuation

scheme asλL

(cid:16) 2πt
t ∈ (cid:104)

T
0, T
2

(cid:17)
(cid:17)

a = 1 − α sin
λR
a = 0,

(cid:16)− h−y

(cid:17)

exp

;

,

d0
a = 1 − α sin
λR

t ∈ (cid:104)
(cid:16) 2πt

0, T
2

(cid:17)

T

(cid:17)

;

exp

λL
a = 0,

(cid:16)− h−y

(cid:17)

,

d0

t ∈(cid:104)T
t ∈(cid:104)T

2 ,T

2 ,T

(cid:105)
(cid:105)

,

.

(4.2)

In the above, d0 controls the steepness of the decay, and is chosen as h/3 in all simulations in this
work. Typical resultant bent shapes are illustrated in Fig. 4.1(a) by a few snapshots of the midline
position (grey dashed lines). Therefore, continuously applying periodic contractions on both sides
lead to an undulatory motion.

To resolve the FESI of soft robotic swimmers, we implement the above active strain model using
a ﬁctitious domain (FD) method [35, 104, 58]. The key idea of the FD method is that the interior
of the solid is assumed to be ﬁlled with a ﬁctitious ﬂuid that is constrained to move at the same
velocity with the solid by a pseudo body force (via a Lagrange multiplier). The advantage of the
FD method is that the ﬂow ﬁelds can be solved eﬃciently with a ﬁxed Cartesian grid. Suppose that
a deformable body of density ρs is immersed in the incompressible Newtonian ﬂuid of viscosity
µ and density ρ f . Let Ω denote the entire computational domain containing both solid and ﬂuid
domains, and S(t) represent the solid domain. It is noted that tracking swimming objects in the
ﬁxed coordinates requires using a very large computational domain, which makes computation
expensive. Alternatively, here we employ an instantaneous inertial frame Ω that co-moves with the
swimmer at a certain reference speed U [56]. Then the dimensionless FD governing equations in
the weak form become

(cid:19)

Þ

S

λ · v f dx,

(4.3)

(4.4)

Þ

(cid:18) ∂u f

∂t

Ω

(cid:19)

+ ˆu f · ∇u f

· v f dx +

−pI +

1
Re

(∇u f )T

: ∇v f dx =

(cid:18)

Þ
Þ

Ω

Ω

q∇ · u f dx = 0,

50

where ˆu f = u f − U. The dimensionless governing equations for neo-Hookean solid material
(τ = G(B − I)) are solved in the absolute references as

· vsdx +

S
dx = −

(∇vs)T : [λ0 ln JI + G(B − I)]dx
λ · vsdx,

(4.5)

(4.6)

Þ

(cid:20)
(cid:18)

S
(∇vs)T :

Þ

S

−

(cid:18) dus

(cid:19)(cid:21)

− Fr

g
g

dt
[∇u f + (∇u f )T]

(cid:19)

(ρr − 1)
1
−pI +
Re

Þ

Þ

S

Þ

S

(u f − us) · ζ dx = 0

e = (F· F−1

a )·(F· F−1

Equations (4.3)-(4.6) represent the ﬂuid momentum equation, the ﬂuid continuity equation, the
solid momentum equation, and the velocity constraint in the solid domain, respectively. In these
equations, u f is the ﬂuid velocity, us the solid velocity, p the ﬂuid pressure, and λ the pseudo
body-force (i.e., the Lagrange multiplier). The variables v f , vs, q and ζ are the corresponding
variations, respectively. For the original FD method for the passive deformation model, B = F · FT
is the left Cauchy-Green deformation tensor, here F being the deformation gradient tensor deﬁned
as: F = ∂x/∂X, in which x and X are the current and reference conﬁgurations of the solid,
respectively. J is the determinant of F, and J =1 for the incompressible solid. In contrast, for the
active strain model studied, B = Fe · FT
a )T, where Fe is the elastic deformation
tensor which causes the elastic stress, and Fa is an input deformation tensor without generating
elastic stress. In addition, g denotes the gravitational acceleration. The following characteristic
scales are used for the non-dimensionlization scheme: Lc for length, Uc for velocity, Lc/Uc for
c/Lc for Lagrange multiplier λ. The following dimensionless
time, ρ f U2
, material parameters λ0 = λ0
control parameters are also introduced: density ratio ρr = ρs
ρ f U2
ρ f
c
. Here λ0 is related
to the compressibility property of the material and G represents the shear modulus. It should be
mentioned that for soft swimmers, the characteristic velocity is deﬁned as Uc = fcLc for a given
reference actuation frequency fc. In addition, in order to enforce the incompressibility constraint
of the solid material, i.e., J = 1, we adopt a penalty function-like approximation by setting a large
enough value for λ0 (104 ∼ 105).

, Froude number Fr = gLc
U2
c

c for pressure p, and ρ f U2

and G = G
ρ f U2
c

, Reynolds number Re =

ρ f UcLc

µ

51

4.2 Controller Design and Optimization

4.2.1 Control Scheme

Figure 4.2: Schematic illustrating the control scheme based on the translational and rotational
motion of a soft swimmer in relation to a moving target xt.

The primary goal for our control scheme is to make the soft swimmer eﬀectively follow a
moving target point. Such situations could arise in applications such as schooling (i.e., leader-
follower conﬁguration) or simply following a pre-determined desired trajectory. By its 2D nature
our soft swimmer (elastic beam) is limited to forward swimming and turning, based on which we
propose the following actuation scheme:

αR(t) = α0(t) +

2 , αL(t) = α0(t) − β(t)
β(t)
2 .

(4.7)

Here αR and αL correspond to the peak actuation strength in the right and left muscles, respectively,
viewed in the body frame of the beam oriented from the tail toward the head. This diﬀers from
previous descriptions where both sides used a single value of constant α. Here α0 is the mean
contraction strength, which is assumed to be limited to the range 0 ≤ α0 ≤ αmax.

Turning can be achieved by biasing the contraction strength between the two actuated sides.
The soft swimmer tends to turn in the direction of the stronger actuation, providing a means to
control the beam’s heading while swimming. Testing showed that this mechanism allows the beam
to turn even while it is at rest with an averaged swimming speed of zero. The bias that leads to
turning is denoted by β and is limited to the range −βmax ≤ β ≤ βmax. Based on the signs in (4.7),
a positive value of β causes the beam to turn to the right and vice versa. Note that αR and αL are

52

xcomxtrhrswimmertargetθ limited to the same range as α0 due to the physical constraints of the actuator, but the presence of β
has the potential to push their values outside of this envelope. To prevent this while still achieving
the desired bias, the portion of the bias outside the envelope is shifted to the opposite side. For
2 ≥ αmax, then αR = αmax and αL = αmax − β. This is suﬃcient assuming
example, if αR = α0 + β
that βmax ≤ αmax and in practice we typically set βmax = 1

2 αmax.

To complete the control scheme it is necessary to connect our control variables, α0 and β, to the
physical state of the swimmer. We elect to use proportional-integral-derivative (PID) controllers as
they have been widely used across disciplines. A general PID controller takes the following form:

Þ t

0

f(t) = kpe(t) + ki

e(t(cid:48))dt(cid:48) + kd

de(t)
dt

(4.8)

where f(t) is the control variable, e(t) is the error corresponding to the diﬀerence between the
target and current state, and the constants kp, ki, kd are tuning parameters. The actuation strength
α0 is calculated based on the error in position, which is deﬁned as:

e1(t) = |r|sgn(r · v) .

(4.9)
In the above, the vector v represents the moving target’s velocity, and r = xt − xcom connects
the swimmer’s center of mass (CoM) to the target position with its magnitude determining the
magnitude of error (see schematic in Fig. 4.2). On its own this error magnitude would lead to the
error being always positive which can potentially promote overshoot and lead the integral term to
blow up. To prevent this we determine the sign of the error by sgn(r · v) (here “·” represents the
vector inner product). If the target is moving away from the swimmer, r · v will be positive driving
the swimmer to accelerate. Conversely, if the target is moving toward the swimmer, r · v will be
negative causing the swimmer to slow down and potentially stop. If the controller returns a value
that is out of bounds, the control variable is set to the respective bound and integration is stopped
to prevent windup of the integral term. The overall result is a follower type scheme where the
swimmer tends to remain behind the target and overshoot behavior is suppressed.

Special attention must be given to the derivative term due to the noise in the error. This is a
common problem in PID controllers as noise can produce large and drastically varying values of

53

the derivative that can lead to undesirable behavior. The undulatory motion of the swimmer leads
to a signiﬁcant amount of sway which appears in the error as an oscillation driven at the forcing
frequency. We ﬁlter out this oscillation by calculating the derivative based on the moving-average
of the error, with the window size for averaging equal to that of the period of the forcing, T, as
shown in (4.10). In particular, the following algorithm is used to approximate the derivative term:

Þ t

t−T

de1(t)
dt

≈ ¯e1(t) − ¯e1(t − ∆t)

∆t

¯e1(t) =

,

1
T

e1(t(cid:48))dt(cid:48)

(4.10)

Directional control is achieved by determining β via a separate PID controller, where the error
is based on the angle between the target vector, r, and the swimmer’s time-averaged heading vector,
¯rh. Due to the highly ﬂexible nature of the swimmer, it is diﬃcult to clearly deﬁne a body frame.
As such we choose to deﬁne our orientation in terms of a chord line running from the tail to the head
of the swimmer, rh. The large-amplitude motion of the head and tail means that this vector will
vary signiﬁcantly across a full swimming stroke, limiting its value as an instantaneous measure.
However, time-averaging like that done in (4.10) provides a stable measure of the swimmer’s
Þ t
heading. In sum, the bias control is based on the error e2(t) = θ, the angle between the vectors r
t−T rh(t(cid:48))dt(cid:48). Here rh = xhead − xtail represents the instantaneous direction vector
and ¯rh = 1
T
that connects the head and tail. From the deﬁnition, since ¯rh is already smoothed by the averaging,
no special treatment is necessary for the derivative term.

4.2.2 Multi-Objective Optimization

Multi-objective optimization makes it possible to view the landscape of potentially optimal conﬁg-
urations. This allows for better understanding of the trade-oﬀs between various objectives, leading
to smarter design decisions. For example, in a swimmer a small decrease in speed may lead to
sharp increase in eﬃciency or vice versa. The possibility of such a trade-oﬀ would not be apparent
using single objective methods. The results of multi-objective optimization appear as a Pareto
front, which is a surface in the objective space containing all solutions where each solution cannot
be further improved in one objective without negatively aﬀecting another objective(s).

54

Our work is focused on employing high-ﬁdelity FSI simulations to investigate soft swimming
mechanisms. This limits the use of reduced-order models, making evolutionary methods the most
promising option for optimization.

Locomotion traditionally lends itself to two major design objectives, swimming speed and
eﬃciency. However, our fully coupled model allows us to consider other important design criteria
such as those related to the actuation of the swimmer. Soft materials tend to have limited lifecycles
and can experience premature fatigue failure compared to traditional actuators, if cyclic loading
conditions are not accounted for. Thus the design of such soft swimmer robots may best be done
with the consideration of 3 or more objectives.

We elect to optimize the soft swimmer by employing a multi-objective evolutionary method,
speciﬁcally the U-NSGA-III method as developed by Seada and Deb [88]. NSGA-III [21] is a
recent extension of the widely used NSGA-II [22] that incorporates the idea of predeﬁned reference
vectors, increasing the possible number of objectives from 2 to 3+. The weakness of NSGA-III is
that this change also makes it ineﬀective in the case of 2 or fewer objectives, meaning it is not a
replacement for NSGA-II but instead a complementary method. U-NSGA-III solves this problem
by extending NSGA-III in such a way that one or two-objective optimization becomes possible,
making it a well rounded method capable of any number of objectives.

The soft swimmer is driven by two PID controllers with three tuning parameters each, leading
to a total of six parameters. We optimize these parameters in terms of three objectives: combined
average error ( ¯e), cost of transport (Cot), and average strain energy ( ¯E), which are elaborated below.
The combined average error is the sum of the individually averaged errors as

¯e = | ¯e1| + | ¯e2|

Þ t
0 |ei(t)|dt, i = 1,2. This
in which we track the entire evolution history during time t via | ¯ei| = 1
provides a measure of tracking performance. The energy eﬃciency of the swimmer is evaluated
through the measure of cost of transport,

(4.11)

t

Cot =

E f
d ,

55

(4.12)

Table 4.1: Summary of the variables (PID tuning parameters) and objectives considered in the
optimization.

¯e
Cot
¯E

Variables Objectives
kp,position
ki,position
kd,position
kp,heading
ki,heading
kd,heading

where E f =Þ t
power P, and d =Þ t

0 Pdt is the total energy imparted to the ﬂuid calculated by integrating the imparted
0 |(cid:104)u(t(cid:48))(cid:105)|dt(cid:48) is the path length determined from the averaged motion of the

Þ udV (here we calculate

center of mass, in which we deﬁne (cid:104)u(cid:105) = 1
T
the integral for the solid velocity u in the Lagrangian domain of size V, and (cid:104)(cid:105) denotes a spatial
average). Finally, the average total elastic strain energy is deﬁned as

V

Þ t
t−T (cid:104)u(t(cid:48))(cid:105) dt(cid:48) and (cid:104)u(cid:105) = 1
Þ t

where we deﬁne E =Þ G (tr(B) − 2) dV for an incompressible neo-Hookean material [5]. This

(4.13)

E(t(cid:48))dt(cid:48)

0

¯E =

1
t

,

relates the the force and energy requirements of the actuators as well as the potential for fatigue due
to cyclic loading.

4.3 Results and Discussion

4.3.1 Undulatory Swimming

We have observed that after the initial transient, the soft swimmer quickly reaches a quasi-steady
forward swimming with periodic undulatory body motion. In Fig. 4.3(a), we show a typical ﬂow
map of the swimmer when choosing Re = 500, G = 1000, δ = 0.4, T = 2, and symmetric
actuation at α = 0.15 (namely, α0 = 0.15, β = 0 in Eq.(4.7)). A pair of vortices of opposite signs is
generated during each tail upward/downward beating, leading to a reverse von Kármán vortex street.
In Fig. 4.3(b), we show that the beam midline deformation extracted from the simulation during

56

Figure 4.3: Free-swimming of a 2D elastic beam at Re = 500 and T = 2.0. (a) Instantaneous
snapshot of vorticity ﬁeld. (b) Envelope trajectories of the beam midline during one actuation
period: (top) numerical results; (bottom) ﬁtted results by Eq. (4.14). (c) Mean values of U, E, and
Cot as function of contraction strength α.

each time period (top), which can be approximately ﬁtted by a weighted sine function (bottom):

κ(x, t) = (ax2 + bx + c) sin[2π (k x + ωt) + φ] .

(4.14)

In this case, we obtain the coeﬃcients a = −1.27, b = 0.92, c = 0.07, k = 0.4, ω = 1.0, and
φ = 0.5. Upon comparison, there is good match between the simulation result and the ﬁtted model
in the middle section while the front and tail sections show slight discrepancies. Overall, we see
that our simple three-parameter (α0, β,T) actuation scheme with muscle contraction is capable of
reproducing various types of traveling-wave-like swimming gaits that typically require a larger
parameter space.

For the forward swimming cases (β = 0), we have evaluated the mean values of swimming
speed U, elastic energy E, and cost of transport Cot as a function of the amplitude α0 of muscle
contraction when ﬁxing Re = 500 and T = 2. As shown in Fig. 4.3(c), the system exhibits largely
monotonic behaviors for the three examined variables when the actuation strength α is selected
between 0.05 and 0.25, with an exception for Cot, which shows a local peak around α = 0.08. For
very small values of α (α < 0.05), the soft swimmer does not show signiﬁcant directional motion;
when α > 0.25, we ﬁnd that the average swimming speed (and hence also Cot) decreases as α
increases, suggesting that large body deformation does not necessarily lead to eﬃcient swimming.
In this study, we focus on the monotonic regime to facilitate controller design.

57

Figure 4.4: Simulation results for the soft swimmer at α0 = 0.2, Re = 500, and G = 1000: (a)
Center of mass (CoM) paths at diﬀerent values of β; (b) Turning curvature R−1 as a function of β.

As discussed previously, specifying the empirical traveling-wave form of Eq. (4.14) only applies
to the scenario of forward swimming. To achieve more elaborate maneuvering such as turning, we
introduce a bias when applying contractions on the two sides of the swimming body. In Fig. 4.4,
we demonstrate the eﬀectiveness of the actuation scheme for turning with ﬁxed actuation strength
α0 (α0 = 0.2) but increasing β. In Fig. 4.4(a), the swimmer’s instantaneous CoM positions clearly
show that it approximately follows a circular trajectory, the radius of which decreases (tighter turn)
as the bias value β increases. Interestingly, the results in Fig. 4.4(b) suggest an almost linear relation
between the bias value of β and the trajectory curvature R−1 (reciprocal of the radius).

4.3.2 Optimization and Analysis

4.3.2.1 Trajectory Comparison

To make sure the controller is eﬀective, we optimize it by accounting for diﬀerent scenarios of the
moving target’s characteristic velocity and the Reynolds number Re as deﬁned at the end of Section
4.1. Note that we deﬁne Reynolds number in terms of the actuation frequency, not a physical
velocity, so by varying Re we are implying a change in the ﬂuid properties, namely the viscosity. In
these studies we employ two diﬀerent types of target trajectories, a straight path as characterized in
Eq. (4.15) and a sinusoidal path as expressed in Eq. (4.16). The swimmer starts from at rest facing

58

XY012345012345β=0.1β=0.075β=0.05β=0.025(a)RβR-100.040.080.1200.511.52(b)along the x axis. We describe the straight trajectory of the target with a constant velocity v as

xt =

yt = y(0)

xt = vt + x(0)
xt = vt + x(0)

xt =

yt = A sin ωt + y(0)

(4.15)

(4.16)

For the sinusoidal trajectory, the moving target has a constant velocity in the x direction, while its
y position follows a sinusoidal function of time:

Here we ﬁx the amplitude, A = 2, and angular frequency, ω = 0.05π. The inclusion of sinusoidal
motion means that the swimmer has to constantly adjust its velocity and direction, which demon-
strates its ability to follow complex paths. All cases are run with a maximum (dimensionless)
simulation time of t = 80, and as a precaution, the forward velocity controller is not activated until
t = 2 and the heading controller is not activated until t = 4. Before these times, α0 = αmax and
β = 0. The idea is to allow a smooth start and prevent potentially large values of β from negatively
aﬀecting the starting behavior of the swimmer, by ﬁrst allowing it to achieve a signiﬁcant forward
velocity. However, testing has shown both restrictions to be unnecessary if βmax is suﬃciently
smaller than αmax.

In Fig. 4.5, we show the performance of a soft swimmer when subjected to diﬀerent PID
controllers by comparing its CoM trajectories with the speciﬁed straight (panel(a)) and sinusoidal
(panel(b)) paths (black lines). Here we select two optimal solutions, corresponding to the case of
lowest error ¯e (min( ¯e)) and the case of a moderate error ( ¯e = 4.5) but lowest elastic energy (min( ¯E))
at Re = 500. For the minimum error examples in both panels, the swimmer quickly converges to
the target point’s position and is able to follow it closely. Note that the target is a moving point
and not simply a designated path, meaning that a close match to the target trajectory implies the
swimmer must be matching both the desired heading and velocity. Variation in the velocity can be
discerned by observing the variation in the wavelength of the undulatory oscillations in the CoM
trajectory. For a ﬁxed driving frequency, the swimmer’s velocity will correspond to the wavelength

59

Figure 4.5: CoM trajectories of a soft swimmer tracking a moving target that is moving along both
the (a) straight and (b) sinusoidal paths (black lines), when the swimmer is subjected to the PID
control. The (red) and (blue) lines correspond to the controllers that have the lowest and a
moderate error ¯e (but lowest elastic energy ¯E), respectively. The black arrows mark the target’s
moving directions.

of this oscillation. At the start the wavelength is long as the swimmer attempts to reach the target
as quickly as possible. The wavelength then shortens as the swimmer approaches the target and
matches its velocity. In the sinusoidal case, the wavelength is shorter near the peaks and longer in
between corresponding to the velocity of the sine function which varies by the cosine. In contrast,
for the moderate error case, we see that the swimmer appears to lag behind, and eventually fails to
catch up to the target, although the small amplitudes (relative to the best ¯e cases) indicates lower
contractile strengths which implies less deformation of the body and thus lower elastic strain energy.

4.3.2.2 Pareto Fronts of Optimized Solutions

Here we show global views of the obtained optimal solutions. Figure 4.6 compares the resulting
Pareto fronts (a set of nondominated optimal solutions if no objective can be improved without
sacriﬁcing at least one other objective) in the space of objectives when the target velocity v takes
diﬀerent values while the Reynolds number Re is kept at 500, for the cases of straight trajectory ((a)-
(d)) and sinusoidal trajectory ((e)-(h)), respectively. The obtained optimal solutions for these two

60

(a)xt = vt + x(0)yt = y(0)min(e)e = 4.5(b)xt = vt + x(0)yt = A sin(ω t) + y(0)min(e)e = 4.5diﬀerent paths are approximately allocated in the same regimes, exhibiting very similar patterns.
Figure 4.7 shows the results when the Reynolds number is varied while the velocity is ﬁxed at 0.2.
All datasets shown correspond to a total of 50 generations of evolution with a child population
of 64. The swimmer’s body is neutrally buoyant with δ = 0.5 and non-dimensional parameters
G = 1000 and λ0 = 104. It is worthwhile to mention that some portions of the objective space,
and the individual solutions (PID parameters) that fall within, may be undesirable. For example,
there are scenarios when the swimmer is barely moving and possibly not swimming at all. This
“trivial” state can be interpreted using the three objectives, where the ideal optimal value of strain
energy, ¯E = 0, coincides with the trivial case of zero actuation, αR = αL = 0. From the optimizer’s
prospective, this is a valid solution because our objective functions do not in any way represent
this as a negative behavior. The objective functions could possibly be modiﬁed to avoid this trivial
solution by adding a penalty function or cutoﬀ, but such changes are undesirable since they could
potentially add constraints or biases on the optimizer in unpredictable ways. Thus we stick to the
original form of objective functions and remove trivial solutions in the end as part of the ﬁnal
decision-making process.

A 2D view of the Pareto front on the projected ( ¯E, ¯e) plane can clearly reveal the eﬀect of
varying v on the resulting optimal solutions, as seen in Fig. 4.6(c) and (g). The distributions
are clearly separate with the proﬁle shifting up and right as v increases, both as expected, since
increasing v indicates the increase in tracking diﬃculty and the requirement of higher actuation
output. In comparison, the eﬀect on Cot as seen in Figs. 4.6 (b), (d), (f), and (h) appears to be
minimal with the proﬁles being nearly independent of v.

Varying Re proves to have a signiﬁcantly diﬀerent eﬀect with dependency apparent primarily
in Cot and ¯E as can be seen in Fig. 4.7. From Figs. 4.7(b), (d), (f), and (h), we see that as Re
increases, Cot drops. From Figs. 4.7(c) and (g), the elastic strain energy also drops with increasing
Re. As Re increases the viscous forces acting on the swimmer decrease, leading to both a lower
Cot and actuation strength required to drive the swimmer.
In all cases, for the regions where
the error ¯e is low, there is an inverse relationship between Cot and the error. The error proves

61

Figure 4.6: Pareto fronts in the space of tracking error ( ¯e), cost of transport (Cot), and elastic
strain energy ( ¯E) for various target speeds v with Re = 500. Blue-square: v = 0.1. Red-triangle: v
= 0.2. Green-circle: v = 0.3. Straight path: (a-d). Sinusoidal path: (e-h). The leftmost panels (a)
and (e) are 3D views of the objective spaces. The panels to their right are 2D projections of the
respective distribution that visualize the trends between each pair of objectives. (b) and (f): ¯e vs
Cot. (c) and (g): ¯e vs ¯E. (d) and (h): Cot vs ¯E.

to be mainly dominated by the positional error, e1, which in turn is logically inversely related to
the swimmer’s velocity. From this we can conclude that Cot increases with increasing swimmer
velocity as expected due to ﬂuid drag.

4.3.2.3 Objective Correlations and Parameter Trends

In Fig. 4.8 we further compare the Pearson correlation coeﬃcients of the objectives for each data set.
The most apparent result is that in all cases the coeﬃcient relating ¯e and ¯E is approximately equal to
−1. This means that they are correlated strongly with an inversely proportional relationship, which
is conﬁrmed by observing panels (c) and (g) in Figs. 4.6 and 4.7 where the population in each
dataset has a highly linear distribution. Physically this indicates that lower error requires stronger
actuation. The strength of this correlation also implies that these objectives are not independent in
our system. For example, it is not possible to ﬁnd an optimal solution that decreases both ¯e and ¯E
simultaneously.

62

Figure 4.7: Pareto fronts in the space of tracking error ( ¯e), cost of transport (Cot), and elastic
strain energy ( ¯E) for various Re with v = 0.2. Blue-square: Re = 500. Red-triangle: Re = 1000.
Green-circle: Re = 1500. Yellow-diamond: Re = 2000. Straight path: (a-d). Sinusoidal path:
(e-h). The leftmost panels (a) and (e) are 3D views of the objective spaces. The panels to their
right are 2D projections of the respective distribution that visualize the trends between each pair
of objectives. (b) and (f): ¯e vs Cot. (c) and (g): ¯e vs ¯E. (d) and (h): Cot vs ¯E.

The Pearson coeﬃcients relating ¯e to Cot as well as relating Cot to ¯E do not show similarly clear
patterns. Qualitatively, in the regions of lower ¯e there appear to be noticeable and consistent trends;
however, when viewed as a whole, these fall apart. Thus without deﬁning a criterion to trim the
data set, a solid conclusion cannot be made. Interestingly it can be noted that coeﬃcients relating
¯e to Cot and those relating Cot to ¯E appear to be equal in magnitude and inverted in sign. This can
be explained by the objectives in each dataset falling roughly in a plane such that the distribution
seen in panels (b) and (f) (corresponding to ¯e to Cot) is the mirror of that seen in panels (d) and(h)
(corresponding to Cot to ¯E).

Furthermore, we seek more quantitative understandings of how the controller parameters vary
when moving through the objective space. Considering the fact that the mapping between the
objective and parameter space is neither necessarily smooth nor continuous, we emphasize that
regions or even individuals in the Pareto optimal objective space must be “classiﬁed” into diﬀerent
groups that needs to be treated separately. For our study there is an obvious choice for classiﬁcation,

63

Figure 4.8: Comparison of the Pearson correlation coeﬃcients (ρ) of the objectives for varying v
(left) and varying Re (right). Blue-square = ¯e to Cot. Red-triangle = ¯e to ¯E. Yellow-circle = Cot to
¯E.

swimmers versus non-swimmers. As noted previously, there are eﬀectively trivial solutions when
¯E is near zero and the error is correspondingly high. This occurs when the the contraction strength
is suﬃciently small such that the soft swimmer barely moves. By excluding these high error
individuals (above roughly half the maximum error), we are left with parameter regions that may
exhibit qualitatively clear trends in the objective space.

Fig. 4.9 shows the trends of kp,position versus error for all cases used in the optimization study at
diﬀerent values of v and Re. When near the minimum error, kp,position appears to reach maximum
values, which corresponds to the optimal solutions near the upper bound of contraction strength α0.
For fairly large error values, kp,position becomes negligible, corresponding to those trivial solutions
with low ¯E. From the PID controller formulation in Eq. (4.8), α0 ∝ kp,positione1, suggesting that
when α0 is properly bounded, larger controller gain kp,position results in reduced error.

Besides the kp,position − error relations, however, we did not ﬁnd other clear trends when we
sweep through the parameter space and connect the tuning parameters to the objectives. The data
appear to be highly noisy, which is primarily due to the stochastic nature of evolutionary methods.
We have observed that often times, they are able to converge rapidly to a certain optimum region,
and then slowly approach the optimal solution(s) within via a stochastic search. Nevertheless,
when comparing the obtained optimal tuning parameters for diﬀerent target path (i.e., straight vs.

64

sinusoidal) as well as across Reynolds numbers (e.g., see Fig. 4.9), we have found that their values
all fall into approximately the same ranges, indicating that these solutions indeed characterize the
global optima. In other words, the proposed optimization approach seems to have produced PID
controllers with near-optimal behavior across varying target behaviors and ﬂuid environments.

Figure 4.9: The trends of kp,position vs. ¯e. (a-b) Cases of varying v. (c-d) Cases of varying Re.
(Left) panels are from the straight path and (right) from the sinusoidal path. The black dashed
lines delineate the borders between the meaningful and trivial solutions.

4.4 Concluding Remarks

In this work we presented a new computational framework for optimizing feedback control of
soft robotic ﬁsh enabled by muscle-like actuation. A novel ﬁctitious domain/active-strain method
[58] was employed to handle the resultant ﬂuid-structure interaction when local contractile active
strains – emulating muscle actuation – are imposed on the soft swimmer. In particular, using a 2D

65

elastic beam with ﬁnite width as a model for the soft swimmer, we demonstrated versatile swimming
motions (forward swimming and turning at diﬀerent velocities and radii, for example) by varying the
strength and bias of the active strains imposed on both sides of the swimmer. Two PID controllers
are adopted for adjusting the actuation strength (α0) and bias (β), respectively, for the swimmer to
track a moving target. With this setup, multi-objective optimization of the controller parameters
was conducted using the U-NSGA-III algorithm, where three objectives were considered, including
the tracking error, cost of transport, and elastic strain energy. The PID controller parameters were
optimized using a series of moving target trajectories, involving diﬀerent speeds and Reynolds
numbers for two types of paths (straight and sinusoidal). By analyzing the resulting Pareto fronts,
the study revealed the trends of how the obtained values of objectives are inﬂuenced by the target
trajectory conditions.
It also showed the correlation and trade-oﬀs among the objectives under
feedback control, as well as how the optimized controller parameters are related to the objectives
and trajectory conditions. Overall, this work has shown the feasibility of using high-ﬁdelity CFD to
design and optimize feedback control for soft robotic ﬁsh subject to multiple, potentially conﬂicting
objectives.

There are a few directions for future work. First, we will consider the extension of the presented
framework to the case where the soft swimmer has a 3D geometry and a more general form
of distributed actuation. Second, while PID controllers were used as an example of model-free
controller in this work, we plan to exploit the computation data and system ID techniques to
obtain the dynamic model of a soft swimmer, based on which model-based controllers (such as
linear quadratic regulator (LQR) and model-predictive control [8]) will be developed and evaluated.
Ultimately, it is our plan to experimentally validate the modeling, control, and optimization methods
with a soft robotic ﬁsh prototype actuated by smart materials (e.g., shape memory alloys or dielectric
elastomer actuators).

66

CHAPTER 5

FUTURE WORK

5.1 Body Shape Optimization for Muscle Driven Undulatory Swimming

In Chapter 4 we demonstrate the eﬀectiveness of coupling our simulation techniques to multi-
objective evolutionary optimizers for design purposes. The question now is how far can we push
the envelope and how many factors can be included in an optimization. Most existing studies
are limited to single objective methods and/or simpliﬁed swimmer models. To build oﬀ of these
studies we attempt a multi-objective optimization of the body shape of a soft, artiﬁcial muscle
driven swimmer using our high ﬁdelity simulation method. Body shape is of course important
from a hydrodynamic perspective but also because it directly inﬂuences the gait of soft swimmers.
How well or even if a soft robot can swim depends on the mechanical response of the body which
is determined by it’s shape and actuation scheme. So to start we elect to optimize the body of a 2D
robotic ﬁsh with a simple actuation scheme similar to that used previous in the plate ﬁsh.

Three objectives were chosen to optimize: swimming speed (max), swimming eﬃciency
(max), and average elastic energy (min). The swimming speed is a self explanatory measure of the
eﬀectiveness of the swimmer. The swimming eﬃciency is deﬁned by the ratio of the kinetic energy
to the energy transferred to the ﬂuid over a given period. This provides a measure of how well the
swimmer translates expended energy into useful locomotion. Lastly, the average elastic energy of
the body relates to the strength of the required actuators (muscles), overall energy expenditure, and
material fatigue.

To vary the shape, the swimmers body is deﬁned via multiple points connected by an interpo-
lating Catmull-Rom spline with bilateral symmetry. Fig. 5.1 shows a randomly sampled swimmer
from an initial population. Actuation is limited to an active region of ﬁxed length and strength,
the position of which is also evolved. As in the plate ﬁsh, the actuation is approximated to be at
the surface via exponential decay into the body. The local depth and direction of the actuation at a

67

point inside the body is determined from the linear surface mesh.

Figure 5.1: An example of a randomly generated swimmer body. The blue line is the body surface
and the dots are interpolation points for the Catmull-Rom spline. Muscle actuation is limited to
within the active region.

Figure 5.2: Sample objective space population distribution for generation 100 of a variable body
shape optimization.

Fig. 5.2 shows a sample population from a variable body shape optimization as was just deﬁned.
By nature it is diﬃcult to visualize the objective space of problems with 3 or more objectives in a
2D format. Here the surface correlating to the population is concave into the image. When fully
converged, each individual represents a potentially optimal conﬁguration. Convergence is diﬃcult to
deﬁne and usually involves simply running the optimization until the population suﬃciently settles.
This does not guarantee that the true Pareto front has been found but for complex systems even a
partial optimization can be useful when the alternative is often none at all. Fig. 5.3 is an attempt

68

at visualizing the convergence by focusing only on the extremes, the best value present for each
objective in a given population. This ignores the intermediate individuals meaning convergence
here does not indicate overall convergence. However, since overall convergence does require the
extremes to be converged it provides a necessary but not suﬃcient convergence condition.

Figure 5.3: Sample convergence history of the extreme (most optimal) values for each objective in
a variable body shape optimization.

The stair step pattern in the convergence history is characteristic of evolutionary optimization.
Each jump corresponds to a random variation that lead to positive change. Plateaus occur where
the variations produce nonviable or negative change. As the generation count increases the plateaus
tend to grow longer since the more converged the solution is, the more likely a random variation will
be harmful. Thus demonstrating how evolutionary methods can quickly ﬁnd the general location of
an optimum but converge to an exact value slowly. The signiﬁcant changes in the objective values
from beginning to end seemed to indicate that the optimization was able to progress meaningfully.
Fig. 5.4 shows the individuals corresponding to the extreme objective values sampled at
progressive generations. Analyzing the ﬁnal generation (gen. 100) it is apparent that the fastest (a)
and lowest average elastic energy (c) individuals have evolved toward a slender body reminiscent
of an eel or tadpole with a streamlined and highly ﬂexible body. The fastest individual resembles
the plate ﬁsh except with a pronounced head that we hypothesize reduces head motion which is a

69

Figure 5.4: Sample conﬁgurations over progressive generations: (a) fastest swimming speed, (b)
highest swimming eﬃciency, and (c) lowest average elastic energy.

major source of drag. In the case of elastic energy, a slender body reduces the bending stiﬀness of
the body and thus the elastic energy stored in the body for a given deformation. The high eﬃciency
individual is the exception with a thick body that seems to limit the amplitude of the propulsive
paddling motion. The exact advantage is not clear but the less chaotic motion of the body may
reduce the presence of drag generating ﬂuid structures such as leading edge vortices.

The current trajectory for work in this area is focused on enabling and optimizing more compli-
cated actuation schemes. The simplicity of the actuation scheme allowed us to qualitatively verify
the eﬀectiveness of the optimization by repeatedly producing the tadpole shape that we see in nature
and with similar performance to our undulating plate. Ideally, we would like to simultaneously
optimize the body shape along with a complex, full body actuation scheme akin to existing studies
[99]. The major issue with evolving the actuation pattern is that for the case of elastic energy, and
to some degree eﬃciency, the optimizer will tend toward the trivial solutions where the actuation
strength is near zero everywhere, such that the swimmer is eﬀectively turned oﬀ. These swimmers
clearly do not achieve any useful swimming velocity which defeats the overall purpose of the study.
Thus a major goal is to ﬁnd a way to promote a minimum level of swimming speed that neither
overly suppresses exploration nor requires modiﬁcation of the evolutionary algorithms.

70

CHAPTER 6

CONCLUSIONS

We set out on this project with the intent to develop a computational framework that would
help drive the growth of soft robotics research and technology in ways existing simulations and
experiments could not yet achieve. Previously existing and commonly used simulation techniques
for FSI generally contain some form of limitation that prevents them from capturing the full, rich
physics of swimming soft robots. Often times these limitations stem from assumptions and reduced
models that greatly simplify the problem at the cost of an acceptable loss of physical accuracy and
generality for their originally intended purposes. For muscle driven, fully soft, swimming robots
the list of valid assumptions becomes much smaller and so to achieve a viable level of physical
accuracy we are forced to look to more complete, more general, and more expensive methods.

In this work we presented a holistic computational framework for the design and optimization
of soft robotic ﬁsh powered by artiﬁcial muscle along with complementary investigations of phe-
nomena related to ﬂuid/elastic-structure interaction. The active strain method when coupled with
the distributed Lagrange multiplier/ﬁctitious domain method allows for high ﬁdelity, fully coupled
simulations of soft robots that employ muscle driven undulating, jetting, and paddling propul-
sion. By coupling these simulations with a multi-objective evolutionary optimization method we
demonstrated that it is possible to design and optimize maneuverable, feedback controlled muscle
driven soft swimming robots. We have also investigated the eﬀect of solid viscoelasticity on the
behavior of instabilities at a ﬂuid-solid interface and noted key regions where the eﬀects are most
pronounced.

We believe we have achieved our original goal and hope that this work can useful in promoting
the ﬁeld of fully soft swimming robots by providing tools and insights that can be used to investigate
designs and materials that are still maturing, too expensive, or too complex for large scale studies.
We believe it also has the potential to act as a reference for developing new reduced order and
machine learning based models where experimental data is currently of limited supply.

71

BIBLIOGRAPHY

72

BIBLIOGRAPHY

[1] S. Alben, L. Miller, and J. Peng. Eﬃcient kinematics for jet-propelled swimming. J. Fluid

Mech., 733:100, 2013.

[2] D. Ambrosi and S. Pezzuto. Active stress vs. active strain in mechanobiology: constitutive

issues. Journal of Elasticity, 107(2):199–212, 2012.

[3] C. Antoci, M. Gallati, and S. Sibilla. Numerical simulation of ﬂuid–structure interaction

by sph. Comput. Struct., 85:879–890, 2007.

[4] N. W. Bartlett, M. T. Tolley, J. T. Overvelde, J. C. Weaver, B. Mosadegh, K. Bertoldi,
G. M. Whitesides, and R. J. Wood. A 3D-printed, functionally graded soft robot powered
by combustion. Science, 349:161–165, 2015.

[5] T. Belytschko, W. K. Liu, B. Moran, and K. Elkhodary. Nonlinear Finite Elements for

Continua and Structures. John Wiley & Sons, 2013.

[6] R. Bhardwaj and R. Mittal. Benchmarking a coupled immersed-boundary-ﬁnite-element

solver for large-scale ﬂow-induced deformation. AIAA Journal, 50:1638–1642, 2011.

[7]

J. C. Case, E. L. White, and R. K. Kramer. Sensor enabled closed-loop bending control of
soft beams. Smart Materials and Structures, 25(4):045018, 2016.

[8] M. Castano and X. Tan. Model predictive control-based path following for tail-actuated
robotic ﬁsh. Journal of Dynamic Systems, Measurement and Control, 141(7):071012, 2019.

[9] S. Chandrasekhar. Hydrodynamic and Hydromagnetic Stability. Oxford University Press,

1961.

[10] D. Chen, S. Cai, Z. Suo, and R. Hayward. Surface energy as a barrier to creasing of
elastomer ﬁlms: an elastic analogy to classical nucleation. Phys. Rev. Lett., 109(3):038001,
2012.

[11] K. Chen. Elastic instability of the interface in couette ﬂow of viscoelastic liquids. J. Non-

Newtonian Fluid Mech, 40:261–267, 1991.

[12] Q. Chen, Q. Zhong, M. Qi, and X. Wang. Comparison of vortex identiﬁcation criteria for

planar velocity ﬁelds in wall turbulence. Physics of Fluids, 27(8):085101, 2015.

[13] S. Chen and G. Doolen. Lattice boltzmann method for ﬂuid ﬂows. Annu. Rev. Fluid Mech.,

30:329, 1998.

73

[14] S. Chester and L. Anand. A coupled theory of ﬂuid permeation and large deformations
for elastomeric materials. Journal of the Mechanics and Physics of Solids, 58:1879–1906,
2010.

[15] C. Christianson, N. N. Goldberg, D. D. Deheyn, S. Cai, and M. T. Tolley. Translucent soft
robots driven by frameless ﬂuid electrode dielectric elastomer actuators. Science Robotics,
3(17):eaat1893, 2018.

[16] S. Colin and J. Costello. Morphology, swimming performance and propulsive mode of six

co-occurring hydromedusae. J. Exp. Biol., 205:427, 2002.

[17] S. Colin, J. Costello, and E. Klos. IN SITU swimming and feeding behavior of eight

co-occurring hydromedusae. Mar. Ecol. Prog. Ser., 253:305, 2003.

[18]

[19]

[20]

J. Costello, S. Colin, and J. Dabiri. Medusan morphospace: phylogenetic constraints, biome-
chanical solutions, and ecological consequences. Invertebr. Biol., 127:265, 2008.

J. Dabiri, S. Colin, J. Costello, and M. Gharib. Flow patterns generated by oblate medusan
jellyﬁsh: ﬁeld measurements and laboratory analyses. J. Exp. Biol., 208:1257, 2005.

J. Dabiri, S. Colin, K. Katija, and J. Costello. A wake-based correlate of swimming perfor-
mance and foraging behavior in seven co-occurring jellyﬁsh species. J. Exp. Biol., 213:1217,
2010.

[21] K. Deb and H. Jain. An evolutionary many-objective optimization algorithm using reference-
point-based nondominated sorting approach, Part I: solving problems with box constraints.
IEEE Transactions on Evolutionary Computation, 18:577–601, 2014.

[22] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan. A fast and elitist multiobjective genetic

algorithm: nsga-ii. IEEE Transactions on Evolutionary Computation, 6:182–197, 2002.

[23] D. Drotman, S. Jadhav, M. Karimi, P. Dezonia, and M. T. Tolley. 3D printed soft actuators
for a legged robot capable of navigating unstructured terrain. In 2017 IEEE International
Conference on Robotics and Automation (ICRA), pages 5532–5538. IEEE, 2017.

[24] F. Fang, K. Ho, L. Ristroph, and M. Shelley. A computational model of the ﬂight dynamics

and aerodynamics of a jellyﬁsh-like ﬂying machine. J. Fluid Mech., 819:621, 2017.

[25] L. Fauci and A. McDonald. Sperm motility in the presence of boundaries. Bulletin of

Mathematical Biology, 57:670–699, 1995.

[26] C. Gambini, B. Abou, A. Ponton, and A. J. M. Cornelissen. Micro- and macrorheology of

jellyﬁsh extracellular matrix. Biophys. J., 102:1–9, 2012.

74

[27] T. Gao and H. H. Hu. Deformation of elastic particles in viscous shear ﬂow. J. Comput.

Phys., 228:2132–2151, 2009.

[28] T. Gao, H. H. Hu, and P. Ponte Castañeda. Rheology of a suspension of elastic particles in

a viscous shear ﬂow. J. Fluid Mech., 687:209–237, 2011.

[29] T. Gao, H. H. Hu, and P. Ponte Castañeda. Shape dynamics and rheology of soft elastic

particles in a shear ﬂow. Phys. Rev. Lett., 108:058302, 2012.

[30] T. Gao, Y. Tseng, and X. Lu. An improved hybrid cartesian/immersed boundary method

for ﬂuid-solid ﬂows. Int. J. Numer. Meth. Fluids, 55:1189–1211, 2007.

[31] M. Gazzola, A. A. Tchieu, D. Alexeev, A. de Brauer, and P. Koumoutsakos. Learning to
school in the presence of hydrodynamic interactions. Journal of Fluid Mechanics, 789:726–
749, 2016.

[32] M. Gazzola, W. M. Van Rees, and P. Koumoutsakos. C-start: optimal start of larval ﬁsh.

Journal of Fluid Mechanics, 698:5–18, 2012.

[33] V. Gkanis and S. Kumar. Instability of creeping couette ﬂow past a neo-hookean solid.

Phys. Fluids, 15:2864–2471, 2003.

[34] V. Gkanis and S. Kumar. Stability of pressure-driven creeping ﬂows in channels lined with

a nonlinear elastic solid. J. Fluid Mech., 524:357–375, 2005.

[35] R. Glowinski, T.-W. Pan, T. I. Hesla, and D. D. Joseph. A distributed lagrange multi-
plier/ﬁctitious domain method for particulate ﬂows. International Journal of Multiphase
Flow, 25:755–794, 1999.

[36] C. Hamlet, L. J. Fauci, and E. D. Tytell. The eﬀect of intrinsic muscular nonlinearities
on the energetics of locomotion in a computational model of an anguilliform swimmer.
Journal of Theoretical Biology, 385:119–129, 2015.

[37] X. He and L. Luo. Lattice boltzmann model for the incompressible navier–stokes equation.

J. Stat. Phys., 88:927, 1997.

[38] G. Hersch and L. Miller. Reynolds number limits for jet propulsion: a numerical study of

simpliﬁed jellyﬁsh. J. Theor. Biol., 285:84, 2011.

[39] A. V. Hill. The heat of shortening and the dynamic constants of muscle. Proceedings of the

Royal Society of London. Series B-Biological Sciences, 126(843):136–195, 1938.

[40] H. H. Hu, M. Y. Zhu, and N. Patankar. Direct numerical simulations of ﬂuid solid systems

using the arbitrary lagrangian eulerian technique. J. Comput. Phys., 169:427–462, 2001.

75

[41] H. Huang and X.-Y. Lu. An ellipsoidal particle in tube poiseuille ﬂow. Journal of Fluid

Mechanics, 822:664–688, 2017.

[42] H. Huang, X. Yang, and X.-y. Lu. Sedimentation of an ellipsoidal particle in narrow tubes.

Physics of Fluids, 26(5):053302, 2014.

[43] W.-X. Huang and H. J. Sung. Three-dimensional simulation of a ﬂapping ﬂag in a uniform

ﬂow. Journal of Fluid Mechanics, 653:301–336, 2010.

[44] H. Jain and K. Deb. An evolutionary many-objective optimization algorithm using reference-
point-based nondominated sorting approach, part ii: handling constraints and extending to
an adaptive approach. IEEE Transactions on Evolutionary Computation, 18:602–622, 2014.

[45] D. D. Joseph. Fluid Dynamics of Viscoelastic Liquids. Applied mathematical sciences, Vol.

84, Springer Verlag, New York, 1990.

[46] V. Kalro and T. Tezduyar. A parallel 3d computational method for ﬂuid-structure interactions

in parachute systems. Comput. Meth. Appl. Mech. Eng., 190:321, 2000.

[47] Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin. Inferring the structure
and dynamics of interactions in schooling ﬁsh. Proceedings of the National Academy of
Sciences, 108(46):18720–18725, 2011.

[48] S. Kern and P. Koumoutsakos. Simulations of optimized anguilliform swimming. Journal

of Experimental Biology, 209(24):4841–4857, 2006.

[49]

J. Kim, D. Kim, and H. Choi. An immersed-boundary ﬁnite-volume method for simulations
of ﬂow in complex geometries. Journal of Computational Physics, 171:132–150, 2001.

[50] S.-W. Kim, J.-G. Lee, S. An, M. Cho, and K.-J. Cho. A large-stroke shape memory
alloy spring actuator using double-coil conﬁguration. Smart Materials and Structures,
24(9):095014, 2015.

[51] V. Kumaran, G. H. Fredrickson, and P. Pincus. Flow-induced instability at the interface

between a ﬂuid and a gel at low reynolds number. J. Phys. II France, 4:893–911, 1994.

[52] V. Kumaran and R. Muralikrishnan. Spontaneous growth of ﬂuctuations in the viscous ﬂow

of a ﬂuid past a soft interface. Phys. Rev. Lett., 84:3310–3313, 2000.

[53] C. Lai, Z. Zheng, E. Dressaire, G. Ramon, H. Huppert, and H. Stone. Elastic relaxation of

ﬂuid-driven cracks and the resulting backﬂow. Phys. Rev. Lett., 117:268001, 2016.

[54] E. H. Lee and D. T. Liu. Finite-strain elastic-plastic theory particularly for plane wave

analysis. Journal of Applied Physics, 38:19–27, 1967.

76

[55] E. H. Lee. Elastic-plastic deformation at ﬁnite strains. Journal of Applied Mechanics,

36(1):1–6, 1969.

[56] L. Li, S. Sherwin, and P. Bearman. A moving frame of reference algorithm for ﬂuid/structure
interaction of rotating and translating bodies. International Journal for Numerical Methods
in Fluids, 38:187–206, 2002.

[57] M. J. Lighthill. Note on the swimming of slender ﬁsh. Journal of Fluid Mechanics, 9:305–

317, 1960.

[58] Z. Lin, A. Hess, Z. Yu, S. Cai, and T. Gao. A ﬂuid-structure interaction study of soft
robotic swimmer using a ﬁctitious domain/active-strain method. Journal of Computational
Physics, 376:1138–1155, 2019.

[59] C. Macosko. Rheology: Principles, Measurements and Applicatfons. VCH Publishers, New

York, 1994.

[60] C. Majidi. Soft robotics: a perspective—current trends and prospects for the future. Soft

Robotics, 1:5–11, 2014.

[61]

J. Martins, E. Pires, R. Salvado, and P. Dinis. A numerical model of passive and active
behavior of skeletal muscles. Computer methods in applied mechanics and engineering,
151(3-4):419–433, 1998.

[62] A. Masud and T. Hughes. A space-time galerkin/least-squares ﬁnite element formulation
of the navier-stokes equations for moving domain problems. Comput. Methods Appl. Mech.
Eng., 146:91–126, 1997.

[63] N. G. McCrum, C. P. Buckley, and C. B. Bucknall. Principles of Polymer Engineering.

Oxford University Press, 1997.

[64] L. Mihai and A. Goriely. Positive or negative poynting eﬀect? the role of adscititious
inequalities in hyperelastic materials. In Proc. R. Soc. Lond. A, volume 467 of number 2136,
pages 3633–3646, 2011.

[65] A. Miriyev, K. Stack, and H. Lipson. Soft material for soft actuators. Nature Communica-

tions, 8(1):596, 2017.

[66] R. Mittal, H. Dong, M. Bozkurttas, F. Najjar, A. Vargas, and A. von Loebbecke. A ver-
satile sharp interface immersed boundary method for incompressible ﬂows with complex
boundaries. Journal of Computational Physics, 227(10):4825–4852, 2008.

[67] R. Mittal and G. Iaccarino. Immersed boundary methods. Annu. Rev. Fluid Mech., 37:239–

261, 2005.

77

[68] S. Mora, T. Phou, J.-M. Fromental, L. M. Pismen, and Y. Pomeau. Capillarity driven

instability of a soft solid. Phys. Rev. Lett, 105:214301, 2010.

[69] R. Muralikrishnan and V. Kumaran. Experimental study of the instability of the viscous

ﬂow past a ﬂexible surface. Phys. Fluids, 14:775–780, 2002.

[70]

J. Nawroth, H. Lee, A. Feinberg, C. Ripplinger, M. McCain, A. Grosberg, J. Dabiri, and
K. Parker. A tissue-engineered jellyﬁsh with biomimetic propulsion. Nature Biotechnology,
30:792–797, 2012.

[71] R. W. Ogden. Nonlinear Elastic Deformations. Dover, 1984.

[72] V. Palmre, J. J. Hubbard, M. Fleming, D. Pugal, S. Kim, K. J. Kim, and K. K. Leang. An
ipmc-enabled bio-inspired bending/twisting ﬁn for underwater applications. Smart Materi-
als and Structures, 22(1):014003, 2012.

[73] S. Park, C. Chang, W. Huang, and H. Sung. Simulation of swimming oblate jellyﬁsh with

a paddling-based locomotion. J. Fluid Mech., 748:731, 2014.

[74] S. Park, M. Gazzola, K. Park, S. Park, V. Di Santo, E. Blevins, J. Lind, P. Campbell, S.
Dauth, A. Capulli, F. Pasqualini, S. Ahn, A. Cho, H. Yuan, B. Maoz, R. Vijaykumar, J.
Choi, K. Deisseroth, G. Lauder, L. Mahadevan, and K. Parker. Phototactic guidance of a
tissue-engineered soft-robotic ray. Science, 353:158–162, 2016.

[75] S. Park, B. Kim, J. Lee, W. Huang, and H. Sung. Dynamics of prolate jellyﬁsh with a

jet-based locomotion. J. Fluids Struct., 57:331, 2015.

[76]

J. K. Parrish, S. V. Viscido, and D. Grunbaum. Self-organized ﬁsh schools: an examination
of emergent properties. The Biological Bulletin, 202(3):296–305, 2002.

[77] C. S. Peskin. The immersed boundary method. Acta Numer., 10:479–517, 2002.

[78]

J. Poynting. On pressure perpendicular to the shear planes in ﬁnite pure shears, and on the
lengthening of loaded wires when twisted. Proc. R. Soc. Lond. A, 82(557):546–559, 1909.

[79] E. M. Purcell. Life at low reynolds number. Am. J. Phys., 45:3–11, 1977.

[80] Y. Renardy. Stability of the interface in two-layer couette ﬂow of upper convected maxwell

liquids. J. Non-Newtonian Fluid Mech., 28:99–115, 1988.

[81] L. Ristroph and S. Childress. Stable hovering of a jellyﬁsh-like ﬂying machine. J. R. Soc.

Interface, 11:20130992, 2014.

78

[82] D. E. Rivera, M. Morari, and S. Skogestad. Internal model control: PID controller design.
Industrial & Engineering Chemistry Process Design and Development, 25(1):252–265,
1986.

[83] M. Al-Rubaiai, T. Pinto, D. Torres, N. Sepulveda, and X. Tan. Characterization of a 3d-
printed conductive pla material with electrically controlled stiﬀness. In ASME 2017 Con-
ference on Smart Materials, Adaptive Structures and Intelligent Systems, V001T01A003–
V001T01A003. American Society of Mechanical Engineers, 2017.

[84] D. Rus and M. T. Tolley. Design, fabrication and control of soft robots. Nature, 521:467–

475, 2015.

[85] M. Sahin and K. Mohseni. An arbitrary lagrangian-eulerian formulation for the numerical
simulation of ﬂow patterns generated by the hydromedusa aequorea victoria. J. Comput.
Phys., 228:4588, 2009.

[86] M. Sahin, K. Mohseni, and S. Colin. The numerical comparison of ﬂow patterns and
propulsive performances for the hydromedusae sarsia tubulosa and aequorea victoria. J.
Exp. Biol., 212:2656, 2009.

[87] B. Saintyves, T. Jules, T. Salez, and L. Mahadevan. Self-sustained lift and low friction via

soft lubrication. Proc. Natl. Acad. Sci., 113(21):5847–5849, 2016.

[88] H. Seada and K. Deb. A uniﬁed evolutionary optimization procedure for single, multiple,

and many objectives. IEEE Transactions on Evolutionary Computation, 20:358, 2016.

[89] W. Shan, S. Diller, A. Tutcuoglu, and C. Majidi. Rigidity-tuning conductive elastomer.

Smart Materials and Structures, 24(6):065001, 2015.

[90]

J. Skotheim and L. Mahadevan. Soft lubrication. Phys. Rev. Lett., 92(24):245509, 2004.

[91] Z. Suo. Theory of dielectric elastomers. Acta Mechanica Solida Sinica, 23:549–578, 2010.

[92] P. L. Tallec and J. Mouro. Fluid structure interaction with large structural displacements.

Comput. Methods Appl. Mech. Eng., 190:3039–3067, 2001.

[93]

J. Teran, L. Fauci, and M. Shelley. Viscoelastic ﬂuid response can increase the speed and
eﬃciency of a free swimmer. Physical Review Letters, 104:038101, 2010.

[94] F.-B. Tian, H. Dai, H. Luo, J. F. Doyle, and B. Rousseau. Fluid–structure interaction
involving large deformations: 3d simulations and applications to biological systems. Journal
of computational physics, 258:451–469, 2014.

79

[95] H. Tian, Z. Wang, Y. Chen, J. Shao, T. Gao, and S. Cai. Polydopamine-coated main-chain
liquid crystal elastomer as optically driven artiﬁcial muscle. ACS Applied Materials &
Interfaces, 10(9):8307–8316, 2018.

[96] S. Trabia, Z. Olsen, and K. J. Kim. Searching for a new ionomer for 3d printable ionic
polymer–metal composites: aquivion as a candidate. Smart Materials and Structures,
26(11):115029, 2017.

[97] S. Turek and J. Hron. Proposal for numerical benchmarking of ﬂuid-structure interaction
between an elastic object and laminar incompressible ﬂow. In Fluid-structure interaction,
pages 371–385. Springer, 2006.

[98] E. D. Tytell, M. C. Leftwich, C.-Y. Hsu, B. E. Griﬃth, A. H. Cohen, A. J. Smits, C. Hamlet,
and L. J. Fauci. Role of body stiﬀness in undulatory swimming: insights from robotic and
computational models. Physical Review Fluids, 1(7):073202, 2016.

[99] W. M. van Rees, M. Gazzola, and P. Koumoutsakos. Optimal morphokinematics for undu-
latory swimmers at intermediate reynolds numbers. Journal of Fluid Mechanics, 775:178–
188, 2015.

[100] T. Williams, G. Bowtell, and N. Curtin. Predicting force generation by lamprey muscle
during applied sinusoidal movement using a simple dynamic model. J. Exp. Biol., 201:869,
1998.

[101]

J.-Z. Wu, H.-Y. Ma, and M.-D. Zhou. Vorticity and Vortex Dynamics. Springer, 2006.

[102] C. Xuan and J. Biggins. Finite-wavelength surface-tension-driven instabilities in soft
solids, including instability in a cylindrical channel through an elastic solid. Phys. Rev.
E, 94(2):023107, 2016.

[103] X. Yang, X. Zhang, Z. Li, and G.-W. He. A smoothing technique for discrete delta functions
with application to immersed boundary method in moving boundary simulations. Journal
of Computational Physics, 228(20):7821–7836, 2009.

[104] Z. Yu. A DLM/FD method for ﬂuid/ﬂexible-body interactions. Journal of Computational

Physics, 207:1–27, 2005.

[105] Z. Yu, Y. Wang, and X. Shao. Numerical simulation of the ﬂapping of a three-dimensional

ﬂexible plate in uniform ﬂow. J. Sound Vib., 331:4448–4463, 2012.

[106] Z. S. Yu and X. M. Shao. A three-dimensional ﬁctitious domain method for the simulation

of ﬂuid-structure interactions. Journal of Hydrodynamics Ser B, 22(5):178–183, 2010.

[107] Z. Yu and X. Shao. A direct-forcing ﬁctitious domain method for particulate ﬂows. Journal

of Computational Physics, 227:292–314, 2007.

80

[108]

J. Zhang, S. Childress, A. Libchaber, and M. Shelley. Flexible ﬁlaments in a ﬂowing soap
ﬁlm as a model for one-dimensional ﬂags in a two-dimensional wind. Nature, 408:835–838,
2000.

[109]

J. Zhou, R. J. Adrian, S. Balachandar, and T. Kendall. Mechanisms for generating coherent
packets of hairpin vortices in channel ﬂow. Journal of ﬂuid mechanics, 387:353–396, 1999.

81